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