1 package Math::BigInt::Calc;
5 # use warnings; # dont use warnings for older Perls
11 # Package to store unsigned big integers in decimal and do math with them
13 # Internally the numbers are stored in an array with at least 1 element, no
14 # leading zero parts (except the first) and in base 1eX where X is determined
15 # automatically at loading time to be the maximum possible value
18 # - fully remove funky $# stuff (maybe)
20 # USE_MUL: due to problems on certain os (os390, posix-bc) "* 1e-5" is used
21 # instead of "/ 1e5" at some places, (marked with USE_MUL). Other platforms
22 # BS2000, some Crays need USE_DIV instead.
23 # The BEGIN block is used to determine which of the two variants gives the
26 # Beware of things like:
27 # $i = $i * $y + $car; $car = int($i / $MBASE); $i = $i % $MBASE;
28 # This works on x86, but fails on ARM (SA1100, iPAQ) due to whoknows what
29 # reasons. So, use this instead (slower, but correct):
30 # $i = $i * $y + $car; $car = int($i / $MBASE); $i -= $MBASE * $car;
32 ##############################################################################
33 # global constants, flags and accessory
35 # constants for easier life
37 my ($MBASE,$BASE,$RBASE,$BASE_LEN,$MAX_VAL,$BASE_LEN2,$BASE_LEN_SMALL);
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 be the testsuite, set is used only by the BEGIN block below
50 # find whether we can use mul or div or none in mul()/div()
51 # (in last case reduce BASE_LEN_SMALL)
52 $BASE_LEN_SMALL = $b+1;
54 while (--$BASE_LEN_SMALL > 5)
56 $MBASE = int("1e".$BASE_LEN_SMALL);
57 $RBASE = abs('1e-'.$BASE_LEN_SMALL); # see USE_MUL
59 $caught += 1 if (int($MBASE * $RBASE) != 1); # should be 1
60 $caught += 2 if (int($MBASE / $MBASE) != 1); # should be 1
63 # BASE_LEN is used for anything else than mul()/div()
64 $BASE_LEN = $BASE_LEN_SMALL;
65 $BASE_LEN = shift if (defined $_[0]); # one more arg?
66 $BASE = int("1e".$BASE_LEN);
68 $BASE_LEN2 = int($BASE_LEN_SMALL / 2); # for mul shortcut
69 $MBASE = int("1e".$BASE_LEN_SMALL);
70 $RBASE = abs('1e-'.$BASE_LEN_SMALL); # see USE_MUL
73 #print "BASE_LEN: $BASE_LEN MAX_VAL: $MAX_VAL BASE: $BASE RBASE: $RBASE ";
74 #print "BASE_LEN_SMALL: $BASE_LEN_SMALL MBASE: $MBASE\n";
79 # $caught & 1 != 0 => cannot use MUL
80 # $caught & 2 != 0 => cannot use DIV
81 # The parens around ($caught & 1) were important, indeed, if we would use
85 # print "# use mul\n";
86 # must USE_MUL since we cannot use DIV
87 *{_mul} = \&_mul_use_mul;
88 *{_div} = \&_div_use_mul;
92 # print "# use div\n";
94 *{_mul} = \&_mul_use_div;
95 *{_div} = \&_div_use_div;
98 return $BASE_LEN unless wantarray;
99 return ($BASE_LEN, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN_SMALL, $MAX_VAL);
104 # from Daniel Pfeiffer: determine largest group of digits that is precisely
105 # multipliable with itself plus carry
106 # Test now changed to expect the proper pattern, not a result off by 1 or 2
107 my ($e, $num) = 3; # lowest value we will use is 3+1-1 = 3
110 $num = ('9' x ++$e) + 0;
112 } while ("$num" =~ /9{$e}0{$e}/); # must be a certain pattern
113 $e--; # last test failed, so retract one step
114 # the limits below brush the problems with the test above under the rug:
115 # the test should be able to find the proper $e automatically
116 $e = 5 if $^O =~ /^uts/; # UTS get's some special treatment
117 $e = 5 if $^O =~ /^unicos/; # unicos is also problematic (6 seems to work
118 # there, but we play safe)
119 $e = 5 if $] < 5.006; # cap, for older Perls
120 $e = 7 if $e > 7; # cap, for VMS, OS/390 and other 64 bit systems
121 # 8 fails inside random testsuite, so take 7
123 # determine how many digits fit into an integer and can be safely added
124 # together plus carry w/o causing an overflow
128 ############################################################################
129 # the next block is no longer important
131 ## this below detects 15 on a 64 bit system, because after that it becomes
132 ## 1e16 and not 1000000 :/ I can make it detect 18, but then I get a lot of
133 ## test failures. Ugh! (Tomake detect 18: uncomment lines marked with *)
135 #my $bi = 5; # approx. 16 bit
136 #$num = int('9' x $bi);
138 ## while ( ($num+$num+1) eq '1' . '9' x $bi) # *
139 #while ( int($num+$num+1) eq '1' . '9' x $bi)
141 # $bi++; $num = int('9' x $bi);
142 # # $bi++; $num *= 10; $num += 9; # *
144 #$bi--; # back off one step
145 # by setting them equal, we ignore the findings and use the default
146 # one-size-fits-all approach from former versions
147 my $bi = $e; # XXX, this should work always
149 __PACKAGE__->_base_len($e,$bi); # set and store
151 # find out how many bits _and, _or and _xor can take (old default = 16)
152 # I don't think anybody has yet 128 bit scalars, so let's play safe.
153 local $^W = 0; # don't warn about 'nonportable number'
154 $AND_BITS = 15; $XOR_BITS = 15; $OR_BITS = 15;
156 # find max bits, we will not go higher than numberofbits that fit into $BASE
157 # to make _and etc simpler (and faster for smaller, slower for large numbers)
159 while (2 ** $max < $BASE) { $max++; }
162 $max = 16 if $] < 5.006; # older Perls might not take >16 too well
167 $x = oct('0b' . '1' x $AND_BITS); $y = $x & $x;
168 $z = (2 ** $AND_BITS) - 1;
169 } while ($AND_BITS < $max && $x == $z && $y == $x);
170 $AND_BITS --; # retreat one step
173 $x = oct('0b' . '1' x $XOR_BITS); $y = $x ^ 0;
174 $z = (2 ** $XOR_BITS) - 1;
175 } while ($XOR_BITS < $max && $x == $z && $y == $x);
176 $XOR_BITS --; # retreat one step
179 $x = oct('0b' . '1' x $OR_BITS); $y = $x | $x;
180 $z = (2 ** $OR_BITS) - 1;
181 } while ($OR_BITS < $max && $x == $z && $y == $x);
182 $OR_BITS --; # retreat one step
186 ###############################################################################
190 # (ref to string) return ref to num_array
191 # Convert a number from string format (without sign) to internal base
192 # 1ex format. Assumes normalized value as input.
194 my $il = length($$d)-1;
196 # < BASE_LEN due len-1 above
197 return [ int($$d) ] if $il < $BASE_LEN; # shortcut for short numbers
199 # this leaves '00000' instead of int 0 and will be corrected after any op
200 [ reverse(unpack("a" . ($il % $BASE_LEN+1)
201 . ("a$BASE_LEN" x ($il / $BASE_LEN)), $$d)) ];
206 $AND_MASK = __PACKAGE__->_new( \( 2 ** $AND_BITS ));
207 $XOR_MASK = __PACKAGE__->_new( \( 2 ** $XOR_BITS ));
208 $OR_MASK = __PACKAGE__->_new( \( 2 ** $OR_BITS ));
225 # create a two (used internally for shifting)
235 # catch and throw away
238 ##############################################################################
239 # convert back to string and number
243 # (ref to BINT) return num_str
244 # Convert number from internal base 100000 format to string format.
245 # internal format is always normalized (no leading zeros, "-0" => "+0")
249 my $l = scalar @$ar; # number of parts
250 return $nan if $l < 1; # should not happen
252 # handle first one different to strip leading zeros from it (there are no
253 # leading zero parts in internal representation)
254 $l --; $ret .= int($ar->[$l]); $l--;
255 # Interestingly, the pre-padd method uses more time
256 # the old grep variant takes longer (14 vs. 10 sec)
257 my $z = '0' x ($BASE_LEN-1);
260 $ret .= substr($z.$ar->[$l],-$BASE_LEN); # fastest way I could think of
268 # Make a number (scalar int/float) from a BigInt object
270 return $x->[0] if scalar @$x == 1; # below $BASE
275 $num += $fac*$_; $fac *= $BASE;
280 ##############################################################################
285 # (ref to int_num_array, ref to int_num_array)
286 # routine to add two base 1eX numbers
287 # stolen from Knuth Vol 2 Algorithm A pg 231
288 # there are separate routines to add and sub as per Knuth pg 233
289 # This routine clobbers up array x, but not y.
293 return $x if (@$y == 1) && $y->[0] == 0; # $x + 0 => $x
294 if ((@$x == 1) && $x->[0] == 0) # 0 + $y => $y->copy
296 # twice as slow as $x = [ @$y ], but necc. to retain $x as ref :(
297 @$x = @$y; return $x;
300 # for each in Y, add Y to X and carry. If after that, something is left in
301 # X, foreach in X add carry to X and then return X, carry
302 # Trades one "$j++" for having to shift arrays
303 my $i; my $car = 0; my $j = 0;
306 $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
311 $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
318 # (ref to int_num_array, ref to int_num_array)
319 # Add 1 to $x, modify $x in place
324 return $x if (($i += 1) < $BASE); # early out
325 $i = 0; # overflow, next
327 push @$x,1 if ($x->[-1] == 0); # last overflowed, so extend
333 # (ref to int_num_array, ref to int_num_array)
334 # Sub 1 from $x, modify $x in place
337 my $MAX = $BASE-1; # since MAX_VAL based on MBASE
340 last if (($i -= 1) >= 0); # early out
341 $i = $MAX; # underflow, next
343 pop @$x if $x->[-1] == 0 && @$x > 1; # last underflowed (but leave 0)
349 # (ref to int_num_array, ref to int_num_array, swap)
350 # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
351 # subtract Y from X by modifying x in place
352 my ($c,$sx,$sy,$s) = @_;
354 my $car = 0; my $i; my $j = 0;
360 last unless defined $sy->[$j] || $car;
361 $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
363 # might leave leading zeros, so fix that
364 return __strip_zeros($sx);
366 #print "case 1 (swap)\n";
369 # we can't do an early out if $x is < than $y, since we
370 # need to copy the high chunks from $y. Found by Bob Mathews.
371 #last unless defined $sy->[$j] || $car;
373 if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
376 # might leave leading zeros, so fix that
382 # (ref to int_num_array, ref to int_num_array)
383 # multiply two numbers in internal representation
384 # modifies first arg, second need not be different from first
385 my ($c,$xv,$yv) = @_;
389 # shortcut for two very short numbers (improved by Nathan Zook)
390 # works also if xv and yv are the same reference, and handles also $x == 0
393 if (($xv->[0] *= $yv->[0]) >= $MBASE)
395 $xv->[0] = $xv->[0] - ($xv->[1] = int($xv->[0] * $RBASE)) * $MBASE;
405 # multiply a large number a by a single element one, so speed up
406 my $y = $yv->[0]; my $car = 0;
409 $i = $i * $y + $car; $car = int($i * $RBASE); $i -= $car * $MBASE;
411 push @$xv, $car if $car != 0;
414 # shortcut for result $x == 0 => result = 0
415 return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
417 # since multiplying $x with $x fails, make copy in this case
418 $yv = [@$xv] if $xv == $yv; # same references?
420 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
429 # $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
431 # $prod - ($car = int($prod * RBASE)) * $MBASE; # see USE_MUL
433 # $prod[$cty] += $car if $car; # need really to check for 0?
437 # looping through this if $xi == 0 is silly - so optimize it away!
438 $xi = (shift @prod || 0), next if $xi == 0;
441 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
442 ## this is actually a tad slower
443 ## $prod = $prod[$cty]; $prod += ($car + $xi * $yi); # no ||0 here
445 $prod - ($car = int($prod * $RBASE)) * $MBASE; # see USE_MUL
447 $prod[$cty] += $car if $car; # need really to check for 0?
448 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
457 # (ref to int_num_array, ref to int_num_array)
458 # multiply two numbers in internal representation
459 # modifies first arg, second need not be different from first
460 my ($c,$xv,$yv) = @_;
464 # shortcut for two small numbers, also handles $x == 0
467 # shortcut for two very short numbers (improved by Nathan Zook)
468 # works also if xv and yv are the same reference, and handles also $x == 0
469 if (($xv->[0] *= $yv->[0]) >= $MBASE)
472 $xv->[0] - ($xv->[1] = int($xv->[0] / $MBASE)) * $MBASE;
482 # multiply a large number a by a single element one, so speed up
483 my $y = $yv->[0]; my $car = 0;
486 $i = $i * $y + $car; $car = int($i / $MBASE); $i -= $car * $MBASE;
488 push @$xv, $car if $car != 0;
491 # shortcut for result $x == 0 => result = 0
492 return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
494 # since multiplying $x with $x fails, make copy in this case
495 $yv = [@$xv] if $xv == $yv; # same references?
497 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
501 # looping through this if $xi == 0 is silly - so optimize it away!
502 $xi = (shift @prod || 0), next if $xi == 0;
505 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
507 $prod - ($car = int($prod / $MBASE)) * $MBASE;
509 $prod[$cty] += $car if $car; # need really to check for 0?
510 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
519 # ref to array, ref to array, modify first array and return remainder if
522 # see comments in _div_use_div() for more explanations
524 my ($c,$x,$yorg) = @_;
526 # the general div algorithmn here is about O(N*N) and thus quite slow, so
527 # we first check for some special cases and use shortcuts to handle them.
529 # This works, because we store the numbers in a chunked format where each
530 # element contains 5..7 digits (depending on system).
532 # if both numbers have only one element:
533 if (@$x == 1 && @$yorg == 1)
535 # shortcut, $yorg and $x are two small numbers
538 my $r = [ $x->[0] % $yorg->[0] ];
539 $x->[0] = int($x->[0] / $yorg->[0]);
544 $x->[0] = int($x->[0] / $yorg->[0]);
549 # if x has more than one, but y has only one element:
553 $rem = _mod($c,[ @$x ],$yorg) if wantarray;
555 # shortcut, $y is < $BASE
556 my $j = scalar @$x; my $r = 0;
557 my $y = $yorg->[0]; my $b;
560 $b = $r * $MBASE + $x->[$j];
561 $x->[$j] = int($b/$y);
564 pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
565 return ($x,$rem) if wantarray;
569 # now x and y have more than one element
571 # check whether y has more elements than x, if yet, the result will be 0
575 $rem = [@$x] if wantarray; # make copy
576 splice (@$x,1); # keep ref to original array
577 $x->[0] = 0; # set to 0
578 return ($x,$rem) if wantarray; # including remainder?
579 return $x; # only x, which is [0] now
581 # check whether the numbers have the same number of elements, in that case
582 # the result will fit into one element and can be computed efficiently
586 # if $yorg has more digits than $x (it's leading element is longer than
587 # the one from $x), the result will also be 0:
588 if (length(int($yorg->[-1])) > length(int($x->[-1])))
590 $rem = [@$x] if wantarray; # make copy
591 splice (@$x,1); # keep ref to org array
592 $x->[0] = 0; # set to 0
593 return ($x,$rem) if wantarray; # including remainder?
596 # now calculate $x / $yorg
597 if (length(int($yorg->[-1])) == length(int($x->[-1])))
599 # same length, so make full compare, and if equal, return 1
600 # hm, same lengths, but same contents? So we need to check all parts:
601 my $a = 0; my $j = scalar @$x - 1;
602 # manual way (abort if unequal, good for early ne)
605 last if ($a = $x->[$j] - $yorg->[$j]); $j--;
607 # $a contains the result of the compare between X and Y
608 # a < 0: x < y, a == 0 => x == y, a > 0: x > y
613 $rem = [ 0 ]; # a = 0 => x == y => rem 1
614 $rem = [@$x] if $a != 0; # a < 0 => x < y => rem = x
616 splice(@$x,1); # keep single element
617 $x->[0] = 0; # if $a < 0
623 return ($x,$rem) if wantarray;
626 # $x >= $y, proceed normally
632 my $y = [ @$yorg ]; # always make copy to preserve
634 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
636 $car = $bar = $prd = 0;
637 if (($dd = int($MBASE/($y->[-1]+1))) != 1)
641 $xi = $xi * $dd + $car;
642 $xi -= ($car = int($xi * $RBASE)) * $MBASE; # see USE_MUL
644 push(@$x, $car); $car = 0;
647 $yi = $yi * $dd + $car;
648 $yi -= ($car = int($yi * $RBASE)) * $MBASE; # see USE_MUL
655 @q = (); ($v2,$v1) = @$y[-2,-1];
659 ($u2,$u1,$u0) = @$x[-3..-1];
661 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
663 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
664 --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
667 ($car, $bar) = (0,0);
668 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
670 $prd = $q * $y->[$yi] + $car;
671 $prd -= ($car = int($prd * $RBASE)) * $MBASE; # see USE_MUL
672 $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
674 if ($x->[-1] < $car + $bar)
677 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
680 if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $MBASE));
693 for $xi (reverse @$x)
695 $prd = $car * $MBASE + $xi;
696 $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
717 # ref to array, ref to array, modify first array and return remainder if
719 my ($c,$x,$yorg) = @_;
721 # the general div algorithmn here is about O(N*N) and thus quite slow, so
722 # we first check for some special cases and use shortcuts to handle them.
724 # This works, because we store the numbers in a chunked format where each
725 # element contains 5..7 digits (depending on system).
727 # if both numbers have only one element:
728 if (@$x == 1 && @$yorg == 1)
730 # shortcut, $yorg and $x are two small numbers
733 my $r = [ $x->[0] % $yorg->[0] ];
734 $x->[0] = int($x->[0] / $yorg->[0]);
739 $x->[0] = int($x->[0] / $yorg->[0]);
743 # if x has more than one, but y has only one element:
747 $rem = _mod($c,[ @$x ],$yorg) if wantarray;
749 # shortcut, $y is < $BASE
750 my $j = scalar @$x; my $r = 0;
751 my $y = $yorg->[0]; my $b;
754 $b = $r * $MBASE + $x->[$j];
755 $x->[$j] = int($b/$y);
758 pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
759 return ($x,$rem) if wantarray;
762 # now x and y have more than one element
764 # check whether y has more elements than x, if yet, the result will be 0
768 $rem = [@$x] if wantarray; # make copy
769 splice (@$x,1); # keep ref to original array
770 $x->[0] = 0; # set to 0
771 return ($x,$rem) if wantarray; # including remainder?
772 return $x; # only x, which is [0] now
774 # check whether the numbers have the same number of elements, in that case
775 # the result will fit into one element and can be computed efficiently
779 # if $yorg has more digits than $x (it's leading element is longer than
780 # the one from $x), the result will also be 0:
781 if (length(int($yorg->[-1])) > length(int($x->[-1])))
783 $rem = [@$x] if wantarray; # make copy
784 splice (@$x,1); # keep ref to org array
785 $x->[0] = 0; # set to 0
786 return ($x,$rem) if wantarray; # including remainder?
789 # now calculate $x / $yorg
791 if (length(int($yorg->[-1])) == length(int($x->[-1])))
793 # same length, so make full compare, and if equal, return 1
794 # hm, same lengths, but same contents? So we need to check all parts:
795 my $a = 0; my $j = scalar @$x - 1;
796 # manual way (abort if unequal, good for early ne)
799 last if ($a = $x->[$j] - $yorg->[$j]); $j--;
801 # $a contains the result of the compare between X and Y
802 # a < 0: x < y, a == 0 => x == y, a > 0: x > y
807 $rem = [ 0 ]; # a = 0 => x == y => rem 1
808 $rem = [@$x] if $a != 0; # a < 0 => x < y => rem = x
810 splice(@$x,1); # keep single element
811 $x->[0] = 0; # if $a < 0
817 return ($x,$rem) if wantarray;
820 # $x >= $y, so proceed normally
826 my $y = [ @$yorg ]; # always make copy to preserve
828 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
830 $car = $bar = $prd = 0;
831 if (($dd = int($MBASE/($y->[-1]+1))) != 1)
835 $xi = $xi * $dd + $car;
836 $xi -= ($car = int($xi / $MBASE)) * $MBASE;
838 push(@$x, $car); $car = 0;
841 $yi = $yi * $dd + $car;
842 $yi -= ($car = int($yi / $MBASE)) * $MBASE;
850 # @q will accumulate the final result, $q contains the current computed
851 # part of the final result
853 @q = (); ($v2,$v1) = @$y[-2,-1];
857 ($u2,$u1,$u0) = @$x[-3..-1];
859 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
861 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
862 --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
865 ($car, $bar) = (0,0);
866 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
868 $prd = $q * $y->[$yi] + $car;
869 $prd -= ($car = int($prd / $MBASE)) * $MBASE;
870 $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
872 if ($x->[-1] < $car + $bar)
875 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
878 if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $MBASE));
882 pop(@$x); unshift(@q, $q);
890 for $xi (reverse @$x)
892 $prd = $car * $MBASE + $xi;
893 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
912 ##############################################################################
917 # internal absolute post-normalized compare (ignore signs)
918 # ref to array, ref to array, return <0, 0, >0
919 # arrays must have at least one entry; this is not checked for
920 my ($c,$cx,$cy) = @_;
922 # shortcut for short numbers
923 return (($cx->[0] <=> $cy->[0]) <=> 0)
924 if scalar @$cx == scalar @$cy && scalar @$cx == 1;
926 # fast comp based on number of array elements (aka pseudo-length)
927 my $lxy = (scalar @$cx - scalar @$cy)
928 # or length of first element if same number of elements (aka difference 0)
930 # need int() here because sometimes the last element is '00018' vs '18'
931 (length(int($cx->[-1])) - length(int($cy->[-1])));
932 return -1 if $lxy < 0; # already differs, ret
933 return 1 if $lxy > 0; # ditto
935 # manual way (abort if unequal, good for early ne)
936 my $a; my $j = scalar @$cx;
939 last if ($a = $cx->[$j] - $cy->[$j]);
946 # compute number of digits
948 # int() because add/sub sometimes leaves strings (like '00005') instead of
949 # '5' in this place, thus causing length() to report wrong length
952 (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
957 # return the nth digit, negative values count backward
958 # zero is rightmost, so _digit(123,0) will give 3
961 my $len = _len('',$x);
963 $n = $len+$n if $n < 0; # -1 last, -2 second-to-last
964 $n = abs($n); # if negative was too big
965 $len--; $n = $len if $n > $len; # n to big?
967 my $elem = int($n / $BASE_LEN); # which array element
968 my $digit = $n % $BASE_LEN; # which digit in this element
969 $elem = '0000'.@$x[$elem]; # get element padded with 0's
970 substr($elem,-$digit-1,1);
975 # return amount of trailing zeros in decimal
976 # check each array elem in _m for having 0 at end as long as elem == 0
977 # Upon finding a elem != 0, stop
979 my $zeros = 0; my $elem;
984 $elem = "$e"; # preserve x
985 $elem =~ s/.*?(0*$)/$1/; # strip anything not zero
986 $zeros *= $BASE_LEN; # elems * 5
987 $zeros += length($elem); # count trailing zeros
990 $zeros ++; # real else branch: 50% slower!
995 ##############################################################################
1000 # return true if arg (BINT or num_str) is zero (array '+', '0')
1003 (((scalar @$x == 1) && ($x->[0] == 0))) <=> 0;
1008 # return true if arg (BINT or num_str) is even
1010 (!($x->[0] & 1)) <=> 0;
1015 # return true if arg (BINT or num_str) is even
1018 (($x->[0] & 1)) <=> 0;
1023 # return true if arg (BINT or num_str) is one (array '+', '1')
1026 (scalar @$x == 1) && ($x->[0] == 1) <=> 0;
1031 # internal normalization function that strips leading zeros from the array
1032 # args: ref to array
1035 my $cnt = scalar @$s; # get count of parts
1037 push @$s,0 if $i < 0; # div might return empty results, so fix it
1039 return $s if @$s == 1; # early out
1041 #print "strip: cnt $cnt i $i\n";
1042 # '0', '3', '4', '0', '0',
1047 # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
1048 # >= 1: skip first part (this can be zero)
1049 while ($i > 0) { last if $s->[$i] != 0; $i--; }
1050 $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
1054 ###############################################################################
1055 # check routine to test internal state for corruptions
1059 # used by the test suite
1062 return "$x is not a reference" if !ref($x);
1064 # are all parts are valid?
1065 my $i = 0; my $j = scalar @$x; my ($e,$try);
1068 $e = $x->[$i]; $e = 'undef' unless defined $e;
1069 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
1070 last if $e !~ /^[+]?[0-9]+$/;
1071 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
1072 last if "$e" !~ /^[+]?[0-9]+$/;
1073 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
1074 last if '' . "$e" !~ /^[+]?[0-9]+$/;
1075 $try = ' < 0 || >= $BASE; '."($x, $e)";
1076 last if $e <0 || $e >= $BASE;
1077 # this test is disabled, since new/bnorm and certain ops (like early out
1078 # in add/sub) are allowed/expected to leave '00000' in some elements
1079 #$try = '=~ /^00+/; '."($x, $e)";
1080 #last if $e =~ /^00+/;
1083 return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
1088 ###############################################################################
1089 ###############################################################################
1090 # some optional routines to make BigInt faster
1094 # if possible, use mod shortcut
1095 my ($c,$x,$yo) = @_;
1097 # slow way since $y to big
1098 if (scalar @$yo > 1)
1100 my ($xo,$rem) = _div($c,$x,$yo);
1105 # both are single element arrays
1106 if (scalar @$x == 1)
1112 # @y is a single element, but @x has more than one element
1116 # when BASE % Y == 0 then (B * BASE) % Y == 0
1117 # (B * BASE) % $y + A % Y => A % Y
1118 # so need to consider only last element: O(1)
1123 # else need to go through all elements: O(N), but loop is a bit simplified
1127 $r = ($r + $_) % $y; # not much faster, but heh...
1128 #$r += $_ % $y; $r %= $y;
1135 # else need to go through all elements: O(N)
1136 my $r = 0; my $bm = 1;
1139 $r = ($_ * $bm + $r) % $y;
1140 $bm = ($bm * $b) % $y;
1142 #$r += ($_ % $y) * $bm;
1150 splice (@$x,1); # keep one element of $x
1154 ##############################################################################
1159 my ($c,$x,$y,$n) = @_;
1163 $n = _new($c,\$n); return _div($c,$x, _pow($c,$n,$y));
1166 # shortcut (faster) for shifting by 10)
1167 # multiples of $BASE_LEN
1168 my $dst = 0; # destination
1169 my $src = _num($c,$y); # as normal int
1170 my $xlen = (@$x-1)*$BASE_LEN+length(int($x->[-1])); # len of x in digits
1171 if ($src > $xlen or ($src == $xlen and ! defined $x->[1]))
1173 # 12345 67890 shifted right by more than 10 digits => 0
1174 splice (@$x,1); # leave only one element
1175 $x->[0] = 0; # set to zero
1178 my $rem = $src % $BASE_LEN; # remainder to shift
1179 $src = int($src / $BASE_LEN); # source
1182 splice (@$x,0,$src); # even faster, 38.4 => 39.3
1186 my $len = scalar @$x - $src; # elems to go
1187 my $vd; my $z = '0'x $BASE_LEN;
1188 $x->[scalar @$x] = 0; # avoid || 0 test inside loop
1191 $vd = $z.$x->[$src];
1192 $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
1194 $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
1195 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1196 $x->[$dst] = int($vd);
1199 splice (@$x,$dst) if $dst > 0; # kill left-over array elems
1200 pop @$x if $x->[-1] == 0 && @$x > 1; # kill last element if 0
1207 my ($c,$x,$y,$n) = @_;
1211 $n = _new($c,\$n); return _mul($c,$x, _pow($c,$n,$y));
1214 # shortcut (faster) for shifting by 10) since we are in base 10eX
1215 # multiples of $BASE_LEN:
1216 my $src = scalar @$x; # source
1217 my $len = _num($c,$y); # shift-len as normal int
1218 my $rem = $len % $BASE_LEN; # remainder to shift
1219 my $dst = $src + int($len/$BASE_LEN); # destination
1220 my $vd; # further speedup
1221 $x->[$src] = 0; # avoid first ||0 for speed
1222 my $z = '0' x $BASE_LEN;
1225 $vd = $x->[$src]; $vd = $z.$vd;
1226 $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
1227 $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
1228 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1229 $x->[$dst] = int($vd);
1232 # set lowest parts to 0
1233 while ($dst >= 0) { $x->[$dst--] = 0; }
1234 # fix spurios last zero element
1235 splice @$x,-1 if $x->[-1] == 0;
1242 # ref to array, ref to array, return ref to array
1243 my ($c,$cx,$cy) = @_;
1245 if (scalar @$cy == 1 && $cy->[0] == 0)
1247 splice (@$cx,1); $cx->[0] = 1; # y == 0 => x => 1
1250 if ((scalar @$cx == 1 && $cx->[0] == 1) || # x == 1
1251 (scalar @$cy == 1 && $cy->[0] == 1)) # or y == 1
1255 if (scalar @$cx == 1 && $cx->[0] == 0)
1257 splice (@$cx,1); $cx->[0] = 0; # 0 ** y => 0 (if not y <= 0)
1263 my $y_bin = ${_as_bin($c,$cy)}; $y_bin =~ s/^0b//;
1264 my $len = length($y_bin);
1267 _mul($c,$pow2,$cx) if substr($y_bin,$len,1) eq '1'; # is odd?
1278 # ref to array, return ref to array
1281 if ((@$cx == 1) && ($cx->[0] <= 2))
1283 $cx->[0] ||= 1; # 0 => 1, 1 => 1, 2 => 2
1287 # go forward until $base is exceeded
1288 # limit is either $x steps (steps == 100 means a result always too high) or
1290 my $steps = 100; $steps = $cx->[0] if @$cx == 1;
1291 my $r = 2; my $cf = 3; my $step = 2; my $last = $r;
1292 while ($r*$cf < $BASE && $step < $steps)
1294 $last = $r; $r *= $cf++; $step++;
1296 if ((@$cx == 1) && $step == $cx->[0])
1298 # completely done, so keep reference to $x and return
1303 # now we must do the left over steps
1304 my $n; # steps still to do
1305 if (scalar @$cx == 1)
1314 $cx->[0] = $last; splice (@$cx,1); # keep ref to $x
1315 my $zero_elements = 0;
1317 # do left-over steps fit into a scalar?
1318 if (ref $n eq 'ARRAY')
1320 # No, so use slower inc() & cmp()
1322 while (_acmp($step,$n) <= 0)
1324 # as soon as the last element of $cx is 0, we split it up and remember
1325 # how many zeors we got so far. The reason is that n! will accumulate
1326 # zeros at the end rather fast.
1329 $zero_elements ++; shift @$cx;
1331 _mul($c,$cx,$step); _inc($c,$step);
1336 # Yes, so we can speed it up slightly
1339 # When the last element of $cx is 0, we split it up and remember
1340 # how many we got so far. The reason is that n! will accumulate
1341 # zeros at the end rather fast.
1344 $zero_elements ++; shift @$cx;
1346 _mul($c,$cx,[$step]); $step++;
1349 # multiply in the zeros again
1350 while ($zero_elements-- > 0)
1354 $cx; # return result
1359 # calculate integer log of $x to base $base
1360 # ref to array, ref to array - return ref to array
1361 my ($c,$x,$base) = @_;
1364 return if (scalar @$x == 1 && $x->[0] == 0);
1365 # BASE 0 or 1 => NaN
1366 return if (scalar @$base == 1 && $base->[0] < 2);
1367 my $cmp = _acmp($c,$x,$base); # X == BASE => 1
1370 splice (@$x,1); $x->[0] = 1;
1376 splice (@$x,1); $x->[0] = 0;
1380 # this trial multiplication is very fast, even for large counts (like for
1381 # 2 ** 1024, since this still requires only 1024 very fast steps
1382 # (multiplication of a large number by a very small number is very fast))
1383 my $x_org = _copy($c,$x); # preserve x
1384 splice(@$x,1); $x->[0] = 1; # keep ref to $x
1386 my $trial = _copy($c,$base);
1388 # XXX TODO this only works if $base has only one element
1389 if (scalar @$base == 1)
1391 # compute int ( length_in_base_10(X) / ( log(base) / log(10) ) )
1392 my $len = _len($c,$x_org);
1393 my $res = int($len / (log($base->[0]) / log(10))) || 1; # avoid $res == 0
1396 $trial = _pow ($c, _copy($c, $base), $x);
1397 my $a = _acmp($x,$trial,$x_org);
1398 return ($x,1) if $a == 0;
1399 # we now know that $res is too small
1402 _mul($c,$trial,$base); _add($c, $x, [1]);
1407 _div($c,$trial,$base); _sub($c, $x, [1]);
1409 # did we now get the right result?
1410 $a = _acmp($x,$trial,$x_org);
1411 return ($x,1) if $a == 0; # yes, exactly
1415 _div($c,$trial,$base); _sub($c, $x, [1]);
1419 # simple loop that increments $x by two in each step, possible overstepping
1420 # the real result by one
1423 my $base_mul = _mul($c, _copy($c,$base), $base);
1425 while (($a = _acmp($x,$trial,$x_org)) < 0)
1427 _mul($c,$trial,$base_mul); _add($c, $x, [2]);
1433 # overstepped the result
1435 _div($c,$trial,$base);
1436 $a = _acmp($x,$trial,$x_org);
1441 $exact = 0 if $a != 0;
1444 ($x,$exact); # return result
1448 use constant DEBUG => 0;
1450 sub steps { $steps };
1454 # square-root of $x in place
1455 # Compute a guess of the result (by rule of thumb), then improve it via
1459 if (scalar @$x == 1)
1461 # fit's into one Perl scalar, so result can be computed directly
1462 $x->[0] = int(sqrt($x->[0]));
1465 my $y = _copy($c,$x);
1466 # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
1467 # since our guess will "grow"
1468 my $l = int((_len($c,$x)-1) / 2);
1470 my $lastelem = $x->[-1]; # for guess
1471 my $elems = scalar @$x - 1;
1472 # not enough digits, but could have more?
1473 if ((length($lastelem) <= 3) && ($elems > 1))
1475 # right-align with zero pad
1476 my $len = length($lastelem) & 1;
1477 print "$lastelem => " if DEBUG;
1478 $lastelem .= substr($x->[-2] . '0' x $BASE_LEN,0,$BASE_LEN);
1479 # former odd => make odd again, or former even to even again
1480 $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
1481 print "$lastelem\n" if DEBUG;
1484 # construct $x (instead of _lsft($c,$x,$l,10)
1485 my $r = $l % $BASE_LEN; # 10000 00000 00000 00000 ($BASE_LEN=5)
1486 $l = int($l / $BASE_LEN);
1487 print "l = $l " if DEBUG;
1489 splice @$x,$l; # keep ref($x), but modify it
1491 # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
1493 # 14400 00000 => sqrt(14400) => guess first digits to be 120
1494 # 144000 000000 => sqrt(144000) => guess 379
1496 print "$lastelem (elems $elems) => " if DEBUG;
1497 $lastelem = $lastelem / 10 if ($elems & 1 == 1); # odd or even?
1498 my $g = sqrt($lastelem); $g =~ s/\.//; # 2.345 => 2345
1499 $r -= 1 if $elems & 1 == 0; # 70 => 7
1501 # padd with zeros if result is too short
1502 $x->[$l--] = int(substr($g . '0' x $r,0,$r+1));
1503 print "now ",$x->[-1] if DEBUG;
1504 print " would have been ", int('1' . '0' x $r),"\n" if DEBUG;
1506 # If @$x > 1, we could compute the second elem of the guess, too, to create
1507 # an even better guess. Not implemented yet. Does it improve performance?
1508 $x->[$l--] = 0 while ($l >= 0); # all other digits of guess are zero
1510 print "start x= ",${_str($c,$x)},"\n" if DEBUG;
1513 my $lastlast = _zero();
1514 $steps = 0 if DEBUG;
1515 while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
1518 $lastlast = _copy($c,$last);
1519 $last = _copy($c,$x);
1520 _add($c,$x, _div($c,_copy($c,$y),$x));
1522 print " x= ",${_str($c,$x)},"\n" if DEBUG;
1524 print "\nsteps in sqrt: $steps, " if DEBUG;
1525 _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0; # overshot?
1526 print " final ",$x->[-1],"\n" if DEBUG;
1532 # take n'th root of $x in place (n >= 3)
1535 if (scalar @$x == 1)
1539 # result will always be smaller than 2 so trunc to 1 at once
1544 # fit's into one Perl scalar, so result can be computed directly
1545 # cannot use int() here, because it rounds wrongly (try
1546 # (81 ** 3) ** (1/3) to see what I mean)
1547 #$x->[0] = int( $x->[0] ** (1 / $n->[0]) );
1548 # round to 8 digits, then truncate result to integer
1549 $x->[0] = int ( sprintf ("%.8f", $x->[0] ** (1 / $n->[0]) ) );
1554 # we know now that X is more than one element long
1556 # if $n is a power of two, we can repeatedly take sqrt($X) and find the
1557 # proper result, because sqrt(sqrt($x)) == root($x,4)
1558 my $b = _as_bin($c,$n);
1559 if ($$b =~ /0b1(0+)$/)
1561 my $count = CORE::length($1); # 0b100 => len('00') => 2
1562 my $cnt = $count; # counter for loop
1563 unshift (@$x, 0); # add one element, together with one
1564 # more below in the loop this makes 2
1567 # 'inflate' $X by adding one element, basically computing
1568 # $x * $BASE * $BASE. This gives us more $BASE_LEN digits for result
1569 # since len(sqrt($X)) approx == len($x) / 2.
1571 # calculate sqrt($x), $x is now one element to big, again. In the next
1572 # round we make that two, again.
1575 # $x is now one element to big, so truncate result by removing it
1580 # trial computation by starting with 2,4,8,16 etc until we overstep
1584 # while still to do more than X steps
1588 while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
1590 _mul ($c, $step, [2]);
1591 _add ($c, $trial, $step);
1595 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) == 0)
1597 @$x = @$trial; # make copy while preserving ref to $x
1600 # overstepped, so go back on step
1601 _sub($c, $trial, $step);
1602 } while (scalar @$step > 1 || $step->[0] > 128);
1606 # add two, because $trial cannot be exactly the result (otherwise we would
1607 # alrady have found it)
1608 _add($c, $trial, $step);
1610 # and now add more and more (2,4,6,8,10 etc)
1611 while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
1613 _add ($c, $trial, $step);
1616 # hit not exactly? (overstepped)
1617 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
1622 # hit not exactly? (overstepped)
1623 # 80 too small, 81 slightly too big, 82 too big
1624 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
1629 @$x = @$trial; # make copy while preserving ref to $x
1635 ##############################################################################
1642 # the shortcut makes equal, large numbers _really_ fast, and makes only a
1643 # very small performance drop for small numbers (e.g. something with less
1644 # than 32 bit) Since we optimize for large numbers, this is enabled.
1645 return $x if _acmp($c,$x,$y) == 0; # shortcut
1647 my $m = _one(); my ($xr,$yr);
1648 my $mask = $AND_MASK;
1651 my $y1 = _copy($c,$y); # make copy
1655 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1657 ($x1, $xr) = _div($c,$x1,$mask);
1658 ($y1, $yr) = _div($c,$y1,$mask);
1660 # make ints() from $xr, $yr
1661 # this is when the AND_BITS are greater tahn $BASE and is slower for
1662 # small (<256 bits) numbers, but faster for large numbers. Disabled
1663 # due to KISS principle
1665 # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1666 # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1667 # _add($c,$x, _mul($c, _new( $c, \($xrr & $yrr) ), $m) );
1669 # 0+ due to '&' doesn't work in strings
1670 _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
1680 return _zero() if _acmp($c,$x,$y) == 0; # shortcut (see -and)
1682 my $m = _one(); my ($xr,$yr);
1683 my $mask = $XOR_MASK;
1686 my $y1 = _copy($c,$y); # make copy
1690 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1692 ($x1, $xr) = _div($c,$x1,$mask);
1693 ($y1, $yr) = _div($c,$y1,$mask);
1694 # make ints() from $xr, $yr (see _and())
1695 #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1696 #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1697 #_add($c,$x, _mul($c, _new( $c, \($xrr ^ $yrr) ), $m) );
1699 # 0+ due to '^' doesn't work in strings
1700 _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
1703 # the loop stops when the shorter of the two numbers is exhausted
1704 # the remainder of the longer one will survive bit-by-bit, so we simple
1705 # multiply-add it in
1706 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1707 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1716 return $x if _acmp($c,$x,$y) == 0; # shortcut (see _and)
1718 my $m = _one(); my ($xr,$yr);
1719 my $mask = $OR_MASK;
1722 my $y1 = _copy($c,$y); # make copy
1726 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1728 ($x1, $xr) = _div($c,$x1,$mask);
1729 ($y1, $yr) = _div($c,$y1,$mask);
1730 # make ints() from $xr, $yr (see _and())
1731 # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1732 # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1733 # _add($c,$x, _mul($c, _new( $c, \($xrr | $yrr) ), $m) );
1735 # 0+ due to '|' doesn't work in strings
1736 _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
1739 # the loop stops when the shorter of the two numbers is exhausted
1740 # the remainder of the longer one will survive bit-by-bit, so we simple
1741 # multiply-add it in
1742 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1743 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1750 # convert a decimal number to hex (ref to array, return ref to string)
1753 # fit's into one element (handle also 0x0 case)
1756 my $t = sprintf("0x%x",$x->[0]);
1760 my $x1 = _copy($c,$x);
1763 my ($xr, $h, $x10000);
1766 $x10000 = [ 0x10000 ]; $h = 'h4';
1770 $x10000 = [ 0x1000 ]; $h = 'h3';
1772 # while (! _is_zero($c,$x1))
1773 while (@$x1 != 1 || $x1->[0] != 0) # _is_zero()
1775 ($x1, $xr) = _div($c,$x1,$x10000);
1776 $es .= unpack($h,pack('v',$xr->[0])); # XXX TODO: why pack('v',...)?
1779 $es =~ s/^[0]+//; # strip leading zeros
1786 # convert a decimal number to bin (ref to array, return ref to string)
1789 # fit's into one element (and Perl recent enough), handle also 0b0 case
1790 # handle zero case for older Perls
1791 if ($] <= 5.005 && @$x == 1 && $x->[0] == 0)
1793 my $t = '0b0'; return \$t;
1795 if (@$x == 1 && $] >= 5.006)
1797 my $t = sprintf("0b%b",$x->[0]);
1800 my $x1 = _copy($c,$x);
1803 my ($xr, $b, $x10000);
1806 $x10000 = [ 0x10000 ]; $b = 'b16';
1810 $x10000 = [ 0x1000 ]; $b = 'b12';
1812 # while (! _is_zero($c,$x1))
1813 while (!(@$x1 == 1 && $x1->[0] == 0)) # _is_zero()
1815 ($x1, $xr) = _div($c,$x1,$x10000);
1816 $es .= unpack($b,pack('v',$xr->[0])); # XXX TODO: why pack('v',...)?
1817 # $es .= unpack($b,$xr->[0]);
1820 $es =~ s/^[0]+//; # strip leading zeros
1827 # convert a hex number to decimal (ref to string, return ref to array)
1831 my $m = [ 0x10000 ]; # 16 bit at a time
1834 my $len = length($$hs)-2;
1835 $len = int($len/4); # 4-digit parts, w/o '0x'
1836 my $val; my $i = -4;
1839 $val = substr($$hs,$i,4);
1840 $val =~ s/^[+-]?0x// if $len == 0; # for last part only because
1841 $val = hex($val); # hex does not like wrong chars
1843 _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1844 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
1851 # convert a hex number to decimal (ref to string, return ref to array)
1854 # instead of converting X (8) bit at a time, it is faster to "convert" the
1855 # number to hex, and then call _from_hex.
1858 $hs =~ s/^[+-]?0b//; # remove sign and 0b
1859 my $l = length($hs); # bits
1860 $hs = '0' x (8-($l % 8)) . $hs if ($l % 8) != 0; # padd left side w/ 0
1861 my $h = unpack('H*', pack ('B*', $hs)); # repack as hex
1863 $c->_from_hex(\('0x'.$h));
1866 ##############################################################################
1867 # special modulus functions
1874 my $u = _zero($c); my $u1 = _one($c);
1875 my $a = _copy($c,$y); my $b = _copy($c,$x);
1877 # Euclid's Algorithm for bgcd(), only that we calc bgcd() ($a) and the
1878 # result ($u) at the same time. See comments in BigInt for why this works.
1880 ($a, $q, $b) = ($b, _div($c,$a,$b)); # step 1
1882 while (!_is_zero($c,$b))
1884 my $t = _add($c, # step 2:
1885 _mul($c,_copy($c,$u1), $q) , # t = u1 * q
1887 $u = $u1; # u = u1, u1 = t
1890 ($a, $q, $b) = ($b, _div($c,$a,$b)); # step 1
1893 # if the gcd is not 1, then return NaN
1894 return (undef,undef) unless _is_one($c,$a);
1896 $sign = $sign == 1 ? '+' : '-';
1902 # modulus of power ($x ** $y) % $z
1903 my ($c,$num,$exp,$mod) = @_;
1905 # in the trivial case,
1906 if (_is_one($c,$mod))
1908 splice @$num,0,1; $num->[0] = 0;
1911 if ((scalar @$num == 1) && (($num->[0] == 0) || ($num->[0] == 1)))
1917 # $num = _mod($c,$num,$mod); # this does not make it faster
1919 my $acc = _copy($c,$num); my $t = _one();
1921 my $expbin = ${_as_bin($c,$exp)}; $expbin =~ s/^0b//;
1922 my $len = length($expbin);
1925 if ( substr($expbin,$len,1) eq '1') # is_odd
1928 $t = _mod($c,$t,$mod);
1931 $acc = _mod($c,$acc,$mod);
1937 ##############################################################################
1938 ##############################################################################
1945 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
1949 Provides support for big integer calculations. Not intended to be used by other
1950 modules. Other modules which sport the same functions can also be used to support
1951 Math::BigInt, like Math::BigInt::GMP or Math::BigInt::Pari.
1955 In order to allow for multiple big integer libraries, Math::BigInt was
1956 rewritten to use library modules for core math routines. Any module which
1957 follows the same API as this can be used instead by using the following:
1959 use Math::BigInt lib => 'libname';
1961 'libname' is either the long name ('Math::BigInt::Pari'), or only the short
1962 version like 'Pari'.
1968 The following functions MUST be defined in order to support the use by
1971 _new(string) return ref to new object from ref to decimal string
1972 _zero() return a new object with value 0
1973 _one() return a new object with value 1
1975 _str(obj) return ref to a string representing the object
1976 _num(obj) returns a Perl integer/floating point number
1977 NOTE: because of Perl numeric notation defaults,
1978 the _num'ified obj may lose accuracy due to
1979 machine-dependend floating point size limitations
1981 _add(obj,obj) Simple addition of two objects
1982 _mul(obj,obj) Multiplication of two objects
1983 _div(obj,obj) Division of the 1st object by the 2nd
1984 In list context, returns (result,remainder).
1985 NOTE: this is integer math, so no
1986 fractional part will be returned.
1987 The second operand will be not be 0, so no need to
1989 _sub(obj,obj) Simple subtraction of 1 object from another
1990 a third, optional parameter indicates that the params
1991 are swapped. In this case, the first param needs to
1992 be preserved, while you can destroy the second.
1993 sub (x,y,1) => return x - y and keep x intact!
1994 _dec(obj) decrement object by one (input is garant. to be > 0)
1995 _inc(obj) increment object by one
1998 _acmp(obj,obj) <=> operator for objects (return -1, 0 or 1)
2000 _len(obj) returns count of the decimal digits of the object
2001 _digit(obj,n) returns the n'th decimal digit of object
2003 _is_one(obj) return true if argument is +1
2004 _is_zero(obj) return true if argument is 0
2005 _is_even(obj) return true if argument is even (0,2,4,6..)
2006 _is_odd(obj) return true if argument is odd (1,3,5,7..)
2008 _copy return a ref to a true copy of the object
2010 _check(obj) check whether internal representation is still intact
2011 return 0 for ok, otherwise error message as string
2013 The following functions are optional, and can be defined if the underlying lib
2014 has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
2015 slow) fallback routines to emulate these:
2017 _from_hex(str) return ref to new object from ref to hexadecimal string
2018 _from_bin(str) return ref to new object from ref to binary string
2020 _as_hex(str) return ref to scalar string containing the value as
2021 unsigned hex string, with the '0x' prepended.
2022 Leading zeros must be stripped.
2023 _as_bin(str) Like as_hex, only as binary string containing only
2024 zeros and ones. Leading zeros must be stripped and a
2025 '0b' must be prepended.
2027 _rsft(obj,N,B) shift object in base B by N 'digits' right
2028 For unsupported bases B, return undef to signal failure
2029 _lsft(obj,N,B) shift object in base B by N 'digits' left
2030 For unsupported bases B, return undef to signal failure
2032 _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
2033 Note: XOR, AND and OR pad with zeros if size mismatches
2034 _and(obj1,obj2) AND (bit-wise) object 1 with object 2
2035 _or(obj1,obj2) OR (bit-wise) object 1 with object 2
2041 _mod(obj,obj) Return remainder of div of the 1st by the 2nd object
2042 _sqrt(obj) return the square root of object (truncated to int)
2043 _root(obj) return the n'th (n >= 3) root of obj (truncated to int)
2044 _fac(obj) return factorial of object 1 (1*2*3*4..)
2045 _pow(obj,obj) return object 1 to the power of object 2
2046 return undef for NaN
2047 _gcd(obj,obj) return Greatest Common Divisor of two objects
2049 _zeros(obj) return number of trailing decimal zeros
2050 _modinv return inverse modulus
2051 _modpow return modulus of power ($x ** $y) % $z
2052 _log_int(X,N) calculate integer log() of X in base N
2053 X >= 0, N >= 0 (return undef for NaN)
2054 returns (RESULT, EXACT) where EXACT is:
2055 1 : result is exactly RESULT
2056 0 : result was truncated to RESULT
2057 undef : unknown whether result is exactly RESULT
2059 Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
2062 So the library needs only to deal with unsigned big integers. Testing of input
2063 parameter validity is done by the caller, so you need not worry about
2064 underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by zero or similar
2067 The first parameter can be modified, that includes the possibility that you
2068 return a reference to a completely different object instead. Although keeping
2069 the reference and just changing it's contents is prefered over creating and
2070 returning a different reference.
2072 Return values are always references to objects, strings, or true/false for
2073 comparisation routines.
2075 Exceptions are C<_lsft()> and C<_rsft()>, which return undef if they can not
2076 shift the argument. This is used to delegate shifting of bases different than
2077 the one you can support back to Math::BigInt, which will use some generic code
2078 to calculate the result.
2080 =head1 WRAP YOUR OWN
2082 If you want to port your own favourite c-lib for big numbers to the
2083 Math::BigInt interface, you can take any of the already existing modules as
2084 a rough guideline. You should really wrap up the latest BigInt and BigFloat
2085 testsuites with your module, and replace in them any of the following:
2091 use Math::BigInt lib => 'yourlib';
2093 This way you ensure that your library really works 100% within Math::BigInt.
2097 This program is free software; you may redistribute it and/or modify it under
2098 the same terms as Perl itself.
2102 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
2104 Seperated from BigInt and shaped API with the help of John Peacock.
2105 Fixed, sped-up and enhanced by Tels http://bloodgate.com 2001-2003.
2109 L<Math::BigInt>, L<Math::BigFloat>, L<Math::BigInt::BitVect>,
2110 L<Math::BigInt::GMP>, L<Math::BigInt::FastCalc> and L<Math::BigInt::Pari>.