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