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