1 package Math::BigInt::Calc;
5 # use warnings; # dont use warnings for older Perls
8 use vars qw/@ISA $VERSION/;
13 # Package to store unsigned big integers in decimal and do math with them
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
20 # - fully remove funky $# stuff (maybe)
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
28 ##############################################################################
29 # global constants, flags and accessory
31 # constants for easier life
33 my ($MBASE,$BASE,$RBASE,$BASE_LEN,$MAX_VAL,$BASE_LEN2,$BASE_LEN_SMALL);
34 my ($AND_BITS,$XOR_BITS,$OR_BITS);
35 my ($AND_MASK,$XOR_MASK,$OR_MASK);
40 # set/get the BASE_LEN and assorted other, connected values
41 # used only be the testsuite, set is used only by the BEGIN block below
47 # find whether we can use mul or div or none in mul()/div()
48 # (in last case reduce BASE_LEN_SMALL)
49 $BASE_LEN_SMALL = $b+1;
51 while (--$BASE_LEN_SMALL > 5)
53 $MBASE = int("1e".$BASE_LEN_SMALL);
54 $RBASE = abs('1e-'.$BASE_LEN_SMALL); # see USE_MUL
56 $caught += 1 if (int($MBASE * $RBASE) != 1); # should be 1
57 $caught += 2 if (int($MBASE / $MBASE) != 1); # should be 1
60 # BASE_LEN is used for anything else than mul()/div()
61 $BASE_LEN = $BASE_LEN_SMALL;
62 $BASE_LEN = shift if (defined $_[0]); # one more arg?
63 $BASE = int("1e".$BASE_LEN);
65 $BASE_LEN2 = int($BASE_LEN_SMALL / 2); # for mul shortcut
66 $MBASE = int("1e".$BASE_LEN_SMALL);
67 $RBASE = abs('1e-'.$BASE_LEN_SMALL); # see USE_MUL
70 $LEN_CONVERT = 1 if $BASE_LEN_SMALL != $BASE_LEN;
72 #print "BASE_LEN: $BASE_LEN MAX_VAL: $MAX_VAL BASE: $BASE RBASE: $RBASE ";
73 #print "BASE_LEN_SMALL: $BASE_LEN_SMALL MBASE: $MBASE\n";
78 *{_mul} = \&_mul_use_mul;
79 *{_div} = \&_div_use_mul;
81 else # $caught must be 2, since it can't be 1 nor 3
84 *{_mul} = \&_mul_use_div;
85 *{_div} = \&_div_use_div;
88 return $BASE_LEN unless wantarray;
89 return ($BASE_LEN, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN_SMALL, $MAX_VAL);
94 # from Daniel Pfeiffer: determine largest group of digits that is precisely
95 # multipliable with itself plus carry
96 # Test now changed to expect the proper pattern, not a result off by 1 or 2
97 my ($e, $num) = 3; # lowest value we will use is 3+1-1 = 3
100 $num = ('9' x ++$e) + 0;
102 } while ("$num" =~ /9{$e}0{$e}/); # must be a certain pattern
103 $e--; # last test failed, so retract one step
104 # the limits below brush the problems with the test above under the rug:
105 # the test should be able to find the proper $e automatically
106 $e = 5 if $^O =~ /^uts/; # UTS get's some special treatment
107 $e = 5 if $^O =~ /^unicos/; # unicos is also problematic (6 seems to work
108 # there, but we play safe)
109 $e = 8 if $e > 8; # cap, for VMS, OS/390 and other 64 bit systems
111 # determine how many digits fit into an integer and can be safely added
112 # together plus carry w/o causing an overflow
114 # this below detects 15 on a 64 bit system, because after that it becomes
115 # 1e16 and not 1000000 :/ I can make it detect 18, but then I get a lot of
116 # test failures. Ugh! (Tomake detect 18: uncomment lines marked with *)
118 my $bi = 5; # approx. 16 bit
119 $num = int('9' x $bi);
121 # while ( ($num+$num+1) eq '1' . '9' x $bi) # *
122 while ( int($num+$num+1) eq '1' . '9' x $bi)
124 $bi++; $num = int('9' x $bi);
125 # $bi++; $num *= 10; $num += 9; # *
127 $bi--; # back off one step
128 # by setting them equal, we ignore the findings and use the default
129 # one-size-fits-all approach from former versions
130 $bi = $e; # XXX, this should work always
132 __PACKAGE__->_base_len($e,$bi); # set and store
134 # find out how many bits _and, _or and _xor can take (old default = 16)
135 # I don't think anybody has yet 128 bit scalars, so let's play safe.
136 local $^W = 0; # don't warn about 'nonportable number'
137 $AND_BITS = 15; $XOR_BITS = 15; $OR_BITS = 15;
139 # find max bits, we will not go higher than numberofbits that fit into $BASE
140 # to make _and etc simpler (and faster for smaller, slower for large numbers)
142 while (2 ** $max < $BASE) { $max++; }
146 $x = oct('0b' . '1' x $AND_BITS); $y = $x & $x;
147 $z = (2 ** $AND_BITS) - 1;
148 } while ($AND_BITS < $max && $x == $z && $y == $x);
149 $AND_BITS --; # retreat one step
152 $x = oct('0b' . '1' x $XOR_BITS); $y = $x ^ 0;
153 $z = (2 ** $XOR_BITS) - 1;
154 } while ($XOR_BITS < $max && $x == $z && $y == $x);
155 $XOR_BITS --; # retreat one step
158 $x = oct('0b' . '1' x $OR_BITS); $y = $x | $x;
159 $z = (2 ** $OR_BITS) - 1;
160 } while ($OR_BITS < $max && $x == $z && $y == $x);
161 $OR_BITS --; # retreat one step
165 ##############################################################################
166 # convert between the "small" and the "large" representation
170 # take an array in base $BASE_LEN_SMALL and convert it in-place to $BASE_LEN
173 # print "_to_large $BASE_LEN_SMALL => $BASE_LEN\n";
175 return $x if $LEN_CONVERT == 0 || # nothing to converconvertor
176 @$x == 1; # only one element => early out
178 # 12345 67890 12345 67890 contents
180 # 123456 7890123 4567890 contents
183 # my @d; my $str = '';
184 # my $z = '0' x $BASE_LEN_SMALL;
187 # # ... . 04321 . 000321
188 # $str = substr($z.$_,-$BASE_LEN_SMALL,$BASE_LEN_SMALL) . $str;
189 # if (length($str) > $BASE_LEN)
191 # push @d, substr($str,-$BASE_LEN,$BASE_LEN); # extract one piece
192 # substr($str,-$BASE_LEN,$BASE_LEN) = ''; # remove it
195 # push @d, $str if $str !~ /^0*$/; # extract last piece
197 # $x->[-1] = int($x->[-1]); # strip leading zero
201 my $l = scalar @$x; # number of parts
202 $l --; $ret .= int($x->[$l]); $l--;
203 my $z = '0' x ($BASE_LEN_SMALL-1);
206 $ret .= substr($z.$x->[$l],-$BASE_LEN_SMALL);
209 my $str = _new($c,\$ret); # make array
210 @$x = @$str; # clobber contents of $x
211 $x->[-1] = int($x->[-1]); # strip leading zero
216 # take an array in base $BASE_LEN and convert it in-place to $BASE_LEN_SMALL
219 return $x if $LEN_CONVERT == 0; # nothing to do
220 return $x if @$x == 1 && length(int($x->[0])) <= $BASE_LEN_SMALL;
223 my $il = length($$d)-1;
224 ## this leaves '00000' instead of int 0 and will be corrected after any op
225 # clobber contents of $x
226 @$x = reverse(unpack("a" . ($il % $BASE_LEN_SMALL+1)
227 . ("a$BASE_LEN_SMALL" x ($il / $BASE_LEN_SMALL)), $$d));
229 $x->[-1] = int($x->[-1]); # strip leading zero
232 ###############################################################################
236 # (ref to string) return ref to num_array
237 # Convert a number from string format (without sign) to internal base
238 # 1ex format. Assumes normalized value as input.
240 my $il = length($$d)-1;
241 # this leaves '00000' instead of int 0 and will be corrected after any op
242 [ reverse(unpack("a" . ($il % $BASE_LEN+1)
243 . ("a$BASE_LEN" x ($il / $BASE_LEN)), $$d)) ];
248 $AND_MASK = __PACKAGE__->_new( \( 2 ** $AND_BITS ));
249 $XOR_MASK = __PACKAGE__->_new( \( 2 ** $XOR_BITS ));
250 $OR_MASK = __PACKAGE__->_new( \( 2 ** $OR_BITS ));
267 # create a two (for _pow)
276 # catch and throw away
279 ##############################################################################
280 # convert back to string and number
284 # (ref to BINT) return num_str
285 # Convert number from internal base 100000 format to string format.
286 # internal format is always normalized (no leading zeros, "-0" => "+0")
290 my $l = scalar @$ar; # number of parts
291 return $nan if $l < 1; # should not happen
293 # handle first one different to strip leading zeros from it (there are no
294 # leading zero parts in internal representation)
295 $l --; $ret .= int($ar->[$l]); $l--;
296 # Interestingly, the pre-padd method uses more time
297 # the old grep variant takes longer (14 to 10 sec)
298 my $z = '0' x ($BASE_LEN-1);
301 $ret .= substr($z.$ar->[$l],-$BASE_LEN); # fastest way I could think of
309 # Make a number (scalar int/float) from a BigInt object
311 return $x->[0] if scalar @$x == 1; # below $BASE
316 $num += $fac*$_; $fac *= $BASE;
321 ##############################################################################
326 # (ref to int_num_array, ref to int_num_array)
327 # routine to add two base 1eX numbers
328 # stolen from Knuth Vol 2 Algorithm A pg 231
329 # there are separate routines to add and sub as per Knuth pg 233
330 # This routine clobbers up array x, but not y.
334 return $x if (@$y == 1) && $y->[0] == 0; # $x + 0 => $x
335 if ((@$x == 1) && $x->[0] == 0) # 0 + $y => $y->copy
337 # twice as slow as $x = [ @$y ], but necc. to retain $x as ref :(
338 @$x = @$y; return $x;
341 # for each in Y, add Y to X and carry. If after that, something is left in
342 # X, foreach in X add carry to X and then return X, carry
343 # Trades one "$j++" for having to shift arrays, $j could be made integer
344 # but this would impose a limit to number-length of 2**32.
345 my $i; my $car = 0; my $j = 0;
348 $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
353 $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
360 # (ref to int_num_array, ref to int_num_array)
361 # routine to add 1 to a base 1eX numbers
362 # This routine clobbers up array x, but not y.
367 return $x if (($i += 1) < $BASE); # early out
368 $i = 0; # overflow, next
370 push @$x,1 if ($x->[-1] == 0); # last overflowed, so extend
376 # (ref to int_num_array, ref to int_num_array)
377 # routine to add 1 to a base 1eX numbers
378 # This routine clobbers up array x, but not y.
381 my $MAX = $BASE-1; # since MAX_VAL based on MBASE
384 last if (($i -= 1) >= 0); # early out
385 $i = $MAX; # overflow, next
387 pop @$x if $x->[-1] == 0 && @$x > 1; # last overflowed (but leave 0)
393 # (ref to int_num_array, ref to int_num_array, swap)
394 # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
395 # subtract Y from X by modifying x in place
396 my ($c,$sx,$sy,$s) = @_;
398 my $car = 0; my $i; my $j = 0;
404 last unless defined $sy->[$j] || $car;
405 $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
407 # might leave leading zeros, so fix that
408 return __strip_zeros($sx);
410 #print "case 1 (swap)\n";
413 # we can't do an early out if $x is than $y, since we
414 # need to copy the high chunks from $y. Found by Bob Mathews.
415 #last unless defined $sy->[$j] || $car;
417 if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
420 # might leave leading zeros, so fix that
426 # compute $x ** 2 or $x * $x in-place and return $x
429 # From: Handbook of Applied Cryptography by A. Menezes, P. van Oorschot and
430 # S. Vanstone., Chapter 14
432 #14.16 Algorithm Multiple-precision squaring
433 #INPUT: positive integer x = (xt 1 xt 2 ... x1 x0)b.
434 #OUTPUT: x * x = x ** 2 in radix b representation.
435 #1. For i from 0 to (2t - 1) do: wi <- 0.
436 #2. For i from 0 to (t - 1) do the following:
437 # 2.1 (uv)b w2i + xi * xi, w2i v, c u.
438 # 2.2 For j from (i + 1)to (t - 1) do the following:
439 # (uv)b <- wi+j + 2*xj * xi + c, wi+j <- v, c <- u.
441 #3. Return((w2t-1 w2t-2 ... w1 w0)b).
443 # # Note: That description is crap. Half of the symbols are not explained or
444 # # used with out beeing set.
445 # my $t = scalar @$x; # count
447 # for ($i = 0; $i < $t; $i++)
449 # $x->[$i] = $x->[$i*2] + $x[$i]*$x[$i];
450 # $x->[$i*2] = $x[$i]; $c = $x[$i];
451 # for ($j = $i+1; $j < $t; $j++)
453 # $x->[$i] = $x->[$i+$j] + 2 * $x->[$i] * $x->[$j];
454 # $x->[$i+$j] = $x[$j]; $c = $x[$i];
456 # $x->[$i+$t] = $x[$i];
463 # (ref to int_num_array, ref to int_num_array)
464 # multiply two numbers in internal representation
465 # modifies first arg, second need not be different from first
466 my ($c,$xv,$yv) = @_;
468 # shortcut for two very short numbers (improved by Nathan Zook)
469 # works also if xv and yv are the same reference
470 if ((@$xv == 1) && (@$yv == 1))
472 if (($xv->[0] *= $yv->[0]) >= $MBASE)
474 $xv->[0] = $xv->[0] - ($xv->[1] = int($xv->[0] * $RBASE)) * $MBASE;
478 # shortcut for result == 0
479 if ( ((@$xv == 1) && ($xv->[0] == 0)) ||
480 ((@$yv == 1) && ($yv->[0] == 0)) )
486 # since multiplying $x with $x fails, make copy in this case
487 $yv = [@$xv] if "$xv" eq "$yv"; # same references?
488 # since multiplying $x with $x would fail here, use the faster squaring
489 # return _square($c,$xv) if "$xv" eq "$yv"; # same reference?
491 if ($LEN_CONVERT != 0)
493 $c->_to_small($xv); $c->_to_small($yv);
496 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
505 # $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
507 # $prod - ($car = int($prod * RBASE)) * $MBASE; # see USE_MUL
509 # $prod[$cty] += $car if $car; # need really to check for 0?
513 # looping through this if $xi == 0 is silly - so optimize it away!
514 $xi = (shift @prod || 0), next if $xi == 0;
517 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
518 ## this is actually a tad slower
519 ## $prod = $prod[$cty]; $prod += ($car + $xi * $yi); # no ||0 here
521 $prod - ($car = int($prod * $RBASE)) * $MBASE; # see USE_MUL
523 $prod[$cty] += $car if $car; # need really to check for 0?
524 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
527 if ($LEN_CONVERT != 0)
541 # (ref to int_num_array, ref to int_num_array)
542 # multiply two numbers in internal representation
543 # modifies first arg, second need not be different from first
544 my ($c,$xv,$yv) = @_;
546 # shortcut for two very short numbers (improved by Nathan Zook)
547 # works also if xv and yv are the same reference
548 if ((@$xv == 1) && (@$yv == 1))
550 if (($xv->[0] *= $yv->[0]) >= $MBASE)
553 $xv->[0] - ($xv->[1] = int($xv->[0] / $MBASE)) * $MBASE;
557 # shortcut for result == 0
558 if ( ((@$xv == 1) && ($xv->[0] == 0)) ||
559 ((@$yv == 1) && ($yv->[0] == 0)) )
566 # since multiplying $x with $x fails, make copy in this case
567 $yv = [@$xv] if "$xv" eq "$yv"; # same references?
568 # since multiplying $x with $x would fail here, use the faster squaring
569 # return _square($c,$xv) if "$xv" eq "$yv"; # same reference?
571 if ($LEN_CONVERT != 0)
573 $c->_to_small($xv); $c->_to_small($yv);
576 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
580 # looping through this if $xi == 0 is silly - so optimize it away!
581 $xi = (shift @prod || 0), next if $xi == 0;
584 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
586 $prod - ($car = int($prod / $MBASE)) * $MBASE;
588 $prod[$cty] += $car if $car; # need really to check for 0?
589 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
592 if ($LEN_CONVERT != 0)
606 # ref to array, ref to array, modify first array and return remainder if
608 my ($c,$x,$yorg) = @_;
610 if (@$x == 1 && @$yorg == 1)
612 # shortcut, $yorg and $x are two small numbers
615 my $r = [ $x->[0] % $yorg->[0] ];
616 $x->[0] = int($x->[0] / $yorg->[0]);
621 $x->[0] = int($x->[0] / $yorg->[0]);
628 $rem = _mod($c,[ @$x ],$yorg) if wantarray;
630 # shortcut, $y is < $BASE
631 my $j = scalar @$x; my $r = 0;
632 my $y = $yorg->[0]; my $b;
635 $b = $r * $MBASE + $x->[$j];
636 $x->[$j] = int($b/$y);
639 pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
640 return ($x,$rem) if wantarray;
645 if ($LEN_CONVERT != 0)
647 $c->_to_small($x); $c->_to_small($y);
650 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
652 $car = $bar = $prd = 0;
653 if (($dd = int($MBASE/($y->[-1]+1))) != 1)
657 $xi = $xi * $dd + $car;
658 $xi -= ($car = int($xi * $RBASE)) * $MBASE; # see USE_MUL
660 push(@$x, $car); $car = 0;
663 $yi = $yi * $dd + $car;
664 $yi -= ($car = int($yi * $RBASE)) * $MBASE; # see USE_MUL
671 @q = (); ($v2,$v1) = @$y[-2,-1];
675 ($u2,$u1,$u0) = @$x[-3..-1];
677 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
679 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
680 --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
683 ($car, $bar) = (0,0);
684 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
686 $prd = $q * $y->[$yi] + $car;
687 $prd -= ($car = int($prd * $RBASE)) * $MBASE; # see USE_MUL
688 $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
690 if ($x->[-1] < $car + $bar)
693 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
696 if ($car = (($x->[$xi] += $y->[$yi] + $car) > $MBASE));
700 pop(@$x); unshift(@q, $q);
708 for $xi (reverse @$x)
710 $prd = $car * $MBASE + $xi;
711 $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
721 if ($LEN_CONVERT != 0)
723 $c->_to_large($x); $c->_to_large($d);
733 if ($LEN_CONVERT != 0)
746 # ref to array, ref to array, modify first array and return remainder if
748 my ($c,$x,$yorg) = @_;
750 if (@$x == 1 && @$yorg == 1)
752 # shortcut, $yorg and $x are two small numbers
755 my $r = [ $x->[0] % $yorg->[0] ];
756 $x->[0] = int($x->[0] / $yorg->[0]);
761 $x->[0] = int($x->[0] / $yorg->[0]);
768 $rem = _mod($c,[ @$x ],$yorg) if wantarray;
770 # shortcut, $y is < $BASE
771 my $j = scalar @$x; my $r = 0;
772 my $y = $yorg->[0]; my $b;
775 $b = $r * $MBASE + $x->[$j];
776 $x->[$j] = int($b/$y);
779 pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
780 return ($x,$rem) if wantarray;
785 if ($LEN_CONVERT != 0)
787 $c->_to_small($x); $c->_to_small($y);
790 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
792 $car = $bar = $prd = 0;
793 if (($dd = int($MBASE/($y->[-1]+1))) != 1)
797 $xi = $xi * $dd + $car;
798 $xi -= ($car = int($xi / $MBASE)) * $MBASE;
800 push(@$x, $car); $car = 0;
803 $yi = $yi * $dd + $car;
804 $yi -= ($car = int($yi / $MBASE)) * $MBASE;
811 @q = (); ($v2,$v1) = @$y[-2,-1];
815 ($u2,$u1,$u0) = @$x[-3..-1];
817 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
819 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
820 --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
823 ($car, $bar) = (0,0);
824 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
826 $prd = $q * $y->[$yi] + $car;
827 $prd -= ($car = int($prd / $MBASE)) * $MBASE;
828 $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
830 if ($x->[-1] < $car + $bar)
833 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
836 if ($car = (($x->[$xi] += $y->[$yi] + $car) > $MBASE));
840 pop(@$x); unshift(@q, $q);
848 for $xi (reverse @$x)
850 $prd = $car * $MBASE + $xi;
851 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
861 if ($LEN_CONVERT != 0)
863 $c->_to_large($x); $c->_to_large($d);
873 if ($LEN_CONVERT != 0)
884 ##############################################################################
889 # internal absolute post-normalized compare (ignore signs)
890 # ref to array, ref to array, return <0, 0, >0
891 # arrays must have at least one entry; this is not checked for
893 my ($c,$cx,$cy) = @_;
895 # fast comp based on array elements
896 my $lxy = scalar @$cx - scalar @$cy;
897 return -1 if $lxy < 0; # already differs, ret
898 return 1 if $lxy > 0; # ditto
900 # now calculate length based on digits, not parts
901 $lxy = _len($c,$cx) - _len($c,$cy); # difference
902 return -1 if $lxy < 0;
903 return 1 if $lxy > 0;
905 # hm, same lengths, but same contents?
907 # first way takes 5.49 sec instead of 4.87, but has the early out advantage
908 # so grep is slightly faster, but more inflexible. hm. $_ instead of $k
909 # yields 5.6 instead of 5.5 sec huh?
910 # manual way (abort if unequal, good for early ne)
911 my $j = scalar @$cx - 1;
914 last if ($a = $cx->[$j] - $cy->[$j]); $j--;
916 # my $j = scalar @$cx;
919 # last if ($a = $cx->[$j] - $cy->[$j]);
925 # while it early aborts, it is even slower than the manual variant
926 #grep { return $a if ($a = $_ - $cy->[$i++]); } @$cx;
927 # grep way, go trough all (bad for early ne)
928 #grep { $a = $_ - $cy->[$i++]; } @$cx;
934 # compute number of digits in bigint, minus the sign
936 # int() because add/sub sometimes leaves strings (like '00005') instead of
937 # '5' in this place, thus causing length() to report wrong length
940 return (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
945 # return the nth digit, negative values count backward
946 # zero is rightmost, so _digit(123,0) will give 3
949 my $len = _len('',$x);
951 $n = $len+$n if $n < 0; # -1 last, -2 second-to-last
952 $n = abs($n); # if negative was too big
953 $len--; $n = $len if $n > $len; # n to big?
955 my $elem = int($n / $BASE_LEN); # which array element
956 my $digit = $n % $BASE_LEN; # which digit in this element
957 $elem = '0000'.@$x[$elem]; # get element padded with 0's
958 return substr($elem,-$digit-1,1);
963 # return amount of trailing zeros in decimal
964 # check each array elem in _m for having 0 at end as long as elem == 0
965 # Upon finding a elem != 0, stop
967 my $zeros = 0; my $elem;
972 $elem = "$e"; # preserve x
973 $elem =~ s/.*?(0*$)/$1/; # strip anything not zero
974 $zeros *= $BASE_LEN; # elems * 5
975 $zeros += length($elem); # count trailing zeros
978 $zeros ++; # real else branch: 50% slower!
983 ##############################################################################
988 # return true if arg (BINT or num_str) is zero (array '+', '0')
991 (((scalar @$x == 1) && ($x->[0] == 0))) <=> 0;
996 # return true if arg (BINT or num_str) is even
998 (!($x->[0] & 1)) <=> 0;
1003 # return true if arg (BINT or num_str) is even
1006 (($x->[0] & 1)) <=> 0;
1011 # return true if arg (BINT or num_str) is one (array '+', '1')
1014 (scalar @$x == 1) && ($x->[0] == 1) <=> 0;
1019 # internal normalization function that strips leading zeros from the array
1020 # args: ref to array
1023 my $cnt = scalar @$s; # get count of parts
1025 push @$s,0 if $i < 0; # div might return empty results, so fix it
1027 return $s if @$s == 1; # early out
1029 #print "strip: cnt $cnt i $i\n";
1030 # '0', '3', '4', '0', '0',
1035 # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
1036 # >= 1: skip first part (this can be zero)
1037 while ($i > 0) { last if $s->[$i] != 0; $i--; }
1038 $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
1042 ###############################################################################
1043 # check routine to test internal state of corruptions
1047 # used by the test suite
1050 return "$x is not a reference" if !ref($x);
1052 # are all parts are valid?
1053 my $i = 0; my $j = scalar @$x; my ($e,$try);
1056 $e = $x->[$i]; $e = 'undef' unless defined $e;
1057 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
1058 last if $e !~ /^[+]?[0-9]+$/;
1059 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
1060 last if "$e" !~ /^[+]?[0-9]+$/;
1061 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
1062 last if '' . "$e" !~ /^[+]?[0-9]+$/;
1063 $try = ' < 0 || >= $BASE; '."($x, $e)";
1064 last if $e <0 || $e >= $BASE;
1065 # this test is disabled, since new/bnorm and certain ops (like early out
1066 # in add/sub) are allowed/expected to leave '00000' in some elements
1067 #$try = '=~ /^00+/; '."($x, $e)";
1068 #last if $e =~ /^00+/;
1071 return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
1076 ###############################################################################
1077 ###############################################################################
1078 # some optional routines to make BigInt faster
1082 # if possible, use mod shortcut
1083 my ($c,$x,$yo) = @_;
1085 # slow way since $y to big
1086 if (scalar @$yo > 1)
1088 my ($xo,$rem) = _div($c,$x,$yo);
1092 # both are single element arrays
1093 if (scalar @$x == 1)
1099 # @y is single element, but @x has more than one
1103 # when BASE % Y == 0 then (B * BASE) % Y == 0
1104 # (B * BASE) % $y + A % Y => A % Y
1105 # so need to consider only last element: O(1)
1110 # else need to go trough all elements: O(N), but loop is a bit simplified
1114 $r = ($r + $_) % $y; # not much faster, but heh...
1115 #$r += $_ % $y; $r %= $y;
1122 # else need to go trough all elements: O(N)
1123 my $r = 0; my $bm = 1;
1126 $r = ($_ * $bm + $r) % $y;
1127 $bm = ($bm * $b) % $y;
1129 #$r += ($_ % $y) * $bm;
1141 ##############################################################################
1146 my ($c,$x,$y,$n) = @_;
1150 $n = _new($c,\$n); return _div($c,$x, _pow($c,$n,$y));
1153 # shortcut (faster) for shifting by 10)
1154 # multiples of $BASE_LEN
1155 my $dst = 0; # destination
1156 my $src = _num($c,$y); # as normal int
1157 my $rem = $src % $BASE_LEN; # remainder to shift
1158 $src = int($src / $BASE_LEN); # source
1161 splice (@$x,0,$src); # even faster, 38.4 => 39.3
1165 my $len = scalar @$x - $src; # elems to go
1166 my $vd; my $z = '0'x $BASE_LEN;
1167 $x->[scalar @$x] = 0; # avoid || 0 test inside loop
1170 $vd = $z.$x->[$src];
1171 $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
1173 $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
1174 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1175 $x->[$dst] = int($vd);
1178 splice (@$x,$dst) if $dst > 0; # kill left-over array elems
1179 pop @$x if $x->[-1] == 0 && @$x > 1; # kill last element if 0
1186 my ($c,$x,$y,$n) = @_;
1190 $n = _new($c,\$n); return _mul($c,$x, _pow($c,$n,$y));
1193 # shortcut (faster) for shifting by 10) since we are in base 10eX
1194 # multiples of $BASE_LEN:
1195 my $src = scalar @$x; # source
1196 my $len = _num($c,$y); # shift-len as normal int
1197 my $rem = $len % $BASE_LEN; # remainder to shift
1198 my $dst = $src + int($len/$BASE_LEN); # destination
1199 my $vd; # further speedup
1200 $x->[$src] = 0; # avoid first ||0 for speed
1201 my $z = '0' x $BASE_LEN;
1204 $vd = $x->[$src]; $vd = $z.$vd;
1205 $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
1206 $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
1207 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1208 $x->[$dst] = int($vd);
1211 # set lowest parts to 0
1212 while ($dst >= 0) { $x->[$dst--] = 0; }
1213 # fix spurios last zero element
1214 splice @$x,-1 if $x->[-1] == 0;
1221 # ref to array, ref to array, return ref to array
1222 my ($c,$cx,$cy) = @_;
1226 my $y1 = _copy($c,$cy);
1227 while (!_is_one($c,$y1))
1229 _mul($c,$pow2,$cx) if _is_odd($c,$y1);
1233 _mul($c,$cx,$pow2) unless _is_one($c,$pow2);
1240 # ref to array, return ref to array
1243 if ((@$cx == 1) && ($cx->[0] <= 2))
1245 $cx->[0] = 1 * ($cx->[0]||1); # 0,1 => 1, 2 => 2
1249 # go forward until $base is exceeded
1250 # limit is either $x or $base (x == 100 means as result too high)
1251 my $steps = 100; $steps = $cx->[0] if @$cx == 1;
1252 my $r = 2; my $cf = 3; my $step = 1; my $last = $r;
1253 while ($r < $BASE && $step < $steps)
1255 $last = $r; $r *= $cf++; $step++;
1257 if ((@$cx == 1) && ($step == $cx->[0]))
1263 my $n = _copy($c,$cx);
1267 while (!(@$n == 1 && $n->[0] == $step))
1269 _mul($c,$cx,$n); _dec($c,$n);
1274 use constant DEBUG => 0;
1278 sub steps { $steps };
1283 # ref to array, return ref to array
1286 if (scalar @$x == 1)
1288 # fit's into one Perl scalar
1289 $x->[0] = int(sqrt($x->[0]));
1292 my $y = _copy($c,$x);
1293 # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
1294 # since our guess will "grow"
1295 my $l = int((_len($c,$x)-1) / 2);
1297 my $lastelem = $x->[-1]; # for guess
1298 my $elems = scalar @$x - 1;
1299 # not enough digits, but could have more?
1300 if ((length($lastelem) <= 3) && ($elems > 1))
1302 # right-align with zero pad
1303 my $len = length($lastelem) & 1;
1304 print "$lastelem => " if DEBUG;
1305 $lastelem .= substr($x->[-2] . '0' x $BASE_LEN,0,$BASE_LEN);
1306 # former odd => make odd again, or former even to even again
1307 $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
1308 print "$lastelem\n" if DEBUG;
1311 # construct $x (instead of _lsft($c,$x,$l,10)
1312 my $r = $l % $BASE_LEN; # 10000 00000 00000 00000 ($BASE_LEN=5)
1313 $l = int($l / $BASE_LEN);
1314 print "l = $l " if DEBUG;
1316 splice @$x,$l; # keep ref($x), but modify it
1318 # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
1320 # 14400 00000 => sqrt(14400) => 120
1321 # 144000 000000 => sqrt(144000) => 379
1323 # $x->[$l--] = int('1' . '0' x $r); # old way of guessing
1324 print "$lastelem (elems $elems) => " if DEBUG;
1325 $lastelem = $lastelem / 10 if ($elems & 1 == 1); # odd or even?
1326 my $g = sqrt($lastelem); $g =~ s/\.//; # 2.345 => 2345
1327 $r -= 1 if $elems & 1 == 0; # 70 => 7
1329 # padd with zeros if result is too short
1330 $x->[$l--] = int(substr($g . '0' x $r,0,$r+1));
1331 print "now ",$x->[-1] if DEBUG;
1332 print " would have been ", int('1' . '0' x $r),"\n" if DEBUG;
1334 # If @$x > 1, we could compute the second elem of the guess, too, to create
1335 # an even better guess. Not implemented yet.
1336 $x->[$l--] = 0 while ($l >= 0); # all other digits of guess are zero
1338 print "start x= ",${_str($c,$x)},"\n" if DEBUG;
1341 my $lastlast = _zero();
1342 $steps = 0 if DEBUG;
1343 while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
1346 $lastlast = _copy($c,$last);
1347 $last = _copy($c,$x);
1348 _add($c,$x, _div($c,_copy($c,$y),$x));
1350 print " x= ",${_str($c,$x)},"\n" if DEBUG;
1352 print "\nsteps in sqrt: $steps, " if DEBUG;
1353 _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0; # overshot?
1354 print " final ",$x->[-1],"\n" if DEBUG;
1358 ##############################################################################
1365 # the shortcut makes equal, large numbers _really_ fast, and makes only a
1366 # very small performance drop for small numbers (e.g. something with less
1367 # than 32 bit) Since we optimize for large numbers, this is enabled.
1368 return $x if _acmp($c,$x,$y) == 0; # shortcut
1370 my $m = _one(); my ($xr,$yr);
1371 my $mask = $AND_MASK;
1374 my $y1 = _copy($c,$y); # make copy
1378 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1380 ($x1, $xr) = _div($c,$x1,$mask);
1381 ($y1, $yr) = _div($c,$y1,$mask);
1383 # make ints() from $xr, $yr
1384 # this is when the AND_BITS are greater tahn $BASE and is slower for
1385 # small (<256 bits) numbers, but faster for large numbers. Disabled
1386 # due to KISS principle
1388 # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1389 # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1390 # _add($c,$x, _mul($c, _new( $c, \($xrr & $yrr) ), $m) );
1392 # 0+ due to '&' doesn't work in strings
1393 _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
1403 return _zero() if _acmp($c,$x,$y) == 0; # shortcut (see -and)
1405 my $m = _one(); my ($xr,$yr);
1406 my $mask = $XOR_MASK;
1409 my $y1 = _copy($c,$y); # make copy
1413 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1415 ($x1, $xr) = _div($c,$x1,$mask);
1416 ($y1, $yr) = _div($c,$y1,$mask);
1417 # make ints() from $xr, $yr (see _and())
1418 #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1419 #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1420 #_add($c,$x, _mul($c, _new( $c, \($xrr ^ $yrr) ), $m) );
1422 # 0+ due to '^' doesn't work in strings
1423 _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
1426 # the loop stops when the shorter of the two numbers is exhausted
1427 # the remainder of the longer one will survive bit-by-bit, so we simple
1428 # multiply-add it in
1429 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1430 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1439 return $x if _acmp($c,$x,$y) == 0; # shortcut (see _and)
1441 my $m = _one(); my ($xr,$yr);
1442 my $mask = $OR_MASK;
1445 my $y1 = _copy($c,$y); # make copy
1449 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1451 ($x1, $xr) = _div($c,$x1,$mask);
1452 ($y1, $yr) = _div($c,$y1,$mask);
1453 # make ints() from $xr, $yr (see _and())
1454 # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1455 # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1456 # _add($c,$x, _mul($c, _new( $c, \($xrr | $yrr) ), $m) );
1458 # 0+ due to '|' doesn't work in strings
1459 _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
1462 # the loop stops when the shorter of the two numbers is exhausted
1463 # the remainder of the longer one will survive bit-by-bit, so we simple
1464 # multiply-add it in
1465 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1466 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1473 # convert a decimal number to hex (ref to array, return ref to string)
1476 my $x1 = _copy($c,$x);
1480 my $x10000 = [ 0x10000 ];
1481 while (! _is_zero($c,$x1))
1483 ($x1, $xr) = _div($c,$x1,$x10000);
1484 $es .= unpack('h4',pack('v',$xr->[0]));
1487 $es =~ s/^[0]+//; # strip leading zeros
1494 # convert a decimal number to bin (ref to array, return ref to string)
1497 my $x1 = _copy($c,$x);
1501 my $x10000 = [ 0x10000 ];
1502 while (! _is_zero($c,$x1))
1504 ($x1, $xr) = _div($c,$x1,$x10000);
1505 $es .= unpack('b16',pack('v',$xr->[0]));
1508 $es =~ s/^[0]+//; # strip leading zeros
1515 # convert a hex number to decimal (ref to string, return ref to array)
1519 my $m = [ 0x10000 ]; # 16 bit at a time
1522 my $len = length($$hs)-2;
1523 $len = int($len/4); # 4-digit parts, w/o '0x'
1524 my $val; my $i = -4;
1527 $val = substr($$hs,$i,4);
1528 $val =~ s/^[+-]?0x// if $len == 0; # for last part only because
1529 $val = hex($val); # hex does not like wrong chars
1531 _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1532 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
1539 # convert a hex number to decimal (ref to string, return ref to array)
1542 # instead of converting 8 bit at a time, it is faster to convert the
1543 # number to hex, and then call _from_hex.
1546 $hs =~ s/^[+-]?0b//; # remove sign and 0b
1547 my $l = length($hs); # bits
1548 $hs = '0' x (8-($l % 8)) . $hs if ($l % 8) != 0; # padd left side w/ 0
1549 my $h = unpack('H*', pack ('B*', $hs)); # repack as hex
1550 return $c->_from_hex(\('0x'.$h));
1553 my $m = [ 0x100 ]; # 8 bit at a time
1556 my $len = length($$bs)-2;
1557 $len = int($len/8); # 4-digit parts, w/o '0x'
1558 my $val; my $i = -8;
1561 $val = substr($$bs,$i,8);
1562 $val =~ s/^[+-]?0b// if $len == 0; # for last part only
1564 $val = ord(pack('B8',substr('00000000'.$val,-8,8)));
1567 _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1568 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
1573 ##############################################################################
1574 ##############################################################################
1581 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
1585 Provides support for big integer calculations. Not intended to be used by other
1586 modules (except Math::BigInt::Cached). Other modules which sport the same
1587 functions can also be used to support Math::Bigint, like Math::BigInt::Pari.
1591 In order to allow for multiple big integer libraries, Math::BigInt was
1592 rewritten to use library modules for core math routines. Any module which
1593 follows the same API as this can be used instead by using the following:
1595 use Math::BigInt lib => 'libname';
1597 'libname' is either the long name ('Math::BigInt::Pari'), or only the short
1598 version like 'Pari'.
1602 The following functions MUST be defined in order to support the use by
1605 _new(string) return ref to new object from ref to decimal string
1606 _zero() return a new object with value 0
1607 _one() return a new object with value 1
1609 _str(obj) return ref to a string representing the object
1610 _num(obj) returns a Perl integer/floating point number
1611 NOTE: because of Perl numeric notation defaults,
1612 the _num'ified obj may lose accuracy due to
1613 machine-dependend floating point size limitations
1615 _add(obj,obj) Simple addition of two objects
1616 _mul(obj,obj) Multiplication of two objects
1617 _div(obj,obj) Division of the 1st object by the 2nd
1618 In list context, returns (result,remainder).
1619 NOTE: this is integer math, so no
1620 fractional part will be returned.
1621 _sub(obj,obj) Simple subtraction of 1 object from another
1622 a third, optional parameter indicates that the params
1623 are swapped. In this case, the first param needs to
1624 be preserved, while you can destroy the second.
1625 sub (x,y,1) => return x - y and keep x intact!
1626 _dec(obj) decrement object by one (input is garant. to be > 0)
1627 _inc(obj) increment object by one
1630 _acmp(obj,obj) <=> operator for objects (return -1, 0 or 1)
1632 _len(obj) returns count of the decimal digits of the object
1633 _digit(obj,n) returns the n'th decimal digit of object
1635 _is_one(obj) return true if argument is +1
1636 _is_zero(obj) return true if argument is 0
1637 _is_even(obj) return true if argument is even (0,2,4,6..)
1638 _is_odd(obj) return true if argument is odd (1,3,5,7..)
1640 _copy return a ref to a true copy of the object
1642 _check(obj) check whether internal representation is still intact
1643 return 0 for ok, otherwise error message as string
1645 The following functions are optional, and can be defined if the underlying lib
1646 has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
1647 slow) fallback routines to emulate these:
1649 _from_hex(str) return ref to new object from ref to hexadecimal string
1650 _from_bin(str) return ref to new object from ref to binary string
1652 _as_hex(str) return ref to scalar string containing the value as
1653 unsigned hex string, with the '0x' prepended.
1654 Leading zeros must be stripped.
1655 _as_bin(str) Like as_hex, only as binary string containing only
1656 zeros and ones. Leading zeros must be stripped and a
1657 '0b' must be prepended.
1659 _rsft(obj,N,B) shift object in base B by N 'digits' right
1660 For unsupported bases B, return undef to signal failure
1661 _lsft(obj,N,B) shift object in base B by N 'digits' left
1662 For unsupported bases B, return undef to signal failure
1664 _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
1665 Note: XOR, AND and OR pad with zeros if size mismatches
1666 _and(obj1,obj2) AND (bit-wise) object 1 with object 2
1667 _or(obj1,obj2) OR (bit-wise) object 1 with object 2
1669 _mod(obj,obj) Return remainder of div of the 1st by the 2nd object
1670 _sqrt(obj) return the square root of object (truncate to int)
1671 _fac(obj) return factorial of object 1 (1*2*3*4..)
1672 _pow(obj,obj) return object 1 to the power of object 2
1673 _gcd(obj,obj) return Greatest Common Divisor of two objects
1675 _zeros(obj) return number of trailing decimal zeros
1677 Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
1680 Testing of input parameter validity is done by the caller, so you need not
1681 worry about underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by
1682 zero or similar cases.
1684 The first parameter can be modified, that includes the possibility that you
1685 return a reference to a completely different object instead. Although keeping
1686 the reference and just changing it's contents is prefered over creating and
1687 returning a different reference.
1689 Return values are always references to objects or strings. Exceptions are
1690 C<_lsft()> and C<_rsft()>, which return undef if they can not shift the
1691 argument. This is used to delegate shifting of bases different than the one
1692 you can support back to Math::BigInt, which will use some generic code to
1693 calculate the result.
1695 =head1 WRAP YOUR OWN
1697 If you want to port your own favourite c-lib for big numbers to the
1698 Math::BigInt interface, you can take any of the already existing modules as
1699 a rough guideline. You should really wrap up the latest BigInt and BigFloat
1700 testsuites with your module, and replace in them any of the following:
1706 use Math::BigInt lib => 'yourlib';
1708 This way you ensure that your library really works 100% within Math::BigInt.
1712 This program is free software; you may redistribute it and/or modify it under
1713 the same terms as Perl itself.
1717 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
1719 Seperated from BigInt and shaped API with the help of John Peacock.
1723 L<Math::BigInt>, L<Math::BigFloat>, L<Math::BigInt::BitVect>,
1724 L<Math::BigInt::GMP>, L<Math::BigInt::Cached> and L<Math::BigInt::Pari>.