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 ($BASE,$RBASE,$BASE_LEN,$MAX_VAL,$BASE_LEN2);
34 my ($AND_BITS,$XOR_BITS,$OR_BITS);
35 my ($AND_MASK,$XOR_MASK,$OR_MASK);
39 # set/get the BASE_LEN and assorted other, connected values
40 # used only be the testsuite, set is used only by the BEGIN block below
46 $b = 5 if $^O =~ /^uts/; # UTS needs 5, because 6 and 7 break
49 while (--$BASE_LEN > 5)
51 $BASE = int("1e".$BASE_LEN);
52 $RBASE = abs('1e-'.$BASE_LEN); # see USE_MUL
54 $caught += 1 if (int($BASE * $RBASE) != 1); # should be 1
55 $caught += 2 if (int($BASE / $BASE) != 1); # should be 1
56 # print "caught $caught\n";
59 $BASE = int("1e".$BASE_LEN);
60 $RBASE = abs('1e-'.$BASE_LEN); # see USE_MUL
62 $BASE_LEN2 = int($BASE_LEN / 2); # for mul shortcut
63 # print "BASE_LEN: $BASE_LEN MAX_VAL: $MAX_VAL BASE: $BASE RBASE: $RBASE\n";
68 *{_mul} = \&_mul_use_mul;
69 *{_div} = \&_div_use_mul;
71 else # $caught must be 2, since it can't be 1 nor 3
74 *{_mul} = \&_mul_use_div;
75 *{_div} = \&_div_use_div;
80 return ($BASE_LEN, $AND_BITS, $XOR_BITS, $OR_BITS);
87 # from Daniel Pfeiffer: determine largest group of digits that is precisely
88 # multipliable with itself plus carry
89 # Test now changed to expect the proper pattern, not a result off by 1 or 2
90 my ($e, $num) = 3; # lowest value we will use is 3+1-1 = 3
93 $num = ('9' x ++$e) + 0;
96 } while ("$num" =~ /9{$e}0{$e}/); # must be a certain pattern
97 $e--; # last test failed, so retract one step
98 # the limits below brush the problems with the test above under the rug:
99 # the test should be able to find the proper $e automatically
100 $e = 5 if $^O =~ /^uts/; # UTS get's some special treatment
101 $e = 5 if $^O =~ /^unicos/; # unicos is also problematic (6 seems to work
102 # there, but we play safe)
103 $e = 8 if $e > 8; # cap, for VMS, OS/390 and other 64 bit systems
105 __PACKAGE__->_base_len($e); # set and store
107 # find out how many bits _and, _or and _xor can take (old default = 16)
108 # I don't think anybody has yet 128 bit scalars, so let's play safe.
110 local $^W = 0; # don't warn about 'nonportable number'
111 $AND_BITS = 15; $XOR_BITS = 15; $OR_BITS = 15;
113 # find max bits, we will not go higher than numberofbits that fit into $BASE
114 # to make _and etc simpler (and faster for smaller, slower for large numbers)
116 while (2 ** $max < $BASE) { $max++; }
120 $x = oct('0b' . '1' x $AND_BITS); $y = $x & $x;
121 $z = (2 ** $AND_BITS) - 1;
122 } while ($AND_BITS < $max && $x == $z && $y == $x);
123 $AND_BITS --; # retreat one step
126 $x = oct('0b' . '1' x $XOR_BITS); $y = $x ^ 0;
127 $z = (2 ** $XOR_BITS) - 1;
128 } while ($XOR_BITS < $max && $x == $z && $y == $x);
129 $XOR_BITS --; # retreat one step
132 $x = oct('0b' . '1' x $OR_BITS); $y = $x | $x;
133 $z = (2 ** $OR_BITS) - 1;
134 } while ($OR_BITS < $max && $x == $z && $y == $x);
135 $OR_BITS --; # retreat one step
137 # print "AND $AND_BITS XOR $XOR_BITS OR $OR_BITS\n";
140 ##############################################################################
141 # create objects from various representations
145 # (ref to string) return ref to num_array
146 # Convert a number from string format to internal base 100000 format.
147 # Assumes normalized value as input.
149 my $il = CORE::length($$d)-1;
150 # these leaves '00000' instead of int 0 and will be corrected after any op
151 return [ reverse(unpack("a" . ($il % $BASE_LEN+1)
152 . ("a$BASE_LEN" x ($il / $BASE_LEN)), $$d)) ];
157 $AND_MASK = __PACKAGE__->_new( \( 2 ** $AND_BITS ));
158 $XOR_MASK = __PACKAGE__->_new( \( 2 ** $XOR_BITS ));
159 $OR_MASK = __PACKAGE__->_new( \( 2 ** $OR_BITS ));
176 # create a two (for _pow)
185 # catch and throw away
188 ##############################################################################
189 # convert back to string and number
193 # (ref to BINT) return num_str
194 # Convert number from internal base 100000 format to string format.
195 # internal format is always normalized (no leading zeros, "-0" => "+0")
198 my $l = scalar @$ar; # number of parts
199 return $nan if $l < 1; # should not happen
200 # handle first one different to strip leading zeros from it (there are no
201 # leading zero parts in internal representation)
202 $l --; $ret .= $ar->[$l]; $l--;
203 # Interestingly, the pre-padd method uses more time
204 # the old grep variant takes longer (14 to 10 sec)
205 my $z = '0' x ($BASE_LEN-1);
208 $ret .= substr($z.$ar->[$l],-$BASE_LEN); # fastest way I could think of
216 # Make a number (scalar int/float) from a BigInt object
218 return $x->[0] if scalar @$x == 1; # below $BASE
223 $num += $fac*$_; $fac *= $BASE;
228 ##############################################################################
233 # (ref to int_num_array, ref to int_num_array)
234 # routine to add two base 1eX numbers
235 # stolen from Knuth Vol 2 Algorithm A pg 231
236 # there are separate routines to add and sub as per Knuth pg 233
237 # This routine clobbers up array x, but not y.
241 # for each in Y, add Y to X and carry. If after that, something is left in
242 # X, foreach in X add carry to X and then return X, carry
243 # Trades one "$j++" for having to shift arrays, $j could be made integer
244 # but this would impose a limit to number-length of 2**32.
245 my $i; my $car = 0; my $j = 0;
248 $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
253 $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
260 # (ref to int_num_array, ref to int_num_array)
261 # routine to add 1 to a base 1eX numbers
262 # This routine clobbers up array x, but not y.
267 return $x if (($i += 1) < $BASE); # early out
270 if ($x->[-1] == 0) # last overflowed
279 # (ref to int_num_array, ref to int_num_array)
280 # routine to add 1 to a base 1eX numbers
281 # This routine clobbers up array x, but not y.
286 last if (($i -= 1) >= 0); # early out
289 pop @$x if $x->[-1] == 0 && @$x > 1; # last overflowed (but leave 0)
295 # (ref to int_num_array, ref to int_num_array)
296 # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
297 # subtract Y from X (X is always greater/equal!) by modifying x in place
298 my ($c,$sx,$sy,$s) = @_;
300 my $car = 0; my $i; my $j = 0;
306 last unless defined $sy->[$j] || $car;
307 $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
309 # might leave leading zeros, so fix that
310 return __strip_zeros($sx);
312 #print "case 1 (swap)\n";
315 last unless defined $sy->[$j] || $car;
317 if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
320 # might leave leading zeros, so fix that
326 # (BINT, BINT) return nothing
327 # multiply two numbers in internal representation
328 # modifies first arg, second need not be different from first
329 my ($c,$xv,$yv) = @_;
331 # shortcut for two very short numbers
332 # +0 since part maybe string '00001' from new()
333 if ((@$xv == 1) && (@$yv == 1)
334 && (length($xv->[0]+0) <= $BASE_LEN2)
335 && (length($yv->[0]+0) <= $BASE_LEN2))
337 $xv->[0] *= $yv->[0];
341 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
342 # since multiplying $x with $x fails, make copy in this case
343 $yv = [@$xv] if "$xv" eq "$yv"; # same references?
351 # $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
353 # $prod - ($car = int($prod * RBASE)) * $BASE; # see USE_MUL
355 # $prod[$cty] += $car if $car; # need really to check for 0?
359 # looping through this if $xi == 0 is silly - so optimize it away!
360 $xi = (shift @prod || 0), next if $xi == 0;
363 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
364 ## this is actually a tad slower
365 ## $prod = $prod[$cty]; $prod += ($car + $xi * $yi); # no ||0 here
367 $prod - ($car = int($prod * $RBASE)) * $BASE; # see USE_MUL
369 $prod[$cty] += $car if $car; # need really to check for 0?
370 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
378 # (BINT, BINT) return nothing
379 # multiply two numbers in internal representation
380 # modifies first arg, second need not be different from first
381 my ($c,$xv,$yv) = @_;
383 # shortcut for two very short numbers
384 # +0 since part maybe string '00001' from new()
385 if ((@$xv == 1) && (@$yv == 1)
386 && (length($xv->[0]+0) <= $BASE_LEN2)
387 && (length($yv->[0]+0) <= $BASE_LEN2))
389 $xv->[0] *= $yv->[0];
393 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
394 # since multiplying $x with $x fails, make copy in this case
395 $yv = [@$xv] if "$xv" eq "$yv"; # same references?
399 # looping through this if $xi == 0 is silly - so optimize it away!
400 $xi = (shift @prod || 0), next if $xi == 0;
403 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
405 $prod - ($car = int($prod / $BASE)) * $BASE;
407 $prod[$cty] += $car if $car; # need really to check for 0?
408 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
416 # ref to array, ref to array, modify first array and return remainder if
418 my ($c,$x,$yorg) = @_;
419 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1);
421 my (@d,$tmp,$q,$u2,$u1,$u0);
423 $car = $bar = $prd = 0;
426 if (($dd = int($BASE/($y->[-1]+1))) != 1)
430 $xi = $xi * $dd + $car;
431 $xi -= ($car = int($xi * $RBASE)) * $BASE; # see USE_MUL
433 push(@$x, $car); $car = 0;
436 $yi = $yi * $dd + $car;
437 $yi -= ($car = int($yi * $RBASE)) * $BASE; # see USE_MUL
444 @q = (); ($v2,$v1) = @$y[-2,-1];
448 ($u2,$u1,$u0) = @$x[-3..-1];
450 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
452 # $q = (($u0 == $v1) ? 99999 : int(($u0*$BASE+$u1)/$v1));
453 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
454 --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
457 ($car, $bar) = (0,0);
458 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
460 $prd = $q * $y->[$yi] + $car;
461 $prd -= ($car = int($prd * $RBASE)) * $BASE; # see USE_MUL
462 $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
464 if ($x->[-1] < $car + $bar)
467 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
470 if ($car = (($x->[$xi] += $y->[$yi] + $car) > $BASE));
474 pop(@$x); unshift(@q, $q);
482 for $xi (reverse @$x)
484 $prd = $car * $BASE + $xi;
485 $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
507 # ref to array, ref to array, modify first array and return remainder if
509 my ($c,$x,$yorg) = @_;
510 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1);
512 my (@d,$tmp,$q,$u2,$u1,$u0);
514 $car = $bar = $prd = 0;
517 if (($dd = int($BASE/($y->[-1]+1))) != 1)
521 $xi = $xi * $dd + $car;
522 $xi -= ($car = int($xi / $BASE)) * $BASE;
524 push(@$x, $car); $car = 0;
527 $yi = $yi * $dd + $car;
528 $yi -= ($car = int($yi / $BASE)) * $BASE;
535 @q = (); ($v2,$v1) = @$y[-2,-1];
539 ($u2,$u1,$u0) = @$x[-3..-1];
541 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
543 # $q = (($u0 == $v1) ? 99999 : int(($u0*$BASE+$u1)/$v1));
544 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
545 --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
548 ($car, $bar) = (0,0);
549 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
551 $prd = $q * $y->[$yi] + $car;
552 $prd -= ($car = int($prd / $BASE)) * $BASE;
553 $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
555 if ($x->[-1] < $car + $bar)
558 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
561 if ($car = (($x->[$xi] += $y->[$yi] + $car) > $BASE));
565 pop(@$x); unshift(@q, $q);
573 for $xi (reverse @$x)
575 $prd = $car * $BASE + $xi;
576 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
593 ##############################################################################
598 # internal absolute post-normalized compare (ignore signs)
599 # ref to array, ref to array, return <0, 0, >0
600 # arrays must have at least one entry; this is not checked for
602 my ($c,$cx,$cy) = @_;
604 # fat comp based on array elements
605 my $lxy = scalar @$cx - scalar @$cy;
606 return -1 if $lxy < 0; # already differs, ret
607 return 1 if $lxy > 0; # ditto
609 # now calculate length based on digits, not parts
610 $lxy = _len($c,$cx) - _len($c,$cy); # difference
611 return -1 if $lxy < 0;
612 return 1 if $lxy > 0;
614 # hm, same lengths, but same contents?
616 # first way takes 5.49 sec instead of 4.87, but has the early out advantage
617 # so grep is slightly faster, but more inflexible. hm. $_ instead of $k
618 # yields 5.6 instead of 5.5 sec huh?
619 # manual way (abort if unequal, good for early ne)
620 my $j = scalar @$cx - 1;
623 last if ($a = $cx->[$j] - $cy->[$j]); $j--;
628 # while it early aborts, it is even slower than the manual variant
629 #grep { return $a if ($a = $_ - $cy->[$i++]); } @$cx;
630 # grep way, go trough all (bad for early ne)
631 #grep { $a = $_ - $cy->[$i++]; } @$cx;
637 # compute number of digits in bigint, minus the sign
639 # int() because add/sub sometimes leaves strings (like '00005') instead of
640 # '5' in this place, thus causing length() to report wrong length
643 return (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
648 # return the nth digit, negative values count backward
649 # zero is rightmost, so _digit(123,0) will give 3
652 my $len = _len('',$x);
654 $n = $len+$n if $n < 0; # -1 last, -2 second-to-last
655 $n = abs($n); # if negative was too big
656 $len--; $n = $len if $n > $len; # n to big?
658 my $elem = int($n / $BASE_LEN); # which array element
659 my $digit = $n % $BASE_LEN; # which digit in this element
660 $elem = '0000'.@$x[$elem]; # get element padded with 0's
661 return substr($elem,-$digit-1,1);
666 # return amount of trailing zeros in decimal
667 # check each array elem in _m for having 0 at end as long as elem == 0
668 # Upon finding a elem != 0, stop
670 my $zeros = 0; my $elem;
675 $elem = "$e"; # preserve x
676 $elem =~ s/.*?(0*$)/$1/; # strip anything not zero
677 $zeros *= $BASE_LEN; # elems * 5
678 $zeros += CORE::length($elem); # count trailing zeros
681 $zeros ++; # real else branch: 50% slower!
686 ##############################################################################
691 # return true if arg (BINT or num_str) is zero (array '+', '0')
693 return (((scalar @$x == 1) && ($x->[0] == 0))) <=> 0;
698 # return true if arg (BINT or num_str) is even
700 return (!($x->[0] & 1)) <=> 0;
705 # return true if arg (BINT or num_str) is even
707 return (($x->[0] & 1)) <=> 0;
712 # return true if arg (BINT or num_str) is one (array '+', '1')
714 return (scalar @$x == 1) && ($x->[0] == 1) <=> 0;
719 # internal normalization function that strips leading zeros from the array
723 my $cnt = scalar @$s; # get count of parts
725 push @$s,0 if $i < 0; # div might return empty results, so fix it
727 #print "strip: cnt $cnt i $i\n";
728 # '0', '3', '4', '0', '0',
733 # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
734 # >= 1: skip first part (this can be zero)
735 while ($i > 0) { last if $s->[$i] != 0; $i--; }
736 $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
740 ###############################################################################
741 # check routine to test internal state of corruptions
745 # used by the test suite
748 return "$x is not a reference" if !ref($x);
750 # are all parts are valid?
751 my $i = 0; my $j = scalar @$x; my ($e,$try);
754 $e = $x->[$i]; $e = 'undef' unless defined $e;
755 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
756 last if $e !~ /^[+]?[0-9]+$/;
757 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
758 last if "$e" !~ /^[+]?[0-9]+$/;
759 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
760 last if '' . "$e" !~ /^[+]?[0-9]+$/;
761 $try = ' < 0 || >= $BASE; '."($x, $e)";
762 last if $e <0 || $e >= $BASE;
763 # this test is disabled, since new/bnorm and certain ops (like early out
764 # in add/sub) are allowed/expected to leave '00000' in some elements
765 #$try = '=~ /^00+/; '."($x, $e)";
766 #last if $e =~ /^00+/;
769 return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
774 ###############################################################################
775 ###############################################################################
776 # some optional routines to make BigInt faster
780 # if possible, use mod shortcut
783 # slow way since $y to big
786 my ($xo,$rem) = _div($c,$x,$yo);
790 # both are single element arrays
797 # @y is single element, but @x has more than one
801 # when BASE % Y == 0 then (B * BASE) % Y == 0
802 # (B * BASE) % $y + A % Y => A % Y
803 # so need to consider only last element: O(1)
808 # else need to go trough all elements: O(N), but loop is a bit simplified
820 # else need to go trough all elements: O(N)
821 my $r = 0; my $bm = 1;
824 $r += ($_ % $y) * $bm;
836 ##############################################################################
841 my ($c,$x,$y,$n) = @_;
845 return; # we cant do this here, due to now _pow, so signal failure
849 # shortcut (faster) for shifting by 10)
850 # multiples of $BASE_LEN
851 my $dst = 0; # destination
852 my $src = _num($c,$y); # as normal int
853 my $rem = $src % $BASE_LEN; # remainder to shift
854 $src = int($src / $BASE_LEN); # source
857 splice (@$x,0,$src); # even faster, 38.4 => 39.3
861 my $len = scalar @$x - $src; # elems to go
862 my $vd; my $z = '0'x $BASE_LEN;
863 $x->[scalar @$x] = 0; # avoid || 0 test inside loop
867 $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
869 $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
870 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
871 $x->[$dst] = int($vd);
874 splice (@$x,$dst) if $dst > 0; # kill left-over array elems
875 pop @$x if $x->[-1] == 0; # kill last element if 0
883 my ($c,$x,$y,$n) = @_;
887 return; # we cant do this here, due to now _pow, so signal failure
891 # shortcut (faster) for shifting by 10) since we are in base 10eX
892 # multiples of $BASE_LEN:
893 my $src = scalar @$x; # source
894 my $len = _num($c,$y); # shift-len as normal int
895 my $rem = $len % $BASE_LEN; # remainder to shift
896 my $dst = $src + int($len/$BASE_LEN); # destination
897 my $vd; # further speedup
898 $x->[$src] = 0; # avoid first ||0 for speed
899 my $z = '0' x $BASE_LEN;
902 $vd = $x->[$src]; $vd = $z.$vd;
903 $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
904 $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
905 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
906 $x->[$dst] = int($vd);
909 # set lowest parts to 0
910 while ($dst >= 0) { $x->[$dst--] = 0; }
911 # fix spurios last zero element
912 splice @$x,-1 if $x->[-1] == 0;
920 # ref to array, ref to array, return ref to array
921 my ($c,$cx,$cy) = @_;
925 my $y1 = _copy($c,$cy);
926 while (!_is_one($c,$y1))
928 _mul($c,$pow2,$cx) if _is_odd($c,$y1);
932 _mul($c,$cx,$pow2) unless _is_one($c,$pow2);
939 # ref to array, return ref to array
944 # fit's into one Perl scalar
945 $x->[0] = int(sqrt($x->[0]));
948 my $y = _copy($c,$x);
949 my $l = [ _len($c,$x) / 2 ];
951 splice @$x,0; $x->[0] = 1; # keep ref($x), but modify it
957 my $lastlast = _zero();
958 while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
960 $lastlast = _copy($c,$last);
961 $last = _copy($c,$x);
962 _add($c,$x, _div($c,_copy($c,$y),$x));
965 _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0; # overshot?
969 ##############################################################################
976 # the shortcut makes equal, large numbers _really_ fast, and makes only a
977 # very small performance drop for small numbers (e.g. something with less
978 # than 32 bit) Since we optimize for large numbers, this is enabled.
979 return $x if _acmp($c,$x,$y) == 0; # shortcut
981 my $m = _one(); my ($xr,$yr);
982 my $mask = $AND_MASK;
985 my $y1 = _copy($c,$y); # make copy
989 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
991 ($x1, $xr) = _div($c,$x1,$mask);
992 ($y1, $yr) = _div($c,$y1,$mask);
994 # make ints() from $xr, $yr
995 # this is when the AND_BITS are greater tahn $BASE and is slower for
996 # small (<256 bits) numbers, but faster for large numbers. Disabled
997 # due to KISS principle
999 # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1000 # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1001 # _add($c,$x, _mul($c, _new( $c, \($xrr & $yrr) ), $m) );
1003 _add($c,$x, _mul($c, [ $xr->[0] & $yr->[0] ], $m) );
1013 return _zero() if _acmp($c,$x,$y) == 0; # shortcut (see -and)
1015 my $m = _one(); my ($xr,$yr);
1016 my $mask = $XOR_MASK;
1019 my $y1 = _copy($c,$y); # make copy
1023 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1025 ($x1, $xr) = _div($c,$x1,$mask);
1026 ($y1, $yr) = _div($c,$y1,$mask);
1027 # make ints() from $xr, $yr (see _and())
1028 #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1029 #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1030 #_add($c,$x, _mul($c, _new( $c, \($xrr ^ $yrr) ), $m) );
1032 _add($c,$x, _mul($c, [ $xr->[0] ^ $yr->[0] ], $m) );
1035 # the loop stops when the shorter of the two numbers is exhausted
1036 # the remainder of the longer one will survive bit-by-bit, so we simple
1037 # multiply-add it in
1038 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1039 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1048 return $x if _acmp($c,$x,$y) == 0; # shortcut (see _and)
1050 my $m = _one(); my ($xr,$yr);
1051 my $mask = $OR_MASK;
1054 my $y1 = _copy($c,$y); # make copy
1058 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1060 ($x1, $xr) = _div($c,$x1,$mask);
1061 ($y1, $yr) = _div($c,$y1,$mask);
1062 # make ints() from $xr, $yr (see _and())
1063 # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1064 # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1065 # _add($c,$x, _mul($c, _new( $c, \($xrr | $yrr) ), $m) );
1067 _add($c,$x, _mul($c, [ $xr->[0] | $yr->[0] ], $m) );
1070 # the loop stops when the shorter of the two numbers is exhausted
1071 # the remainder of the longer one will survive bit-by-bit, so we simple
1072 # multiply-add it in
1073 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1074 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1081 # convert a hex number to decimal (ref to string, return ref to array)
1085 my $m = [ 0x10000 ]; # 16 bit at a time
1088 my $len = CORE::length($$hs)-2;
1089 $len = int($len/4); # 4-digit parts, w/o '0x'
1090 my $val; my $i = -4;
1093 $val = substr($$hs,$i,4);
1094 $val =~ s/^[+-]?0x// if $len == 0; # for last part only because
1095 $val = hex($val); # hex does not like wrong chars
1097 _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1098 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
1105 # convert a hex number to decimal (ref to string, return ref to array)
1109 my $m = [ 0x100 ]; # 8 bit at a time
1112 my $len = CORE::length($$bs)-2;
1113 $len = int($len/8); # 4-digit parts, w/o '0x'
1114 my $val; my $i = -8;
1117 $val = substr($$bs,$i,8);
1118 $val =~ s/^[+-]?0b// if $len == 0; # for last part only
1120 #$val = oct('0b'.$val); # does not work on Perl prior to 5.6.0
1121 # $val = ('0' x (8-CORE::length($val))).$val if CORE::length($val) < 8;
1122 $val = ord(pack('B8',substr('00000000'.$val,-8,8)));
1125 _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1126 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
1131 ##############################################################################
1132 ##############################################################################
1139 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
1143 Provides support for big integer calculations. Not intended to be used by other
1144 modules (except Math::BigInt::Cached). Other modules which sport the same
1145 functions can also be used to support Math::Bigint, like Math::BigInt::Pari.
1149 In order to allow for multiple big integer libraries, Math::BigInt was
1150 rewritten to use library modules for core math routines. Any module which
1151 follows the same API as this can be used instead by using the following:
1153 use Math::BigInt lib => 'libname';
1155 'libname' is either the long name ('Math::BigInt::Pari'), or only the short
1156 version like 'Pari'.
1160 The following functions MUST be defined in order to support the use by
1163 _new(string) return ref to new object from ref to decimal string
1164 _zero() return a new object with value 0
1165 _one() return a new object with value 1
1167 _str(obj) return ref to a string representing the object
1168 _num(obj) returns a Perl integer/floating point number
1169 NOTE: because of Perl numeric notation defaults,
1170 the _num'ified obj may lose accuracy due to
1171 machine-dependend floating point size limitations
1173 _add(obj,obj) Simple addition of two objects
1174 _mul(obj,obj) Multiplication of two objects
1175 _div(obj,obj) Division of the 1st object by the 2nd
1176 In list context, returns (result,remainder).
1177 NOTE: this is integer math, so no
1178 fractional part will be returned.
1179 _sub(obj,obj) Simple subtraction of 1 object from another
1180 a third, optional parameter indicates that the params
1181 are swapped. In this case, the first param needs to
1182 be preserved, while you can destroy the second.
1183 sub (x,y,1) => return x - y and keep x intact!
1184 _dec(obj) decrement object by one (input is garant. to be > 0)
1185 _inc(obj) increment object by one
1188 _acmp(obj,obj) <=> operator for objects (return -1, 0 or 1)
1190 _len(obj) returns count of the decimal digits of the object
1191 _digit(obj,n) returns the n'th decimal digit of object
1193 _is_one(obj) return true if argument is +1
1194 _is_zero(obj) return true if argument is 0
1195 _is_even(obj) return true if argument is even (0,2,4,6..)
1196 _is_odd(obj) return true if argument is odd (1,3,5,7..)
1198 _copy return a ref to a true copy of the object
1200 _check(obj) check whether internal representation is still intact
1201 return 0 for ok, otherwise error message as string
1203 The following functions are optional, and can be defined if the underlying lib
1204 has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
1205 slow) fallback routines to emulate these:
1207 _from_hex(str) return ref to new object from ref to hexadecimal string
1208 _from_bin(str) return ref to new object from ref to binary string
1210 _as_hex(str) return ref to scalar string containing the value as
1211 unsigned hex string, with the '0x' prepended.
1212 Leading zeros must be stripped.
1213 _as_bin(str) Like as_hex, only as binary string containing only
1214 zeros and ones. Leading zeros must be stripped and a
1215 '0b' must be prepended.
1217 _rsft(obj,N,B) shift object in base B by N 'digits' right
1218 For unsupported bases B, return undef to signal failure
1219 _lsft(obj,N,B) shift object in base B by N 'digits' left
1220 For unsupported bases B, return undef to signal failure
1222 _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
1223 Note: XOR, AND and OR pad with zeros if size mismatches
1224 _and(obj1,obj2) AND (bit-wise) object 1 with object 2
1225 _or(obj1,obj2) OR (bit-wise) object 1 with object 2
1227 _mod(obj,obj) Return remainder of div of the 1st by the 2nd object
1228 _sqrt(obj) return the square root of object (truncate to int)
1229 _pow(obj,obj) return object 1 to the power of object 2
1230 _gcd(obj,obj) return Greatest Common Divisor of two objects
1232 _zeros(obj) return number of trailing decimal zeros
1234 Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
1237 Testing of input parameter validity is done by the caller, so you need not
1238 worry about underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by
1239 zero or similar cases.
1241 The first parameter can be modified, that includes the possibility that you
1242 return a reference to a completely different object instead. Although keeping
1243 the reference and just changing it's contents is prefered over creating and
1244 returning a different reference.
1246 Return values are always references to objects or strings. Exceptions are
1247 C<_lsft()> and C<_rsft()>, which return undef if they can not shift the
1248 argument. This is used to delegate shifting of bases different than the one
1249 you can support back to Math::BigInt, which will use some generic code to
1250 calculate the result.
1252 =head1 WRAP YOUR OWN
1254 If you want to port your own favourite c-lib for big numbers to the
1255 Math::BigInt interface, you can take any of the already existing modules as
1256 a rough guideline. You should really wrap up the latest BigInt and BigFloat
1257 testsuites with your module, and replace in them any of the following:
1263 use Math::BigInt lib => 'yourlib';
1265 This way you ensure that your library really works 100% within Math::BigInt.
1269 This program is free software; you may redistribute it and/or modify it under
1270 the same terms as Perl itself.
1274 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
1276 Seperated from BigInt and shaped API with the help of John Peacock.
1280 L<Math::BigInt>, L<Math::BigFloat>, L<Math::BigInt::BitVect>,
1281 L<Math::BigInt::GMP>, L<Math::BigInt::Cached> and L<Math::BigInt::Pari>.