lib/Math/BigInt/Calc.pm

   1 package Math::BigInt::Calc;
   2
   3 use 5.005;
   4 use strict;
   5 # use warnings; # dont use warnings for older Perls
   6
   7 require Exporter;
   8 use vars qw/@ISA $VERSION/;
   9 @ISA = qw(Exporter);
  10
  11 $VERSION = '0.14';
  12
  13 # Package to store unsigned big integers in decimal and do math with them
  14
  15 # Internally the numbers are stored in an array with at least 1 element, no
  16 # leading zero parts (except the first) and in base 100000
  17
  18 # todo:
  19 # - fully remove funky $# stuff (maybe)
  20
  21 # USE_MUL: due to problems on certain os (os390, posix-bc) "* 1e-5" is used
  22 # instead of "/ 1e5" at some places, (marked with USE_MUL). Other platforms
  23 # BS2000, some Crays need USE_DIV instead.
  24 # The BEGIN block is used to determine which of the two variants gives the
  25 # correct result.
  26
  27 ##############################################################################
  28 # global constants, flags and accessory
  29
  30 # constants for easier life
  31 my $nan = 'NaN';
  32 my ($BASE,$RBASE,$BASE_LEN,$MAX_VAL);
  33
  34 sub _base_len
  35   {
  36   # set/get the BASE_LEN and assorted other, connected values
  37   # used only be the testsuite, set is used only by the BEGIN block below
  38   my $b = shift;
  39   if (defined $b)
  40     {
  41     $b = 8 if $b > 8;                   # cap, for VMS, OS/390 and other 64 bit
  42     $BASE_LEN = $b;
  43     $BASE = int("1e".$BASE_LEN);
  44     $RBASE = abs('1e-'.$BASE_LEN);      # see USE_MUL
  45     $MAX_VAL = $BASE-1;
  46     # print "BASE_LEN: $BASE_LEN MAX_VAL: $MAX_VAL\n";
  47     # print "int: ",int($BASE * $RBASE),"\n";
  48     if (int($BASE * $RBASE) == 0)               # should be 1
  49       {
  50       # must USE_MUL
  51       *{_mul} = \&_mul_use_mul;
  52       *{_div} = \&_div_use_mul;
  53       }
  54     else
  55       {
  56       # can USE_DIV instead
  57       *{_mul} = \&_mul_use_div;
  58       *{_div} = \&_div_use_div;
  59       }
  60     }
  61   $BASE_LEN;
  62   }
  63
  64 BEGIN
  65   {
  66   # from Daniel Pfeiffer: determine largest group of digits that is precisely
  67   # multipliable with itself plus carry
  68   # Test now changed to expect the proper pattern, not a result off by 1 or 2
  69   my ($e, $num) = 3;    # lowest value we will use is 3+1-1 = 3
  70   do
  71     {
  72     $num = ('9' x ++$e) + 0;
  73     $num *= $num + 1;
  74     # print "$num $e\n";
  75     } while ("$num" =~ /9{$e}0{$e}/);           # must be a certain pattern
  76   # last test failed, so retract one step:
  77   _base_len($e-1);
  78   }
  79
  80 ##############################################################################
  81 # create objects from various representations
  82
  83 sub _new
  84   {
  85   # (string) return ref to num_array
  86   # Convert a number from string format to internal base 100000 format.
  87   # Assumes normalized value as input.
  88   my $d = $_[1];
  89   # print "_new $d $$d\n";
  90   my $il = CORE::length($$d)-1;
  91   # these leaves '00000' instead of int 0 and will be corrected after any op
  92   return [ reverse(unpack("a" . ($il % $BASE_LEN+1)
  93     . ("a$BASE_LEN" x ($il / $BASE_LEN)), $$d)) ];
  94   }
  95
  96 sub _zero
  97   {
  98   # create a zero
  99   return [ 0 ];
 100   }
 101
 102 sub _one
 103   {
 104   # create a one
 105   return [ 1 ];
 106   }
 107
 108 sub _copy
 109   {
 110   return [ @{$_[1]} ];
 111   }
 112
 113 # catch and throw away
 114 sub import { }
 115
 116 ##############################################################################
 117 # convert back to string and number
 118
 119 sub _str
 120   {
 121   # (ref to BINT) return num_str
 122   # Convert number from internal base 100000 format to string format.
 123   # internal format is always normalized (no leading zeros, "-0" => "+0")
 124   my $ar = $_[1];
 125   my $ret = "";
 126   my $l = scalar @$ar;         # number of parts
 127   return $nan if $l < 1;       # should not happen
 128   # handle first one different to strip leading zeros from it (there are no
 129   # leading zero parts in internal representation)
 130   $l --; $ret .= $ar->[$l]; $l--;
 131   # Interestingly, the pre-padd method uses more time
 132   # the old grep variant takes longer (14 to 10 sec)
 133   my $z = '0' x ($BASE_LEN-1);
 134   while ($l >= 0)
 135     {
 136     $ret .= substr($z.$ar->[$l],-$BASE_LEN); # fastest way I could think of
 137     $l--;
 138     }
 139   return \$ret;
 140   }
 141
 142 sub _num
 143   {
 144   # Make a number (scalar int/float) from a BigInt object
 145   my $x = $_[1];
 146   return $x->[0] if scalar @$x == 1;  # below $BASE
 147   my $fac = 1;
 148   my $num = 0;
 149   foreach (@$x)
 150     {
 151     $num += $fac*$_; $fac *= $BASE;
 152     }
 153   return $num;
 154   }
 155
 156 ##############################################################################
 157 # actual math code
 158
 159 sub _add
 160   {
 161   # (ref to int_num_array, ref to int_num_array)
 162   # routine to add two base 1eX numbers
 163   # stolen from Knuth Vol 2 Algorithm A pg 231
 164   # there are separate routines to add and sub as per Knuth pg 233
 165   # This routine clobbers up array x, but not y.
 166
 167   my ($c,$x,$y) = @_;
 168
 169   # for each in Y, add Y to X and carry. If after that, something is left in
 170   # X, foreach in X add carry to X and then return X, carry
 171   # Trades one "$j++" for having to shift arrays, $j could be made integer
 172   # but this would impose a limit to number-length of 2**32.
 173   my $i; my $car = 0; my $j = 0;
 174   for $i (@$y)
 175     {
 176     $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
 177     $j++;
 178     }
 179   while ($car != 0)
 180     {
 181     $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
 182     }
 183   return $x;
 184   }
 185
 186 sub _inc
 187   {
 188   # (ref to int_num_array, ref to int_num_array)
 189   # routine to add 1 to a base 1eX numbers
 190   # This routine clobbers up array x, but not y.
 191   my ($c,$x) = @_;
 192
 193   for my $i (@$x)
 194     {
 195     return $x if (($i += 1) < $BASE);           # early out
 196     $i -= $BASE;
 197     }
 198   if ($x->[-1] == 0)                            # last overflowed
 199     {
 200     push @$x,1;                                 # extend
 201     }
 202   return $x;
 203   }
 204
 205 sub _dec
 206   {
 207   # (ref to int_num_array, ref to int_num_array)
 208   # routine to add 1 to a base 1eX numbers
 209   # This routine clobbers up array x, but not y.
 210   my ($c,$x) = @_;
 211
 212   for my $i (@$x)
 213     {
 214     last if (($i -= 1) >= 0);                   # early out
 215     $i = $MAX_VAL;
 216     }
 217   pop @$x if $x->[-1] == 0 && @$x > 1;          # last overflowed (but leave 0)
 218   return $x;
 219   }
 220
 221 sub _sub
 222   {
 223   # (ref to int_num_array, ref to int_num_array)
 224   # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
 225   # subtract Y from X (X is always greater/equal!) by modifying x in place
 226   my ($c,$sx,$sy,$s) = @_;
 227
 228   my $car = 0; my $i; my $j = 0;
 229   if (!$s)
 230     {
 231     #print "case 2\n";
 232     for $i (@$sx)
 233       {
 234       last unless defined $sy->[$j] || $car;
 235       #print "x: $i y: $sy->[$j] c: $car\n";
 236       $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
 237       #print "x: $i y: $sy->[$j-1] c: $car\n";
 238       }
 239     # might leave leading zeros, so fix that
 240     __strip_zeros($sx);
 241     return $sx;
 242     }
 243   else
 244     {
 245     #print "case 1 (swap)\n";
 246     for $i (@$sx)
 247       {
 248       last unless defined $sy->[$j] || $car;
 249       #print "$sy->[$j] $i $car => $sx->[$j]\n";
 250       $sy->[$j] += $BASE
 251        if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
 252       #print "$sy->[$j] $i $car => $sy->[$j]\n";
 253       $j++;
 254       }
 255     # might leave leading zeros, so fix that
 256     __strip_zeros($sy);
 257     return $sy;
 258     }
 259   }
 260
 261 sub _mul_use_mul
 262   {
 263   # (BINT, BINT) return nothing
 264   # multiply two numbers in internal representation
 265   # modifies first arg, second need not be different from first
 266   my ($c,$xv,$yv) = @_;
 267
 268   my @prod = (); my ($prod,$car,$cty,$xi,$yi);
 269   # since multiplying $x with $x fails, make copy in this case
 270   $yv = [@$xv] if "$xv" eq "$yv";       # same references?
 271   for $xi (@$xv)
 272     {
 273     $car = 0; $cty = 0;
 274
 275     # slow variant
 276 #    for $yi (@$yv)
 277 #      {
 278 #      $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
 279 #      $prod[$cty++] =
 280 #       $prod - ($car = int($prod * RBASE)) * $BASE;  # see USE_MUL
 281 #      }
 282 #    $prod[$cty] += $car if $car; # need really to check for 0?
 283 #    $xi = shift @prod;
 284
 285     # faster variant
 286     # looping through this if $xi == 0 is silly - so optimize it away!
 287     $xi = (shift @prod || 0), next if $xi == 0;
 288     for $yi (@$yv)
 289       {
 290       $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
 291 ##     this is actually a tad slower
 292 ##        $prod = $prod[$cty]; $prod += ($car + $xi * $yi);     # no ||0 here
 293       $prod[$cty++] =
 294        $prod - ($car = int($prod * $RBASE)) * $BASE;  # see USE_MUL
 295       }
 296     $prod[$cty] += $car if $car; # need really to check for 0?
 297     $xi = shift @prod;
 298     }
 299   push @$xv, @prod;
 300   __strip_zeros($xv);
 301   # normalize (handled last to save check for $y->is_zero()
 302   return $xv;
 303   }
 304
 305 sub _mul_use_div
 306   {
 307   # (BINT, BINT) return nothing
 308   # multiply two numbers in internal representation
 309   # modifies first arg, second need not be different from first
 310   my ($c,$xv,$yv) = @_;
 311
 312   my @prod = (); my ($prod,$car,$cty,$xi,$yi);
 313   # since multiplying $x with $x fails, make copy in this case
 314   $yv = [@$xv] if "$xv" eq "$yv";       # same references?
 315   for $xi (@$xv)
 316     {
 317     $car = 0; $cty = 0;
 318     # looping through this if $xi == 0 is silly - so optimize it away!
 319     $xi = (shift @prod || 0), next if $xi == 0;
 320     for $yi (@$yv)
 321       {
 322       $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
 323       $prod[$cty++] =
 324        $prod - ($car = int($prod / $BASE)) * $BASE;
 325       }
 326     $prod[$cty] += $car if $car; # need really to check for 0?
 327     $xi = shift @prod;
 328     }
 329   push @$xv, @prod;
 330   __strip_zeros($xv);
 331   # normalize (handled last to save check for $y->is_zero()
 332   return $xv;
 333   }
 334
 335 sub _div_use_mul
 336   {
 337   # ref to array, ref to array, modify first array and return remainder if
 338   # in list context
 339   # no longer handles sign
 340   my ($c,$x,$yorg) = @_;
 341   my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1);
 342
 343   my (@d,$tmp,$q,$u2,$u1,$u0);
 344
 345   $car = $bar = $prd = 0;
 346
 347   my $y = [ @$yorg ];
 348   if (($dd = int($BASE/($y->[-1]+1))) != 1)
 349     {
 350     for $xi (@$x)
 351       {
 352       $xi = $xi * $dd + $car;
 353       $xi -= ($car = int($xi * $RBASE)) * $BASE;        # see USE_MUL
 354       }
 355     push(@$x, $car); $car = 0;
 356     for $yi (@$y)
 357       {
 358       $yi = $yi * $dd + $car;
 359       $yi -= ($car = int($yi * $RBASE)) * $BASE;        # see USE_MUL
 360       }
 361     }
 362   else
 363     {
 364     push(@$x, 0);
 365     }
 366   @q = (); ($v2,$v1) = @$y[-2,-1];
 367   $v2 = 0 unless $v2;
 368   while ($#$x > $#$y)
 369     {
 370     ($u2,$u1,$u0) = @$x[-3..-1];
 371     $u2 = 0 unless $u2;
 372     #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
 373     # if $v1 == 0;
 374     # $q = (($u0 == $v1) ? 99999 : int(($u0*$BASE+$u1)/$v1));
 375      $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
 376     --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
 377     if ($q)
 378       {
 379       ($car, $bar) = (0,0);
 380       for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
 381         {
 382         $prd = $q * $y->[$yi] + $car;
 383         $prd -= ($car = int($prd * $RBASE)) * $BASE;    # see USE_MUL
 384         $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
 385         }
 386       if ($x->[-1] < $car + $bar)
 387         {
 388         $car = 0; --$q;
 389         for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
 390           {
 391           $x->[$xi] -= $BASE
 392            if ($car = (($x->[$xi] += $y->[$yi] + $car) > $BASE));
 393           }
 394         }
 395       }
 396       pop(@$x); unshift(@q, $q);
 397     }
 398   if (wantarray)
 399     {
 400     @d = ();
 401     if ($dd != 1)
 402       {
 403       $car = 0;
 404       for $xi (reverse @$x)
 405         {
 406         $prd = $car * $BASE + $xi;
 407         $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
 408         unshift(@d, $tmp);
 409         }
 410       }
 411     else
 412       {
 413       @d = @$x;
 414       }
 415     @$x = @q;
 416     __strip_zeros($x);
 417     __strip_zeros(\@d);
 418     return ($x,\@d);
 419     }
 420   @$x = @q;
 421   __strip_zeros($x);
 422   return $x;
 423   }
 424
 425 sub _div_use_div
 426   {
 427   # ref to array, ref to array, modify first array and return remainder if
 428   # in list context
 429   # no longer handles sign
 430   my ($c,$x,$yorg) = @_;
 431   my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1);
 432
 433   my (@d,$tmp,$q,$u2,$u1,$u0);
 434
 435   $car = $bar = $prd = 0;
 436
 437   my $y = [ @$yorg ];
 438   if (($dd = int($BASE/($y->[-1]+1))) != 1)
 439     {
 440     for $xi (@$x)
 441       {
 442       $xi = $xi * $dd + $car;
 443       $xi -= ($car = int($xi / $BASE)) * $BASE;
 444       }
 445     push(@$x, $car); $car = 0;
 446     for $yi (@$y)
 447       {
 448       $yi = $yi * $dd + $car;
 449       $yi -= ($car = int($yi / $BASE)) * $BASE;
 450       }
 451     }
 452   else
 453     {
 454     push(@$x, 0);
 455     }
 456   @q = (); ($v2,$v1) = @$y[-2,-1];
 457   $v2 = 0 unless $v2;
 458   while ($#$x > $#$y)
 459     {
 460     ($u2,$u1,$u0) = @$x[-3..-1];
 461     $u2 = 0 unless $u2;
 462     #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
 463     # if $v1 == 0;
 464     # $q = (($u0 == $v1) ? 99999 : int(($u0*$BASE+$u1)/$v1));
 465      $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
 466     --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
 467     if ($q)
 468       {
 469       ($car, $bar) = (0,0);
 470       for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
 471         {
 472         $prd = $q * $y->[$yi] + $car;
 473         $prd -= ($car = int($prd / $BASE)) * $BASE;
 474         $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
 475         }
 476       if ($x->[-1] < $car + $bar)
 477         {
 478         $car = 0; --$q;
 479         for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
 480           {
 481           $x->[$xi] -= $BASE
 482            if ($car = (($x->[$xi] += $y->[$yi] + $car) > $BASE));
 483           }
 484         }
 485       }
 486       pop(@$x); unshift(@q, $q);
 487     }
 488   if (wantarray)
 489     {
 490     @d = ();
 491     if ($dd != 1)
 492       {
 493       $car = 0;
 494       for $xi (reverse @$x)
 495         {
 496         $prd = $car * $BASE + $xi;
 497         $car = $prd - ($tmp = int($prd / $dd)) * $dd;
 498         unshift(@d, $tmp);
 499         }
 500       }
 501     else
 502       {
 503       @d = @$x;
 504       }
 505     @$x = @q;
 506     __strip_zeros($x);
 507     __strip_zeros(\@d);
 508     return ($x,\@d);
 509     }
 510   @$x = @q;
 511   __strip_zeros($x);
 512   return $x;
 513   }
 514
 515 sub _mod
 516   {
 517   # if possible, use mod shortcut
 518   my ($c,$x,$yo) = @_;
 519
 520   # slow way since $y to big
 521   if (scalar @$yo > 1)
 522     {
 523     my ($xo,$rem) = _div($c,$x,$yo);
 524     return $rem;
 525     }
 526   my $y = $yo->[0];
 527   # both are single element
 528   if (scalar @$x == 1)
 529     {
 530     $x->[0] %= $y;
 531     return $x;
 532     }
 533
 534   my $b = $BASE % $y;
 535   if ($b == 0)
 536     {
 537     # when BASE % Y == 0 then (B * BASE) % Y == 0
 538     # (B * BASE) % $y + A % Y => A % Y
 539     # so need to consider only last element: O(1)
 540     $x->[0] %= $y;
 541     }
 542   else
 543     {
 544     # else need to go trough all elemens: O(N)
 545     # XXX not ready yet
 546     my ($xo,$rem) = _div($c,$x,$yo);
 547     return $rem;
 548
 549 #    my $i = 0; my $r = 1;
 550 #    print "Multi: ";
 551 #    foreach (@$x)
 552 #      {
 553 #      print "$_ $r $b $y\n";
 554 #      print "\$_ % \$y = ",$_ % $y,"\n";
 555 #      print "\$_ % \$y * \$b = ",($_ % $y) * $b,"\n";
 556 #      $r += ($_ % $y) * $b;
 557 #      print "$r $b $y =>";
 558 #      $r %= $y if $r > $y;
 559 #      print " $r\n";
 560 #      }
 561 #    $x->[0] = $r;
 562     }
 563   splice (@$x,1);
 564   return $x;
 565   }
 566
 567 ##############################################################################
 568 # shifts
 569
 570 sub _rsft
 571   {
 572   my ($c,$x,$y,$n) = @_;
 573
 574   if ($n != 10)
 575     {
 576     return;     # we cant do this here, due to now _pow, so signal failure
 577     }
 578   else
 579     {
 580     # shortcut (faster) for shifting by 10)
 581     # multiples of $BASE_LEN
 582     my $dst = 0;                                # destination
 583     my $src = _num($c,$y);                      # as normal int
 584     my $rem = $src % $BASE_LEN;                 # remainder to shift
 585     $src = int($src / $BASE_LEN);               # source
 586     if ($rem == 0)
 587       {
 588       splice (@$x,0,$src);                      # even faster, 38.4 => 39.3
 589       }
 590     else
 591       {
 592       my $len = scalar @$x - $src;              # elems to go
 593       my $vd; my $z = '0'x $BASE_LEN;
 594       $x->[scalar @$x] = 0;                     # avoid || 0 test inside loop
 595       while ($dst < $len)
 596         {
 597         $vd = $z.$x->[$src];
 598         #print "$dst $src '$vd' ";
 599         $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
 600         #print "'$vd' ";
 601         $src++;
 602         $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
 603         #print "'$vd1' ";
 604         #print "'$vd'\n";
 605         $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
 606         $x->[$dst] = int($vd);
 607         $dst++;
 608         }
 609       splice (@$x,$dst) if $dst > 0;            # kill left-over array elems
 610       pop @$x if $x->[-1] == 0;                 # kill last element if 0
 611       } # else rem == 0
 612     }
 613   $x;
 614   }
 615
 616 sub _lsft
 617   {
 618   my ($c,$x,$y,$n) = @_;
 619
 620   if ($n != 10)
 621     {
 622     return;     # we cant do this here, due to now _pow, so signal failure
 623     }
 624   else
 625     {
 626     # shortcut (faster) for shifting by 10) since we are in base 10eX
 627     # multiples of $BASE_LEN:
 628     my $src = scalar @$x;                       # source
 629     my $len = _num($c,$y);                      # shift-len as normal int
 630     my $rem = $len % $BASE_LEN;                 # remainder to shift
 631     my $dst = $src + int($len/$BASE_LEN);       # destination
 632     my $vd;                                     # further speedup
 633     #print "src $src:",$x->[$src]||0," dst $dst:",$v->[$dst]||0," rem $rem\n";
 634     $x->[$src] = 0;                             # avoid first ||0 for speed
 635     my $z = '0' x $BASE_LEN;
 636     while ($src >= 0)
 637       {
 638       $vd = $x->[$src]; $vd = $z.$vd;
 639       #print "s $src d $dst '$vd' ";
 640       $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
 641       #print "'$vd' ";
 642       $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
 643       #print "'$vd' ";
 644       $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
 645       #print "'$vd'\n";
 646       $x->[$dst] = int($vd);
 647       $dst--; $src--;
 648       }
 649     # set lowest parts to 0
 650     while ($dst >= 0) { $x->[$dst--] = 0; }
 651     # fix spurios last zero element
 652     splice @$x,-1 if $x->[-1] == 0;
 653     #print "elems: "; my $i = 0;
 654     #foreach (reverse @$v) { print "$i $_ "; $i++; } print "\n";
 655     }
 656   $x;
 657   }
 658
 659 ##############################################################################
 660 # testing
 661
 662 sub _acmp
 663   {
 664   # internal absolute post-normalized compare (ignore signs)
 665   # ref to array, ref to array, return <0, 0, >0
 666   # arrays must have at least one entry; this is not checked for
 667
 668   my ($c,$cx, $cy) = @_;
 669
 670   #print "$cx $cy\n";
 671   my ($i,$a,$x,$y,$k);
 672   # calculate length based on digits, not parts
 673   $x = _len('',$cx); $y = _len('',$cy);
 674   # print "length: ",($x-$y),"\n";
 675   my $lxy = $x - $y;                            # if different in length
 676   return -1 if $lxy < 0;
 677   return 1 if $lxy > 0;
 678   #print "full compare\n";
 679   $i = 0; $a = 0;
 680   # first way takes 5.49 sec instead of 4.87, but has the early out advantage
 681   # so grep is slightly faster, but more inflexible. hm. $_ instead of $k
 682   # yields 5.6 instead of 5.5 sec huh?
 683   # manual way (abort if unequal, good for early ne)
 684   my $j = scalar @$cx - 1;
 685   while ($j >= 0)
 686    {
 687    # print "$cx->[$j] $cy->[$j] $a",$cx->[$j]-$cy->[$j],"\n";
 688    last if ($a = $cx->[$j] - $cy->[$j]); $j--;
 689    }
 690   return 1 if $a > 0;
 691   return -1 if $a < 0;
 692   return 0;                                     # equal
 693   # while it early aborts, it is even slower than the manual variant
 694   #grep { return $a if ($a = $_ - $cy->[$i++]); } @$cx;
 695   # grep way, go trough all (bad for early ne)
 696   #grep { $a = $_ - $cy->[$i++]; } @$cx;
 697   #return $a;
 698   }
 699
 700 sub _len
 701   {
 702   # compute number of digits in bigint, minus the sign
 703   # int() because add/sub sometimes leaves strings (like '00005') instead of
 704   # int ('5') in this place, thus causing length() to report wrong length
 705   my $cx = $_[1];
 706
 707   return (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
 708   }
 709
 710 sub _digit
 711   {
 712   # return the nth digit, negative values count backward
 713   # zero is rightmost, so _digit(123,0) will give 3
 714   my ($c,$x,$n) = @_;
 715
 716   my $len = _len('',$x);
 717
 718   $n = $len+$n if $n < 0;               # -1 last, -2 second-to-last
 719   $n = abs($n);                         # if negative was too big
 720   $len--; $n = $len if $n > $len;       # n to big?
 721
 722   my $elem = int($n / $BASE_LEN);       # which array element
 723   my $digit = $n % $BASE_LEN;           # which digit in this element
 724   $elem = '0000'.@$x[$elem];            # get element padded with 0's
 725   return substr($elem,-$digit-1,1);
 726   }
 727
 728 sub _zeros
 729   {
 730   # return amount of trailing zeros in decimal
 731   # check each array elem in _m for having 0 at end as long as elem == 0
 732   # Upon finding a elem != 0, stop
 733   my $x = $_[1];
 734   my $zeros = 0; my $elem;
 735   foreach my $e (@$x)
 736     {
 737     if ($e != 0)
 738       {
 739       $elem = "$e";                             # preserve x
 740       $elem =~ s/.*?(0*$)/$1/;                  # strip anything not zero
 741       $zeros *= $BASE_LEN;                      # elems * 5
 742       $zeros += CORE::length($elem);            # count trailing zeros
 743       last;                                     # early out
 744       }
 745     $zeros ++;                                  # real else branch: 50% slower!
 746     }
 747   return $zeros;
 748   }
 749
 750 ##############################################################################
 751 # _is_* routines
 752
 753 sub _is_zero
 754   {
 755   # return true if arg (BINT or num_str) is zero (array '+', '0')
 756   my $x = $_[1];
 757   return (((scalar @$x == 1) && ($x->[0] == 0))) <=> 0;
 758   }
 759
 760 sub _is_even
 761   {
 762   # return true if arg (BINT or num_str) is even
 763   my $x = $_[1];
 764   return (!($x->[0] & 1)) <=> 0;
 765   }
 766
 767 sub _is_odd
 768   {
 769   # return true if arg (BINT or num_str) is even
 770   my $x = $_[1];
 771   return (($x->[0] & 1)) <=> 0;
 772   }
 773
 774 sub _is_one
 775   {
 776   # return true if arg (BINT or num_str) is one (array '+', '1')
 777   my $x = $_[1];
 778   return (scalar @$x == 1) && ($x->[0] == 1) <=> 0;
 779   }
 780
 781 sub __strip_zeros
 782   {
 783   # internal normalization function that strips leading zeros from the array
 784   # args: ref to array
 785   my $s = shift;
 786
 787   my $cnt = scalar @$s; # get count of parts
 788   my $i = $cnt-1;
 789   #print "strip: cnt $cnt i $i\n";
 790   # '0', '3', '4', '0', '0',
 791   #  0    1    2    3    4
 792   # cnt = 5, i = 4
 793   # i = 4
 794   # i = 3
 795   # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
 796   # >= 1: skip first part (this can be zero)
 797   while ($i > 0) { last if $s->[$i] != 0; $i--; }
 798   $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
 799   return $s;
 800   }
 801
 802 ###############################################################################
 803 # check routine to test internal state of corruptions
 804
 805 sub _check
 806   {
 807   # used by the test suite
 808   my $x = $_[1];
 809
 810   return "$x is not a reference" if !ref($x);
 811
 812   # are all parts are valid?
 813   my $i = 0; my $j = scalar @$x; my ($e,$try);
 814   while ($i < $j)
 815     {
 816     $e = $x->[$i]; $e = 'undef' unless defined $e;
 817     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
 818     last if $e !~ /^[+]?[0-9]+$/;
 819     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
 820     last if "$e" !~ /^[+]?[0-9]+$/;
 821     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
 822     last if '' . "$e" !~ /^[+]?[0-9]+$/;
 823     $try = ' < 0 || >= $BASE; '."($x, $e)";
 824     last if $e <0 || $e >= $BASE;
 825     # this test is disabled, since new/bnorm and certain ops (like early out
 826     # in add/sub) are allowed/expected to leave '00000' in some elements
 827     #$try = '=~ /^00+/; '."($x, $e)";
 828     #last if $e =~ /^00+/;
 829     $i++;
 830     }
 831   return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
 832   return 0;
 833   }
 834
 835 1;
 836 __END__
 837
 838 =head1 NAME
 839
 840 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
 841
 842 =head1 SYNOPSIS
 843
 844 Provides support for big integer calculations. Not intended to be used by other
 845 modules (except Math::BigInt::Cached). Other modules which sport the same
 846 functions can also be used to support Math::Bigint, like Math::BigInt::Pari.
 847
 848 =head1 DESCRIPTION
 849
 850 In order to allow for multiple big integer libraries, Math::BigInt
 851 was rewritten to use library modules for core math routines. Any
 852 module which follows the same API as this can be used instead by
 853 using the following call:
 854
 855         use Math::BigInt lib => 'libname';
 856
 857 =head1 EXPORT
 858
 859 The following functions MUST be defined in order to support
 860 the use by Math::BigInt:
 861
 862         _new(string)    return ref to new object from ref to decimal string
 863         _zero()         return a new object with value 0
 864         _one()          return a new object with value 1
 865
 866         _str(obj)       return ref to a string representing the object
 867         _num(obj)       returns a Perl integer/floating point number
 868                         NOTE: because of Perl numeric notation defaults,
 869                         the _num'ified obj may lose accuracy due to
 870                         machine-dependend floating point size limitations
 871
 872         _add(obj,obj)   Simple addition of two objects
 873         _mul(obj,obj)   Multiplication of two objects
 874         _div(obj,obj)   Division of the 1st object by the 2nd
 875                         In list context, returns (result,remainder).
 876                         NOTE: this is integer math, so no
 877                         fractional part will be returned.
 878         _sub(obj,obj)   Simple subtraction of 1 object from another
 879                         a third, optional parameter indicates that the params
 880                         are swapped. In this case, the first param needs to
 881                         be preserved, while you can destroy the second.
 882                         sub (x,y,1) => return x - y and keep x intact!
 883         _dec(obj)       decrement object by one (input is garant. to be > 0)
 884         _inc(obj)       increment object by one
 885
 886
 887         _acmp(obj,obj)  <=> operator for objects (return -1, 0 or 1)
 888
 889         _len(obj)       returns count of the decimal digits of the object
 890         _digit(obj,n)   returns the n'th decimal digit of object
 891
 892         _is_one(obj)    return true if argument is +1
 893         _is_zero(obj)   return true if argument is 0
 894         _is_even(obj)   return true if argument is even (0,2,4,6..)
 895         _is_odd(obj)    return true if argument is odd (1,3,5,7..)
 896
 897         _copy           return a ref to a true copy of the object
 898
 899         _check(obj)     check whether internal representation is still intact
 900                         return 0 for ok, otherwise error message as string
 901
 902 The following functions are optional, and can be defined if the underlying lib
 903 has a fast way to do them. If undefined, Math::BigInt will use a pure, but
 904 slow, Perl way as fallback to emulate these:
 905
 906         _from_hex(str)  return ref to new object from ref to hexadecimal string
 907         _from_bin(str)  return ref to new object from ref to binary string
 908
 909         _as_hex(str)    return ref to scalar string containing the value as
 910                         unsigned hex string, with the '0x' prepended.
 911                         Leading zeros must be stripped.
 912         _as_bin(str)    Like as_hex, only as binary string containing only
 913                         zeros and ones. Leading zeros must be stripped and a
 914                         '0b' must be prepended.
 915
 916         _rsft(obj,N,B)  shift object in base B by N 'digits' right
 917                         For unsupported bases B, return undef to signal failure
 918         _lsft(obj,N,B)  shift object in base B by N 'digits' left
 919                         For unsupported bases B, return undef to signal failure
 920
 921         _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
 922                         Note: XOR, AND and OR pad with zeros if size mismatches
 923         _and(obj1,obj2) AND (bit-wise) object 1 with object 2
 924         _or(obj1,obj2)  OR (bit-wise) object 1 with object 2
 925
 926         _mod(obj,obj)   Return remainder of div of the 1st by the 2nd object
 927         _sqrt(obj)      return the square root of object
 928         _pow(obj,obj)   return object 1 to the power of object 2
 929         _gcd(obj,obj)   return Greatest Common Divisor of two objects
 930
 931         _zeros(obj)     return number of trailing decimal zeros
 932
 933 Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
 934 or '0b1101').
 935
 936 Testing of input parameter validity is done by the caller, so you need not
 937 worry about underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by
 938 zero or similar cases.
 939
 940 The first parameter can be modified, that includes the possibility that you
 941 return a reference to a completely different object instead. Although keeping
 942 the reference and just changing it's contents is prefered over creating and
 943 returning a different reference.
 944
 945 Return values are always references to objects or strings. Exceptions are
 946 C<_lsft()> and C<_rsft()>, which return undef if they can not shift the
 947 argument. This is used to delegate shifting of bases different than 10 back
 948 to Math::BigInt, which will use some generic code to calculate the result.
 949
 950 =head1 WRAP YOUR OWN
 951
 952 If you want to port your own favourite c-lib for big numbers to the
 953 Math::BigInt interface, you can take any of the already existing modules as
 954 a rough guideline. You should really wrap up the latest BigInt and BigFloat
 955 testsuites with your module, and replace in them any of the following:
 956
 957         use Math::BigInt;
 958
 959 by this:
 960
 961         use Math::BigInt lib => 'yourlib';
 962
 963 This way you ensure that your library really works 100% within Math::BigInt.
 964
 965 =head1 LICENSE
 966
 967 This program is free software; you may redistribute it and/or modify it under
 968 the same terms as Perl itself.
 969
 970 =head1 AUTHORS
 971
 972 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
 973 in late 2000, 2001.
 974 Seperated from BigInt and shaped API with the help of John Peacock.
 975
 976 =head1 SEE ALSO
 977
 978 L<Math::BigInt>, L<Math::BigFloat>, L<Math::BigInt::BitVect>,
 979 L<Math::BigInt::GMP>, L<Math::BigInt::Cached> and L<Math::BigInt::Pari>.
 980
 981 =cut