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, $BASE, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN, $MAX_VAL,);
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, $BASE, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN, $MAX_VAL);
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 odd
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! 1*2*3*4 * 5*6*7 5 * 6 * 7 6 7
1540 # ( - ) = --------- = --------------- = --------- = 5 * - * -
1541 # ( 3 ) (7-3)! 3! 1*2*3*4 * 1*2*3 1 * 2 * 3 2 3
1543 if (!_is_zero($c,$k))
1545 my $x = _copy($c,$n);
1548 my $f = _copy($c,$n); _inc($c,$f); # n = 5, f = 6, d = 2
1550 while (_acmp($c,$f,$x) <= 0) # f <= n ?
1552 # n = (n * f / d) == 5 * 6 / 2
1553 $n = _mul($c,$n,$f); $n = _div($c,$n,$d);
1555 _inc($c,$f); _inc($c,$d);
1560 # keep ref to $n and set it to 1
1561 splice (@$n,1); $n->[0] = 1;
1580 # ref to array, return ref to array
1583 if ((@$cx == 1) && ($cx->[0] <= 7))
1585 $cx->[0] = $factorials[$cx->[0]]; # 0 => 1, 1 => 1, 2 => 2 etc.
1589 if ((@$cx == 1) && # we do this only if $x >= 12 and $x <= 7000
1590 ($cx->[0] >= 12 && $cx->[0] < 7000))
1593 # Calculate (k-j) * (k-j+1) ... k .. (k+j-1) * (k + j)
1594 # See http://blogten.blogspot.com/2007/01/calculating-n.html
1595 # The above series can be expressed as factors:
1596 # k * k - (j - i) * 2
1597 # We cache k*k, and calculate (j * j) as the sum of the first j odd integers
1599 # This will not work when N exceeds the storage of a Perl scalar, however,
1600 # in this case the algorithm would be way to slow to terminate, anyway.
1602 # As soon as the last element of $cx is 0, we split it up and remember
1603 # how many zeors we got so far. The reason is that n! will accumulate
1604 # zeros at the end rather fast.
1605 my $zero_elements = 0;
1607 # If n is even, set n = n -1
1608 my $k = _num($c,$cx); my $even = 1;
1613 # set k to the center point
1615 # print "k $k even: $even\n";
1616 # now calculate k * k
1618 my $odd = 1; my $sum = 1;
1620 # keep reference to x
1621 my $new_x = _new($c, $k * $even);
1625 $zero_elements ++; shift @$cx;
1627 # print STDERR "x = ", _str($c,$cx),"\n";
1628 my $BASE2 = int(sqrt($BASE))-1;
1632 my $m = ($k2 - $sum); $odd += 2; $sum += $odd; $j++;
1633 while ($j <= $i && ($m < $BASE2) && (($k2 - $sum) < $BASE2))
1636 $odd += 2; $sum += $odd; $j++;
1637 # print STDERR "\n k2 $k2 m $m sum $sum odd $odd\n"; sleep(1);
1645 _mul($c,$cx,$c->_new($m));
1649 $zero_elements ++; shift @$cx;
1651 # print STDERR "Calculate $k2 - $sum = $m (x = ", _str($c,$cx),")\n";
1653 # multiply in the zeros again
1654 unshift @$cx, (0) x $zero_elements;
1658 # go forward until $base is exceeded
1659 # limit is either $x steps (steps == 100 means a result always too high) or
1661 my $steps = 100; $steps = $cx->[0] if @$cx == 1;
1662 my $r = 2; my $cf = 3; my $step = 2; my $last = $r;
1663 while ($r*$cf < $BASE && $step < $steps)
1665 $last = $r; $r *= $cf++; $step++;
1667 if ((@$cx == 1) && $step == $cx->[0])
1669 # completely done, so keep reference to $x and return
1674 # now we must do the left over steps
1675 my $n; # steps still to do
1676 if (scalar @$cx == 1)
1685 # Set $cx to the last result below $BASE (but keep ref to $x)
1686 $cx->[0] = $last; splice (@$cx,1);
1687 # As soon as the last element of $cx is 0, we split it up and remember
1688 # how many zeors we got so far. The reason is that n! will accumulate
1689 # zeros at the end rather fast.
1690 my $zero_elements = 0;
1692 # do left-over steps fit into a scalar?
1693 if (ref $n eq 'ARRAY')
1695 # No, so use slower inc() & cmp()
1696 # ($n is at least $BASE here)
1697 my $base_2 = int(sqrt($BASE)) - 1;
1698 #print STDERR "base_2: $base_2\n";
1699 while ($step < $base_2)
1703 $zero_elements ++; shift @$cx;
1705 my $b = $step * ($step + 1); $step += 2;
1709 while (_acmp($c,$step,$n) <= 0)
1713 $zero_elements ++; shift @$cx;
1715 _mul($c,$cx,$step); _inc($c,$step);
1720 # Yes, so we can speed it up slightly
1722 # print "# left over steps $n\n";
1724 my $base_4 = int(sqrt(sqrt($BASE))) - 2;
1725 #print STDERR "base_4: $base_4\n";
1727 while ($step < $n4 && $step < $base_4)
1731 $zero_elements ++; shift @$cx;
1733 my $b = $step * ($step + 1); $step += 2; $b *= $step * ($step + 1); $step += 2;
1736 my $base_2 = int(sqrt($BASE)) - 1;
1738 #print STDERR "base_2: $base_2\n";
1739 while ($step < $n2 && $step < $base_2)
1743 $zero_elements ++; shift @$cx;
1745 my $b = $step * ($step + 1); $step += 2;
1748 # do what's left over
1751 _mul($c,$cx,[$step]); $step++;
1754 $zero_elements ++; shift @$cx;
1758 # multiply in the zeros again
1759 unshift @$cx, (0) x $zero_elements;
1760 $cx; # return result
1763 #############################################################################
1767 # calculate integer log of $x to base $base
1768 # ref to array, ref to array - return ref to array
1769 my ($c,$x,$base) = @_;
1772 return if (scalar @$x == 1 && $x->[0] == 0);
1773 # BASE 0 or 1 => NaN
1774 return if (scalar @$base == 1 && $base->[0] < 2);
1775 my $cmp = _acmp($c,$x,$base); # X == BASE => 1
1778 splice (@$x,1); $x->[0] = 1;
1784 splice (@$x,1); $x->[0] = 0;
1788 my $x_org = _copy($c,$x); # preserve x
1789 splice(@$x,1); $x->[0] = 1; # keep ref to $x
1791 # Compute a guess for the result based on:
1792 # $guess = int ( length_in_base_10(X) / ( log(base) / log(10) ) )
1793 my $len = _len($c,$x_org);
1794 my $log = log($base->[-1]) / log(10);
1796 # for each additional element in $base, we add $BASE_LEN to the result,
1797 # based on the observation that log($BASE,10) is BASE_LEN and
1798 # log(x*y) == log(x) + log(y):
1799 $log += ((scalar @$base)-1) * $BASE_LEN;
1801 # calculate now a guess based on the values obtained above:
1802 my $res = int($len / $log);
1805 my $trial = _pow ($c, _copy($c, $base), $x);
1806 my $a = _acmp($c,$trial,$x_org);
1808 # print STDERR "# trial ", _str($c,$x)," was: $a (0 = exact, -1 too small, +1 too big)\n";
1810 # found an exact result?
1811 return ($x,1) if $a == 0;
1816 _div($c,$trial,$base); _dec($c, $x);
1817 while (($a = _acmp($c,$trial,$x_org)) > 0)
1819 # print STDERR "# big _log_int at ", _str($c,$x), "\n";
1820 _div($c,$trial,$base); _dec($c, $x);
1822 # result is now exact (a == 0), or too small (a < 0)
1823 return ($x, $a == 0 ? 1 : 0);
1826 # else: result was to small
1827 _mul($c,$trial,$base);
1829 # did we now get the right result?
1830 $a = _acmp($c,$trial,$x_org);
1832 if ($a == 0) # yes, exactly
1837 return ($x,0) if $a > 0;
1839 # Result still too small (we should come here only if the estimate above
1840 # was very off base):
1842 # Now let the normal trial run obtain the real result
1843 # Simple loop that increments $x by 2 in each step, possible overstepping
1846 my $base_mul = _mul($c, _copy($c,$base), $base); # $base * $base
1848 while (($a = _acmp($c,$trial,$x_org)) < 0)
1850 # print STDERR "# small _log_int at ", _str($c,$x), "\n";
1851 _mul($c,$trial,$base_mul); _add($c, $x, [2]);
1857 # overstepped the result
1859 _div($c,$trial,$base);
1860 $a = _acmp($c,$trial,$x_org);
1865 $exact = 0 if $a != 0; # a = -1 => not exact result, a = 0 => exact
1868 ($x,$exact); # return result
1872 use constant DEBUG => 0;
1874 sub steps { $steps };
1878 # square-root of $x in place
1879 # Compute a guess of the result (by rule of thumb), then improve it via
1883 if (scalar @$x == 1)
1885 # fits into one Perl scalar, so result can be computed directly
1886 $x->[0] = int(sqrt($x->[0]));
1889 my $y = _copy($c,$x);
1890 # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
1891 # since our guess will "grow"
1892 my $l = int((_len($c,$x)-1) / 2);
1894 my $lastelem = $x->[-1]; # for guess
1895 my $elems = scalar @$x - 1;
1896 # not enough digits, but could have more?
1897 if ((length($lastelem) <= 3) && ($elems > 1))
1899 # right-align with zero pad
1900 my $len = length($lastelem) & 1;
1901 print "$lastelem => " if DEBUG;
1902 $lastelem .= substr($x->[-2] . '0' x $BASE_LEN,0,$BASE_LEN);
1903 # former odd => make odd again, or former even to even again
1904 $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
1905 print "$lastelem\n" if DEBUG;
1908 # construct $x (instead of _lsft($c,$x,$l,10)
1909 my $r = $l % $BASE_LEN; # 10000 00000 00000 00000 ($BASE_LEN=5)
1910 $l = int($l / $BASE_LEN);
1911 print "l = $l " if DEBUG;
1913 splice @$x,$l; # keep ref($x), but modify it
1915 # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
1917 # 14400 00000 => sqrt(14400) => guess first digits to be 120
1918 # 144000 000000 => sqrt(144000) => guess 379
1920 print "$lastelem (elems $elems) => " if DEBUG;
1921 $lastelem = $lastelem / 10 if ($elems & 1 == 1); # odd or even?
1922 my $g = sqrt($lastelem); $g =~ s/\.//; # 2.345 => 2345
1923 $r -= 1 if $elems & 1 == 0; # 70 => 7
1925 # padd with zeros if result is too short
1926 $x->[$l--] = int(substr($g . '0' x $r,0,$r+1));
1927 print "now ",$x->[-1] if DEBUG;
1928 print " would have been ", int('1' . '0' x $r),"\n" if DEBUG;
1930 # If @$x > 1, we could compute the second elem of the guess, too, to create
1931 # an even better guess. Not implemented yet. Does it improve performance?
1932 $x->[$l--] = 0 while ($l >= 0); # all other digits of guess are zero
1934 print "start x= ",_str($c,$x),"\n" if DEBUG;
1937 my $lastlast = _zero();
1938 $steps = 0 if DEBUG;
1939 while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
1942 $lastlast = _copy($c,$last);
1943 $last = _copy($c,$x);
1944 _add($c,$x, _div($c,_copy($c,$y),$x));
1946 print " x= ",_str($c,$x),"\n" if DEBUG;
1948 print "\nsteps in sqrt: $steps, " if DEBUG;
1949 _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0; # overshot?
1950 print " final ",$x->[-1],"\n" if DEBUG;
1956 # take n'th root of $x in place (n >= 3)
1959 if (scalar @$x == 1)
1963 # result will always be smaller than 2 so trunc to 1 at once
1968 # fits into one Perl scalar, so result can be computed directly
1969 # cannot use int() here, because it rounds wrongly (try
1970 # (81 ** 3) ** (1/3) to see what I mean)
1971 #$x->[0] = int( $x->[0] ** (1 / $n->[0]) );
1972 # round to 8 digits, then truncate result to integer
1973 $x->[0] = int ( sprintf ("%.8f", $x->[0] ** (1 / $n->[0]) ) );
1978 # we know now that X is more than one element long
1980 # if $n is a power of two, we can repeatedly take sqrt($X) and find the
1981 # proper result, because sqrt(sqrt($x)) == root($x,4)
1982 my $b = _as_bin($c,$n);
1983 if ($b =~ /0b1(0+)$/)
1985 my $count = CORE::length($1); # 0b100 => len('00') => 2
1986 my $cnt = $count; # counter for loop
1987 unshift (@$x, 0); # add one element, together with one
1988 # more below in the loop this makes 2
1991 # 'inflate' $X by adding one element, basically computing
1992 # $x * $BASE * $BASE. This gives us more $BASE_LEN digits for result
1993 # since len(sqrt($X)) approx == len($x) / 2.
1995 # calculate sqrt($x), $x is now one element to big, again. In the next
1996 # round we make that two, again.
1999 # $x is now one element to big, so truncate result by removing it
2004 # trial computation by starting with 2,4,8,16 etc until we overstep
2008 # while still to do more than X steps
2012 while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
2014 _mul ($c, $step, [2]);
2015 _add ($c, $trial, $step);
2019 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) == 0)
2021 @$x = @$trial; # make copy while preserving ref to $x
2024 # overstepped, so go back on step
2025 _sub($c, $trial, $step);
2026 } while (scalar @$step > 1 || $step->[0] > 128);
2030 # add two, because $trial cannot be exactly the result (otherwise we would
2031 # alrady have found it)
2032 _add($c, $trial, $step);
2034 # and now add more and more (2,4,6,8,10 etc)
2035 while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
2037 _add ($c, $trial, $step);
2040 # hit not exactly? (overstepped)
2041 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
2046 # hit not exactly? (overstepped)
2047 # 80 too small, 81 slightly too big, 82 too big
2048 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
2053 @$x = @$trial; # make copy while preserving ref to $x
2059 ##############################################################################
2066 # the shortcut makes equal, large numbers _really_ fast, and makes only a
2067 # very small performance drop for small numbers (e.g. something with less
2068 # than 32 bit) Since we optimize for large numbers, this is enabled.
2069 return $x if _acmp($c,$x,$y) == 0; # shortcut
2071 my $m = _one(); my ($xr,$yr);
2072 my $mask = $AND_MASK;
2075 my $y1 = _copy($c,$y); # make copy
2079 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2081 ($x1, $xr) = _div($c,$x1,$mask);
2082 ($y1, $yr) = _div($c,$y1,$mask);
2084 # make ints() from $xr, $yr
2085 # this is when the AND_BITS are greater than $BASE and is slower for
2086 # small (<256 bits) numbers, but faster for large numbers. Disabled
2087 # due to KISS principle
2089 # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2090 # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2091 # _add($c,$x, _mul($c, _new( $c, ($xrr & $yrr) ), $m) );
2093 # 0+ due to '&' doesn't work in strings
2094 _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
2104 return _zero() if _acmp($c,$x,$y) == 0; # shortcut (see -and)
2106 my $m = _one(); my ($xr,$yr);
2107 my $mask = $XOR_MASK;
2110 my $y1 = _copy($c,$y); # make copy
2114 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2116 ($x1, $xr) = _div($c,$x1,$mask);
2117 ($y1, $yr) = _div($c,$y1,$mask);
2118 # make ints() from $xr, $yr (see _and())
2119 #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2120 #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2121 #_add($c,$x, _mul($c, _new( $c, ($xrr ^ $yrr) ), $m) );
2123 # 0+ due to '^' doesn't work in strings
2124 _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
2127 # the loop stops when the shorter of the two numbers is exhausted
2128 # the remainder of the longer one will survive bit-by-bit, so we simple
2129 # multiply-add it in
2130 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
2131 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
2140 return $x if _acmp($c,$x,$y) == 0; # shortcut (see _and)
2142 my $m = _one(); my ($xr,$yr);
2143 my $mask = $OR_MASK;
2146 my $y1 = _copy($c,$y); # make copy
2150 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2152 ($x1, $xr) = _div($c,$x1,$mask);
2153 ($y1, $yr) = _div($c,$y1,$mask);
2154 # make ints() from $xr, $yr (see _and())
2155 # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2156 # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2157 # _add($c,$x, _mul($c, _new( $c, ($xrr | $yrr) ), $m) );
2159 # 0+ due to '|' doesn't work in strings
2160 _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
2163 # the loop stops when the shorter of the two numbers is exhausted
2164 # the remainder of the longer one will survive bit-by-bit, so we simple
2165 # multiply-add it in
2166 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
2167 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
2174 # convert a decimal number to hex (ref to array, return ref to string)
2177 # fits into one element (handle also 0x0 case)
2178 return sprintf("0x%x",$x->[0]) if @$x == 1;
2180 my $x1 = _copy($c,$x);
2183 my ($xr, $h, $x10000);
2186 $x10000 = [ 0x10000 ]; $h = 'h4';
2190 $x10000 = [ 0x1000 ]; $h = 'h3';
2192 while (@$x1 != 1 || $x1->[0] != 0) # _is_zero()
2194 ($x1, $xr) = _div($c,$x1,$x10000);
2195 $es .= unpack($h,pack('V',$xr->[0]));
2198 $es =~ s/^[0]+//; # strip leading zeros
2199 '0x' . $es; # return result prepended with 0x
2204 # convert a decimal number to bin (ref to array, return ref to string)
2207 # fits into one element (and Perl recent enough), handle also 0b0 case
2208 # handle zero case for older Perls
2209 if ($] <= 5.005 && @$x == 1 && $x->[0] == 0)
2211 my $t = '0b0'; return $t;
2213 if (@$x == 1 && $] >= 5.006)
2215 my $t = sprintf("0b%b",$x->[0]);
2218 my $x1 = _copy($c,$x);
2221 my ($xr, $b, $x10000);
2224 $x10000 = [ 0x10000 ]; $b = 'b16';
2228 $x10000 = [ 0x1000 ]; $b = 'b12';
2230 while (!(@$x1 == 1 && $x1->[0] == 0)) # _is_zero()
2232 ($x1, $xr) = _div($c,$x1,$x10000);
2233 $es .= unpack($b,pack('v',$xr->[0]));
2236 $es =~ s/^[0]+//; # strip leading zeros
2237 '0b' . $es; # return result prepended with 0b
2242 # convert a decimal number to octal (ref to array, return ref to string)
2245 # fits into one element (handle also 0 case)
2246 return sprintf("0%o",$x->[0]) if @$x == 1;
2248 my $x1 = _copy($c,$x);
2252 my $x1000 = [ 0100000 ];
2253 while (@$x1 != 1 || $x1->[0] != 0) # _is_zero()
2255 ($x1, $xr) = _div($c,$x1,$x1000);
2256 $es .= reverse sprintf("%05o", $xr->[0]);
2259 $es =~ s/^[0]+//; # strip leading zeros
2260 '0' . $es; # return result prepended with 0
2265 # convert a octal number to decimal (string, return ref to array)
2268 # for older Perls, play safe
2269 my $m = [ 0100000 ];
2270 my $d = 5; # 5 digits at a time
2275 my $len = int( (length($os)-1)/$d ); # $d digit parts, w/o the '0'
2276 my $val; my $i = -$d;
2279 $val = substr($os,$i,$d); # get oct digits
2280 $val = CORE::oct($val);
2282 my $adder = [ $val ];
2283 _add ($c, $x, _mul ($c, $adder, $mul ) ) if $val != 0;
2284 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
2291 # convert a hex number to decimal (string, return ref to array)
2294 my $m = _new($c, 0x10000000); # 28 bit at a time (<32 bit!)
2295 my $d = 7; # 7 digits at a time
2298 # for older Perls, play safe
2299 $m = [ 0x10000 ]; # 16 bit at a time (<32 bit!)
2300 $d = 4; # 4 digits at a time
2306 my $len = int( (length($hs)-2)/$d ); # $d digit parts, w/o the '0x'
2307 my $val; my $i = -$d;
2310 $val = substr($hs,$i,$d); # get hex digits
2311 $val =~ s/^0x// if $len == 0; # for last part only because
2312 $val = CORE::hex($val); # hex does not like wrong chars
2314 my $adder = [ $val ];
2315 # if the resulting number was to big to fit into one element, create a
2316 # two-element version (bug found by Mark Lakata - Thanx!)
2317 if (CORE::length($val) > $BASE_LEN)
2319 $adder = _new($c,$val);
2321 _add ($c, $x, _mul ($c, $adder, $mul ) ) if $val != 0;
2322 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
2329 # convert a hex number to decimal (string, return ref to array)
2332 # instead of converting X (8) bit at a time, it is faster to "convert" the
2333 # number to hex, and then call _from_hex.
2336 $hs =~ s/^[+-]?0b//; # remove sign and 0b
2337 my $l = length($hs); # bits
2338 $hs = '0' x (8-($l % 8)) . $hs if ($l % 8) != 0; # padd left side w/ 0
2339 my $h = '0x' . unpack('H*', pack ('B*', $hs)); # repack as hex
2344 ##############################################################################
2345 # special modulus functions
2352 my $u = _zero($c); my $u1 = _one($c);
2353 my $a = _copy($c,$y); my $b = _copy($c,$x);
2355 # Euclid's Algorithm for bgcd(), only that we calc bgcd() ($a) and the
2356 # result ($u) at the same time. See comments in BigInt for why this works.
2358 ($a, $q, $b) = ($b, _div($c,$a,$b)); # step 1
2360 while (!_is_zero($c,$b))
2362 my $t = _add($c, # step 2:
2363 _mul($c,_copy($c,$u1), $q) , # t = u1 * q
2365 $u = $u1; # u = u1, u1 = t
2368 ($a, $q, $b) = ($b, _div($c,$a,$b)); # step 1
2371 # if the gcd is not 1, then return NaN
2372 return (undef,undef) unless _is_one($c,$a);
2374 ($u1, $sign == 1 ? '+' : '-');
2379 # modulus of power ($x ** $y) % $z
2380 my ($c,$num,$exp,$mod) = @_;
2382 # in the trivial case,
2383 if (_is_one($c,$mod))
2385 splice @$num,0,1; $num->[0] = 0;
2388 if ((scalar @$num == 1) && (($num->[0] == 0) || ($num->[0] == 1)))
2394 # $num = _mod($c,$num,$mod); # this does not make it faster
2396 my $acc = _copy($c,$num); my $t = _one();
2398 my $expbin = _as_bin($c,$exp); $expbin =~ s/^0b//;
2399 my $len = length($expbin);
2402 if ( substr($expbin,$len,1) eq '1') # is_odd
2405 $t = _mod($c,$t,$mod);
2408 $acc = _mod($c,$acc,$mod);
2416 # greatest common divisor
2419 while ( (scalar @$y != 1) || ($y->[0] != 0) ) # while ($y != 0)
2421 my $t = _copy($c,$y);
2422 $y = _mod($c, $x, $y);
2428 ##############################################################################
2429 ##############################################################################
2436 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
2440 Provides support for big integer calculations. Not intended to be used by other
2441 modules. Other modules which sport the same functions can also be used to support
2442 Math::BigInt, like Math::BigInt::GMP or Math::BigInt::Pari.
2446 In order to allow for multiple big integer libraries, Math::BigInt was
2447 rewritten to use library modules for core math routines. Any module which
2448 follows the same API as this can be used instead by using the following:
2450 use Math::BigInt lib => 'libname';
2452 'libname' is either the long name ('Math::BigInt::Pari'), or only the short
2453 version like 'Pari'.
2459 The following functions MUST be defined in order to support the use by
2460 Math::BigInt v1.70 or later:
2462 api_version() return API version, 1 for v1.70, 2 for v1.83
2463 _new(string) return ref to new object from ref to decimal string
2464 _zero() return a new object with value 0
2465 _one() return a new object with value 1
2466 _two() return a new object with value 2
2467 _ten() return a new object with value 10
2469 _str(obj) return ref to a string representing the object
2470 _num(obj) returns a Perl integer/floating point number
2471 NOTE: because of Perl numeric notation defaults,
2472 the _num'ified obj may lose accuracy due to
2473 machine-dependent floating point size limitations
2475 _add(obj,obj) Simple addition of two objects
2476 _mul(obj,obj) Multiplication of two objects
2477 _div(obj,obj) Division of the 1st object by the 2nd
2478 In list context, returns (result,remainder).
2479 NOTE: this is integer math, so no
2480 fractional part will be returned.
2481 The second operand will be not be 0, so no need to
2483 _sub(obj,obj) Simple subtraction of 1 object from another
2484 a third, optional parameter indicates that the params
2485 are swapped. In this case, the first param needs to
2486 be preserved, while you can destroy the second.
2487 sub (x,y,1) => return x - y and keep x intact!
2488 _dec(obj) decrement object by one (input is guaranteed to be > 0)
2489 _inc(obj) increment object by one
2492 _acmp(obj,obj) <=> operator for objects (return -1, 0 or 1)
2494 _len(obj) returns count of the decimal digits of the object
2495 _digit(obj,n) returns the n'th decimal digit of object
2497 _is_one(obj) return true if argument is 1
2498 _is_two(obj) return true if argument is 2
2499 _is_ten(obj) return true if argument is 10
2500 _is_zero(obj) return true if argument is 0
2501 _is_even(obj) return true if argument is even (0,2,4,6..)
2502 _is_odd(obj) return true if argument is odd (1,3,5,7..)
2504 _copy return a ref to a true copy of the object
2506 _check(obj) check whether internal representation is still intact
2507 return 0 for ok, otherwise error message as string
2509 _from_hex(str) return new object from a hexadecimal string
2510 _from_bin(str) return new object from a binary string
2511 _from_oct(str) return new object from an octal string
2513 _as_hex(str) return string containing the value as
2514 unsigned hex string, with the '0x' prepended.
2515 Leading zeros must be stripped.
2516 _as_bin(str) Like as_hex, only as binary string containing only
2517 zeros and ones. Leading zeros must be stripped and a
2518 '0b' must be prepended.
2520 _rsft(obj,N,B) shift object in base B by N 'digits' right
2521 _lsft(obj,N,B) shift object in base B by N 'digits' left
2523 _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
2524 Note: XOR, AND and OR pad with zeros if size mismatches
2525 _and(obj1,obj2) AND (bit-wise) object 1 with object 2
2526 _or(obj1,obj2) OR (bit-wise) object 1 with object 2
2528 _mod(obj1,obj2) Return remainder of div of the 1st by the 2nd object
2529 _sqrt(obj) return the square root of object (truncated to int)
2530 _root(obj) return the n'th (n >= 3) root of obj (truncated to int)
2531 _fac(obj) return factorial of object 1 (1*2*3*4..)
2532 _pow(obj1,obj2) return object 1 to the power of object 2
2533 return undef for NaN
2534 _zeros(obj) return number of trailing decimal zeros
2535 _modinv return inverse modulus
2536 _modpow return modulus of power ($x ** $y) % $z
2537 _log_int(X,N) calculate integer log() of X in base N
2538 X >= 0, N >= 0 (return undef for NaN)
2539 returns (RESULT, EXACT) where EXACT is:
2540 1 : result is exactly RESULT
2541 0 : result was truncated to RESULT
2542 undef : unknown whether result is exactly RESULT
2543 _gcd(obj,obj) return Greatest Common Divisor of two objects
2545 The following functions are REQUIRED for an api_version of 2 or greater:
2547 _1ex($x) create the number 1Ex where x >= 0
2548 _alen(obj) returns approximate count of the decimal digits of the
2549 object. This estimate MUST always be greater or equal
2550 to what _len() returns.
2551 _nok(n,k) calculate n over k (binomial coefficient)
2553 The following functions are optional, and can be defined if the underlying lib
2554 has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
2555 slow) fallback routines to emulate these:
2561 Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
2564 So the library needs only to deal with unsigned big integers. Testing of input
2565 parameter validity is done by the caller, so you need not worry about
2566 underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by zero or similar
2569 The first parameter can be modified, that includes the possibility that you
2570 return a reference to a completely different object instead. Although keeping
2571 the reference and just changing its contents is preferred over creating and
2572 returning a different reference.
2574 Return values are always references to objects, strings, or true/false for
2575 comparison routines.
2577 =head1 WRAP YOUR OWN
2579 If you want to port your own favourite c-lib for big numbers to the
2580 Math::BigInt interface, you can take any of the already existing modules as
2581 a rough guideline. You should really wrap up the latest BigInt and BigFloat
2582 testsuites with your module, and replace in them any of the following:
2588 use Math::BigInt lib => 'yourlib';
2590 This way you ensure that your library really works 100% within Math::BigInt.
2594 This program is free software; you may redistribute it and/or modify it under
2595 the same terms as Perl itself.
2599 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
2601 Separated from BigInt and shaped API with the help of John Peacock.
2603 Fixed, speed-up, streamlined and enhanced by Tels 2001 - 2007.
2607 L<Math::BigInt>, L<Math::BigFloat>,
2608 L<Math::BigInt::GMP>, L<Math::BigInt::FastCalc> and L<Math::BigInt::Pari>.