1 package Math::BigInt::Calc;
5 # use warnings; # dont use warnings for older Perls
9 use vars qw/ @ISA @EXPORT $VERSION/;
13 _add _mul _div _mod _sub
19 _check _zero _one _copy _zeros
24 # Package to store unsigned big integers in decimal and do math with them
26 # Internally the numbers are stored in an array with at least 1 element, no
27 # leading zero parts (except the first) and in base 100000
30 # - fully remove funky $# stuff (maybe)
31 # - use integer; vs 1e7 as base
33 # USE_MUL: due to problems on certain os (os390, posix-bc) "* 1e-5" is used
34 # instead of "/ 1e5" at some places, (marked with USE_MUL). But instead of
35 # using the reverse only on problematic machines, I used it everytime to avoid
36 # the costly comparisons. This _should_ work everywhere. Thanx Peter Prymmer
38 ##############################################################################
39 # global constants, flags and accessory
41 # constants for easier life
45 my $BASE = int("1e".$BASE_LEN); # var for trying to change it to 1e7
46 my $RBASE = abs('1e-'.$BASE_LEN); # see USE_MUL
50 # Daniel Pfeiffer: determine largest group of digits that is precisely
51 # multipliable with itself plus carry
54 $num = ('9' x ++$e) + 0;
56 } until ($num == $num - 1 or $num - 1 == $num - 2);
58 $BASE = int("1e".$BASE_LEN);
59 $RBASE = abs('1e-'.$BASE_LEN); # see USE_MUL
62 # for quering and setting, to debug/benchmark things
69 $BASE = int("1e".$BASE_LEN);
70 $RBASE = abs('1e-'.$BASE_LEN); # see USE_MUL
75 ##############################################################################
76 # create objects from various representations
80 # (string) return ref to num_array
81 # Convert a number from string format to internal base 100000 format.
82 # Assumes normalized value as input.
84 # print "_new $d $$d\n";
85 my $il = CORE::length($$d)-1;
86 # these leaves '00000' instead of int 0 and will be corrected after any op
87 return [ reverse(unpack("a" . ($il % $BASE_LEN+1)
88 . ("a$BASE_LEN" x ($il / $BASE_LEN)), $$d)) ];
108 ##############################################################################
109 # convert back to string and number
113 # (ref to BINT) return num_str
114 # Convert number from internal base 100000 format to string format.
115 # internal format is always normalized (no leading zeros, "-0" => "+0")
118 my $l = scalar @$ar; # number of parts
119 return $nan if $l < 1; # should not happen
120 # handle first one different to strip leading zeros from it (there are no
121 # leading zero parts in internal representation)
122 $l --; $ret .= $ar->[$l]; $l--;
123 # Interestingly, the pre-padd method uses more time
124 # the old grep variant takes longer (14 to 10 sec)
125 my $z = '0' x ($BASE_LEN-1);
128 $ret .= substr($z.$ar->[$l],-$BASE_LEN); # fastest way I could think of
136 # Make a number (scalar int/float) from a BigInt object
138 return $x->[0] if scalar @$x == 1; # below $BASE
143 $num += $fac*$_; $fac *= $BASE;
148 ##############################################################################
153 # (ref to int_num_array, ref to int_num_array)
154 # routine to add two base 1eX numbers
155 # stolen from Knuth Vol 2 Algorithm A pg 231
156 # there are separate routines to add and sub as per Knuth pg 233
157 # This routine clobbers up array x, but not y.
161 # for each in Y, add Y to X and carry. If after that, something is left in
162 # X, foreach in X add carry to X and then return X, carry
163 # Trades one "$j++" for having to shift arrays, $j could be made integer
164 # but this would impose a limit to number-length of 2**32.
165 my $i; my $car = 0; my $j = 0;
169 if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
174 $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
181 # (ref to int_num_array, ref to int_num_array)
182 # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
183 # subtract Y from X (X is always greater/equal!) by modifying x in place
184 my ($c,$sx,$sy,$s) = @_;
186 my $car = 0; my $i; my $j = 0;
192 last unless defined $sy->[$j] || $car;
193 #print "x: $i y: $sy->[$j] c: $car\n";
194 $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
195 #print "x: $i y: $sy->[$j-1] c: $car\n";
197 # might leave leading zeros, so fix that
203 #print "case 1 (swap)\n";
206 last unless defined $sy->[$j] || $car;
207 #print "$sy->[$j] $i $car => $sx->[$j]\n";
209 if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
210 #print "$sy->[$j] $i $car => $sy->[$j]\n";
213 # might leave leading zeros, so fix that
221 # (BINT, BINT) return nothing
222 # multiply two numbers in internal representation
223 # modifies first arg, second need not be different from first
224 my ($c,$xv,$yv) = @_;
226 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
227 # since multiplying $x with $x fails, make copy in this case
228 $yv = [@$xv] if "$xv" eq "$yv"; # same references?
236 # $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
238 # $prod - ($car = int($prod * RBASE)) * $BASE; # see USE_MUL
240 # $prod[$cty] += $car if $car; # need really to check for 0?
244 # looping through this if $xi == 0 is silly - so optimize it away!
245 $xi = (shift @prod || 0), next if $xi == 0;
248 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
249 ## this is actually a tad slower
250 ## $prod = $prod[$cty]; $prod += ($car + $xi * $yi); # no ||0 here
252 $prod - ($car = int($prod * $RBASE)) * $BASE; # see USE_MUL
254 $prod[$cty] += $car if $car; # need really to check for 0?
259 # normalize (handled last to save check for $y->is_zero()
265 # ref to array, ref to array, modify first array and return remainder if
267 # no longer handles sign
268 my ($c,$x,$yorg) = @_;
269 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1);
271 my (@d,$tmp,$q,$u2,$u1,$u0);
273 $car = $bar = $prd = 0;
276 if (($dd = int($BASE/($y->[-1]+1))) != 1)
280 $xi = $xi * $dd + $car;
281 $xi -= ($car = int($xi * $RBASE)) * $BASE; # see USE_MUL
283 push(@$x, $car); $car = 0;
286 $yi = $yi * $dd + $car;
287 $yi -= ($car = int($yi * $RBASE)) * $BASE; # see USE_MUL
294 @q = (); ($v2,$v1) = @$y[-2,-1];
298 ($u2,$u1,$u0) = @$x[-3..-1];
300 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
302 $q = (($u0 == $v1) ? 99999 : int(($u0*$BASE+$u1)/$v1));
303 --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
306 ($car, $bar) = (0,0);
307 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
309 $prd = $q * $y->[$yi] + $car;
310 $prd -= ($car = int($prd * $RBASE)) * $BASE; # see USE_MUL
311 $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
313 if ($x->[-1] < $car + $bar)
316 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
319 if ($car = (($x->[$xi] += $y->[$yi] + $car) > $BASE));
323 pop(@$x); unshift(@q, $q);
331 for $xi (reverse @$x)
333 $prd = $car * $BASE + $xi;
334 $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
352 ##############################################################################
357 my ($c,$x,$y,$n) = @_;
361 return; # we cant do this here, due to now _pow, so signal failure
365 # shortcut (faster) for shifting by 10)
366 # multiples of $BASE_LEN
367 my $dst = 0; # destination
368 my $src = _num($c,$y); # as normal int
369 my $rem = $src % $BASE_LEN; # reminder to shift
370 $src = int($src / $BASE_LEN); # source
373 splice (@$x,0,$src); # even faster, 38.4 => 39.3
377 my $len = scalar @$x - $src; # elems to go
378 my $vd; my $z = '0'x $BASE_LEN;
379 $x->[scalar @$x] = 0; # avoid || 0 test inside loop
383 #print "$dst $src '$vd' ";
384 $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
387 $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
390 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
391 $x->[$dst] = int($vd);
394 splice (@$x,$dst) if $dst > 0; # kill left-over array elems
395 pop @$x if $x->[-1] == 0; # kill last element if 0
403 my ($c,$x,$y,$n) = @_;
407 return; # we cant do this here, due to now _pow, so signal failure
411 # shortcut (faster) for shifting by 10) since we are in base 10eX
412 # multiples of $BASE_LEN:
413 my $src = scalar @$x; # source
414 my $len = _num($c,$y); # shift-len as normal int
415 my $rem = $len % $BASE_LEN; # reminder to shift
416 my $dst = $src + int($len/$BASE_LEN); # destination
417 my $vd; # further speedup
418 #print "src $src:",$x->[$src]||0," dst $dst:",$v->[$dst]||0," rem $rem\n";
419 $x->[$src] = 0; # avoid first ||0 for speed
420 my $z = '0' x $BASE_LEN;
423 $vd = $x->[$src]; $vd = $z.$vd;
424 #print "s $src d $dst '$vd' ";
425 $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
427 $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
429 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
431 $x->[$dst] = int($vd);
434 # set lowest parts to 0
435 while ($dst >= 0) { $x->[$dst--] = 0; }
436 # fix spurios last zero element
437 splice @$x,-1 if $x->[-1] == 0;
438 #print "elems: "; my $i = 0;
439 #foreach (reverse @$v) { print "$i $_ "; $i++; } print "\n";
444 ##############################################################################
449 # internal absolute post-normalized compare (ignore signs)
450 # ref to array, ref to array, return <0, 0, >0
451 # arrays must have at least one entry; this is not checked for
453 my ($c,$cx, $cy) = @_;
457 # calculate length based on digits, not parts
458 $x = _len('',$cx); $y = _len('',$cy);
459 # print "length: ",($x-$y),"\n";
460 my $lxy = $x - $y; # if different in length
461 return -1 if $lxy < 0;
462 return 1 if $lxy > 0;
463 #print "full compare\n";
465 # first way takes 5.49 sec instead of 4.87, but has the early out advantage
466 # so grep is slightly faster, but more inflexible. hm. $_ instead of $k
467 # yields 5.6 instead of 5.5 sec huh?
468 # manual way (abort if unequal, good for early ne)
469 my $j = scalar @$cx - 1;
472 # print "$cx->[$j] $cy->[$j] $a",$cx->[$j]-$cy->[$j],"\n";
473 last if ($a = $cx->[$j] - $cy->[$j]); $j--;
478 # while it early aborts, it is even slower than the manual variant
479 #grep { return $a if ($a = $_ - $cy->[$i++]); } @$cx;
480 # grep way, go trough all (bad for early ne)
481 #grep { $a = $_ - $cy->[$i++]; } @$cx;
487 # computer number of digits in bigint, minus the sign
488 # int() because add/sub sometimes leaves strings (like '00005') instead of
489 # int ('5') in this place, causing length to fail
492 return (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
497 # return the nth digit, negative values count backward
498 # zero is rightmost, so _digit(123,0) will give 3
501 my $len = _len('',$x);
503 $n = $len+$n if $n < 0; # -1 last, -2 second-to-last
504 $n = abs($n); # if negative was too big
505 $len--; $n = $len if $n > $len; # n to big?
507 my $elem = int($n / $BASE_LEN); # which array element
508 my $digit = $n % $BASE_LEN; # which digit in this element
509 $elem = '0000'.@$x[$elem]; # get element padded with 0's
510 return substr($elem,-$digit-1,1);
515 # return amount of trailing zeros in decimal
516 # check each array elem in _m for having 0 at end as long as elem == 0
517 # Upon finding a elem != 0, stop
519 my $zeros = 0; my $elem;
524 $elem = "$e"; # preserve x
525 $elem =~ s/.*?(0*$)/$1/; # strip anything not zero
526 $zeros *= $BASE_LEN; # elems * 5
527 $zeros += CORE::length($elem); # count trailing zeros
530 $zeros ++; # real else branch: 50% slower!
535 ##############################################################################
540 # return true if arg (BINT or num_str) is zero (array '+', '0')
542 return (((scalar @$x == 1) && ($x->[0] == 0))) <=> 0;
547 # return true if arg (BINT or num_str) is even
549 return (!($x->[0] & 1)) <=> 0;
554 # return true if arg (BINT or num_str) is even
556 return (($x->[0] & 1)) <=> 0;
561 # return true if arg (BINT or num_str) is one (array '+', '1')
563 return (scalar @$x == 1) && ($x->[0] == 1) <=> 0;
568 # internal normalization function that strips leading zeros from the array
572 my $cnt = scalar @$s; # get count of parts
574 #print "strip: cnt $cnt i $i\n";
575 # '0', '3', '4', '0', '0',
580 # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
581 # >= 1: skip first part (this can be zero)
582 while ($i > 0) { last if $s->[$i] != 0; $i--; }
583 $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
587 ###############################################################################
588 # check routine to test internal state of corruptions
592 # no checks yet, pull it out from the test suite
595 return "$x is not a reference" if !ref($x);
597 # are all parts are valid?
598 my $i = 0; my $j = scalar @$x; my ($e,$try);
601 $e = $x->[$i]; $e = 'undef' unless defined $e;
602 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
603 last if $e !~ /^[+]?[0-9]+$/;
604 $try = ' < 0 || >= $BASE; '."($x, $e)";
605 last if $e <0 || $e >= $BASE;
606 # this test is disabled, since new/bnorm and certain ops (like early out
607 # in add/sub) are allowed/expected to leave '00000' in some elements
608 #$try = '=~ /^00+/; '."($x, $e)";
609 #last if $e =~ /^00+/;
612 return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
621 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
625 Provides support for big integer calculations. Not intended
626 to be used by other modules. Other modules which export the
627 same functions can also be used to support Math::Bigint
631 In order to allow for multiple big integer libraries, Math::BigInt
632 was rewritten to use library modules for core math routines. Any
633 module which follows the same API as this can be used instead by
634 using the following call:
636 use Math::BigInt lib => BigNum;
640 The following functions MUST be exported in order to support
641 the use by Math::BigInt:
643 _new(string) return ref to new object from ref to decimal string
644 _zero() return a new object with value 0
645 _one() return a new object with value 1
647 _str(obj) return ref to a string representing the object
648 _num(obj) returns a Perl integer/floating point number
649 NOTE: because of Perl numeric notation defaults,
650 the _num'ified obj may lose accuracy due to
651 machine-dependend floating point size limitations
653 _add(obj,obj) Simple addition of two objects
654 _mul(obj,obj) Multiplication of two objects
655 _div(obj,obj) Division of the 1st object by the 2nd
656 In list context, returns (result,remainder).
657 NOTE: this is integer math, so no
658 fractional part will be returned.
659 _sub(obj,obj) Simple subtraction of 1 object from another
660 a third, optional parameter indicates that the params
661 are swapped. In this case, the first param needs to
662 be preserved, while you can destroy the second.
663 sub (x,y,1) => return x - y and keep x intact!
665 _acmp(obj,obj) <=> operator for objects (return -1, 0 or 1)
667 _len(obj) returns count of the decimal digits of the object
668 _digit(obj,n) returns the n'th decimal digit of object
670 _is_one(obj) return true if argument is +1
671 _is_zero(obj) return true if argument is 0
672 _is_even(obj) return true if argument is even (0,2,4,6..)
673 _is_odd(obj) return true if argument is odd (1,3,5,7..)
675 _copy return a ref to a true copy of the object
677 _check(obj) check whether internal representation is still intact
678 return 0 for ok, otherwise error message as string
680 The following functions are optional, and can be exported if the underlying lib
681 has a fast way to do them. If not defined, Math::BigInt will use a pure, but
682 slow, Perl function as fallback to emulate these:
684 _from_hex(str) return ref to new object from ref to hexadecimal string
685 _from_bin(str) return ref to new object from ref to binary string
687 _rsft(obj,N,B) shift object in base B by N 'digits' right
688 _lsft(obj,N,B) shift object in base B by N 'digits' left
690 _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
691 Mote: XOR, AND and OR pad with zeros if size mismatches
692 _and(obj1,obj2) AND (bit-wise) object 1 with object 2
693 _or(obj1,obj2) OR (bit-wise) object 1 with object 2
695 _sqrt(obj) return the square root of object
696 _pow(obj,obj) return object 1 to the power of object 2
697 _gcd(obj,obj) return Greatest Common Divisor of two objects
699 _zeros(obj) return number of trailing decimal zeros
701 _dec(obj) decrement object by one (input is >= 1)
702 _inc(obj) increment object by one
704 Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
707 Testing of input parameter validity is done by the caller, so you need not
708 worry about underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by
709 zero or similar cases.
711 The first parameter can be modified, that includes the possibility that you
712 return a reference to a completely different object instead. Although keeping
713 the reference the same is prefered.
715 Return values are always references to objects or strings. Exceptions are
716 C<_lsft()> and C<_rsft()>, which return undef if they can not shift the
717 argument. This is used to delegate shifting of bases different than 10 back
718 to BigInt, which will use some generic code to calculate the result.
722 If you want to port your own favourite c-lib for big numbers to the
723 Math::BigInt interface, you can take any of the already existing modules as
724 a rough guideline. You should really wrap up the latest BigInt and BigFloat
725 testsuites with your module, and replace the following line:
731 use Math::BigInt lib => 'yourlib';
733 This way you ensure that your library really works 100% within Math::BigInt.
737 This program is free software; you may redistribute it and/or modify it under
738 the same terms as Perl itself.
742 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
744 Seperated from BigInt and shaped API with the help of John Peacock.
748 L<Math::BigInt>, L<Math::BigFloat>, L<Math::BigInt::BitVect> and
749 L<Math::BigInt::Pari>.