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 100000
19 # - fully remove funky $# stuff (maybe)
21 # USE_MUL: due to problems on certain os (os390, posix-bc) "* 1e-5" is used
22 # instead of "/ 1e5" at some places, (marked with USE_MUL). Other platforms
23 # BS2000, some Crays need USE_DIV instead.
24 # The BEGIN block is used to determine which of the two variants gives the
27 ##############################################################################
28 # global constants, flags and accessory
30 # constants for easier life
33 my ($BASE,$RBASE,$BASE_LEN);
37 # from Daniel Pfeiffer: determine largest group of digits that is precisely
38 # multipliable with itself plus carry
42 $num = ('9' x ++$e) + 0;
44 } until ($num == $num - 1 or $num - 1 == $num - 2);
46 $BASE = int("1e".$BASE_LEN);
47 $RBASE = abs('1e-'.$BASE_LEN); # see USE_MUL
48 if (int($BASE * $RBASE) == 0) # should be 1
51 *{_mul} = \&_mul_use_div;
52 *{_div} = \&_div_use_div;
56 *{_mul} = \&_mul_use_mul;
57 *{_div} = \&_div_use_mul;
61 # for quering and setting, to debug/benchmark things
68 $BASE = int("1e".$BASE_LEN);
69 $RBASE = abs('1e-'.$BASE_LEN); # see USE_MUL
74 ##############################################################################
75 # create objects from various representations
79 # (string) return ref to num_array
80 # Convert a number from string format to internal base 100000 format.
81 # Assumes normalized value as input.
83 # print "_new $d $$d\n";
84 my $il = CORE::length($$d)-1;
85 # these leaves '00000' instead of int 0 and will be corrected after any op
86 return [ reverse(unpack("a" . ($il % $BASE_LEN+1)
87 . ("a$BASE_LEN" x ($il / $BASE_LEN)), $$d)) ];
107 # catch and throw away
110 ##############################################################################
111 # convert back to string and number
115 # (ref to BINT) return num_str
116 # Convert number from internal base 100000 format to string format.
117 # internal format is always normalized (no leading zeros, "-0" => "+0")
120 my $l = scalar @$ar; # number of parts
121 return $nan if $l < 1; # should not happen
122 # handle first one different to strip leading zeros from it (there are no
123 # leading zero parts in internal representation)
124 $l --; $ret .= $ar->[$l]; $l--;
125 # Interestingly, the pre-padd method uses more time
126 # the old grep variant takes longer (14 to 10 sec)
127 my $z = '0' x ($BASE_LEN-1);
130 $ret .= substr($z.$ar->[$l],-$BASE_LEN); # fastest way I could think of
138 # Make a number (scalar int/float) from a BigInt object
140 return $x->[0] if scalar @$x == 1; # below $BASE
145 $num += $fac*$_; $fac *= $BASE;
150 ##############################################################################
155 # (ref to int_num_array, ref to int_num_array)
156 # routine to add two base 1eX numbers
157 # stolen from Knuth Vol 2 Algorithm A pg 231
158 # there are separate routines to add and sub as per Knuth pg 233
159 # This routine clobbers up array x, but not y.
163 # for each in Y, add Y to X and carry. If after that, something is left in
164 # X, foreach in X add carry to X and then return X, carry
165 # Trades one "$j++" for having to shift arrays, $j could be made integer
166 # but this would impose a limit to number-length of 2**32.
167 my $i; my $car = 0; my $j = 0;
171 if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
176 $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
183 # (ref to int_num_array, ref to int_num_array)
184 # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
185 # subtract Y from X (X is always greater/equal!) by modifying x in place
186 my ($c,$sx,$sy,$s) = @_;
188 my $car = 0; my $i; my $j = 0;
194 last unless defined $sy->[$j] || $car;
195 #print "x: $i y: $sy->[$j] c: $car\n";
196 $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
197 #print "x: $i y: $sy->[$j-1] c: $car\n";
199 # might leave leading zeros, so fix that
205 #print "case 1 (swap)\n";
208 last unless defined $sy->[$j] || $car;
209 #print "$sy->[$j] $i $car => $sx->[$j]\n";
211 if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
212 #print "$sy->[$j] $i $car => $sy->[$j]\n";
215 # might leave leading zeros, so fix that
223 # (BINT, BINT) return nothing
224 # multiply two numbers in internal representation
225 # modifies first arg, second need not be different from first
226 my ($c,$xv,$yv) = @_;
228 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
229 # since multiplying $x with $x fails, make copy in this case
230 $yv = [@$xv] if "$xv" eq "$yv"; # same references?
238 # $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
240 # $prod - ($car = int($prod * RBASE)) * $BASE; # see USE_MUL
242 # $prod[$cty] += $car if $car; # need really to check for 0?
246 # looping through this if $xi == 0 is silly - so optimize it away!
247 $xi = (shift @prod || 0), next if $xi == 0;
250 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
251 ## this is actually a tad slower
252 ## $prod = $prod[$cty]; $prod += ($car + $xi * $yi); # no ||0 here
254 $prod - ($car = int($prod * $RBASE)) * $BASE; # see USE_MUL
256 $prod[$cty] += $car if $car; # need really to check for 0?
261 # normalize (handled last to save check for $y->is_zero()
267 # ref to array, ref to array, modify first array and return remainder if
269 # no longer handles sign
270 my ($c,$x,$yorg) = @_;
271 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1);
273 my (@d,$tmp,$q,$u2,$u1,$u0);
275 $car = $bar = $prd = 0;
278 if (($dd = int($BASE/($y->[-1]+1))) != 1)
282 $xi = $xi * $dd + $car;
283 $xi -= ($car = int($xi * $RBASE)) * $BASE; # see USE_MUL
285 push(@$x, $car); $car = 0;
288 $yi = $yi * $dd + $car;
289 $yi -= ($car = int($yi * $RBASE)) * $BASE; # see USE_MUL
296 @q = (); ($v2,$v1) = @$y[-2,-1];
300 ($u2,$u1,$u0) = @$x[-3..-1];
302 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
304 $q = (($u0 == $v1) ? 99999 : int(($u0*$BASE+$u1)/$v1));
305 --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
308 ($car, $bar) = (0,0);
309 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
311 $prd = $q * $y->[$yi] + $car;
312 $prd -= ($car = int($prd * $RBASE)) * $BASE; # see USE_MUL
313 $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
315 if ($x->[-1] < $car + $bar)
318 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
321 if ($car = (($x->[$xi] += $y->[$yi] + $car) > $BASE));
325 pop(@$x); unshift(@q, $q);
333 for $xi (reverse @$x)
335 $prd = $car * $BASE + $xi;
336 $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
354 ##############################################################################
359 my ($c,$x,$y,$n) = @_;
363 return; # we cant do this here, due to now _pow, so signal failure
367 # shortcut (faster) for shifting by 10)
368 # multiples of $BASE_LEN
369 my $dst = 0; # destination
370 my $src = _num($c,$y); # as normal int
371 my $rem = $src % $BASE_LEN; # reminder to shift
372 $src = int($src / $BASE_LEN); # source
375 splice (@$x,0,$src); # even faster, 38.4 => 39.3
379 my $len = scalar @$x - $src; # elems to go
380 my $vd; my $z = '0'x $BASE_LEN;
381 $x->[scalar @$x] = 0; # avoid || 0 test inside loop
385 #print "$dst $src '$vd' ";
386 $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
389 $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
392 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
393 $x->[$dst] = int($vd);
396 splice (@$x,$dst) if $dst > 0; # kill left-over array elems
397 pop @$x if $x->[-1] == 0; # kill last element if 0
405 my ($c,$x,$y,$n) = @_;
409 return; # we cant do this here, due to now _pow, so signal failure
413 # shortcut (faster) for shifting by 10) since we are in base 10eX
414 # multiples of $BASE_LEN:
415 my $src = scalar @$x; # source
416 my $len = _num($c,$y); # shift-len as normal int
417 my $rem = $len % $BASE_LEN; # reminder to shift
418 my $dst = $src + int($len/$BASE_LEN); # destination
419 my $vd; # further speedup
420 #print "src $src:",$x->[$src]||0," dst $dst:",$v->[$dst]||0," rem $rem\n";
421 $x->[$src] = 0; # avoid first ||0 for speed
422 my $z = '0' x $BASE_LEN;
425 $vd = $x->[$src]; $vd = $z.$vd;
426 #print "s $src d $dst '$vd' ";
427 $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
429 $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
431 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
433 $x->[$dst] = int($vd);
436 # set lowest parts to 0
437 while ($dst >= 0) { $x->[$dst--] = 0; }
438 # fix spurios last zero element
439 splice @$x,-1 if $x->[-1] == 0;
440 #print "elems: "; my $i = 0;
441 #foreach (reverse @$v) { print "$i $_ "; $i++; } print "\n";
446 ##############################################################################
451 # internal absolute post-normalized compare (ignore signs)
452 # ref to array, ref to array, return <0, 0, >0
453 # arrays must have at least one entry; this is not checked for
455 my ($c,$cx, $cy) = @_;
459 # calculate length based on digits, not parts
460 $x = _len('',$cx); $y = _len('',$cy);
461 # print "length: ",($x-$y),"\n";
462 my $lxy = $x - $y; # if different in length
463 return -1 if $lxy < 0;
464 return 1 if $lxy > 0;
465 #print "full compare\n";
467 # first way takes 5.49 sec instead of 4.87, but has the early out advantage
468 # so grep is slightly faster, but more inflexible. hm. $_ instead of $k
469 # yields 5.6 instead of 5.5 sec huh?
470 # manual way (abort if unequal, good for early ne)
471 my $j = scalar @$cx - 1;
474 # print "$cx->[$j] $cy->[$j] $a",$cx->[$j]-$cy->[$j],"\n";
475 last if ($a = $cx->[$j] - $cy->[$j]); $j--;
480 # while it early aborts, it is even slower than the manual variant
481 #grep { return $a if ($a = $_ - $cy->[$i++]); } @$cx;
482 # grep way, go trough all (bad for early ne)
483 #grep { $a = $_ - $cy->[$i++]; } @$cx;
489 # computer number of digits in bigint, minus the sign
490 # int() because add/sub sometimes leaves strings (like '00005') instead of
491 # int ('5') in this place, causing length to fail
494 return (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
499 # return the nth digit, negative values count backward
500 # zero is rightmost, so _digit(123,0) will give 3
503 my $len = _len('',$x);
505 $n = $len+$n if $n < 0; # -1 last, -2 second-to-last
506 $n = abs($n); # if negative was too big
507 $len--; $n = $len if $n > $len; # n to big?
509 my $elem = int($n / $BASE_LEN); # which array element
510 my $digit = $n % $BASE_LEN; # which digit in this element
511 $elem = '0000'.@$x[$elem]; # get element padded with 0's
512 return substr($elem,-$digit-1,1);
517 # return amount of trailing zeros in decimal
518 # check each array elem in _m for having 0 at end as long as elem == 0
519 # Upon finding a elem != 0, stop
521 my $zeros = 0; my $elem;
526 $elem = "$e"; # preserve x
527 $elem =~ s/.*?(0*$)/$1/; # strip anything not zero
528 $zeros *= $BASE_LEN; # elems * 5
529 $zeros += CORE::length($elem); # count trailing zeros
532 $zeros ++; # real else branch: 50% slower!
537 ##############################################################################
542 # return true if arg (BINT or num_str) is zero (array '+', '0')
544 return (((scalar @$x == 1) && ($x->[0] == 0))) <=> 0;
549 # return true if arg (BINT or num_str) is even
551 return (!($x->[0] & 1)) <=> 0;
556 # return true if arg (BINT or num_str) is even
558 return (($x->[0] & 1)) <=> 0;
563 # return true if arg (BINT or num_str) is one (array '+', '1')
565 return (scalar @$x == 1) && ($x->[0] == 1) <=> 0;
570 # internal normalization function that strips leading zeros from the array
574 my $cnt = scalar @$s; # get count of parts
576 #print "strip: cnt $cnt i $i\n";
577 # '0', '3', '4', '0', '0',
582 # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
583 # >= 1: skip first part (this can be zero)
584 while ($i > 0) { last if $s->[$i] != 0; $i--; }
585 $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
589 ###############################################################################
590 # check routine to test internal state of corruptions
594 # used by the test suite
597 return "$x is not a reference" if !ref($x);
599 # are all parts are valid?
600 my $i = 0; my $j = scalar @$x; my ($e,$try);
603 $e = $x->[$i]; $e = 'undef' unless defined $e;
604 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
605 last if $e !~ /^[+]?[0-9]+$/;
606 $try = ' < 0 || >= $BASE; '."($x, $e)";
607 last if $e <0 || $e >= $BASE;
608 # this test is disabled, since new/bnorm and certain ops (like early out
609 # in add/sub) are allowed/expected to leave '00000' in some elements
610 #$try = '=~ /^00+/; '."($x, $e)";
611 #last if $e =~ /^00+/;
614 return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
623 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
627 Provides support for big integer calculations. Not intended
628 to be used by other modules. Other modules which export the
629 same functions can also be used to support Math::Bigint
633 In order to allow for multiple big integer libraries, Math::BigInt
634 was rewritten to use library modules for core math routines. Any
635 module which follows the same API as this can be used instead by
636 using the following call:
638 use Math::BigInt lib => BigNum;
642 The following functions MUST be defined in order to support
643 the use by Math::BigInt:
645 _new(string) return ref to new object from ref to decimal string
646 _zero() return a new object with value 0
647 _one() return a new object with value 1
649 _str(obj) return ref to a string representing the object
650 _num(obj) returns a Perl integer/floating point number
651 NOTE: because of Perl numeric notation defaults,
652 the _num'ified obj may lose accuracy due to
653 machine-dependend floating point size limitations
655 _add(obj,obj) Simple addition of two objects
656 _mul(obj,obj) Multiplication of two objects
657 _div(obj,obj) Division of the 1st object by the 2nd
658 In list context, returns (result,remainder).
659 NOTE: this is integer math, so no
660 fractional part will be returned.
661 _sub(obj,obj) Simple subtraction of 1 object from another
662 a third, optional parameter indicates that the params
663 are swapped. In this case, the first param needs to
664 be preserved, while you can destroy the second.
665 sub (x,y,1) => return x - y and keep x intact!
667 _acmp(obj,obj) <=> operator for objects (return -1, 0 or 1)
669 _len(obj) returns count of the decimal digits of the object
670 _digit(obj,n) returns the n'th decimal digit of object
672 _is_one(obj) return true if argument is +1
673 _is_zero(obj) return true if argument is 0
674 _is_even(obj) return true if argument is even (0,2,4,6..)
675 _is_odd(obj) return true if argument is odd (1,3,5,7..)
677 _copy return a ref to a true copy of the object
679 _check(obj) check whether internal representation is still intact
680 return 0 for ok, otherwise error message as string
682 The following functions are optional, and can be defined if the underlying lib
683 has a fast way to do them. If not defined, Math::BigInt will use a pure, but
684 slow, Perl way as fallback to emulate these:
686 _from_hex(str) return ref to new object from ref to hexadecimal string
687 _from_bin(str) return ref to new object from ref to binary string
689 _rsft(obj,N,B) shift object in base B by N 'digits' right
690 _lsft(obj,N,B) shift object in base B by N 'digits' left
692 _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
693 Mote: XOR, AND and OR pad with zeros if size mismatches
694 _and(obj1,obj2) AND (bit-wise) object 1 with object 2
695 _or(obj1,obj2) OR (bit-wise) object 1 with object 2
697 _sqrt(obj) return the square root of object
698 _pow(obj,obj) return object 1 to the power of object 2
699 _gcd(obj,obj) return Greatest Common Divisor of two objects
701 _zeros(obj) return number of trailing decimal zeros
703 _dec(obj) decrement object by one (input is >= 1)
704 _inc(obj) increment object by one
706 Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
709 Testing of input parameter validity is done by the caller, so you need not
710 worry about underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by
711 zero or similar cases.
713 The first parameter can be modified, that includes the possibility that you
714 return a reference to a completely different object instead. Although keeping
715 the reference is prefered over creating and returning a different one.
717 Return values are always references to objects or strings. Exceptions are
718 C<_lsft()> and C<_rsft()>, which return undef if they can not shift the
719 argument. This is used to delegate shifting of bases different than 10 back
720 to BigInt, which will use some generic code to calculate the result.
724 If you want to port your own favourite c-lib for big numbers to the
725 Math::BigInt interface, you can take any of the already existing modules as
726 a rough guideline. You should really wrap up the latest BigInt and BigFloat
727 testsuites with your module, and replace in them any of the following:
733 use Math::BigInt lib => 'yourlib';
735 This way you ensure that your library really works 100% within Math::BigInt.
739 This program is free software; you may redistribute it and/or modify it under
740 the same terms as Perl itself.
744 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
746 Seperated from BigInt and shaped API with the help of John Peacock.
750 L<Math::BigInt>, L<Math::BigFloat>, L<Math::BigInt::BitVect> and
751 L<Math::BigInt::Pari>.