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 (X is always greater/equal!) 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 last unless defined $sy->[$j] || $car;
415 if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
418 # might leave leading zeros, so fix that
424 # compute $x ** 2 or $x * $x in-place and return $x
427 # From: Handbook of Applied Cryptography by A. Menezes, P. van Oorschot and
428 # S. Vanstone., Chapter 14
430 #14.16 Algorithm Multiple-precision squaring
431 #INPUT: positive integer x = (xt 1 xt 2 ... x1 x0)b.
432 #OUTPUT: x * x = x ** 2 in radix b representation.
433 #1. For i from 0 to (2t - 1) do: wi <- 0.
434 #2. For i from 0 to (t - 1) do the following:
435 # 2.1 (uv)b w2i + xi * xi, w2i v, c u.
436 # 2.2 For j from (i + 1)to (t - 1) do the following:
437 # (uv)b <- wi+j + 2*xj * xi + c, wi+j <- v, c <- u.
439 #3. Return((w2t-1 w2t-2 ... w1 w0)b).
441 # # Note: That description is crap. Half of the symbols are not explained or
442 # # used with out beeing set.
443 # my $t = scalar @$x; # count
445 # for ($i = 0; $i < $t; $i++)
447 # $x->[$i] = $x->[$i*2] + $x[$i]*$x[$i];
448 # $x->[$i*2] = $x[$i]; $c = $x[$i];
449 # for ($j = $i+1; $j < $t; $j++)
451 # $x->[$i] = $x->[$i+$j] + 2 * $x->[$i] * $x->[$j];
452 # $x->[$i+$j] = $x[$j]; $c = $x[$i];
454 # $x->[$i+$t] = $x[$i];
461 # (ref to int_num_array, ref to int_num_array)
462 # multiply two numbers in internal representation
463 # modifies first arg, second need not be different from first
464 my ($c,$xv,$yv) = @_;
466 # shortcut for two very short numbers (improved by Nathan Zook)
467 # works also if xv and yv are the same reference
468 if ((@$xv == 1) && (@$yv == 1))
470 if (($xv->[0] *= $yv->[0]) >= $MBASE)
472 $xv->[0] = $xv->[0] - ($xv->[1] = int($xv->[0] * $RBASE)) * $MBASE;
476 # shortcut for result == 0
477 if ( ((@$xv == 1) && ($xv->[0] == 0)) ||
478 ((@$yv == 1) && ($yv->[0] == 0)) )
484 # since multiplying $x with $x fails, make copy in this case
485 $yv = [@$xv] if "$xv" eq "$yv"; # same references?
486 # since multiplying $x with $x would fail here, use the faster squaring
487 # return _square($c,$xv) if "$xv" eq "$yv"; # same reference?
489 if ($LEN_CONVERT != 0)
491 $c->_to_small($xv); $c->_to_small($yv);
494 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
503 # $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
505 # $prod - ($car = int($prod * RBASE)) * $MBASE; # see USE_MUL
507 # $prod[$cty] += $car if $car; # need really to check for 0?
511 # looping through this if $xi == 0 is silly - so optimize it away!
512 $xi = (shift @prod || 0), next if $xi == 0;
515 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
516 ## this is actually a tad slower
517 ## $prod = $prod[$cty]; $prod += ($car + $xi * $yi); # no ||0 here
519 $prod - ($car = int($prod * $RBASE)) * $MBASE; # see USE_MUL
521 $prod[$cty] += $car if $car; # need really to check for 0?
522 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
525 if ($LEN_CONVERT != 0)
539 # (ref to int_num_array, ref to int_num_array)
540 # multiply two numbers in internal representation
541 # modifies first arg, second need not be different from first
542 my ($c,$xv,$yv) = @_;
544 # shortcut for two very short numbers (improved by Nathan Zook)
545 # works also if xv and yv are the same reference
546 if ((@$xv == 1) && (@$yv == 1))
548 if (($xv->[0] *= $yv->[0]) >= $MBASE)
551 $xv->[0] - ($xv->[1] = int($xv->[0] / $MBASE)) * $MBASE;
555 # shortcut for result == 0
556 if ( ((@$xv == 1) && ($xv->[0] == 0)) ||
557 ((@$yv == 1) && ($yv->[0] == 0)) )
564 # since multiplying $x with $x fails, make copy in this case
565 $yv = [@$xv] if "$xv" eq "$yv"; # same references?
566 # since multiplying $x with $x would fail here, use the faster squaring
567 # return _square($c,$xv) if "$xv" eq "$yv"; # same reference?
569 if ($LEN_CONVERT != 0)
571 $c->_to_small($xv); $c->_to_small($yv);
574 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
578 # looping through this if $xi == 0 is silly - so optimize it away!
579 $xi = (shift @prod || 0), next if $xi == 0;
582 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
584 $prod - ($car = int($prod / $MBASE)) * $MBASE;
586 $prod[$cty] += $car if $car; # need really to check for 0?
587 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
590 if ($LEN_CONVERT != 0)
604 # ref to array, ref to array, modify first array and return remainder if
606 my ($c,$x,$yorg) = @_;
608 if (@$x == 1 && @$yorg == 1)
610 # shortcut, $yorg and $x are two small numbers
613 my $r = [ $x->[0] % $yorg->[0] ];
614 $x->[0] = int($x->[0] / $yorg->[0]);
619 $x->[0] = int($x->[0] / $yorg->[0]);
626 $rem = _mod($c,[ @$x ],$yorg) if wantarray;
628 # shortcut, $y is < $BASE
629 my $j = scalar @$x; my $r = 0;
630 my $y = $yorg->[0]; my $b;
633 $b = $r * $MBASE + $x->[$j];
634 $x->[$j] = int($b/$y);
637 pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
638 return ($x,$rem) if wantarray;
643 if ($LEN_CONVERT != 0)
645 $c->_to_small($x); $c->_to_small($y);
648 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
650 $car = $bar = $prd = 0;
651 if (($dd = int($MBASE/($y->[-1]+1))) != 1)
655 $xi = $xi * $dd + $car;
656 $xi -= ($car = int($xi * $RBASE)) * $MBASE; # see USE_MUL
658 push(@$x, $car); $car = 0;
661 $yi = $yi * $dd + $car;
662 $yi -= ($car = int($yi * $RBASE)) * $MBASE; # see USE_MUL
669 @q = (); ($v2,$v1) = @$y[-2,-1];
673 ($u2,$u1,$u0) = @$x[-3..-1];
675 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
677 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
678 --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
681 ($car, $bar) = (0,0);
682 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
684 $prd = $q * $y->[$yi] + $car;
685 $prd -= ($car = int($prd * $RBASE)) * $MBASE; # see USE_MUL
686 $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
688 if ($x->[-1] < $car + $bar)
691 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
694 if ($car = (($x->[$xi] += $y->[$yi] + $car) > $MBASE));
698 pop(@$x); unshift(@q, $q);
706 for $xi (reverse @$x)
708 $prd = $car * $MBASE + $xi;
709 $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
719 if ($LEN_CONVERT != 0)
721 $c->_to_large($x); $c->_to_large($d);
731 if ($LEN_CONVERT != 0)
744 # ref to array, ref to array, modify first array and return remainder if
746 my ($c,$x,$yorg) = @_;
748 if (@$x == 1 && @$yorg == 1)
750 # shortcut, $yorg and $x are two small numbers
753 my $r = [ $x->[0] % $yorg->[0] ];
754 $x->[0] = int($x->[0] / $yorg->[0]);
759 $x->[0] = int($x->[0] / $yorg->[0]);
766 $rem = _mod($c,[ @$x ],$yorg) if wantarray;
768 # shortcut, $y is < $BASE
769 my $j = scalar @$x; my $r = 0;
770 my $y = $yorg->[0]; my $b;
773 $b = $r * $MBASE + $x->[$j];
774 $x->[$j] = int($b/$y);
777 pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
778 return ($x,$rem) if wantarray;
783 if ($LEN_CONVERT != 0)
785 $c->_to_small($x); $c->_to_small($y);
788 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
790 $car = $bar = $prd = 0;
791 if (($dd = int($MBASE/($y->[-1]+1))) != 1)
795 $xi = $xi * $dd + $car;
796 $xi -= ($car = int($xi / $MBASE)) * $MBASE;
798 push(@$x, $car); $car = 0;
801 $yi = $yi * $dd + $car;
802 $yi -= ($car = int($yi / $MBASE)) * $MBASE;
809 @q = (); ($v2,$v1) = @$y[-2,-1];
813 ($u2,$u1,$u0) = @$x[-3..-1];
815 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
817 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
818 --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
821 ($car, $bar) = (0,0);
822 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
824 $prd = $q * $y->[$yi] + $car;
825 $prd -= ($car = int($prd / $MBASE)) * $MBASE;
826 $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
828 if ($x->[-1] < $car + $bar)
831 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
834 if ($car = (($x->[$xi] += $y->[$yi] + $car) > $MBASE));
838 pop(@$x); unshift(@q, $q);
846 for $xi (reverse @$x)
848 $prd = $car * $MBASE + $xi;
849 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
859 if ($LEN_CONVERT != 0)
861 $c->_to_large($x); $c->_to_large($d);
871 if ($LEN_CONVERT != 0)
882 ##############################################################################
887 # internal absolute post-normalized compare (ignore signs)
888 # ref to array, ref to array, return <0, 0, >0
889 # arrays must have at least one entry; this is not checked for
891 my ($c,$cx,$cy) = @_;
893 # fast comp based on array elements
894 my $lxy = scalar @$cx - scalar @$cy;
895 return -1 if $lxy < 0; # already differs, ret
896 return 1 if $lxy > 0; # ditto
898 # now calculate length based on digits, not parts
899 $lxy = _len($c,$cx) - _len($c,$cy); # difference
900 return -1 if $lxy < 0;
901 return 1 if $lxy > 0;
903 # hm, same lengths, but same contents?
905 # first way takes 5.49 sec instead of 4.87, but has the early out advantage
906 # so grep is slightly faster, but more inflexible. hm. $_ instead of $k
907 # yields 5.6 instead of 5.5 sec huh?
908 # manual way (abort if unequal, good for early ne)
909 my $j = scalar @$cx - 1;
912 last if ($a = $cx->[$j] - $cy->[$j]); $j--;
914 # my $j = scalar @$cx;
917 # last if ($a = $cx->[$j] - $cy->[$j]);
923 # while it early aborts, it is even slower than the manual variant
924 #grep { return $a if ($a = $_ - $cy->[$i++]); } @$cx;
925 # grep way, go trough all (bad for early ne)
926 #grep { $a = $_ - $cy->[$i++]; } @$cx;
932 # compute number of digits in bigint, minus the sign
934 # int() because add/sub sometimes leaves strings (like '00005') instead of
935 # '5' in this place, thus causing length() to report wrong length
938 return (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
943 # return the nth digit, negative values count backward
944 # zero is rightmost, so _digit(123,0) will give 3
947 my $len = _len('',$x);
949 $n = $len+$n if $n < 0; # -1 last, -2 second-to-last
950 $n = abs($n); # if negative was too big
951 $len--; $n = $len if $n > $len; # n to big?
953 my $elem = int($n / $BASE_LEN); # which array element
954 my $digit = $n % $BASE_LEN; # which digit in this element
955 $elem = '0000'.@$x[$elem]; # get element padded with 0's
956 return substr($elem,-$digit-1,1);
961 # return amount of trailing zeros in decimal
962 # check each array elem in _m for having 0 at end as long as elem == 0
963 # Upon finding a elem != 0, stop
965 my $zeros = 0; my $elem;
970 $elem = "$e"; # preserve x
971 $elem =~ s/.*?(0*$)/$1/; # strip anything not zero
972 $zeros *= $BASE_LEN; # elems * 5
973 $zeros += length($elem); # count trailing zeros
976 $zeros ++; # real else branch: 50% slower!
981 ##############################################################################
986 # return true if arg (BINT or num_str) is zero (array '+', '0')
989 (((scalar @$x == 1) && ($x->[0] == 0))) <=> 0;
994 # return true if arg (BINT or num_str) is even
996 (!($x->[0] & 1)) <=> 0;
1001 # return true if arg (BINT or num_str) is even
1004 (($x->[0] & 1)) <=> 0;
1009 # return true if arg (BINT or num_str) is one (array '+', '1')
1012 (scalar @$x == 1) && ($x->[0] == 1) <=> 0;
1017 # internal normalization function that strips leading zeros from the array
1018 # args: ref to array
1021 my $cnt = scalar @$s; # get count of parts
1023 push @$s,0 if $i < 0; # div might return empty results, so fix it
1025 return $s if @$s == 1; # early out
1027 #print "strip: cnt $cnt i $i\n";
1028 # '0', '3', '4', '0', '0',
1033 # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
1034 # >= 1: skip first part (this can be zero)
1035 while ($i > 0) { last if $s->[$i] != 0; $i--; }
1036 $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
1040 ###############################################################################
1041 # check routine to test internal state of corruptions
1045 # used by the test suite
1048 return "$x is not a reference" if !ref($x);
1050 # are all parts are valid?
1051 my $i = 0; my $j = scalar @$x; my ($e,$try);
1054 $e = $x->[$i]; $e = 'undef' unless defined $e;
1055 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
1056 last if $e !~ /^[+]?[0-9]+$/;
1057 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
1058 last if "$e" !~ /^[+]?[0-9]+$/;
1059 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
1060 last if '' . "$e" !~ /^[+]?[0-9]+$/;
1061 $try = ' < 0 || >= $BASE; '."($x, $e)";
1062 last if $e <0 || $e >= $BASE;
1063 # this test is disabled, since new/bnorm and certain ops (like early out
1064 # in add/sub) are allowed/expected to leave '00000' in some elements
1065 #$try = '=~ /^00+/; '."($x, $e)";
1066 #last if $e =~ /^00+/;
1069 return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
1074 ###############################################################################
1075 ###############################################################################
1076 # some optional routines to make BigInt faster
1080 # if possible, use mod shortcut
1081 my ($c,$x,$yo) = @_;
1083 # slow way since $y to big
1084 if (scalar @$yo > 1)
1086 my ($xo,$rem) = _div($c,$x,$yo);
1090 # both are single element arrays
1091 if (scalar @$x == 1)
1097 # @y is single element, but @x has more than one
1101 # when BASE % Y == 0 then (B * BASE) % Y == 0
1102 # (B * BASE) % $y + A % Y => A % Y
1103 # so need to consider only last element: O(1)
1108 # else need to go trough all elements: O(N), but loop is a bit simplified
1112 $r = ($r + $_) % $y; # not much faster, but heh...
1113 #$r += $_ % $y; $r %= $y;
1120 # else need to go trough all elements: O(N)
1121 my $r = 0; my $bm = 1;
1124 $r = ($_ * $bm + $r) % $y;
1125 $bm = ($bm * $b) % $y;
1127 #$r += ($_ % $y) * $bm;
1139 ##############################################################################
1144 my ($c,$x,$y,$n) = @_;
1148 $n = _new($c,\$n); return _div($c,$x, _pow($c,$n,$y));
1151 # shortcut (faster) for shifting by 10)
1152 # multiples of $BASE_LEN
1153 my $dst = 0; # destination
1154 my $src = _num($c,$y); # as normal int
1155 my $rem = $src % $BASE_LEN; # remainder to shift
1156 $src = int($src / $BASE_LEN); # source
1159 splice (@$x,0,$src); # even faster, 38.4 => 39.3
1163 my $len = scalar @$x - $src; # elems to go
1164 my $vd; my $z = '0'x $BASE_LEN;
1165 $x->[scalar @$x] = 0; # avoid || 0 test inside loop
1168 $vd = $z.$x->[$src];
1169 $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
1171 $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
1172 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1173 $x->[$dst] = int($vd);
1176 splice (@$x,$dst) if $dst > 0; # kill left-over array elems
1177 pop @$x if $x->[-1] == 0; # kill last element if 0
1184 my ($c,$x,$y,$n) = @_;
1188 $n = _new($c,\$n); return _mul($c,$x, _pow($c,$n,$y));
1191 # shortcut (faster) for shifting by 10) since we are in base 10eX
1192 # multiples of $BASE_LEN:
1193 my $src = scalar @$x; # source
1194 my $len = _num($c,$y); # shift-len as normal int
1195 my $rem = $len % $BASE_LEN; # remainder to shift
1196 my $dst = $src + int($len/$BASE_LEN); # destination
1197 my $vd; # further speedup
1198 $x->[$src] = 0; # avoid first ||0 for speed
1199 my $z = '0' x $BASE_LEN;
1202 $vd = $x->[$src]; $vd = $z.$vd;
1203 $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
1204 $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
1205 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1206 $x->[$dst] = int($vd);
1209 # set lowest parts to 0
1210 while ($dst >= 0) { $x->[$dst--] = 0; }
1211 # fix spurios last zero element
1212 splice @$x,-1 if $x->[-1] == 0;
1219 # ref to array, ref to array, return ref to array
1220 my ($c,$cx,$cy) = @_;
1224 my $y1 = _copy($c,$cy);
1225 while (!_is_one($c,$y1))
1227 _mul($c,$pow2,$cx) if _is_odd($c,$y1);
1231 _mul($c,$cx,$pow2) unless _is_one($c,$pow2);
1238 # ref to array, return ref to array
1241 if ((@$cx == 1) && ($cx->[0] <= 2))
1243 $cx->[0] = 1 * ($cx->[0]||1); # 0,1 => 1, 2 => 2
1247 # go forward until $base is exceeded
1248 # limit is either $x or $base (x == 100 means as result too high)
1249 my $steps = 100; $steps = $cx->[0] if @$cx == 1;
1250 my $r = 2; my $cf = 3; my $step = 1; my $last = $r;
1251 while ($r < $BASE && $step < $steps)
1253 $last = $r; $r *= $cf++; $step++;
1255 if ((@$cx == 1) && ($step == $cx->[0]))
1261 my $n = _copy($c,$cx);
1265 while (!(@$n == 1 && $n->[0] == $step))
1267 _mul($c,$cx,$n); _dec($c,$n);
1272 use constant DEBUG => 0;
1276 sub steps { $steps };
1281 # ref to array, return ref to array
1284 if (scalar @$x == 1)
1286 # fit's into one Perl scalar
1287 $x->[0] = int(sqrt($x->[0]));
1290 my $y = _copy($c,$x);
1291 # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
1292 # since our guess will "grow"
1293 my $l = int((_len($c,$x)-1) / 2);
1295 my $lastelem = $x->[-1]; # for guess
1296 my $elems = scalar @$x - 1;
1297 # not enough digits, but could have more?
1298 if ((length($lastelem) <= 3) && ($elems > 1))
1300 # right-align with zero pad
1301 my $len = length($lastelem) & 1;
1302 print "$lastelem => " if DEBUG;
1303 $lastelem .= substr($x->[-2] . '0' x $BASE_LEN,0,$BASE_LEN);
1304 # former odd => make odd again, or former even to even again
1305 $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
1306 print "$lastelem\n" if DEBUG;
1309 # construct $x (instead of _lsft($c,$x,$l,10)
1310 my $r = $l % $BASE_LEN; # 10000 00000 00000 00000 ($BASE_LEN=5)
1311 $l = int($l / $BASE_LEN);
1312 print "l = $l " if DEBUG;
1314 splice @$x,$l; # keep ref($x), but modify it
1316 # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
1318 # 14400 00000 => sqrt(14400) => 120
1319 # 144000 000000 => sqrt(144000) => 379
1321 # $x->[$l--] = int('1' . '0' x $r); # old way of guessing
1322 print "$lastelem (elems $elems) => " if DEBUG;
1323 $lastelem = $lastelem / 10 if ($elems & 1 == 1); # odd or even?
1324 my $g = sqrt($lastelem); $g =~ s/\.//; # 2.345 => 2345
1325 $r -= 1 if $elems & 1 == 0; # 70 => 7
1327 # padd with zeros if result is too short
1328 $x->[$l--] = int(substr($g . '0' x $r,0,$r+1));
1329 print "now ",$x->[-1] if DEBUG;
1330 print " would have been ", int('1' . '0' x $r),"\n" if DEBUG;
1332 # If @$x > 1, we could compute the second elem of the guess, too, to create
1333 # an even better guess. Not implemented yet.
1334 $x->[$l--] = 0 while ($l >= 0); # all other digits of guess are zero
1336 print "start x= ",${_str($c,$x)},"\n" if DEBUG;
1339 my $lastlast = _zero();
1340 $steps = 0 if DEBUG;
1341 while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
1344 $lastlast = _copy($c,$last);
1345 $last = _copy($c,$x);
1346 _add($c,$x, _div($c,_copy($c,$y),$x));
1348 print " x= ",${_str($c,$x)},"\n" if DEBUG;
1350 print "\nsteps in sqrt: $steps, " if DEBUG;
1351 _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0; # overshot?
1352 print " final ",$x->[-1],"\n" if DEBUG;
1356 ##############################################################################
1363 # the shortcut makes equal, large numbers _really_ fast, and makes only a
1364 # very small performance drop for small numbers (e.g. something with less
1365 # than 32 bit) Since we optimize for large numbers, this is enabled.
1366 return $x if _acmp($c,$x,$y) == 0; # shortcut
1368 my $m = _one(); my ($xr,$yr);
1369 my $mask = $AND_MASK;
1372 my $y1 = _copy($c,$y); # make copy
1376 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1378 ($x1, $xr) = _div($c,$x1,$mask);
1379 ($y1, $yr) = _div($c,$y1,$mask);
1381 # make ints() from $xr, $yr
1382 # this is when the AND_BITS are greater tahn $BASE and is slower for
1383 # small (<256 bits) numbers, but faster for large numbers. Disabled
1384 # due to KISS principle
1386 # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1387 # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1388 # _add($c,$x, _mul($c, _new( $c, \($xrr & $yrr) ), $m) );
1390 # 0+ due to '&' doesn't work in strings
1391 _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
1401 return _zero() if _acmp($c,$x,$y) == 0; # shortcut (see -and)
1403 my $m = _one(); my ($xr,$yr);
1404 my $mask = $XOR_MASK;
1407 my $y1 = _copy($c,$y); # make copy
1411 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1413 ($x1, $xr) = _div($c,$x1,$mask);
1414 ($y1, $yr) = _div($c,$y1,$mask);
1415 # make ints() from $xr, $yr (see _and())
1416 #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1417 #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1418 #_add($c,$x, _mul($c, _new( $c, \($xrr ^ $yrr) ), $m) );
1420 # 0+ due to '^' doesn't work in strings
1421 _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
1424 # the loop stops when the shorter of the two numbers is exhausted
1425 # the remainder of the longer one will survive bit-by-bit, so we simple
1426 # multiply-add it in
1427 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1428 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1437 return $x if _acmp($c,$x,$y) == 0; # shortcut (see _and)
1439 my $m = _one(); my ($xr,$yr);
1440 my $mask = $OR_MASK;
1443 my $y1 = _copy($c,$y); # make copy
1447 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1449 ($x1, $xr) = _div($c,$x1,$mask);
1450 ($y1, $yr) = _div($c,$y1,$mask);
1451 # make ints() from $xr, $yr (see _and())
1452 # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1453 # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1454 # _add($c,$x, _mul($c, _new( $c, \($xrr | $yrr) ), $m) );
1456 # 0+ due to '|' doesn't work in strings
1457 _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
1460 # the loop stops when the shorter of the two numbers is exhausted
1461 # the remainder of the longer one will survive bit-by-bit, so we simple
1462 # multiply-add it in
1463 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1464 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1471 # convert a decimal number to hex (ref to array, return ref to string)
1474 my $x1 = _copy($c,$x);
1478 my $x10000 = [ 0x10000 ];
1479 while (! _is_zero($c,$x1))
1481 ($x1, $xr) = _div($c,$x1,$x10000);
1482 $es .= unpack('h4',pack('v',$xr->[0]));
1485 $es =~ s/^[0]+//; # strip leading zeros
1492 # convert a decimal number to bin (ref to array, return ref to string)
1495 my $x1 = _copy($c,$x);
1499 my $x10000 = [ 0x10000 ];
1500 while (! _is_zero($c,$x1))
1502 ($x1, $xr) = _div($c,$x1,$x10000);
1503 $es .= unpack('b16',pack('v',$xr->[0]));
1506 $es =~ s/^[0]+//; # strip leading zeros
1513 # convert a hex number to decimal (ref to string, return ref to array)
1517 my $m = [ 0x10000 ]; # 16 bit at a time
1520 my $len = length($$hs)-2;
1521 $len = int($len/4); # 4-digit parts, w/o '0x'
1522 my $val; my $i = -4;
1525 $val = substr($$hs,$i,4);
1526 $val =~ s/^[+-]?0x// if $len == 0; # for last part only because
1527 $val = hex($val); # hex does not like wrong chars
1529 _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1530 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
1537 # convert a hex number to decimal (ref to string, return ref to array)
1540 # instead of converting 8 bit at a time, it is faster to convert the
1541 # number to hex, and then call _from_hex.
1544 $hs =~ s/^[+-]?0b//; # remove sign and 0b
1545 my $l = length($hs); # bits
1546 $hs = '0' x (8-($l % 8)) . $hs if ($l % 8) != 0; # padd left side w/ 0
1547 my $h = unpack('H*', pack ('B*', $hs)); # repack as hex
1548 return $c->_from_hex(\('0x'.$h));
1551 my $m = [ 0x100 ]; # 8 bit at a time
1554 my $len = length($$bs)-2;
1555 $len = int($len/8); # 4-digit parts, w/o '0x'
1556 my $val; my $i = -8;
1559 $val = substr($$bs,$i,8);
1560 $val =~ s/^[+-]?0b// if $len == 0; # for last part only
1562 $val = ord(pack('B8',substr('00000000'.$val,-8,8)));
1565 _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1566 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
1571 ##############################################################################
1572 ##############################################################################
1579 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
1583 Provides support for big integer calculations. Not intended to be used by other
1584 modules (except Math::BigInt::Cached). Other modules which sport the same
1585 functions can also be used to support Math::Bigint, like Math::BigInt::Pari.
1589 In order to allow for multiple big integer libraries, Math::BigInt was
1590 rewritten to use library modules for core math routines. Any module which
1591 follows the same API as this can be used instead by using the following:
1593 use Math::BigInt lib => 'libname';
1595 'libname' is either the long name ('Math::BigInt::Pari'), or only the short
1596 version like 'Pari'.
1600 The following functions MUST be defined in order to support the use by
1603 _new(string) return ref to new object from ref to decimal string
1604 _zero() return a new object with value 0
1605 _one() return a new object with value 1
1607 _str(obj) return ref to a string representing the object
1608 _num(obj) returns a Perl integer/floating point number
1609 NOTE: because of Perl numeric notation defaults,
1610 the _num'ified obj may lose accuracy due to
1611 machine-dependend floating point size limitations
1613 _add(obj,obj) Simple addition of two objects
1614 _mul(obj,obj) Multiplication of two objects
1615 _div(obj,obj) Division of the 1st object by the 2nd
1616 In list context, returns (result,remainder).
1617 NOTE: this is integer math, so no
1618 fractional part will be returned.
1619 _sub(obj,obj) Simple subtraction of 1 object from another
1620 a third, optional parameter indicates that the params
1621 are swapped. In this case, the first param needs to
1622 be preserved, while you can destroy the second.
1623 sub (x,y,1) => return x - y and keep x intact!
1624 _dec(obj) decrement object by one (input is garant. to be > 0)
1625 _inc(obj) increment object by one
1628 _acmp(obj,obj) <=> operator for objects (return -1, 0 or 1)
1630 _len(obj) returns count of the decimal digits of the object
1631 _digit(obj,n) returns the n'th decimal digit of object
1633 _is_one(obj) return true if argument is +1
1634 _is_zero(obj) return true if argument is 0
1635 _is_even(obj) return true if argument is even (0,2,4,6..)
1636 _is_odd(obj) return true if argument is odd (1,3,5,7..)
1638 _copy return a ref to a true copy of the object
1640 _check(obj) check whether internal representation is still intact
1641 return 0 for ok, otherwise error message as string
1643 The following functions are optional, and can be defined if the underlying lib
1644 has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
1645 slow) fallback routines to emulate these:
1647 _from_hex(str) return ref to new object from ref to hexadecimal string
1648 _from_bin(str) return ref to new object from ref to binary string
1650 _as_hex(str) return ref to scalar string containing the value as
1651 unsigned hex string, with the '0x' prepended.
1652 Leading zeros must be stripped.
1653 _as_bin(str) Like as_hex, only as binary string containing only
1654 zeros and ones. Leading zeros must be stripped and a
1655 '0b' must be prepended.
1657 _rsft(obj,N,B) shift object in base B by N 'digits' right
1658 For unsupported bases B, return undef to signal failure
1659 _lsft(obj,N,B) shift object in base B by N 'digits' left
1660 For unsupported bases B, return undef to signal failure
1662 _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
1663 Note: XOR, AND and OR pad with zeros if size mismatches
1664 _and(obj1,obj2) AND (bit-wise) object 1 with object 2
1665 _or(obj1,obj2) OR (bit-wise) object 1 with object 2
1667 _mod(obj,obj) Return remainder of div of the 1st by the 2nd object
1668 _sqrt(obj) return the square root of object (truncate to int)
1669 _fac(obj) return factorial of object 1 (1*2*3*4..)
1670 _pow(obj,obj) return object 1 to the power of object 2
1671 _gcd(obj,obj) return Greatest Common Divisor of two objects
1673 _zeros(obj) return number of trailing decimal zeros
1675 Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
1678 Testing of input parameter validity is done by the caller, so you need not
1679 worry about underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by
1680 zero or similar cases.
1682 The first parameter can be modified, that includes the possibility that you
1683 return a reference to a completely different object instead. Although keeping
1684 the reference and just changing it's contents is prefered over creating and
1685 returning a different reference.
1687 Return values are always references to objects or strings. Exceptions are
1688 C<_lsft()> and C<_rsft()>, which return undef if they can not shift the
1689 argument. This is used to delegate shifting of bases different than the one
1690 you can support back to Math::BigInt, which will use some generic code to
1691 calculate the result.
1693 =head1 WRAP YOUR OWN
1695 If you want to port your own favourite c-lib for big numbers to the
1696 Math::BigInt interface, you can take any of the already existing modules as
1697 a rough guideline. You should really wrap up the latest BigInt and BigFloat
1698 testsuites with your module, and replace in them any of the following:
1704 use Math::BigInt lib => 'yourlib';
1706 This way you ensure that your library really works 100% within Math::BigInt.
1710 This program is free software; you may redistribute it and/or modify it under
1711 the same terms as Perl itself.
1715 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
1717 Seperated from BigInt and shaped API with the help of John Peacock.
1721 L<Math::BigInt>, L<Math::BigFloat>, L<Math::BigInt::BitVect>,
1722 L<Math::BigInt::GMP>, L<Math::BigInt::Cached> and L<Math::BigInt::Pari>.