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 to internal base 100000 format.
238 # 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 # for each in Y, add Y to X and carry. If after that, something is left in
335 # X, foreach in X add carry to X and then return X, carry
336 # Trades one "$j++" for having to shift arrays, $j could be made integer
337 # but this would impose a limit to number-length of 2**32.
338 my $i; my $car = 0; my $j = 0;
341 $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
346 $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
353 # (ref to int_num_array, ref to int_num_array)
354 # routine to add 1 to a base 1eX numbers
355 # This routine clobbers up array x, but not y.
360 return $x if (($i += 1) < $BASE); # early out
361 $i = 0; # overflow, next
363 push @$x,1 if ($x->[-1] == 0); # last overflowed, so extend
369 # (ref to int_num_array, ref to int_num_array)
370 # routine to add 1 to a base 1eX numbers
371 # This routine clobbers up array x, but not y.
374 my $MAX = $BASE-1; # since MAX_VAL based on MBASE
377 last if (($i -= 1) >= 0); # early out
378 $i = $MAX; # overflow, next
380 pop @$x if $x->[-1] == 0 && @$x > 1; # last overflowed (but leave 0)
386 # (ref to int_num_array, ref to int_num_array)
387 # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
388 # subtract Y from X (X is always greater/equal!) by modifying x in place
389 my ($c,$sx,$sy,$s) = @_;
391 my $car = 0; my $i; my $j = 0;
397 last unless defined $sy->[$j] || $car;
398 $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
400 # might leave leading zeros, so fix that
401 return __strip_zeros($sx);
403 #print "case 1 (swap)\n";
406 last unless defined $sy->[$j] || $car;
408 if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
411 # might leave leading zeros, so fix that
417 # (BINT, BINT) return nothing
418 # multiply two numbers in internal representation
419 # modifies first arg, second need not be different from first
420 my ($c,$xv,$yv) = @_;
422 # shortcut for two very short numbers
423 # +0 since part maybe string '00001' from new()
424 # works also if xv and yv are the same reference
425 if ((@$xv == 1) && (@$yv == 1)
426 && (length($xv->[0]+0) <= $BASE_LEN2)
427 && (length($yv->[0]+0) <= $BASE_LEN2))
429 $xv->[0] *= $yv->[0];
433 # since multiplying $x with $x fails, make copy in this case
434 $yv = [@$xv] if "$xv" eq "$yv"; # same references?
435 if ($LEN_CONVERT != 0)
437 $c->_to_small($xv); $c->_to_small($yv);
440 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
449 # $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
451 # $prod - ($car = int($prod * RBASE)) * $MBASE; # see USE_MUL
453 # $prod[$cty] += $car if $car; # need really to check for 0?
457 # looping through this if $xi == 0 is silly - so optimize it away!
458 $xi = (shift @prod || 0), next if $xi == 0;
461 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
462 ## this is actually a tad slower
463 ## $prod = $prod[$cty]; $prod += ($car + $xi * $yi); # no ||0 here
465 $prod - ($car = int($prod * $RBASE)) * $MBASE; # see USE_MUL
467 $prod[$cty] += $car if $car; # need really to check for 0?
468 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
471 if ($LEN_CONVERT != 0)
485 # (BINT, BINT) return nothing
486 # multiply two numbers in internal representation
487 # modifies first arg, second need not be different from first
488 my ($c,$xv,$yv) = @_;
490 # shortcut for two very short numbers
491 # +0 since part maybe string '00001' from new()
492 # works also if xv and yv are the same reference
493 if ((@$xv == 1) && (@$yv == 1)
494 && (length($xv->[0]+0) <= $BASE_LEN2)
495 && (length($yv->[0]+0) <= $BASE_LEN2))
497 $xv->[0] *= $yv->[0];
501 # since multiplying $x with $x fails, make copy in this case
502 $yv = [@$xv] if "$xv" eq "$yv"; # same references?
503 if ($LEN_CONVERT != 0)
505 $c->_to_small($xv); $c->_to_small($yv);
508 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
512 # looping through this if $xi == 0 is silly - so optimize it away!
513 $xi = (shift @prod || 0), next if $xi == 0;
516 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
518 $prod - ($car = int($prod / $MBASE)) * $MBASE;
520 $prod[$cty] += $car if $car; # need really to check for 0?
521 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
524 if ($LEN_CONVERT != 0)
538 # ref to array, ref to array, modify first array and return remainder if
540 my ($c,$x,$yorg) = @_;
542 if (@$x == 1 && @$yorg == 1)
544 # shortcut, $y is smaller than $x
547 my $r = [ $x->[0] % $yorg->[0] ];
548 $x->[0] = int($x->[0] / $yorg->[0]);
553 $x->[0] = int($x->[0] / $yorg->[0]);
559 if ($LEN_CONVERT != 0)
561 $c->_to_small($x); $c->_to_small($y);
564 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
566 $car = $bar = $prd = 0;
567 if (($dd = int($MBASE/($y->[-1]+1))) != 1)
571 $xi = $xi * $dd + $car;
572 $xi -= ($car = int($xi * $RBASE)) * $MBASE; # see USE_MUL
574 push(@$x, $car); $car = 0;
577 $yi = $yi * $dd + $car;
578 $yi -= ($car = int($yi * $RBASE)) * $MBASE; # see USE_MUL
585 @q = (); ($v2,$v1) = @$y[-2,-1];
589 ($u2,$u1,$u0) = @$x[-3..-1];
591 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
593 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
594 --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
597 ($car, $bar) = (0,0);
598 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
600 $prd = $q * $y->[$yi] + $car;
601 $prd -= ($car = int($prd * $RBASE)) * $MBASE; # see USE_MUL
602 $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
604 if ($x->[-1] < $car + $bar)
607 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
610 if ($car = (($x->[$xi] += $y->[$yi] + $car) > $MBASE));
614 pop(@$x); unshift(@q, $q);
622 for $xi (reverse @$x)
624 $prd = $car * $MBASE + $xi;
625 $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
635 if ($LEN_CONVERT != 0)
637 $c->_to_large($x); $c->_to_large($d);
647 if ($LEN_CONVERT != 0)
660 # ref to array, ref to array, modify first array and return remainder if
662 my ($c,$x,$yorg) = @_;
664 if (@$x == 1 && @$yorg == 1)
666 # shortcut, $y is smaller than $x
669 my $r = [ $x->[0] % $yorg->[0] ];
670 $x->[0] = int($x->[0] / $yorg->[0]);
675 $x->[0] = int($x->[0] / $yorg->[0]);
681 if ($LEN_CONVERT != 0)
683 $c->_to_small($x); $c->_to_small($y);
686 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
688 $car = $bar = $prd = 0;
689 if (($dd = int($MBASE/($y->[-1]+1))) != 1)
693 $xi = $xi * $dd + $car;
694 $xi -= ($car = int($xi / $MBASE)) * $MBASE;
696 push(@$x, $car); $car = 0;
699 $yi = $yi * $dd + $car;
700 $yi -= ($car = int($yi / $MBASE)) * $MBASE;
707 @q = (); ($v2,$v1) = @$y[-2,-1];
711 ($u2,$u1,$u0) = @$x[-3..-1];
713 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
715 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
716 --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
719 ($car, $bar) = (0,0);
720 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
722 $prd = $q * $y->[$yi] + $car;
723 $prd -= ($car = int($prd / $MBASE)) * $MBASE;
724 $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
726 if ($x->[-1] < $car + $bar)
729 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
732 if ($car = (($x->[$xi] += $y->[$yi] + $car) > $MBASE));
736 pop(@$x); unshift(@q, $q);
744 for $xi (reverse @$x)
746 $prd = $car * $MBASE + $xi;
747 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
757 if ($LEN_CONVERT != 0)
759 $c->_to_large($x); $c->_to_large($d);
769 if ($LEN_CONVERT != 0)
780 ##############################################################################
785 # internal absolute post-normalized compare (ignore signs)
786 # ref to array, ref to array, return <0, 0, >0
787 # arrays must have at least one entry; this is not checked for
789 my ($c,$cx,$cy) = @_;
791 # fast comp based on array elements
792 my $lxy = scalar @$cx - scalar @$cy;
793 return -1 if $lxy < 0; # already differs, ret
794 return 1 if $lxy > 0; # ditto
796 # now calculate length based on digits, not parts
797 $lxy = _len($c,$cx) - _len($c,$cy); # difference
798 return -1 if $lxy < 0;
799 return 1 if $lxy > 0;
801 # hm, same lengths, but same contents?
803 # first way takes 5.49 sec instead of 4.87, but has the early out advantage
804 # so grep is slightly faster, but more inflexible. hm. $_ instead of $k
805 # yields 5.6 instead of 5.5 sec huh?
806 # manual way (abort if unequal, good for early ne)
807 my $j = scalar @$cx - 1;
810 last if ($a = $cx->[$j] - $cy->[$j]); $j--;
816 # while it early aborts, it is even slower than the manual variant
817 #grep { return $a if ($a = $_ - $cy->[$i++]); } @$cx;
818 # grep way, go trough all (bad for early ne)
819 #grep { $a = $_ - $cy->[$i++]; } @$cx;
825 # compute number of digits in bigint, minus the sign
827 # int() because add/sub sometimes leaves strings (like '00005') instead of
828 # '5' in this place, thus causing length() to report wrong length
831 return (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
836 # return the nth digit, negative values count backward
837 # zero is rightmost, so _digit(123,0) will give 3
840 my $len = _len('',$x);
842 $n = $len+$n if $n < 0; # -1 last, -2 second-to-last
843 $n = abs($n); # if negative was too big
844 $len--; $n = $len if $n > $len; # n to big?
846 my $elem = int($n / $BASE_LEN); # which array element
847 my $digit = $n % $BASE_LEN; # which digit in this element
848 $elem = '0000'.@$x[$elem]; # get element padded with 0's
849 return substr($elem,-$digit-1,1);
854 # return amount of trailing zeros in decimal
855 # check each array elem in _m for having 0 at end as long as elem == 0
856 # Upon finding a elem != 0, stop
858 my $zeros = 0; my $elem;
863 $elem = "$e"; # preserve x
864 $elem =~ s/.*?(0*$)/$1/; # strip anything not zero
865 $zeros *= $BASE_LEN; # elems * 5
866 $zeros += length($elem); # count trailing zeros
869 $zeros ++; # real else branch: 50% slower!
874 ##############################################################################
879 # return true if arg (BINT or num_str) is zero (array '+', '0')
882 (((scalar @$x == 1) && ($x->[0] == 0))) <=> 0;
887 # return true if arg (BINT or num_str) is even
889 (!($x->[0] & 1)) <=> 0;
894 # return true if arg (BINT or num_str) is even
897 (($x->[0] & 1)) <=> 0;
902 # return true if arg (BINT or num_str) is one (array '+', '1')
905 (scalar @$x == 1) && ($x->[0] == 1) <=> 0;
910 # internal normalization function that strips leading zeros from the array
914 my $cnt = scalar @$s; # get count of parts
916 push @$s,0 if $i < 0; # div might return empty results, so fix it
918 return $s if @$s == 1; # early out
920 #print "strip: cnt $cnt i $i\n";
921 # '0', '3', '4', '0', '0',
926 # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
927 # >= 1: skip first part (this can be zero)
928 while ($i > 0) { last if $s->[$i] != 0; $i--; }
929 $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
933 ###############################################################################
934 # check routine to test internal state of corruptions
938 # used by the test suite
941 return "$x is not a reference" if !ref($x);
943 # are all parts are valid?
944 my $i = 0; my $j = scalar @$x; my ($e,$try);
947 $e = $x->[$i]; $e = 'undef' unless defined $e;
948 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
949 last if $e !~ /^[+]?[0-9]+$/;
950 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
951 last if "$e" !~ /^[+]?[0-9]+$/;
952 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
953 last if '' . "$e" !~ /^[+]?[0-9]+$/;
954 $try = ' < 0 || >= $BASE; '."($x, $e)";
955 last if $e <0 || $e >= $BASE;
956 # this test is disabled, since new/bnorm and certain ops (like early out
957 # in add/sub) are allowed/expected to leave '00000' in some elements
958 #$try = '=~ /^00+/; '."($x, $e)";
959 #last if $e =~ /^00+/;
962 return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
967 ###############################################################################
968 ###############################################################################
969 # some optional routines to make BigInt faster
973 # if possible, use mod shortcut
976 # slow way since $y to big
979 my ($xo,$rem) = _div($c,$x,$yo);
983 # both are single element arrays
990 # @y is single element, but @x has more than one
994 # when BASE % Y == 0 then (B * BASE) % Y == 0
995 # (B * BASE) % $y + A % Y => A % Y
996 # so need to consider only last element: O(1)
1001 # else need to go trough all elements: O(N), but loop is a bit simplified
1013 # else need to go trough all elements: O(N)
1014 my $r = 0; my $bm = 1;
1017 $r += ($_ % $y) * $bm;
1029 ##############################################################################
1034 my ($c,$x,$y,$n) = @_;
1038 $n = _new($c,\$n); return _div($c,$x, _pow($c,$n,$y));
1041 # shortcut (faster) for shifting by 10)
1042 # multiples of $BASE_LEN
1043 my $dst = 0; # destination
1044 my $src = _num($c,$y); # as normal int
1045 my $rem = $src % $BASE_LEN; # remainder to shift
1046 $src = int($src / $BASE_LEN); # source
1049 splice (@$x,0,$src); # even faster, 38.4 => 39.3
1053 my $len = scalar @$x - $src; # elems to go
1054 my $vd; my $z = '0'x $BASE_LEN;
1055 $x->[scalar @$x] = 0; # avoid || 0 test inside loop
1058 $vd = $z.$x->[$src];
1059 $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
1061 $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
1062 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1063 $x->[$dst] = int($vd);
1066 splice (@$x,$dst) if $dst > 0; # kill left-over array elems
1067 pop @$x if $x->[-1] == 0; # kill last element if 0
1074 my ($c,$x,$y,$n) = @_;
1078 $n = _new($c,\$n); return _mul($c,$x, _pow($c,$n,$y));
1081 # shortcut (faster) for shifting by 10) since we are in base 10eX
1082 # multiples of $BASE_LEN:
1083 my $src = scalar @$x; # source
1084 my $len = _num($c,$y); # shift-len as normal int
1085 my $rem = $len % $BASE_LEN; # remainder to shift
1086 my $dst = $src + int($len/$BASE_LEN); # destination
1087 my $vd; # further speedup
1088 $x->[$src] = 0; # avoid first ||0 for speed
1089 my $z = '0' x $BASE_LEN;
1092 $vd = $x->[$src]; $vd = $z.$vd;
1093 $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
1094 $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
1095 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1096 $x->[$dst] = int($vd);
1099 # set lowest parts to 0
1100 while ($dst >= 0) { $x->[$dst--] = 0; }
1101 # fix spurios last zero element
1102 splice @$x,-1 if $x->[-1] == 0;
1109 # ref to array, ref to array, return ref to array
1110 my ($c,$cx,$cy) = @_;
1114 my $y1 = _copy($c,$cy);
1115 while (!_is_one($c,$y1))
1117 _mul($c,$pow2,$cx) if _is_odd($c,$y1);
1121 _mul($c,$cx,$pow2) unless _is_one($c,$pow2);
1128 # ref to array, return ref to array
1131 if (scalar @$x == 1)
1133 # fit's into one Perl scalar
1134 $x->[0] = int(sqrt($x->[0]));
1137 my $y = _copy($c,$x);
1138 my $l = _len($c,$x) / 2; # hopefully _len/2 is < $BASE
1139 # my $l2 = [ _len($c,$x) / 2 ]; # old way: hopefully _len/2 is < $BASE
1141 splice @$x,0; $x->[0] = 1; # keep ref($x), but modify it
1144 # _lsft($c,$x,$l2,10);
1146 # construct $x (instead of _lsft($c,$x,$l,10)
1147 my $r = $l % $BASE_LEN; # 10000 00000 00000 00000 ($BASE_LEN=5)
1148 $l = int($l / $BASE_LEN);
1149 $x->[$l--] = int('1' . '0' x $r);
1150 $x->[$l--] = 0 while ($l >= 0);
1154 my $lastlast = _zero();
1155 while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
1157 $lastlast = _copy($c,$last);
1158 $last = _copy($c,$x);
1159 _add($c,$x, _div($c,_copy($c,$y),$x));
1162 _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0; # overshot?
1166 ##############################################################################
1173 # the shortcut makes equal, large numbers _really_ fast, and makes only a
1174 # very small performance drop for small numbers (e.g. something with less
1175 # than 32 bit) Since we optimize for large numbers, this is enabled.
1176 return $x if _acmp($c,$x,$y) == 0; # shortcut
1178 my $m = _one(); my ($xr,$yr);
1179 my $mask = $AND_MASK;
1182 my $y1 = _copy($c,$y); # make copy
1186 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1188 ($x1, $xr) = _div($c,$x1,$mask);
1189 ($y1, $yr) = _div($c,$y1,$mask);
1191 # make ints() from $xr, $yr
1192 # this is when the AND_BITS are greater tahn $BASE and is slower for
1193 # small (<256 bits) numbers, but faster for large numbers. Disabled
1194 # due to KISS principle
1196 # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1197 # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1198 # _add($c,$x, _mul($c, _new( $c, \($xrr & $yrr) ), $m) );
1200 # 0+ due to '&' doesn't work in strings
1201 _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
1211 return _zero() if _acmp($c,$x,$y) == 0; # shortcut (see -and)
1213 my $m = _one(); my ($xr,$yr);
1214 my $mask = $XOR_MASK;
1217 my $y1 = _copy($c,$y); # make copy
1221 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1223 ($x1, $xr) = _div($c,$x1,$mask);
1224 ($y1, $yr) = _div($c,$y1,$mask);
1225 # make ints() from $xr, $yr (see _and())
1226 #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1227 #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1228 #_add($c,$x, _mul($c, _new( $c, \($xrr ^ $yrr) ), $m) );
1230 # 0+ due to '^' doesn't work in strings
1231 _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
1234 # the loop stops when the shorter of the two numbers is exhausted
1235 # the remainder of the longer one will survive bit-by-bit, so we simple
1236 # multiply-add it in
1237 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1238 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1247 return $x if _acmp($c,$x,$y) == 0; # shortcut (see _and)
1249 my $m = _one(); my ($xr,$yr);
1250 my $mask = $OR_MASK;
1253 my $y1 = _copy($c,$y); # make copy
1257 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1259 ($x1, $xr) = _div($c,$x1,$mask);
1260 ($y1, $yr) = _div($c,$y1,$mask);
1261 # make ints() from $xr, $yr (see _and())
1262 # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1263 # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1264 # _add($c,$x, _mul($c, _new( $c, \($xrr | $yrr) ), $m) );
1266 # 0+ due to '|' doesn't work in strings
1267 _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
1270 # the loop stops when the shorter of the two numbers is exhausted
1271 # the remainder of the longer one will survive bit-by-bit, so we simple
1272 # multiply-add it in
1273 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1274 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1281 # convert a decimal number to hex (ref to array, return ref to string)
1284 my $x1 = _copy($c,$x);
1288 my $x10000 = [ 0x10000 ];
1289 while (! _is_zero($c,$x1))
1291 ($x1, $xr) = _div($c,$x1,$x10000);
1292 $es .= unpack('h4',pack('v',$xr->[0]));
1295 $es =~ s/^[0]+//; # strip leading zeros
1302 # convert a decimal number to bin (ref to array, return ref to string)
1305 my $x1 = _copy($c,$x);
1309 my $x10000 = [ 0x10000 ];
1310 while (! _is_zero($c,$x1))
1312 ($x1, $xr) = _div($c,$x1,$x10000);
1313 $es .= unpack('b16',pack('v',$xr->[0]));
1316 $es =~ s/^[0]+//; # strip leading zeros
1323 # convert a hex number to decimal (ref to string, return ref to array)
1327 my $m = [ 0x10000 ]; # 16 bit at a time
1330 my $len = length($$hs)-2;
1331 $len = int($len/4); # 4-digit parts, w/o '0x'
1332 my $val; my $i = -4;
1335 $val = substr($$hs,$i,4);
1336 $val =~ s/^[+-]?0x// if $len == 0; # for last part only because
1337 $val = hex($val); # hex does not like wrong chars
1339 _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1340 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
1347 # convert a hex number to decimal (ref to string, return ref to array)
1351 my $m = [ 0x100 ]; # 8 bit at a time
1354 my $len = length($$bs)-2;
1355 $len = int($len/8); # 4-digit parts, w/o '0x'
1356 my $val; my $i = -8;
1359 $val = substr($$bs,$i,8);
1360 $val =~ s/^[+-]?0b// if $len == 0; # for last part only
1362 $val = ord(pack('B8',substr('00000000'.$val,-8,8)));
1365 _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1366 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
1371 ##############################################################################
1372 ##############################################################################
1379 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
1383 Provides support for big integer calculations. Not intended to be used by other
1384 modules (except Math::BigInt::Cached). Other modules which sport the same
1385 functions can also be used to support Math::Bigint, like Math::BigInt::Pari.
1389 In order to allow for multiple big integer libraries, Math::BigInt was
1390 rewritten to use library modules for core math routines. Any module which
1391 follows the same API as this can be used instead by using the following:
1393 use Math::BigInt lib => 'libname';
1395 'libname' is either the long name ('Math::BigInt::Pari'), or only the short
1396 version like 'Pari'.
1400 The following functions MUST be defined in order to support the use by
1403 _new(string) return ref to new object from ref to decimal string
1404 _zero() return a new object with value 0
1405 _one() return a new object with value 1
1407 _str(obj) return ref to a string representing the object
1408 _num(obj) returns a Perl integer/floating point number
1409 NOTE: because of Perl numeric notation defaults,
1410 the _num'ified obj may lose accuracy due to
1411 machine-dependend floating point size limitations
1413 _add(obj,obj) Simple addition of two objects
1414 _mul(obj,obj) Multiplication of two objects
1415 _div(obj,obj) Division of the 1st object by the 2nd
1416 In list context, returns (result,remainder).
1417 NOTE: this is integer math, so no
1418 fractional part will be returned.
1419 _sub(obj,obj) Simple subtraction of 1 object from another
1420 a third, optional parameter indicates that the params
1421 are swapped. In this case, the first param needs to
1422 be preserved, while you can destroy the second.
1423 sub (x,y,1) => return x - y and keep x intact!
1424 _dec(obj) decrement object by one (input is garant. to be > 0)
1425 _inc(obj) increment object by one
1428 _acmp(obj,obj) <=> operator for objects (return -1, 0 or 1)
1430 _len(obj) returns count of the decimal digits of the object
1431 _digit(obj,n) returns the n'th decimal digit of object
1433 _is_one(obj) return true if argument is +1
1434 _is_zero(obj) return true if argument is 0
1435 _is_even(obj) return true if argument is even (0,2,4,6..)
1436 _is_odd(obj) return true if argument is odd (1,3,5,7..)
1438 _copy return a ref to a true copy of the object
1440 _check(obj) check whether internal representation is still intact
1441 return 0 for ok, otherwise error message as string
1443 The following functions are optional, and can be defined if the underlying lib
1444 has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
1445 slow) fallback routines to emulate these:
1447 _from_hex(str) return ref to new object from ref to hexadecimal string
1448 _from_bin(str) return ref to new object from ref to binary string
1450 _as_hex(str) return ref to scalar string containing the value as
1451 unsigned hex string, with the '0x' prepended.
1452 Leading zeros must be stripped.
1453 _as_bin(str) Like as_hex, only as binary string containing only
1454 zeros and ones. Leading zeros must be stripped and a
1455 '0b' must be prepended.
1457 _rsft(obj,N,B) shift object in base B by N 'digits' right
1458 For unsupported bases B, return undef to signal failure
1459 _lsft(obj,N,B) shift object in base B by N 'digits' left
1460 For unsupported bases B, return undef to signal failure
1462 _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
1463 Note: XOR, AND and OR pad with zeros if size mismatches
1464 _and(obj1,obj2) AND (bit-wise) object 1 with object 2
1465 _or(obj1,obj2) OR (bit-wise) object 1 with object 2
1467 _mod(obj,obj) Return remainder of div of the 1st by the 2nd object
1468 _sqrt(obj) return the square root of object (truncate to int)
1469 _pow(obj,obj) return object 1 to the power of object 2
1470 _gcd(obj,obj) return Greatest Common Divisor of two objects
1472 _zeros(obj) return number of trailing decimal zeros
1474 Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
1477 Testing of input parameter validity is done by the caller, so you need not
1478 worry about underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by
1479 zero or similar cases.
1481 The first parameter can be modified, that includes the possibility that you
1482 return a reference to a completely different object instead. Although keeping
1483 the reference and just changing it's contents is prefered over creating and
1484 returning a different reference.
1486 Return values are always references to objects or strings. Exceptions are
1487 C<_lsft()> and C<_rsft()>, which return undef if they can not shift the
1488 argument. This is used to delegate shifting of bases different than the one
1489 you can support back to Math::BigInt, which will use some generic code to
1490 calculate the result.
1492 =head1 WRAP YOUR OWN
1494 If you want to port your own favourite c-lib for big numbers to the
1495 Math::BigInt interface, you can take any of the already existing modules as
1496 a rough guideline. You should really wrap up the latest BigInt and BigFloat
1497 testsuites with your module, and replace in them any of the following:
1503 use Math::BigInt lib => 'yourlib';
1505 This way you ensure that your library really works 100% within Math::BigInt.
1509 This program is free software; you may redistribute it and/or modify it under
1510 the same terms as Perl itself.
1514 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
1516 Seperated from BigInt and shaped API with the help of John Peacock.
1520 L<Math::BigInt>, L<Math::BigFloat>, L<Math::BigInt::BitVect>,
1521 L<Math::BigInt::GMP>, L<Math::BigInt::Cached> and L<Math::BigInt::Pari>.