1 package Math::BigInt::Calc;
5 # use warnings; # dont use warnings for older Perls
9 # Package to store unsigned big integers in decimal and do math with them
11 # Internally the numbers are stored in an array with at least 1 element, no
12 # leading zero parts (except the first) and in base 1eX where X is determined
13 # automatically at loading time to be the maximum possible value
16 # - fully remove funky $# stuff in div() (maybe - that code scares me...)
18 # USE_MUL: due to problems on certain os (os390, posix-bc) "* 1e-5" is used
19 # instead of "/ 1e5" at some places, (marked with USE_MUL). Other platforms
20 # BS2000, some Crays need USE_DIV instead.
21 # The BEGIN block is used to determine which of the two variants gives the
24 # Beware of things like:
25 # $i = $i * $y + $car; $car = int($i / $BASE); $i = $i % $BASE;
26 # This works on x86, but fails on ARM (SA1100, iPAQ) due to whoknows what
27 # reasons. So, use this instead (slower, but correct):
28 # $i = $i * $y + $car; $car = int($i / $BASE); $i -= $BASE * $car;
30 ##############################################################################
31 # global constants, flags and accessory
33 # announce that we are compatible with MBI v1.83 and up
34 sub api_version () { 2; }
36 # constants for easier life
37 my ($BASE,$BASE_LEN,$RBASE,$MAX_VAL);
38 my ($AND_BITS,$XOR_BITS,$OR_BITS);
39 my ($AND_MASK,$XOR_MASK,$OR_MASK);
43 # Set/get the BASE_LEN and assorted other, connected values.
44 # Used only by the testsuite, the set variant is used only by the BEGIN
55 if ($] >= 5.008 && $int && $b > 7)
58 *_mul = \&_mul_use_div_64;
59 *_div = \&_div_use_div_64;
60 $BASE = int("1e".$BASE_LEN);
62 return $BASE_LEN unless wantarray;
63 return ($BASE_LEN, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN, $MAX_VAL, $BASE);
66 # find whether we can use mul or div in mul()/div()
69 while (--$BASE_LEN > 5)
71 $BASE = int("1e".$BASE_LEN);
72 $RBASE = abs('1e-'.$BASE_LEN); # see USE_MUL
74 $caught += 1 if (int($BASE * $RBASE) != 1); # should be 1
75 $caught += 2 if (int($BASE / $BASE) != 1); # should be 1
78 $BASE = int("1e".$BASE_LEN);
79 $RBASE = abs('1e-'.$BASE_LEN); # see USE_MUL
82 # ($caught & 1) != 0 => cannot use MUL
83 # ($caught & 2) != 0 => cannot use DIV
86 # must USE_MUL since we cannot use DIV
87 *_mul = \&_mul_use_mul;
88 *_div = \&_div_use_mul;
93 *_mul = \&_mul_use_div;
94 *_div = \&_div_use_div;
97 return $BASE_LEN unless wantarray;
98 return ($BASE_LEN, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN, $MAX_VAL, $BASE);
103 # (ref to string) return ref to num_array
104 # Convert a number from string format (without sign) to internal base
105 # 1ex format. Assumes normalized value as input.
106 my $il = length($_[1])-1;
108 # < BASE_LEN due len-1 above
109 return [ int($_[1]) ] if $il < $BASE_LEN; # shortcut for short numbers
111 # this leaves '00000' instead of int 0 and will be corrected after any op
112 [ reverse(unpack("a" . ($il % $BASE_LEN+1)
113 . ("a$BASE_LEN" x ($il / $BASE_LEN)), $_[1])) ];
118 # from Daniel Pfeiffer: determine largest group of digits that is precisely
119 # multipliable with itself plus carry
120 # Test now changed to expect the proper pattern, not a result off by 1 or 2
121 my ($e, $num) = 3; # lowest value we will use is 3+1-1 = 3
124 $num = ('9' x ++$e) + 0;
126 } while ("$num" =~ /9{$e}0{$e}/); # must be a certain pattern
127 $e--; # last test failed, so retract one step
128 # the limits below brush the problems with the test above under the rug:
129 # the test should be able to find the proper $e automatically
130 $e = 5 if $^O =~ /^uts/; # UTS get's some special treatment
131 $e = 5 if $^O =~ /^unicos/; # unicos is also problematic (6 seems to work
132 # there, but we play safe)
142 $num = ('9' x ++$e1) + 0;
144 } while ("$num" =~ /9{$e1}0{$e1}/); # must be a certain pattern
145 $e1--; # last test failed, so retract one step
152 __PACKAGE__->_base_len($e,$int); # set and store
155 # find out how many bits _and, _or and _xor can take (old default = 16)
156 # I don't think anybody has yet 128 bit scalars, so let's play safe.
157 local $^W = 0; # don't warn about 'nonportable number'
158 $AND_BITS = 15; $XOR_BITS = 15; $OR_BITS = 15;
160 # find max bits, we will not go higher than numberofbits that fit into $BASE
161 # to make _and etc simpler (and faster for smaller, slower for large numbers)
163 while (2 ** $max < $BASE) { $max++; }
166 $max = 16 if $] < 5.006; # older Perls might not take >16 too well
171 $x = CORE::oct('0b' . '1' x $AND_BITS); $y = $x & $x;
172 $z = (2 ** $AND_BITS) - 1;
173 } while ($AND_BITS < $max && $x == $z && $y == $x);
174 $AND_BITS --; # retreat one step
177 $x = CORE::oct('0b' . '1' x $XOR_BITS); $y = $x ^ 0;
178 $z = (2 ** $XOR_BITS) - 1;
179 } while ($XOR_BITS < $max && $x == $z && $y == $x);
180 $XOR_BITS --; # retreat one step
183 $x = CORE::oct('0b' . '1' x $OR_BITS); $y = $x | $x;
184 $z = (2 ** $OR_BITS) - 1;
185 } while ($OR_BITS < $max && $x == $z && $y == $x);
186 $OR_BITS --; # retreat one step
188 $AND_MASK = __PACKAGE__->_new( ( 2 ** $AND_BITS ));
189 $XOR_MASK = __PACKAGE__->_new( ( 2 ** $XOR_BITS ));
190 $OR_MASK = __PACKAGE__->_new( ( 2 ** $OR_BITS ));
192 # We can compute the approximate lenght no faster than the real length:
196 ###############################################################################
212 # create a two (used internally for shifting)
218 # create a 10 (used internally for shifting)
225 my $rem = $_[1] % $BASE_LEN; # remainder
226 my $parts = $_[1] / $BASE_LEN; # parts
228 # 000000, 000000, 100
229 [ (0) x $parts, '1' . ('0' x $rem) ];
238 # catch and throw away
241 ##############################################################################
242 # convert back to string and number
246 # (ref to BINT) return num_str
247 # Convert number from internal base 100000 format to string format.
248 # internal format is always normalized (no leading zeros, "-0" => "+0")
251 my $l = scalar @$ar; # number of parts
252 if ($l < 1) # should not happen
255 Carp::croak("$_[1] has no elements");
259 # handle first one different to strip leading zeros from it (there are no
260 # leading zero parts in internal representation)
261 $l --; $ret .= int($ar->[$l]); $l--;
262 # Interestingly, the pre-padd method uses more time
263 # the old grep variant takes longer (14 vs. 10 sec)
264 my $z = '0' x ($BASE_LEN-1);
267 $ret .= substr($z.$ar->[$l],-$BASE_LEN); # fastest way I could think of
275 # Make a number (scalar int/float) from a BigInt object
278 return 0+$x->[0] if scalar @$x == 1; # below $BASE
283 $num += $fac*$_; $fac *= $BASE;
288 ##############################################################################
293 # (ref to int_num_array, ref to int_num_array)
294 # routine to add two base 1eX numbers
295 # stolen from Knuth Vol 2 Algorithm A pg 231
296 # there are separate routines to add and sub as per Knuth pg 233
297 # This routine clobbers up array x, but not y.
301 return $x if (@$y == 1) && $y->[0] == 0; # $x + 0 => $x
302 if ((@$x == 1) && $x->[0] == 0) # 0 + $y => $y->copy
304 # twice as slow as $x = [ @$y ], but nec. to retain $x as ref :(
305 @$x = @$y; return $x;
308 # for each in Y, add Y to X and carry. If after that, something is left in
309 # X, foreach in X add carry to X and then return X, carry
310 # Trades one "$j++" for having to shift arrays
311 my $i; my $car = 0; my $j = 0;
314 $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
319 $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
326 # (ref to int_num_array, ref to int_num_array)
327 # Add 1 to $x, modify $x in place
332 return $x if (($i += 1) < $BASE); # early out
333 $i = 0; # overflow, next
335 push @$x,1 if (($x->[-1] || 0) == 0); # last overflowed, so extend
341 # (ref to int_num_array, ref to int_num_array)
342 # Sub 1 from $x, modify $x in place
345 my $MAX = $BASE-1; # since MAX_VAL based on BASE
348 last if (($i -= 1) >= 0); # early out
349 $i = $MAX; # underflow, next
351 pop @$x if $x->[-1] == 0 && @$x > 1; # last underflowed (but leave 0)
357 # (ref to int_num_array, ref to int_num_array, swap)
358 # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
359 # subtract Y from X by modifying x in place
360 my ($c,$sx,$sy,$s) = @_;
362 my $car = 0; my $i; my $j = 0;
367 last unless defined $sy->[$j] || $car;
368 $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
370 # might leave leading zeros, so fix that
371 return __strip_zeros($sx);
375 # we can't do an early out if $x is < than $y, since we
376 # need to copy the high chunks from $y. Found by Bob Mathews.
377 #last unless defined $sy->[$j] || $car;
379 if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
382 # might leave leading zeros, so fix that
388 # (ref to int_num_array, ref to int_num_array)
389 # multiply two numbers in internal representation
390 # modifies first arg, second need not be different from first
391 my ($c,$xv,$yv) = @_;
395 # shortcut for two very short numbers (improved by Nathan Zook)
396 # works also if xv and yv are the same reference, and handles also $x == 0
399 if (($xv->[0] *= $yv->[0]) >= $BASE)
401 $xv->[0] = $xv->[0] - ($xv->[1] = int($xv->[0] * $RBASE)) * $BASE;
411 # multiply a large number a by a single element one, so speed up
412 my $y = $yv->[0]; my $car = 0;
415 $i = $i * $y + $car; $car = int($i * $RBASE); $i -= $car * $BASE;
417 push @$xv, $car if $car != 0;
420 # shortcut for result $x == 0 => result = 0
421 return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
423 # since multiplying $x with $x fails, make copy in this case
424 $yv = [@$xv] if $xv == $yv; # same references?
426 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
435 # $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
437 # $prod - ($car = int($prod * RBASE)) * $BASE; # see USE_MUL
439 # $prod[$cty] += $car if $car; # need really to check for 0?
443 # looping through this if $xi == 0 is silly - so optimize it away!
444 $xi = (shift @prod || 0), next if $xi == 0;
447 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
448 ## this is actually a tad slower
449 ## $prod = $prod[$cty]; $prod += ($car + $xi * $yi); # no ||0 here
451 $prod - ($car = int($prod * $RBASE)) * $BASE; # see USE_MUL
453 $prod[$cty] += $car if $car; # need really to check for 0?
454 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
457 # can't have leading zeros
458 # __strip_zeros($xv);
464 # (ref to int_num_array, ref to int_num_array)
465 # multiply two numbers in internal representation
466 # modifies first arg, second need not be different from first
467 # works for 64 bit integer with "use integer"
468 my ($c,$xv,$yv) = @_;
473 # shortcut for two small numbers, also handles $x == 0
476 # shortcut for two very short numbers (improved by Nathan Zook)
477 # works also if xv and yv are the same reference, and handles also $x == 0
478 if (($xv->[0] *= $yv->[0]) >= $BASE)
481 $xv->[0] - ($xv->[1] = $xv->[0] / $BASE) * $BASE;
491 # multiply a large number a by a single element one, so speed up
492 my $y = $yv->[0]; my $car = 0;
495 #$i = $i * $y + $car; $car = $i / $BASE; $i -= $car * $BASE;
496 $i = $i * $y + $car; $i -= ($car = $i / $BASE) * $BASE;
498 push @$xv, $car if $car != 0;
501 # shortcut for result $x == 0 => result = 0
502 return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
504 # since multiplying $x with $x fails, make copy in this case
505 $yv = [@$xv] if $xv == $yv; # same references?
507 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
511 # looping through this if $xi == 0 is silly - so optimize it away!
512 $xi = (shift @prod || 0), next if $xi == 0;
515 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
516 $prod[$cty++] = $prod - ($car = $prod / $BASE) * $BASE;
518 $prod[$cty] += $car if $car; # need really to check for 0?
519 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
527 # (ref to int_num_array, ref to int_num_array)
528 # multiply two numbers in internal representation
529 # modifies first arg, second need not be different from first
530 my ($c,$xv,$yv) = @_;
534 # shortcut for two small numbers, also handles $x == 0
537 # shortcut for two very short numbers (improved by Nathan Zook)
538 # works also if xv and yv are the same reference, and handles also $x == 0
539 if (($xv->[0] *= $yv->[0]) >= $BASE)
542 $xv->[0] - ($xv->[1] = int($xv->[0] / $BASE)) * $BASE;
552 # multiply a large number a by a single element one, so speed up
553 my $y = $yv->[0]; my $car = 0;
556 $i = $i * $y + $car; $car = int($i / $BASE); $i -= $car * $BASE;
557 # This (together with use integer;) does not work on 32-bit Perls
558 #$i = $i * $y + $car; $i -= ($car = $i / $BASE) * $BASE;
560 push @$xv, $car if $car != 0;
563 # shortcut for result $x == 0 => result = 0
564 return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
566 # since multiplying $x with $x fails, make copy in this case
567 $yv = [@$xv] if $xv == $yv; # same references?
569 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
573 # looping through this if $xi == 0 is silly - so optimize it away!
574 $xi = (shift @prod || 0), next if $xi == 0;
577 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
578 $prod[$cty++] = $prod - ($car = int($prod / $BASE)) * $BASE;
580 $prod[$cty] += $car if $car; # need really to check for 0?
581 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
584 # can't have leading zeros
585 # __strip_zeros($xv);
591 # ref to array, ref to array, modify first array and return remainder if
594 # see comments in _div_use_div() for more explanations
596 my ($c,$x,$yorg) = @_;
598 # the general div algorithmn here is about O(N*N) and thus quite slow, so
599 # we first check for some special cases and use shortcuts to handle them.
601 # This works, because we store the numbers in a chunked format where each
602 # element contains 5..7 digits (depending on system).
604 # if both numbers have only one element:
605 if (@$x == 1 && @$yorg == 1)
607 # shortcut, $yorg and $x are two small numbers
610 my $r = [ $x->[0] % $yorg->[0] ];
611 $x->[0] = int($x->[0] / $yorg->[0]);
616 $x->[0] = int($x->[0] / $yorg->[0]);
621 # if x has more than one, but y has only one element:
625 $rem = _mod($c,[ @$x ],$yorg) if wantarray;
627 # shortcut, $y is < $BASE
628 my $j = scalar @$x; my $r = 0;
629 my $y = $yorg->[0]; my $b;
632 $b = $r * $BASE + $x->[$j];
633 $x->[$j] = int($b/$y);
636 pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
637 return ($x,$rem) if wantarray;
641 # now x and y have more than one element
643 # check whether y has more elements than x, if yet, the result will be 0
647 $rem = [@$x] if wantarray; # make copy
648 splice (@$x,1); # keep ref to original array
649 $x->[0] = 0; # set to 0
650 return ($x,$rem) if wantarray; # including remainder?
651 return $x; # only x, which is [0] now
653 # check whether the numbers have the same number of elements, in that case
654 # the result will fit into one element and can be computed efficiently
658 # if $yorg has more digits than $x (it's leading element is longer than
659 # the one from $x), the result will also be 0:
660 if (length(int($yorg->[-1])) > length(int($x->[-1])))
662 $rem = [@$x] if wantarray; # make copy
663 splice (@$x,1); # keep ref to org array
664 $x->[0] = 0; # set to 0
665 return ($x,$rem) if wantarray; # including remainder?
668 # now calculate $x / $yorg
669 if (length(int($yorg->[-1])) == length(int($x->[-1])))
671 # same length, so make full compare
673 my $a = 0; my $j = scalar @$x - 1;
674 # manual way (abort if unequal, good for early ne)
677 last if ($a = $x->[$j] - $yorg->[$j]); $j--;
679 # $a contains the result of the compare between X and Y
680 # a < 0: x < y, a == 0: x == y, a > 0: x > y
683 $rem = [ 0 ]; # a = 0 => x == y => rem 0
684 $rem = [@$x] if $a != 0; # a < 0 => x < y => rem = x
685 splice(@$x,1); # keep single element
686 $x->[0] = 0; # if $a < 0
687 $x->[0] = 1 if $a == 0; # $x == $y
688 return ($x,$rem) if wantarray;
691 # $x >= $y, so proceed normally
697 my $y = [ @$yorg ]; # always make copy to preserve
699 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
701 $car = $bar = $prd = 0;
702 if (($dd = int($BASE/($y->[-1]+1))) != 1)
706 $xi = $xi * $dd + $car;
707 $xi -= ($car = int($xi * $RBASE)) * $BASE; # see USE_MUL
709 push(@$x, $car); $car = 0;
712 $yi = $yi * $dd + $car;
713 $yi -= ($car = int($yi * $RBASE)) * $BASE; # see USE_MUL
720 @q = (); ($v2,$v1) = @$y[-2,-1];
724 ($u2,$u1,$u0) = @$x[-3..-1];
726 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
728 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
729 --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
732 ($car, $bar) = (0,0);
733 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
735 $prd = $q * $y->[$yi] + $car;
736 $prd -= ($car = int($prd * $RBASE)) * $BASE; # see USE_MUL
737 $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
739 if ($x->[-1] < $car + $bar)
742 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
745 if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
758 for $xi (reverse @$x)
760 $prd = $car * $BASE + $xi;
761 $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
782 # ref to array, ref to array, modify first array and return remainder if
784 # This version works on 64 bit integers
785 my ($c,$x,$yorg) = @_;
788 # the general div algorithmn here is about O(N*N) and thus quite slow, so
789 # we first check for some special cases and use shortcuts to handle them.
791 # This works, because we store the numbers in a chunked format where each
792 # element contains 5..7 digits (depending on system).
794 # if both numbers have only one element:
795 if (@$x == 1 && @$yorg == 1)
797 # shortcut, $yorg and $x are two small numbers
800 my $r = [ $x->[0] % $yorg->[0] ];
801 $x->[0] = int($x->[0] / $yorg->[0]);
806 $x->[0] = int($x->[0] / $yorg->[0]);
810 # if x has more than one, but y has only one element:
814 $rem = _mod($c,[ @$x ],$yorg) if wantarray;
816 # shortcut, $y is < $BASE
817 my $j = scalar @$x; my $r = 0;
818 my $y = $yorg->[0]; my $b;
821 $b = $r * $BASE + $x->[$j];
822 $x->[$j] = int($b/$y);
825 pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
826 return ($x,$rem) if wantarray;
829 # now x and y have more than one element
831 # check whether y has more elements than x, if yet, the result will be 0
835 $rem = [@$x] if wantarray; # make copy
836 splice (@$x,1); # keep ref to original array
837 $x->[0] = 0; # set to 0
838 return ($x,$rem) if wantarray; # including remainder?
839 return $x; # only x, which is [0] now
841 # check whether the numbers have the same number of elements, in that case
842 # the result will fit into one element and can be computed efficiently
846 # if $yorg has more digits than $x (it's leading element is longer than
847 # the one from $x), the result will also be 0:
848 if (length(int($yorg->[-1])) > length(int($x->[-1])))
850 $rem = [@$x] if wantarray; # make copy
851 splice (@$x,1); # keep ref to org array
852 $x->[0] = 0; # set to 0
853 return ($x,$rem) if wantarray; # including remainder?
856 # now calculate $x / $yorg
858 if (length(int($yorg->[-1])) == length(int($x->[-1])))
860 # same length, so make full compare
862 my $a = 0; my $j = scalar @$x - 1;
863 # manual way (abort if unequal, good for early ne)
866 last if ($a = $x->[$j] - $yorg->[$j]); $j--;
868 # $a contains the result of the compare between X and Y
869 # a < 0: x < y, a == 0: x == y, a > 0: x > y
872 $rem = [ 0 ]; # a = 0 => x == y => rem 0
873 $rem = [@$x] if $a != 0; # a < 0 => x < y => rem = x
874 splice(@$x,1); # keep single element
875 $x->[0] = 0; # if $a < 0
876 $x->[0] = 1 if $a == 0; # $x == $y
877 return ($x,$rem) if wantarray; # including remainder?
880 # $x >= $y, so proceed normally
887 my $y = [ @$yorg ]; # always make copy to preserve
889 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
891 $car = $bar = $prd = 0;
892 if (($dd = int($BASE/($y->[-1]+1))) != 1)
896 $xi = $xi * $dd + $car;
897 $xi -= ($car = int($xi / $BASE)) * $BASE;
899 push(@$x, $car); $car = 0;
902 $yi = $yi * $dd + $car;
903 $yi -= ($car = int($yi / $BASE)) * $BASE;
911 # @q will accumulate the final result, $q contains the current computed
912 # part of the final result
914 @q = (); ($v2,$v1) = @$y[-2,-1];
918 ($u2,$u1,$u0) = @$x[-3..-1];
920 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
922 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
923 --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
926 ($car, $bar) = (0,0);
927 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
929 $prd = $q * $y->[$yi] + $car;
930 $prd -= ($car = int($prd / $BASE)) * $BASE;
931 $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
933 if ($x->[-1] < $car + $bar)
936 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
939 if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
943 pop(@$x); unshift(@q, $q);
951 for $xi (reverse @$x)
953 $prd = $car * $BASE + $xi;
954 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
975 # ref to array, ref to array, modify first array and return remainder if
977 my ($c,$x,$yorg) = @_;
979 # the general div algorithmn here is about O(N*N) and thus quite slow, so
980 # we first check for some special cases and use shortcuts to handle them.
982 # This works, because we store the numbers in a chunked format where each
983 # element contains 5..7 digits (depending on system).
985 # if both numbers have only one element:
986 if (@$x == 1 && @$yorg == 1)
988 # shortcut, $yorg and $x are two small numbers
991 my $r = [ $x->[0] % $yorg->[0] ];
992 $x->[0] = int($x->[0] / $yorg->[0]);
997 $x->[0] = int($x->[0] / $yorg->[0]);
1001 # if x has more than one, but y has only one element:
1005 $rem = _mod($c,[ @$x ],$yorg) if wantarray;
1007 # shortcut, $y is < $BASE
1008 my $j = scalar @$x; my $r = 0;
1009 my $y = $yorg->[0]; my $b;
1012 $b = $r * $BASE + $x->[$j];
1013 $x->[$j] = int($b/$y);
1016 pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
1017 return ($x,$rem) if wantarray;
1020 # now x and y have more than one element
1022 # check whether y has more elements than x, if yet, the result will be 0
1026 $rem = [@$x] if wantarray; # make copy
1027 splice (@$x,1); # keep ref to original array
1028 $x->[0] = 0; # set to 0
1029 return ($x,$rem) if wantarray; # including remainder?
1030 return $x; # only x, which is [0] now
1032 # check whether the numbers have the same number of elements, in that case
1033 # the result will fit into one element and can be computed efficiently
1037 # if $yorg has more digits than $x (it's leading element is longer than
1038 # the one from $x), the result will also be 0:
1039 if (length(int($yorg->[-1])) > length(int($x->[-1])))
1041 $rem = [@$x] if wantarray; # make copy
1042 splice (@$x,1); # keep ref to org array
1043 $x->[0] = 0; # set to 0
1044 return ($x,$rem) if wantarray; # including remainder?
1047 # now calculate $x / $yorg
1049 if (length(int($yorg->[-1])) == length(int($x->[-1])))
1051 # same length, so make full compare
1053 my $a = 0; my $j = scalar @$x - 1;
1054 # manual way (abort if unequal, good for early ne)
1057 last if ($a = $x->[$j] - $yorg->[$j]); $j--;
1059 # $a contains the result of the compare between X and Y
1060 # a < 0: x < y, a == 0: x == y, a > 0: x > y
1063 $rem = [ 0 ]; # a = 0 => x == y => rem 0
1064 $rem = [@$x] if $a != 0; # a < 0 => x < y => rem = x
1065 splice(@$x,1); # keep single element
1066 $x->[0] = 0; # if $a < 0
1067 $x->[0] = 1 if $a == 0; # $x == $y
1068 return ($x,$rem) if wantarray; # including remainder?
1071 # $x >= $y, so proceed normally
1078 my $y = [ @$yorg ]; # always make copy to preserve
1080 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
1082 $car = $bar = $prd = 0;
1083 if (($dd = int($BASE/($y->[-1]+1))) != 1)
1087 $xi = $xi * $dd + $car;
1088 $xi -= ($car = int($xi / $BASE)) * $BASE;
1090 push(@$x, $car); $car = 0;
1093 $yi = $yi * $dd + $car;
1094 $yi -= ($car = int($yi / $BASE)) * $BASE;
1102 # @q will accumulate the final result, $q contains the current computed
1103 # part of the final result
1105 @q = (); ($v2,$v1) = @$y[-2,-1];
1109 ($u2,$u1,$u0) = @$x[-3..-1];
1111 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
1113 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
1114 --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
1117 ($car, $bar) = (0,0);
1118 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
1120 $prd = $q * $y->[$yi] + $car;
1121 $prd -= ($car = int($prd / $BASE)) * $BASE;
1122 $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
1124 if ($x->[-1] < $car + $bar)
1127 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
1130 if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
1134 pop(@$x); unshift(@q, $q);
1142 for $xi (reverse @$x)
1144 $prd = $car * $BASE + $xi;
1145 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
1164 ##############################################################################
1169 # internal absolute post-normalized compare (ignore signs)
1170 # ref to array, ref to array, return <0, 0, >0
1171 # arrays must have at least one entry; this is not checked for
1172 my ($c,$cx,$cy) = @_;
1174 # shortcut for short numbers
1175 return (($cx->[0] <=> $cy->[0]) <=> 0)
1176 if scalar @$cx == scalar @$cy && scalar @$cx == 1;
1178 # fast comp based on number of array elements (aka pseudo-length)
1179 my $lxy = (scalar @$cx - scalar @$cy)
1180 # or length of first element if same number of elements (aka difference 0)
1182 # need int() here because sometimes the last element is '00018' vs '18'
1183 (length(int($cx->[-1])) - length(int($cy->[-1])));
1184 return -1 if $lxy < 0; # already differs, ret
1185 return 1 if $lxy > 0; # ditto
1187 # manual way (abort if unequal, good for early ne)
1188 my $a; my $j = scalar @$cx;
1191 last if ($a = $cx->[$j] - $cy->[$j]);
1198 # compute number of digits in base 10
1200 # int() because add/sub sometimes leaves strings (like '00005') instead of
1201 # '5' in this place, thus causing length() to report wrong length
1204 (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
1209 # return the nth digit, negative values count backward
1210 # zero is rightmost, so _digit(123,0) will give 3
1213 my $len = _len('',$x);
1215 $n = $len+$n if $n < 0; # -1 last, -2 second-to-last
1216 $n = abs($n); # if negative was too big
1217 $len--; $n = $len if $n > $len; # n to big?
1219 my $elem = int($n / $BASE_LEN); # which array element
1220 my $digit = $n % $BASE_LEN; # which digit in this element
1221 $elem = '0' x $BASE_LEN . @$x[$elem]; # get element padded with 0's
1222 substr($elem,-$digit-1,1);
1227 # return amount of trailing zeros in decimal
1228 # check each array elem in _m for having 0 at end as long as elem == 0
1229 # Upon finding a elem != 0, stop
1232 return 0 if scalar @$x == 1 && $x->[0] == 0;
1234 my $zeros = 0; my $elem;
1239 $elem = "$e"; # preserve x
1240 $elem =~ s/.*?(0*$)/$1/; # strip anything not zero
1241 $zeros *= $BASE_LEN; # elems * 5
1242 $zeros += length($elem); # count trailing zeros
1245 $zeros ++; # real else branch: 50% slower!
1250 ##############################################################################
1255 # return true if arg is zero
1256 (((scalar @{$_[1]} == 1) && ($_[1]->[0] == 0))) <=> 0;
1261 # return true if arg is even
1262 (!($_[1]->[0] & 1)) <=> 0;
1267 # return true if arg is even
1268 (($_[1]->[0] & 1)) <=> 0;
1273 # return true if arg is one
1274 (scalar @{$_[1]} == 1) && ($_[1]->[0] == 1) <=> 0;
1279 # return true if arg is two
1280 (scalar @{$_[1]} == 1) && ($_[1]->[0] == 2) <=> 0;
1285 # return true if arg is ten
1286 (scalar @{$_[1]} == 1) && ($_[1]->[0] == 10) <=> 0;
1291 # internal normalization function that strips leading zeros from the array
1292 # args: ref to array
1295 my $cnt = scalar @$s; # get count of parts
1297 push @$s,0 if $i < 0; # div might return empty results, so fix it
1299 return $s if @$s == 1; # early out
1301 #print "strip: cnt $cnt i $i\n";
1302 # '0', '3', '4', '0', '0',
1307 # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
1308 # >= 1: skip first part (this can be zero)
1309 while ($i > 0) { last if $s->[$i] != 0; $i--; }
1310 $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
1314 ###############################################################################
1315 # check routine to test internal state for corruptions
1319 # used by the test suite
1322 return "$x is not a reference" if !ref($x);
1324 # are all parts are valid?
1325 my $i = 0; my $j = scalar @$x; my ($e,$try);
1328 $e = $x->[$i]; $e = 'undef' unless defined $e;
1329 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
1330 last if $e !~ /^[+]?[0-9]+$/;
1331 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
1332 last if "$e" !~ /^[+]?[0-9]+$/;
1333 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
1334 last if '' . "$e" !~ /^[+]?[0-9]+$/;
1335 $try = ' < 0 || >= $BASE; '."($x, $e)";
1336 last if $e <0 || $e >= $BASE;
1337 # this test is disabled, since new/bnorm and certain ops (like early out
1338 # in add/sub) are allowed/expected to leave '00000' in some elements
1339 #$try = '=~ /^00+/; '."($x, $e)";
1340 #last if $e =~ /^00+/;
1343 return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
1348 ###############################################################################
1352 # if possible, use mod shortcut
1353 my ($c,$x,$yo) = @_;
1355 # slow way since $y to big
1356 if (scalar @$yo > 1)
1358 my ($xo,$rem) = _div($c,$x,$yo);
1363 # both are single element arrays
1364 if (scalar @$x == 1)
1370 # @y is a single element, but @x has more than one element
1374 # when BASE % Y == 0 then (B * BASE) % Y == 0
1375 # (B * BASE) % $y + A % Y => A % Y
1376 # so need to consider only last element: O(1)
1381 # else need to go through all elements: O(N), but loop is a bit simplified
1385 $r = ($r + $_) % $y; # not much faster, but heh...
1386 #$r += $_ % $y; $r %= $y;
1393 # else need to go through all elements: O(N)
1394 my $r = 0; my $bm = 1;
1397 $r = ($_ * $bm + $r) % $y;
1398 $bm = ($bm * $b) % $y;
1400 #$r += ($_ % $y) * $bm;
1408 splice (@$x,1); # keep one element of $x
1412 ##############################################################################
1417 my ($c,$x,$y,$n) = @_;
1421 $n = _new($c,$n); return _div($c,$x, _pow($c,$n,$y));
1424 # shortcut (faster) for shifting by 10)
1425 # multiples of $BASE_LEN
1426 my $dst = 0; # destination
1427 my $src = _num($c,$y); # as normal int
1428 my $xlen = (@$x-1)*$BASE_LEN+length(int($x->[-1])); # len of x in digits
1429 if ($src >= $xlen or ($src == $xlen and ! defined $x->[1]))
1431 # 12345 67890 shifted right by more than 10 digits => 0
1432 splice (@$x,1); # leave only one element
1433 $x->[0] = 0; # set to zero
1436 my $rem = $src % $BASE_LEN; # remainder to shift
1437 $src = int($src / $BASE_LEN); # source
1440 splice (@$x,0,$src); # even faster, 38.4 => 39.3
1444 my $len = scalar @$x - $src; # elems to go
1445 my $vd; my $z = '0'x $BASE_LEN;
1446 $x->[scalar @$x] = 0; # avoid || 0 test inside loop
1449 $vd = $z.$x->[$src];
1450 $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
1452 $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
1453 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1454 $x->[$dst] = int($vd);
1457 splice (@$x,$dst) if $dst > 0; # kill left-over array elems
1458 pop @$x if $x->[-1] == 0 && @$x > 1; # kill last element if 0
1465 my ($c,$x,$y,$n) = @_;
1469 $n = _new($c,$n); return _mul($c,$x, _pow($c,$n,$y));
1472 # shortcut (faster) for shifting by 10) since we are in base 10eX
1473 # multiples of $BASE_LEN:
1474 my $src = scalar @$x; # source
1475 my $len = _num($c,$y); # shift-len as normal int
1476 my $rem = $len % $BASE_LEN; # remainder to shift
1477 my $dst = $src + int($len/$BASE_LEN); # destination
1478 my $vd; # further speedup
1479 $x->[$src] = 0; # avoid first ||0 for speed
1480 my $z = '0' x $BASE_LEN;
1483 $vd = $x->[$src]; $vd = $z.$vd;
1484 $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
1485 $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
1486 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1487 $x->[$dst] = int($vd);
1490 # set lowest parts to 0
1491 while ($dst >= 0) { $x->[$dst--] = 0; }
1492 # fix spurios last zero element
1493 splice @$x,-1 if $x->[-1] == 0;
1500 # ref to array, ref to array, return ref to array
1501 my ($c,$cx,$cy) = @_;
1503 if (scalar @$cy == 1 && $cy->[0] == 0)
1505 splice (@$cx,1); $cx->[0] = 1; # y == 0 => x => 1
1508 if ((scalar @$cx == 1 && $cx->[0] == 1) || # x == 1
1509 (scalar @$cy == 1 && $cy->[0] == 1)) # or y == 1
1513 if (scalar @$cx == 1 && $cx->[0] == 0)
1515 splice (@$cx,1); $cx->[0] = 0; # 0 ** y => 0 (if not y <= 0)
1521 my $y_bin = _as_bin($c,$cy); $y_bin =~ s/^0b//;
1522 my $len = length($y_bin);
1525 _mul($c,$pow2,$cx) if substr($y_bin,$len,1) eq '1'; # is odd?
1536 # ref to array, return ref to array
1539 # ( 7 ) 7! 7*6*5 * 4*3*2*1 7 * 6 * 5
1540 # ( - ) = --------- = --------------- = ---------
1541 # ( 3 ) 3! (7-3)! 3*2*1 * 4*3*2*1 3 * 2 * 1
1543 # compute n - k + 2 (so we start with 5 in the example above)
1544 my $x = _copy($c,$n);
1547 if (!_is_one($c,$n))
1550 my $f = _copy($c,$n); _inc($c,$f); # n = 5, f = 6, d = 2
1552 while (_acmp($c,$f,$x) <= 0) # f < n ?
1554 # n = (n * f / d) == 5 * 6 / 2 => n == 3
1555 $n = _mul($c,$n,$f); $n = _div($c,$n,$d);
1557 _inc($c,$f); _inc($c,$d);
1562 # keep ref to $n and set it to 1
1563 splice (@$n,1); $n->[0] = 1;
1582 # ref to array, return ref to array
1585 if ((@$cx == 1) && ($cx->[0] <= 7))
1587 $cx->[0] = $factorials[$cx->[0]]; # 0 => 1, 1 => 1, 2 => 2 etc.
1591 if ((@$cx == 1) && # we do this only if $x >= 12 and $x <= 7000
1592 ($cx->[0] >= 12 && $cx->[0] < 7000))
1595 # Calculate (k-j) * (k-j+1) ... k .. (k+j-1) * (k + j)
1596 # See http://blogten.blogspot.com/2007/01/calculating-n.html
1597 # The above series can be expressed as factors:
1598 # k * k - (j - i) * 2
1599 # We cache k*k, and calculate (j * j) as the sum of the first j odd integers
1601 # This will not work when N exceeds the storage of a Perl scalar, however,
1602 # in this case the algorithm would be way to slow to terminate, anyway.
1604 # As soon as the last element of $cx is 0, we split it up and remember
1605 # how many zeors we got so far. The reason is that n! will accumulate
1606 # zeros at the end rather fast.
1607 my $zero_elements = 0;
1609 # If n is even, set n = n -1
1610 my $k = _num($c,$cx); my $even = 1;
1615 # set k to the center point
1617 # print "k $k even: $even\n";
1618 # now calculate k * k
1620 my $odd = 1; my $sum = 1;
1622 # keep reference to x
1623 my $new_x = _new($c, $k * $even);
1627 $zero_elements ++; shift @$cx;
1629 # print STDERR "x = ", _str($c,$cx),"\n";
1630 my $BASE2 = int(sqrt($BASE))-1;
1634 my $m = ($k2 - $sum); $odd += 2; $sum += $odd; $j++;
1635 while ($j <= $i && ($m < $BASE2) && (($k2 - $sum) < $BASE2))
1638 $odd += 2; $sum += $odd; $j++;
1639 # print STDERR "\n k2 $k2 m $m sum $sum odd $odd\n"; sleep(1);
1647 _mul($c,$cx,$c->_new($m));
1651 $zero_elements ++; shift @$cx;
1653 # print STDERR "Calculate $k2 - $sum = $m (x = ", _str($c,$cx),")\n";
1655 # multiply in the zeros again
1656 unshift @$cx, (0) x $zero_elements;
1660 # go forward until $base is exceeded
1661 # limit is either $x steps (steps == 100 means a result always too high) or
1663 my $steps = 100; $steps = $cx->[0] if @$cx == 1;
1664 my $r = 2; my $cf = 3; my $step = 2; my $last = $r;
1665 while ($r*$cf < $BASE && $step < $steps)
1667 $last = $r; $r *= $cf++; $step++;
1669 if ((@$cx == 1) && $step == $cx->[0])
1671 # completely done, so keep reference to $x and return
1676 # now we must do the left over steps
1677 my $n; # steps still to do
1678 if (scalar @$cx == 1)
1687 # Set $cx to the last result below $BASE (but keep ref to $x)
1688 $cx->[0] = $last; splice (@$cx,1);
1689 # As soon as the last element of $cx is 0, we split it up and remember
1690 # how many zeors we got so far. The reason is that n! will accumulate
1691 # zeros at the end rather fast.
1692 my $zero_elements = 0;
1694 # do left-over steps fit into a scalar?
1695 if (ref $n eq 'ARRAY')
1697 # No, so use slower inc() & cmp()
1698 # ($n is at least $BASE here)
1699 my $base_2 = int(sqrt($BASE)) - 1;
1700 #print STDERR "base_2: $base_2\n";
1701 while ($step < $base_2)
1705 $zero_elements ++; shift @$cx;
1707 my $b = $step * ($step + 1); $step += 2;
1711 while (_acmp($c,$step,$n) <= 0)
1715 $zero_elements ++; shift @$cx;
1717 _mul($c,$cx,$step); _inc($c,$step);
1722 # Yes, so we can speed it up slightly
1724 # print "# left over steps $n\n";
1726 my $base_4 = int(sqrt(sqrt($BASE))) - 2;
1727 #print STDERR "base_4: $base_4\n";
1729 while ($step < $n4 && $step < $base_4)
1733 $zero_elements ++; shift @$cx;
1735 my $b = $step * ($step + 1); $step += 2; $b *= $step * ($step + 1); $step += 2;
1738 my $base_2 = int(sqrt($BASE)) - 1;
1740 #print STDERR "base_2: $base_2\n";
1741 while ($step < $n2 && $step < $base_2)
1745 $zero_elements ++; shift @$cx;
1747 my $b = $step * ($step + 1); $step += 2;
1750 # do what's left over
1753 _mul($c,$cx,[$step]); $step++;
1756 $zero_elements ++; shift @$cx;
1760 # multiply in the zeros again
1761 unshift @$cx, (0) x $zero_elements;
1762 $cx; # return result
1765 #############################################################################
1769 # calculate integer log of $x to base $base
1770 # ref to array, ref to array - return ref to array
1771 my ($c,$x,$base) = @_;
1774 return if (scalar @$x == 1 && $x->[0] == 0);
1775 # BASE 0 or 1 => NaN
1776 return if (scalar @$base == 1 && $base->[0] < 2);
1777 my $cmp = _acmp($c,$x,$base); # X == BASE => 1
1780 splice (@$x,1); $x->[0] = 1;
1786 splice (@$x,1); $x->[0] = 0;
1790 my $x_org = _copy($c,$x); # preserve x
1791 splice(@$x,1); $x->[0] = 1; # keep ref to $x
1793 # Compute a guess for the result based on:
1794 # $guess = int ( length_in_base_10(X) / ( log(base) / log(10) ) )
1795 my $len = _len($c,$x_org);
1796 my $log = log($base->[-1]) / log(10);
1798 # for each additional element in $base, we add $BASE_LEN to the result,
1799 # based on the observation that log($BASE,10) is BASE_LEN and
1800 # log(x*y) == log(x) + log(y):
1801 $log += ((scalar @$base)-1) * $BASE_LEN;
1803 # calculate now a guess based on the values obtained above:
1804 my $res = int($len / $log);
1807 my $trial = _pow ($c, _copy($c, $base), $x);
1808 my $a = _acmp($c,$trial,$x_org);
1810 # print STDERR "# trial ", _str($c,$x)," was: $a (0 = exact, -1 too small, +1 too big)\n";
1812 # found an exact result?
1813 return ($x,1) if $a == 0;
1818 _div($c,$trial,$base); _dec($c, $x);
1819 while (($a = _acmp($c,$trial,$x_org)) > 0)
1821 # print STDERR "# big _log_int at ", _str($c,$x), "\n";
1822 _div($c,$trial,$base); _dec($c, $x);
1824 # result is now exact (a == 0), or too small (a < 0)
1825 return ($x, $a == 0 ? 1 : 0);
1828 # else: result was to small
1829 _mul($c,$trial,$base);
1831 # did we now get the right result?
1832 $a = _acmp($c,$trial,$x_org);
1834 if ($a == 0) # yes, exactly
1839 return ($x,0) if $a > 0;
1841 # Result still too small (we should come here only if the estimate above
1842 # was very off base):
1844 # Now let the normal trial run obtain the real result
1845 # Simple loop that increments $x by 2 in each step, possible overstepping
1848 my $base_mul = _mul($c, _copy($c,$base), $base); # $base * $base
1850 while (($a = _acmp($c,$trial,$x_org)) < 0)
1852 # print STDERR "# small _log_int at ", _str($c,$x), "\n";
1853 _mul($c,$trial,$base_mul); _add($c, $x, [2]);
1859 # overstepped the result
1861 _div($c,$trial,$base);
1862 $a = _acmp($c,$trial,$x_org);
1867 $exact = 0 if $a != 0; # a = -1 => not exact result, a = 0 => exact
1870 ($x,$exact); # return result
1874 use constant DEBUG => 0;
1876 sub steps { $steps };
1880 # square-root of $x in place
1881 # Compute a guess of the result (by rule of thumb), then improve it via
1885 if (scalar @$x == 1)
1887 # fits into one Perl scalar, so result can be computed directly
1888 $x->[0] = int(sqrt($x->[0]));
1891 my $y = _copy($c,$x);
1892 # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
1893 # since our guess will "grow"
1894 my $l = int((_len($c,$x)-1) / 2);
1896 my $lastelem = $x->[-1]; # for guess
1897 my $elems = scalar @$x - 1;
1898 # not enough digits, but could have more?
1899 if ((length($lastelem) <= 3) && ($elems > 1))
1901 # right-align with zero pad
1902 my $len = length($lastelem) & 1;
1903 print "$lastelem => " if DEBUG;
1904 $lastelem .= substr($x->[-2] . '0' x $BASE_LEN,0,$BASE_LEN);
1905 # former odd => make odd again, or former even to even again
1906 $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
1907 print "$lastelem\n" if DEBUG;
1910 # construct $x (instead of _lsft($c,$x,$l,10)
1911 my $r = $l % $BASE_LEN; # 10000 00000 00000 00000 ($BASE_LEN=5)
1912 $l = int($l / $BASE_LEN);
1913 print "l = $l " if DEBUG;
1915 splice @$x,$l; # keep ref($x), but modify it
1917 # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
1919 # 14400 00000 => sqrt(14400) => guess first digits to be 120
1920 # 144000 000000 => sqrt(144000) => guess 379
1922 print "$lastelem (elems $elems) => " if DEBUG;
1923 $lastelem = $lastelem / 10 if ($elems & 1 == 1); # odd or even?
1924 my $g = sqrt($lastelem); $g =~ s/\.//; # 2.345 => 2345
1925 $r -= 1 if $elems & 1 == 0; # 70 => 7
1927 # padd with zeros if result is too short
1928 $x->[$l--] = int(substr($g . '0' x $r,0,$r+1));
1929 print "now ",$x->[-1] if DEBUG;
1930 print " would have been ", int('1' . '0' x $r),"\n" if DEBUG;
1932 # If @$x > 1, we could compute the second elem of the guess, too, to create
1933 # an even better guess. Not implemented yet. Does it improve performance?
1934 $x->[$l--] = 0 while ($l >= 0); # all other digits of guess are zero
1936 print "start x= ",_str($c,$x),"\n" if DEBUG;
1939 my $lastlast = _zero();
1940 $steps = 0 if DEBUG;
1941 while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
1944 $lastlast = _copy($c,$last);
1945 $last = _copy($c,$x);
1946 _add($c,$x, _div($c,_copy($c,$y),$x));
1948 print " x= ",_str($c,$x),"\n" if DEBUG;
1950 print "\nsteps in sqrt: $steps, " if DEBUG;
1951 _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0; # overshot?
1952 print " final ",$x->[-1],"\n" if DEBUG;
1958 # take n'th root of $x in place (n >= 3)
1961 if (scalar @$x == 1)
1965 # result will always be smaller than 2 so trunc to 1 at once
1970 # fits into one Perl scalar, so result can be computed directly
1971 # cannot use int() here, because it rounds wrongly (try
1972 # (81 ** 3) ** (1/3) to see what I mean)
1973 #$x->[0] = int( $x->[0] ** (1 / $n->[0]) );
1974 # round to 8 digits, then truncate result to integer
1975 $x->[0] = int ( sprintf ("%.8f", $x->[0] ** (1 / $n->[0]) ) );
1980 # we know now that X is more than one element long
1982 # if $n is a power of two, we can repeatedly take sqrt($X) and find the
1983 # proper result, because sqrt(sqrt($x)) == root($x,4)
1984 my $b = _as_bin($c,$n);
1985 if ($b =~ /0b1(0+)$/)
1987 my $count = CORE::length($1); # 0b100 => len('00') => 2
1988 my $cnt = $count; # counter for loop
1989 unshift (@$x, 0); # add one element, together with one
1990 # more below in the loop this makes 2
1993 # 'inflate' $X by adding one element, basically computing
1994 # $x * $BASE * $BASE. This gives us more $BASE_LEN digits for result
1995 # since len(sqrt($X)) approx == len($x) / 2.
1997 # calculate sqrt($x), $x is now one element to big, again. In the next
1998 # round we make that two, again.
2001 # $x is now one element to big, so truncate result by removing it
2006 # trial computation by starting with 2,4,8,16 etc until we overstep
2010 # while still to do more than X steps
2014 while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
2016 _mul ($c, $step, [2]);
2017 _add ($c, $trial, $step);
2021 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) == 0)
2023 @$x = @$trial; # make copy while preserving ref to $x
2026 # overstepped, so go back on step
2027 _sub($c, $trial, $step);
2028 } while (scalar @$step > 1 || $step->[0] > 128);
2032 # add two, because $trial cannot be exactly the result (otherwise we would
2033 # alrady have found it)
2034 _add($c, $trial, $step);
2036 # and now add more and more (2,4,6,8,10 etc)
2037 while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
2039 _add ($c, $trial, $step);
2042 # hit not exactly? (overstepped)
2043 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
2048 # hit not exactly? (overstepped)
2049 # 80 too small, 81 slightly too big, 82 too big
2050 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
2055 @$x = @$trial; # make copy while preserving ref to $x
2061 ##############################################################################
2068 # the shortcut makes equal, large numbers _really_ fast, and makes only a
2069 # very small performance drop for small numbers (e.g. something with less
2070 # than 32 bit) Since we optimize for large numbers, this is enabled.
2071 return $x if _acmp($c,$x,$y) == 0; # shortcut
2073 my $m = _one(); my ($xr,$yr);
2074 my $mask = $AND_MASK;
2077 my $y1 = _copy($c,$y); # make copy
2081 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2083 ($x1, $xr) = _div($c,$x1,$mask);
2084 ($y1, $yr) = _div($c,$y1,$mask);
2086 # make ints() from $xr, $yr
2087 # this is when the AND_BITS are greater than $BASE and is slower for
2088 # small (<256 bits) numbers, but faster for large numbers. Disabled
2089 # due to KISS principle
2091 # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2092 # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2093 # _add($c,$x, _mul($c, _new( $c, ($xrr & $yrr) ), $m) );
2095 # 0+ due to '&' doesn't work in strings
2096 _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
2106 return _zero() if _acmp($c,$x,$y) == 0; # shortcut (see -and)
2108 my $m = _one(); my ($xr,$yr);
2109 my $mask = $XOR_MASK;
2112 my $y1 = _copy($c,$y); # make copy
2116 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2118 ($x1, $xr) = _div($c,$x1,$mask);
2119 ($y1, $yr) = _div($c,$y1,$mask);
2120 # make ints() from $xr, $yr (see _and())
2121 #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2122 #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2123 #_add($c,$x, _mul($c, _new( $c, ($xrr ^ $yrr) ), $m) );
2125 # 0+ due to '^' doesn't work in strings
2126 _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
2129 # the loop stops when the shorter of the two numbers is exhausted
2130 # the remainder of the longer one will survive bit-by-bit, so we simple
2131 # multiply-add it in
2132 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
2133 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
2142 return $x if _acmp($c,$x,$y) == 0; # shortcut (see _and)
2144 my $m = _one(); my ($xr,$yr);
2145 my $mask = $OR_MASK;
2148 my $y1 = _copy($c,$y); # make copy
2152 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2154 ($x1, $xr) = _div($c,$x1,$mask);
2155 ($y1, $yr) = _div($c,$y1,$mask);
2156 # make ints() from $xr, $yr (see _and())
2157 # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2158 # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2159 # _add($c,$x, _mul($c, _new( $c, ($xrr | $yrr) ), $m) );
2161 # 0+ due to '|' doesn't work in strings
2162 _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
2165 # the loop stops when the shorter of the two numbers is exhausted
2166 # the remainder of the longer one will survive bit-by-bit, so we simple
2167 # multiply-add it in
2168 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
2169 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
2176 # convert a decimal number to hex (ref to array, return ref to string)
2179 # fits into one element (handle also 0x0 case)
2180 return sprintf("0x%x",$x->[0]) if @$x == 1;
2182 my $x1 = _copy($c,$x);
2185 my ($xr, $h, $x10000);
2188 $x10000 = [ 0x10000 ]; $h = 'h4';
2192 $x10000 = [ 0x1000 ]; $h = 'h3';
2194 while (@$x1 != 1 || $x1->[0] != 0) # _is_zero()
2196 ($x1, $xr) = _div($c,$x1,$x10000);
2197 $es .= unpack($h,pack('V',$xr->[0]));
2200 $es =~ s/^[0]+//; # strip leading zeros
2201 '0x' . $es; # return result prepended with 0x
2206 # convert a decimal number to bin (ref to array, return ref to string)
2209 # fits into one element (and Perl recent enough), handle also 0b0 case
2210 # handle zero case for older Perls
2211 if ($] <= 5.005 && @$x == 1 && $x->[0] == 0)
2213 my $t = '0b0'; return $t;
2215 if (@$x == 1 && $] >= 5.006)
2217 my $t = sprintf("0b%b",$x->[0]);
2220 my $x1 = _copy($c,$x);
2223 my ($xr, $b, $x10000);
2226 $x10000 = [ 0x10000 ]; $b = 'b16';
2230 $x10000 = [ 0x1000 ]; $b = 'b12';
2232 while (!(@$x1 == 1 && $x1->[0] == 0)) # _is_zero()
2234 ($x1, $xr) = _div($c,$x1,$x10000);
2235 $es .= unpack($b,pack('v',$xr->[0]));
2238 $es =~ s/^[0]+//; # strip leading zeros
2239 '0b' . $es; # return result prepended with 0b
2244 # convert a decimal number to octal (ref to array, return ref to string)
2247 # fits into one element (handle also 0 case)
2248 return sprintf("0%o",$x->[0]) if @$x == 1;
2250 my $x1 = _copy($c,$x);
2254 my $x1000 = [ 0100000 ];
2255 while (@$x1 != 1 || $x1->[0] != 0) # _is_zero()
2257 ($x1, $xr) = _div($c,$x1,$x1000);
2258 $es .= reverse sprintf("%05o", $xr->[0]);
2261 $es =~ s/^[0]+//; # strip leading zeros
2262 '0' . $es; # return result prepended with 0
2267 # convert a octal number to decimal (string, return ref to array)
2270 # for older Perls, play safe
2271 my $m = [ 0100000 ];
2272 my $d = 5; # 5 digits at a time
2277 my $len = int( (length($os)-1)/$d ); # $d digit parts, w/o the '0'
2278 my $val; my $i = -$d;
2281 $val = substr($os,$i,$d); # get oct digits
2282 $val = CORE::oct($val);
2284 my $adder = [ $val ];
2285 _add ($c, $x, _mul ($c, $adder, $mul ) ) if $val != 0;
2286 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
2293 # convert a hex number to decimal (string, return ref to array)
2296 my $m = _new($c, 0x10000000); # 28 bit at a time (<32 bit!)
2297 my $d = 7; # 7 digits at a time
2300 # for older Perls, play safe
2301 $m = [ 0x10000 ]; # 16 bit at a time (<32 bit!)
2302 $d = 4; # 4 digits at a time
2308 my $len = int( (length($hs)-2)/$d ); # $d digit parts, w/o the '0x'
2309 my $val; my $i = -$d;
2312 $val = substr($hs,$i,$d); # get hex digits
2313 $val =~ s/^0x// if $len == 0; # for last part only because
2314 $val = CORE::hex($val); # hex does not like wrong chars
2316 my $adder = [ $val ];
2317 # if the resulting number was to big to fit into one element, create a
2318 # two-element version (bug found by Mark Lakata - Thanx!)
2319 if (CORE::length($val) > $BASE_LEN)
2321 $adder = _new($c,$val);
2323 _add ($c, $x, _mul ($c, $adder, $mul ) ) if $val != 0;
2324 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
2331 # convert a hex number to decimal (string, return ref to array)
2334 # instead of converting X (8) bit at a time, it is faster to "convert" the
2335 # number to hex, and then call _from_hex.
2338 $hs =~ s/^[+-]?0b//; # remove sign and 0b
2339 my $l = length($hs); # bits
2340 $hs = '0' x (8-($l % 8)) . $hs if ($l % 8) != 0; # padd left side w/ 0
2341 my $h = '0x' . unpack('H*', pack ('B*', $hs)); # repack as hex
2346 ##############################################################################
2347 # special modulus functions
2354 my $u = _zero($c); my $u1 = _one($c);
2355 my $a = _copy($c,$y); my $b = _copy($c,$x);
2357 # Euclid's Algorithm for bgcd(), only that we calc bgcd() ($a) and the
2358 # result ($u) at the same time. See comments in BigInt for why this works.
2360 ($a, $q, $b) = ($b, _div($c,$a,$b)); # step 1
2362 while (!_is_zero($c,$b))
2364 my $t = _add($c, # step 2:
2365 _mul($c,_copy($c,$u1), $q) , # t = u1 * q
2367 $u = $u1; # u = u1, u1 = t
2370 ($a, $q, $b) = ($b, _div($c,$a,$b)); # step 1
2373 # if the gcd is not 1, then return NaN
2374 return (undef,undef) unless _is_one($c,$a);
2376 ($u1, $sign == 1 ? '+' : '-');
2381 # modulus of power ($x ** $y) % $z
2382 my ($c,$num,$exp,$mod) = @_;
2384 # in the trivial case,
2385 if (_is_one($c,$mod))
2387 splice @$num,0,1; $num->[0] = 0;
2390 if ((scalar @$num == 1) && (($num->[0] == 0) || ($num->[0] == 1)))
2396 # $num = _mod($c,$num,$mod); # this does not make it faster
2398 my $acc = _copy($c,$num); my $t = _one();
2400 my $expbin = _as_bin($c,$exp); $expbin =~ s/^0b//;
2401 my $len = length($expbin);
2404 if ( substr($expbin,$len,1) eq '1') # is_odd
2407 $t = _mod($c,$t,$mod);
2410 $acc = _mod($c,$acc,$mod);
2418 # greatest common divisor
2421 while ( (scalar @$y != 1) || ($y->[0] != 0) ) # while ($y != 0)
2423 my $t = _copy($c,$y);
2424 $y = _mod($c, $x, $y);
2430 ##############################################################################
2431 ##############################################################################
2438 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
2442 Provides support for big integer calculations. Not intended to be used by other
2443 modules. Other modules which sport the same functions can also be used to support
2444 Math::BigInt, like Math::BigInt::GMP or Math::BigInt::Pari.
2448 In order to allow for multiple big integer libraries, Math::BigInt was
2449 rewritten to use library modules for core math routines. Any module which
2450 follows the same API as this can be used instead by using the following:
2452 use Math::BigInt lib => 'libname';
2454 'libname' is either the long name ('Math::BigInt::Pari'), or only the short
2455 version like 'Pari'.
2461 The following functions MUST be defined in order to support the use by
2462 Math::BigInt v1.70 or later:
2464 api_version() return API version, 1 for v1.70, 2 for v1.83
2465 _new(string) return ref to new object from ref to decimal string
2466 _zero() return a new object with value 0
2467 _one() return a new object with value 1
2468 _two() return a new object with value 2
2469 _ten() return a new object with value 10
2471 _str(obj) return ref to a string representing the object
2472 _num(obj) returns a Perl integer/floating point number
2473 NOTE: because of Perl numeric notation defaults,
2474 the _num'ified obj may lose accuracy due to
2475 machine-dependent floating point size limitations
2477 _add(obj,obj) Simple addition of two objects
2478 _mul(obj,obj) Multiplication of two objects
2479 _div(obj,obj) Division of the 1st object by the 2nd
2480 In list context, returns (result,remainder).
2481 NOTE: this is integer math, so no
2482 fractional part will be returned.
2483 The second operand will be not be 0, so no need to
2485 _sub(obj,obj) Simple subtraction of 1 object from another
2486 a third, optional parameter indicates that the params
2487 are swapped. In this case, the first param needs to
2488 be preserved, while you can destroy the second.
2489 sub (x,y,1) => return x - y and keep x intact!
2490 _dec(obj) decrement object by one (input is guaranteed to be > 0)
2491 _inc(obj) increment object by one
2494 _acmp(obj,obj) <=> operator for objects (return -1, 0 or 1)
2496 _len(obj) returns count of the decimal digits of the object
2497 _digit(obj,n) returns the n'th decimal digit of object
2499 _is_one(obj) return true if argument is 1
2500 _is_two(obj) return true if argument is 2
2501 _is_ten(obj) return true if argument is 10
2502 _is_zero(obj) return true if argument is 0
2503 _is_even(obj) return true if argument is even (0,2,4,6..)
2504 _is_odd(obj) return true if argument is odd (1,3,5,7..)
2506 _copy return a ref to a true copy of the object
2508 _check(obj) check whether internal representation is still intact
2509 return 0 for ok, otherwise error message as string
2511 _from_hex(str) return new object from a hexadecimal string
2512 _from_bin(str) return new object from a binary string
2513 _from_oct(str) return new object from an octal string
2515 _as_hex(str) return string containing the value as
2516 unsigned hex string, with the '0x' prepended.
2517 Leading zeros must be stripped.
2518 _as_bin(str) Like as_hex, only as binary string containing only
2519 zeros and ones. Leading zeros must be stripped and a
2520 '0b' must be prepended.
2522 _rsft(obj,N,B) shift object in base B by N 'digits' right
2523 _lsft(obj,N,B) shift object in base B by N 'digits' left
2525 _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
2526 Note: XOR, AND and OR pad with zeros if size mismatches
2527 _and(obj1,obj2) AND (bit-wise) object 1 with object 2
2528 _or(obj1,obj2) OR (bit-wise) object 1 with object 2
2530 _mod(obj1,obj2) Return remainder of div of the 1st by the 2nd object
2531 _sqrt(obj) return the square root of object (truncated to int)
2532 _root(obj) return the n'th (n >= 3) root of obj (truncated to int)
2533 _fac(obj) return factorial of object 1 (1*2*3*4..)
2534 _pow(obj1,obj2) return object 1 to the power of object 2
2535 return undef for NaN
2536 _zeros(obj) return number of trailing decimal zeros
2537 _modinv return inverse modulus
2538 _modpow return modulus of power ($x ** $y) % $z
2539 _log_int(X,N) calculate integer log() of X in base N
2540 X >= 0, N >= 0 (return undef for NaN)
2541 returns (RESULT, EXACT) where EXACT is:
2542 1 : result is exactly RESULT
2543 0 : result was truncated to RESULT
2544 undef : unknown whether result is exactly RESULT
2545 _gcd(obj,obj) return Greatest Common Divisor of two objects
2547 The following functions are REQUIRED for an api_version of 2 or greater:
2549 _1ex($x) create the number 1Ex where x >= 0
2550 _alen(obj) returns approximate count of the decimal digits of the
2551 object. This estimate MUST always be greater or equal
2552 to what _len() returns.
2553 _nok(n,k) calculate n over k (binomial coefficient)
2555 The following functions are optional, and can be defined if the underlying lib
2556 has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
2557 slow) fallback routines to emulate these:
2563 Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
2566 So the library needs only to deal with unsigned big integers. Testing of input
2567 parameter validity is done by the caller, so you need not worry about
2568 underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by zero or similar
2571 The first parameter can be modified, that includes the possibility that you
2572 return a reference to a completely different object instead. Although keeping
2573 the reference and just changing its contents is preferred over creating and
2574 returning a different reference.
2576 Return values are always references to objects, strings, or true/false for
2577 comparison routines.
2579 =head1 WRAP YOUR OWN
2581 If you want to port your own favourite c-lib for big numbers to the
2582 Math::BigInt interface, you can take any of the already existing modules as
2583 a rough guideline. You should really wrap up the latest BigInt and BigFloat
2584 testsuites with your module, and replace in them any of the following:
2590 use Math::BigInt lib => 'yourlib';
2592 This way you ensure that your library really works 100% within Math::BigInt.
2596 This program is free software; you may redistribute it and/or modify it under
2597 the same terms as Perl itself.
2601 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
2603 Seperated from BigInt and shaped API with the help of John Peacock.
2605 Fixed, speed-up, streamlined and enhanced by Tels 2001 - 2007.
2609 L<Math::BigInt>, L<Math::BigFloat>,
2610 L<Math::BigInt::GMP>, L<Math::BigInt::FastCalc> and L<Math::BigInt::Pari>.