1 package Math::BigInt::Calc;
5 # use warnings; # do not use warnings for older Perls
7 our $VERSION = '1.999701';
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 who knows 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 length 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 Perl scalar number (int/float) from a BigInt object.
278 return 0 + $x->[0] if scalar @$x == 1; # below $BASE
280 # Start with the most significant element and work towards the least
281 # significant element. Avoid multiplying "inf" (which happens if the number
282 # overflows) with "0" (if there are zero elements in $x) since this gives
283 # "nan" which propagates to the output.
286 for (my $i = $#$x ; $i >= 0 ; --$i) {
293 ##############################################################################
298 # (ref to int_num_array, ref to int_num_array)
299 # routine to add two base 1eX numbers
300 # stolen from Knuth Vol 2 Algorithm A pg 231
301 # there are separate routines to add and sub as per Knuth pg 233
302 # This routine modifies array x, but not y.
306 return $x if (@$y == 1) && $y->[0] == 0; # $x + 0 => $x
307 if ((@$x == 1) && $x->[0] == 0) # 0 + $y => $y->copy
309 # twice as slow as $x = [ @$y ], but nec. to retain $x as ref :(
310 @$x = @$y; return $x;
313 # for each in Y, add Y to X and carry. If after that, something is left in
314 # X, foreach in X add carry to X and then return X, carry
315 # Trades one "$j++" for having to shift arrays
316 my $i; my $car = 0; my $j = 0;
319 $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
324 $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
331 # (ref to int_num_array, ref to int_num_array)
332 # Add 1 to $x, modify $x in place
337 return $x if (($i += 1) < $BASE); # early out
338 $i = 0; # overflow, next
340 push @$x,1 if (($x->[-1] || 0) == 0); # last overflowed, so extend
346 # (ref to int_num_array, ref to int_num_array)
347 # Sub 1 from $x, modify $x in place
350 my $MAX = $BASE-1; # since MAX_VAL based on BASE
353 last if (($i -= 1) >= 0); # early out
354 $i = $MAX; # underflow, next
356 pop @$x if $x->[-1] == 0 && @$x > 1; # last underflowed (but leave 0)
362 # (ref to int_num_array, ref to int_num_array, swap)
363 # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
364 # subtract Y from X by modifying x in place
365 my ($c,$sx,$sy,$s) = @_;
367 my $car = 0; my $i; my $j = 0;
372 last unless defined $sy->[$j] || $car;
373 $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
375 # might leave leading zeros, so fix that
376 return __strip_zeros($sx);
380 # we can't do an early out if $x is < than $y, since we
381 # need to copy the high chunks from $y. Found by Bob Mathews.
382 #last unless defined $sy->[$j] || $car;
384 if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
387 # might leave leading zeros, so fix that
393 # (ref to int_num_array, ref to int_num_array)
394 # multiply two numbers in internal representation
395 # modifies first arg, second need not be different from first
396 my ($c,$xv,$yv) = @_;
400 # shortcut for two very short numbers (improved by Nathan Zook)
401 # works also if xv and yv are the same reference, and handles also $x == 0
404 if (($xv->[0] *= $yv->[0]) >= $BASE)
406 $xv->[0] = $xv->[0] - ($xv->[1] = int($xv->[0] * $RBASE)) * $BASE;
416 # multiply a large number a by a single element one, so speed up
417 my $y = $yv->[0]; my $car = 0;
420 $i = $i * $y + $car; $car = int($i * $RBASE); $i -= $car * $BASE;
422 push @$xv, $car if $car != 0;
425 # shortcut for result $x == 0 => result = 0
426 return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
428 # since multiplying $x with $x fails, make copy in this case
429 $yv = [@$xv] if $xv == $yv; # same references?
431 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
440 # $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
442 # $prod - ($car = int($prod * RBASE)) * $BASE; # see USE_MUL
444 # $prod[$cty] += $car if $car; # need really to check for 0?
448 # looping through this if $xi == 0 is silly - so optimize it away!
449 $xi = (shift @prod || 0), next if $xi == 0;
452 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
453 ## this is actually a tad slower
454 ## $prod = $prod[$cty]; $prod += ($car + $xi * $yi); # no ||0 here
456 $prod - ($car = int($prod * $RBASE)) * $BASE; # see USE_MUL
458 $prod[$cty] += $car if $car; # need really to check for 0?
459 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
462 # can't have leading zeros
463 # __strip_zeros($xv);
469 # (ref to int_num_array, ref to int_num_array)
470 # multiply two numbers in internal representation
471 # modifies first arg, second need not be different from first
472 # works for 64 bit integer with "use integer"
473 my ($c,$xv,$yv) = @_;
478 # shortcut for two small numbers, also handles $x == 0
481 # shortcut for two very short numbers (improved by Nathan Zook)
482 # works also if xv and yv are the same reference, and handles also $x == 0
483 if (($xv->[0] *= $yv->[0]) >= $BASE)
486 $xv->[0] - ($xv->[1] = $xv->[0] / $BASE) * $BASE;
496 # multiply a large number a by a single element one, so speed up
497 my $y = $yv->[0]; my $car = 0;
500 #$i = $i * $y + $car; $car = $i / $BASE; $i -= $car * $BASE;
501 $i = $i * $y + $car; $i -= ($car = $i / $BASE) * $BASE;
503 push @$xv, $car if $car != 0;
506 # shortcut for result $x == 0 => result = 0
507 return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
509 # since multiplying $x with $x fails, make copy in this case
510 $yv = [@$xv] if $xv == $yv; # same references?
512 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
516 # looping through this if $xi == 0 is silly - so optimize it away!
517 $xi = (shift @prod || 0), next if $xi == 0;
520 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
521 $prod[$cty++] = $prod - ($car = $prod / $BASE) * $BASE;
523 $prod[$cty] += $car if $car; # need really to check for 0?
524 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
532 # (ref to int_num_array, ref to int_num_array)
533 # multiply two numbers in internal representation
534 # modifies first arg, second need not be different from first
535 my ($c,$xv,$yv) = @_;
539 # shortcut for two small numbers, also handles $x == 0
542 # shortcut for two very short numbers (improved by Nathan Zook)
543 # works also if xv and yv are the same reference, and handles also $x == 0
544 if (($xv->[0] *= $yv->[0]) >= $BASE)
547 $xv->[0] - ($xv->[1] = int($xv->[0] / $BASE)) * $BASE;
557 # multiply a large number a by a single element one, so speed up
558 my $y = $yv->[0]; my $car = 0;
561 $i = $i * $y + $car; $car = int($i / $BASE); $i -= $car * $BASE;
562 # This (together with use integer;) does not work on 32-bit Perls
563 #$i = $i * $y + $car; $i -= ($car = $i / $BASE) * $BASE;
565 push @$xv, $car if $car != 0;
568 # shortcut for result $x == 0 => result = 0
569 return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
571 # since multiplying $x with $x fails, make copy in this case
572 $yv = [@$xv] if $xv == $yv; # same references?
574 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
578 # looping through this if $xi == 0 is silly - so optimize it away!
579 $xi = (shift @prod || 0), next if $xi == 0;
582 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
583 $prod[$cty++] = $prod - ($car = int($prod / $BASE)) * $BASE;
585 $prod[$cty] += $car if $car; # need really to check for 0?
586 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
589 # can't have leading zeros
590 # __strip_zeros($xv);
596 # ref to array, ref to array, modify first array and return remainder if
599 # see comments in _div_use_div() for more explanations
601 my ($c,$x,$yorg) = @_;
603 # the general div algorithm here is about O(N*N) and thus quite slow, so
604 # we first check for some special cases and use shortcuts to handle them.
606 # This works, because we store the numbers in a chunked format where each
607 # element contains 5..7 digits (depending on system).
609 # if both numbers have only one element:
610 if (@$x == 1 && @$yorg == 1)
612 # shortcut, $yorg and $x are two small numbers
615 my $r = [ $x->[0] % $yorg->[0] ];
616 $x->[0] = int($x->[0] / $yorg->[0]);
621 $x->[0] = int($x->[0] / $yorg->[0]);
626 # if x has more than one, but y has only one element:
630 $rem = _mod($c,[ @$x ],$yorg) if wantarray;
632 # shortcut, $y is < $BASE
633 my $j = scalar @$x; my $r = 0;
634 my $y = $yorg->[0]; my $b;
637 $b = $r * $BASE + $x->[$j];
638 $x->[$j] = int($b/$y);
641 pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
642 return ($x,$rem) if wantarray;
646 # now x and y have more than one element
648 # check whether y has more elements than x, if yet, the result will be 0
652 $rem = [@$x] if wantarray; # make copy
653 splice (@$x,1); # keep ref to original array
654 $x->[0] = 0; # set to 0
655 return ($x,$rem) if wantarray; # including remainder?
656 return $x; # only x, which is [0] now
658 # check whether the numbers have the same number of elements, in that case
659 # the result will fit into one element and can be computed efficiently
663 # if $yorg has more digits than $x (it's leading element is longer than
664 # the one from $x), the result will also be 0:
665 if (length(int($yorg->[-1])) > length(int($x->[-1])))
667 $rem = [@$x] if wantarray; # make copy
668 splice (@$x,1); # keep ref to org array
669 $x->[0] = 0; # set to 0
670 return ($x,$rem) if wantarray; # including remainder?
673 # now calculate $x / $yorg
674 if (length(int($yorg->[-1])) == length(int($x->[-1])))
676 # same length, so make full compare
678 my $a = 0; my $j = scalar @$x - 1;
679 # manual way (abort if unequal, good for early ne)
682 last if ($a = $x->[$j] - $yorg->[$j]); $j--;
684 # $a contains the result of the compare between X and Y
685 # a < 0: x < y, a == 0: x == y, a > 0: x > y
688 $rem = [ 0 ]; # a = 0 => x == y => rem 0
689 $rem = [@$x] if $a != 0; # a < 0 => x < y => rem = x
690 splice(@$x,1); # keep single element
691 $x->[0] = 0; # if $a < 0
692 $x->[0] = 1 if $a == 0; # $x == $y
693 return ($x,$rem) if wantarray;
696 # $x >= $y, so proceed normally
702 my $y = [ @$yorg ]; # always make copy to preserve
704 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
706 $car = $bar = $prd = 0;
707 if (($dd = int($BASE/($y->[-1]+1))) != 1)
711 $xi = $xi * $dd + $car;
712 $xi -= ($car = int($xi * $RBASE)) * $BASE; # see USE_MUL
714 push(@$x, $car); $car = 0;
717 $yi = $yi * $dd + $car;
718 $yi -= ($car = int($yi * $RBASE)) * $BASE; # see USE_MUL
725 @q = (); ($v2,$v1) = @$y[-2,-1];
729 ($u2,$u1,$u0) = @$x[-3..-1];
731 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
733 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
734 --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
737 ($car, $bar) = (0,0);
738 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
740 $prd = $q * $y->[$yi] + $car;
741 $prd -= ($car = int($prd * $RBASE)) * $BASE; # see USE_MUL
742 $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
744 if ($x->[-1] < $car + $bar)
747 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
750 if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
763 for $xi (reverse @$x)
765 $prd = $car * $BASE + $xi;
766 $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
787 # ref to array, ref to array, modify first array and return remainder if
789 # This version works on 64 bit integers
790 my ($c,$x,$yorg) = @_;
793 # the general div algorithm here is about O(N*N) and thus quite slow, so
794 # we first check for some special cases and use shortcuts to handle them.
796 # This works, because we store the numbers in a chunked format where each
797 # element contains 5..7 digits (depending on system).
799 # if both numbers have only one element:
800 if (@$x == 1 && @$yorg == 1)
802 # shortcut, $yorg and $x are two small numbers
805 my $r = [ $x->[0] % $yorg->[0] ];
806 $x->[0] = int($x->[0] / $yorg->[0]);
811 $x->[0] = int($x->[0] / $yorg->[0]);
815 # if x has more than one, but y has only one element:
819 $rem = _mod($c,[ @$x ],$yorg) if wantarray;
821 # shortcut, $y is < $BASE
822 my $j = scalar @$x; my $r = 0;
823 my $y = $yorg->[0]; my $b;
826 $b = $r * $BASE + $x->[$j];
827 $x->[$j] = int($b/$y);
830 pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
831 return ($x,$rem) if wantarray;
834 # now x and y have more than one element
836 # check whether y has more elements than x, if yet, the result will be 0
840 $rem = [@$x] if wantarray; # make copy
841 splice (@$x,1); # keep ref to original array
842 $x->[0] = 0; # set to 0
843 return ($x,$rem) if wantarray; # including remainder?
844 return $x; # only x, which is [0] now
846 # check whether the numbers have the same number of elements, in that case
847 # the result will fit into one element and can be computed efficiently
851 # if $yorg has more digits than $x (it's leading element is longer than
852 # the one from $x), the result will also be 0:
853 if (length(int($yorg->[-1])) > length(int($x->[-1])))
855 $rem = [@$x] if wantarray; # make copy
856 splice (@$x,1); # keep ref to org array
857 $x->[0] = 0; # set to 0
858 return ($x,$rem) if wantarray; # including remainder?
861 # now calculate $x / $yorg
863 if (length(int($yorg->[-1])) == length(int($x->[-1])))
865 # same length, so make full compare
867 my $a = 0; my $j = scalar @$x - 1;
868 # manual way (abort if unequal, good for early ne)
871 last if ($a = $x->[$j] - $yorg->[$j]); $j--;
873 # $a contains the result of the compare between X and Y
874 # a < 0: x < y, a == 0: x == y, a > 0: x > y
877 $rem = [ 0 ]; # a = 0 => x == y => rem 0
878 $rem = [@$x] if $a != 0; # a < 0 => x < y => rem = x
879 splice(@$x,1); # keep single element
880 $x->[0] = 0; # if $a < 0
881 $x->[0] = 1 if $a == 0; # $x == $y
882 return ($x,$rem) if wantarray; # including remainder?
885 # $x >= $y, so proceed normally
892 my $y = [ @$yorg ]; # always make copy to preserve
894 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
896 $car = $bar = $prd = 0;
897 if (($dd = int($BASE/($y->[-1]+1))) != 1)
901 $xi = $xi * $dd + $car;
902 $xi -= ($car = int($xi / $BASE)) * $BASE;
904 push(@$x, $car); $car = 0;
907 $yi = $yi * $dd + $car;
908 $yi -= ($car = int($yi / $BASE)) * $BASE;
916 # @q will accumulate the final result, $q contains the current computed
917 # part of the final result
919 @q = (); ($v2,$v1) = @$y[-2,-1];
923 ($u2,$u1,$u0) = @$x[-3..-1];
925 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
927 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
928 --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
931 ($car, $bar) = (0,0);
932 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
934 $prd = $q * $y->[$yi] + $car;
935 $prd -= ($car = int($prd / $BASE)) * $BASE;
936 $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
938 if ($x->[-1] < $car + $bar)
941 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
944 if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
948 pop(@$x); unshift(@q, $q);
956 for $xi (reverse @$x)
958 $prd = $car * $BASE + $xi;
959 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
980 # ref to array, ref to array, modify first array and return remainder if
982 my ($c,$x,$yorg) = @_;
984 # the general div algorithm here is about O(N*N) and thus quite slow, so
985 # we first check for some special cases and use shortcuts to handle them.
987 # This works, because we store the numbers in a chunked format where each
988 # element contains 5..7 digits (depending on system).
990 # if both numbers have only one element:
991 if (@$x == 1 && @$yorg == 1)
993 # shortcut, $yorg and $x are two small numbers
996 my $r = [ $x->[0] % $yorg->[0] ];
997 $x->[0] = int($x->[0] / $yorg->[0]);
1002 $x->[0] = int($x->[0] / $yorg->[0]);
1006 # if x has more than one, but y has only one element:
1010 $rem = _mod($c,[ @$x ],$yorg) if wantarray;
1012 # shortcut, $y is < $BASE
1013 my $j = scalar @$x; my $r = 0;
1014 my $y = $yorg->[0]; my $b;
1017 $b = $r * $BASE + $x->[$j];
1018 $x->[$j] = int($b/$y);
1021 pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
1022 return ($x,$rem) if wantarray;
1025 # now x and y have more than one element
1027 # check whether y has more elements than x, if yet, the result will be 0
1031 $rem = [@$x] if wantarray; # make copy
1032 splice (@$x,1); # keep ref to original array
1033 $x->[0] = 0; # set to 0
1034 return ($x,$rem) if wantarray; # including remainder?
1035 return $x; # only x, which is [0] now
1037 # check whether the numbers have the same number of elements, in that case
1038 # the result will fit into one element and can be computed efficiently
1042 # if $yorg has more digits than $x (it's leading element is longer than
1043 # the one from $x), the result will also be 0:
1044 if (length(int($yorg->[-1])) > length(int($x->[-1])))
1046 $rem = [@$x] if wantarray; # make copy
1047 splice (@$x,1); # keep ref to org array
1048 $x->[0] = 0; # set to 0
1049 return ($x,$rem) if wantarray; # including remainder?
1052 # now calculate $x / $yorg
1054 if (length(int($yorg->[-1])) == length(int($x->[-1])))
1056 # same length, so make full compare
1058 my $a = 0; my $j = scalar @$x - 1;
1059 # manual way (abort if unequal, good for early ne)
1062 last if ($a = $x->[$j] - $yorg->[$j]); $j--;
1064 # $a contains the result of the compare between X and Y
1065 # a < 0: x < y, a == 0: x == y, a > 0: x > y
1068 $rem = [ 0 ]; # a = 0 => x == y => rem 0
1069 $rem = [@$x] if $a != 0; # a < 0 => x < y => rem = x
1070 splice(@$x,1); # keep single element
1071 $x->[0] = 0; # if $a < 0
1072 $x->[0] = 1 if $a == 0; # $x == $y
1073 return ($x,$rem) if wantarray; # including remainder?
1076 # $x >= $y, so proceed normally
1083 my $y = [ @$yorg ]; # always make copy to preserve
1085 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
1087 $car = $bar = $prd = 0;
1088 if (($dd = int($BASE/($y->[-1]+1))) != 1)
1092 $xi = $xi * $dd + $car;
1093 $xi -= ($car = int($xi / $BASE)) * $BASE;
1095 push(@$x, $car); $car = 0;
1098 $yi = $yi * $dd + $car;
1099 $yi -= ($car = int($yi / $BASE)) * $BASE;
1107 # @q will accumulate the final result, $q contains the current computed
1108 # part of the final result
1110 @q = (); ($v2,$v1) = @$y[-2,-1];
1114 ($u2,$u1,$u0) = @$x[-3..-1];
1116 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
1118 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
1119 --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
1122 ($car, $bar) = (0,0);
1123 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
1125 $prd = $q * $y->[$yi] + $car;
1126 $prd -= ($car = int($prd / $BASE)) * $BASE;
1127 $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
1129 if ($x->[-1] < $car + $bar)
1132 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
1135 if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
1139 pop(@$x); unshift(@q, $q);
1147 for $xi (reverse @$x)
1149 $prd = $car * $BASE + $xi;
1150 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
1169 ##############################################################################
1174 # internal absolute post-normalized compare (ignore signs)
1175 # ref to array, ref to array, return <0, 0, >0
1176 # arrays must have at least one entry; this is not checked for
1177 my ($c,$cx,$cy) = @_;
1179 # shortcut for short numbers
1180 return (($cx->[0] <=> $cy->[0]) <=> 0)
1181 if scalar @$cx == scalar @$cy && scalar @$cx == 1;
1183 # fast comp based on number of array elements (aka pseudo-length)
1184 my $lxy = (scalar @$cx - scalar @$cy)
1185 # or length of first element if same number of elements (aka difference 0)
1187 # need int() here because sometimes the last element is '00018' vs '18'
1188 (length(int($cx->[-1])) - length(int($cy->[-1])));
1189 return -1 if $lxy < 0; # already differs, ret
1190 return 1 if $lxy > 0; # ditto
1192 # manual way (abort if unequal, good for early ne)
1193 my $a; my $j = scalar @$cx;
1196 last if ($a = $cx->[$j] - $cy->[$j]);
1203 # compute number of digits in base 10
1205 # int() because add/sub sometimes leaves strings (like '00005') instead of
1206 # '5' in this place, thus causing length() to report wrong length
1209 (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
1214 # Return the nth digit. Zero is rightmost, so _digit(123,0) gives 3.
1215 # Negative values count from the left, so _digit(123, -1) gives 1.
1218 my $len = _len('',$x);
1220 $n += $len if $n < 0; # -1 last, -2 second-to-last
1221 return "0" if $n < 0 || $n >= $len; # return 0 for digits out of range
1223 my $elem = int($n / $BASE_LEN); # which array element
1224 my $digit = $n % $BASE_LEN; # which digit in this element
1225 substr("$x->[$elem]", -$digit-1, 1);
1230 # return amount of trailing zeros in decimal
1231 # check each array elem in _m for having 0 at end as long as elem == 0
1232 # Upon finding a elem != 0, stop
1235 return 0 if scalar @$x == 1 && $x->[0] == 0;
1237 my $zeros = 0; my $elem;
1242 $elem = "$e"; # preserve x
1243 $elem =~ s/.*?(0*$)/$1/; # strip anything not zero
1244 $zeros *= $BASE_LEN; # elems * 5
1245 $zeros += length($elem); # count trailing zeros
1248 $zeros ++; # real else branch: 50% slower!
1253 ##############################################################################
1258 # return true if arg is zero
1259 (((scalar @{$_[1]} == 1) && ($_[1]->[0] == 0))) <=> 0;
1264 # return true if arg is even
1265 (!($_[1]->[0] & 1)) <=> 0;
1270 # return true if arg is odd
1271 (($_[1]->[0] & 1)) <=> 0;
1276 # return true if arg is one
1277 (scalar @{$_[1]} == 1) && ($_[1]->[0] == 1) <=> 0;
1282 # return true if arg is two
1283 (scalar @{$_[1]} == 1) && ($_[1]->[0] == 2) <=> 0;
1288 # return true if arg is ten
1289 (scalar @{$_[1]} == 1) && ($_[1]->[0] == 10) <=> 0;
1294 # internal normalization function that strips leading zeros from the array
1295 # args: ref to array
1298 my $cnt = scalar @$s; # get count of parts
1300 push @$s,0 if $i < 0; # div might return empty results, so fix it
1302 return $s if @$s == 1; # early out
1304 #print "strip: cnt $cnt i $i\n";
1305 # '0', '3', '4', '0', '0',
1310 # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
1311 # >= 1: skip first part (this can be zero)
1312 while ($i > 0) { last if $s->[$i] != 0; $i--; }
1313 $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
1317 ###############################################################################
1318 # check routine to test internal state for corruptions
1322 # used by the test suite
1325 return "$x is not a reference" if !ref($x);
1327 # are all parts are valid?
1328 my $i = 0; my $j = scalar @$x; my ($e,$try);
1331 $e = $x->[$i]; $e = 'undef' unless defined $e;
1332 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
1333 last if $e !~ /^[+]?[0-9]+$/;
1334 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
1335 last if "$e" !~ /^[+]?[0-9]+$/;
1336 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
1337 last if '' . "$e" !~ /^[+]?[0-9]+$/;
1338 $try = ' < 0 || >= $BASE; '."($x, $e)";
1339 last if $e <0 || $e >= $BASE;
1340 # this test is disabled, since new/bnorm and certain ops (like early out
1341 # in add/sub) are allowed/expected to leave '00000' in some elements
1342 #$try = '=~ /^00+/; '."($x, $e)";
1343 #last if $e =~ /^00+/;
1346 return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
1351 ###############################################################################
1355 # if possible, use mod shortcut
1356 my ($c,$x,$yo) = @_;
1358 # slow way since $y too big
1359 if (scalar @$yo > 1)
1361 my ($xo,$rem) = _div($c,$x,$yo);
1368 # if both are single element arrays
1369 if (scalar @$x == 1)
1375 # if @$x has more than one element, but @$y is a single element
1379 # when BASE % Y == 0 then (B * BASE) % Y == 0
1380 # (B * BASE) % $y + A % Y => A % Y
1381 # so need to consider only last element: O(1)
1386 # else need to go through all elements in @$x: O(N), but loop is a bit
1391 $r = ($r + $_) % $y; # not much faster, but heh...
1392 #$r += $_ % $y; $r %= $y;
1399 # else need to go through all elements in @$x: O(N)
1404 $r = ($_ * $bm + $r) % $y;
1405 $bm = ($bm * $b) % $y;
1407 #$r += ($_ % $y) * $bm;
1415 @$x = $x->[0]; # keep one element of @$x
1419 ##############################################################################
1424 my ($c,$x,$y,$n) = @_;
1428 $n = _new($c,$n); return _div($c,$x, _pow($c,$n,$y));
1431 # shortcut (faster) for shifting by 10)
1432 # multiples of $BASE_LEN
1433 my $dst = 0; # destination
1434 my $src = _num($c,$y); # as normal int
1435 my $xlen = (@$x-1)*$BASE_LEN+length(int($x->[-1])); # len of x in digits
1436 if ($src >= $xlen or ($src == $xlen and ! defined $x->[1]))
1438 # 12345 67890 shifted right by more than 10 digits => 0
1439 splice (@$x,1); # leave only one element
1440 $x->[0] = 0; # set to zero
1443 my $rem = $src % $BASE_LEN; # remainder to shift
1444 $src = int($src / $BASE_LEN); # source
1447 splice (@$x,0,$src); # even faster, 38.4 => 39.3
1451 my $len = scalar @$x - $src; # elems to go
1452 my $vd; my $z = '0'x $BASE_LEN;
1453 $x->[scalar @$x] = 0; # avoid || 0 test inside loop
1456 $vd = $z.$x->[$src];
1457 $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
1459 $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
1460 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1461 $x->[$dst] = int($vd);
1464 splice (@$x,$dst) if $dst > 0; # kill left-over array elems
1465 pop @$x if $x->[-1] == 0 && @$x > 1; # kill last element if 0
1472 my ($c,$x,$y,$n) = @_;
1476 $n = _new($c,$n); return _mul($c,$x, _pow($c,$n,$y));
1479 # shortcut (faster) for shifting by 10) since we are in base 10eX
1480 # multiples of $BASE_LEN:
1481 my $src = scalar @$x; # source
1482 my $len = _num($c,$y); # shift-len as normal int
1483 my $rem = $len % $BASE_LEN; # remainder to shift
1484 my $dst = $src + int($len/$BASE_LEN); # destination
1485 my $vd; # further speedup
1486 $x->[$src] = 0; # avoid first ||0 for speed
1487 my $z = '0' x $BASE_LEN;
1490 $vd = $x->[$src]; $vd = $z.$vd;
1491 $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
1492 $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
1493 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1494 $x->[$dst] = int($vd);
1497 # set lowest parts to 0
1498 while ($dst >= 0) { $x->[$dst--] = 0; }
1499 # fix spurious last zero element
1500 splice @$x,-1 if $x->[-1] == 0;
1507 # ref to array, ref to array, return ref to array
1508 my ($c,$cx,$cy) = @_;
1510 if (scalar @$cy == 1 && $cy->[0] == 0)
1512 splice (@$cx,1); $cx->[0] = 1; # y == 0 => x => 1
1515 if ((scalar @$cx == 1 && $cx->[0] == 1) || # x == 1
1516 (scalar @$cy == 1 && $cy->[0] == 1)) # or y == 1
1520 if (scalar @$cx == 1 && $cx->[0] == 0)
1522 splice (@$cx,1); $cx->[0] = 0; # 0 ** y => 0 (if not y <= 0)
1528 my $y_bin = _as_bin($c,$cy); $y_bin =~ s/^0b//;
1529 my $len = length($y_bin);
1532 _mul($c,$pow2,$cx) if substr($y_bin,$len,1) eq '1'; # is odd?
1541 # Return binomial coefficient (n over k).
1542 # Given refs to arrays, return ref to array.
1543 # First input argument is modified.
1545 my ($c, $n, $k) = @_;
1547 # If k > n/2, or, equivalently, 2*k > n, compute nok(n, k) as
1548 # nok(n, n-k), to minimize the number if iterations in the loop.
1551 my $twok = _mul($c, _two($c), _copy($c, $k)); # 2 * k
1552 if (_acmp($c, $twok, $n) > 0) { # if 2*k > n
1553 $k = _sub($c, _copy($c, $n), $k); # k = n - k
1559 # / 7 \ 7! 1*2*3*4 * 5*6*7 5 * 6 * 7 6 7
1560 # | | = --------- = --------------- = --------- = 5 * - * -
1561 # \ 3 / (7-3)! 3! 1*2*3*4 * 1*2*3 1 * 2 * 3 2 3
1563 if (_is_zero($c, $k)) {
1569 # Make a copy of the original n, since we'll be modifying n in-place.
1571 my $n_orig = _copy($c, $n);
1573 # n = 5, f = 6, d = 2 (cf. example above)
1578 my $f = _copy($c, $n);
1583 # while f <= n (the original n, that is) ...
1585 while (_acmp($c, $f, $n_orig) <= 0) {
1587 # n = (n * f / d) == 5 * 6 / 2 (cf. example above)
1592 # f = 7, d = 3 (cf. example above)
1617 # ref to array, return ref to array
1620 if ((@$cx == 1) && ($cx->[0] <= 7))
1622 $cx->[0] = $factorials[$cx->[0]]; # 0 => 1, 1 => 1, 2 => 2 etc.
1626 if ((@$cx == 1) && # we do this only if $x >= 12 and $x <= 7000
1627 ($cx->[0] >= 12 && $cx->[0] < 7000))
1630 # Calculate (k-j) * (k-j+1) ... k .. (k+j-1) * (k + j)
1631 # See http://blogten.blogspot.com/2007/01/calculating-n.html
1632 # The above series can be expressed as factors:
1633 # k * k - (j - i) * 2
1634 # We cache k*k, and calculate (j * j) as the sum of the first j odd integers
1636 # This will not work when N exceeds the storage of a Perl scalar, however,
1637 # in this case the algorithm would be way to slow to terminate, anyway.
1639 # As soon as the last element of $cx is 0, we split it up and remember
1640 # how many zeors we got so far. The reason is that n! will accumulate
1641 # zeros at the end rather fast.
1642 my $zero_elements = 0;
1644 # If n is even, set n = n -1
1645 my $k = _num($c,$cx); my $even = 1;
1650 # set k to the center point
1652 # print "k $k even: $even\n";
1653 # now calculate k * k
1655 my $odd = 1; my $sum = 1;
1657 # keep reference to x
1658 my $new_x = _new($c, $k * $even);
1662 $zero_elements ++; shift @$cx;
1664 # print STDERR "x = ", _str($c,$cx),"\n";
1665 my $BASE2 = int(sqrt($BASE))-1;
1669 my $m = ($k2 - $sum); $odd += 2; $sum += $odd; $j++;
1670 while ($j <= $i && ($m < $BASE2) && (($k2 - $sum) < $BASE2))
1673 $odd += 2; $sum += $odd; $j++;
1674 # print STDERR "\n k2 $k2 m $m sum $sum odd $odd\n"; sleep(1);
1682 _mul($c,$cx,$c->_new($m));
1686 $zero_elements ++; shift @$cx;
1688 # print STDERR "Calculate $k2 - $sum = $m (x = ", _str($c,$cx),")\n";
1690 # multiply in the zeros again
1691 unshift @$cx, (0) x $zero_elements;
1695 # go forward until $base is exceeded
1696 # limit is either $x steps (steps == 100 means a result always too high) or
1698 my $steps = 100; $steps = $cx->[0] if @$cx == 1;
1699 my $r = 2; my $cf = 3; my $step = 2; my $last = $r;
1700 while ($r*$cf < $BASE && $step < $steps)
1702 $last = $r; $r *= $cf++; $step++;
1704 if ((@$cx == 1) && $step == $cx->[0])
1706 # completely done, so keep reference to $x and return
1711 # now we must do the left over steps
1712 my $n; # steps still to do
1713 if (scalar @$cx == 1)
1722 # Set $cx to the last result below $BASE (but keep ref to $x)
1723 $cx->[0] = $last; splice (@$cx,1);
1724 # As soon as the last element of $cx is 0, we split it up and remember
1725 # how many zeors we got so far. The reason is that n! will accumulate
1726 # zeros at the end rather fast.
1727 my $zero_elements = 0;
1729 # do left-over steps fit into a scalar?
1730 if (ref $n eq 'ARRAY')
1732 # No, so use slower inc() & cmp()
1733 # ($n is at least $BASE here)
1734 my $base_2 = int(sqrt($BASE)) - 1;
1735 #print STDERR "base_2: $base_2\n";
1736 while ($step < $base_2)
1740 $zero_elements ++; shift @$cx;
1742 my $b = $step * ($step + 1); $step += 2;
1746 while (_acmp($c,$step,$n) <= 0)
1750 $zero_elements ++; shift @$cx;
1752 _mul($c,$cx,$step); _inc($c,$step);
1757 # Yes, so we can speed it up slightly
1759 # print "# left over steps $n\n";
1761 my $base_4 = int(sqrt(sqrt($BASE))) - 2;
1762 #print STDERR "base_4: $base_4\n";
1764 while ($step < $n4 && $step < $base_4)
1768 $zero_elements ++; shift @$cx;
1770 my $b = $step * ($step + 1); $step += 2; $b *= $step * ($step + 1); $step += 2;
1773 my $base_2 = int(sqrt($BASE)) - 1;
1775 #print STDERR "base_2: $base_2\n";
1776 while ($step < $n2 && $step < $base_2)
1780 $zero_elements ++; shift @$cx;
1782 my $b = $step * ($step + 1); $step += 2;
1785 # do what's left over
1788 _mul($c,$cx,[$step]); $step++;
1791 $zero_elements ++; shift @$cx;
1795 # multiply in the zeros again
1796 unshift @$cx, (0) x $zero_elements;
1797 $cx; # return result
1800 #############################################################################
1804 # calculate integer log of $x to base $base
1805 # ref to array, ref to array - return ref to array
1806 my ($c,$x,$base) = @_;
1809 return if (scalar @$x == 1 && $x->[0] == 0);
1810 # BASE 0 or 1 => NaN
1811 return if (scalar @$base == 1 && $base->[0] < 2);
1812 my $cmp = _acmp($c,$x,$base); # X == BASE => 1
1815 splice (@$x,1); $x->[0] = 1;
1821 splice (@$x,1); $x->[0] = 0;
1825 my $x_org = _copy($c,$x); # preserve x
1826 splice(@$x,1); $x->[0] = 1; # keep ref to $x
1828 # Compute a guess for the result based on:
1829 # $guess = int ( length_in_base_10(X) / ( log(base) / log(10) ) )
1830 my $len = _len($c,$x_org);
1831 my $log = log($base->[-1]) / log(10);
1833 # for each additional element in $base, we add $BASE_LEN to the result,
1834 # based on the observation that log($BASE,10) is BASE_LEN and
1835 # log(x*y) == log(x) + log(y):
1836 $log += ((scalar @$base)-1) * $BASE_LEN;
1838 # calculate now a guess based on the values obtained above:
1839 my $res = int($len / $log);
1842 my $trial = _pow ($c, _copy($c, $base), $x);
1843 my $a = _acmp($c,$trial,$x_org);
1845 # print STDERR "# trial ", _str($c,$x)," was: $a (0 = exact, -1 too small, +1 too big)\n";
1847 # found an exact result?
1848 return ($x,1) if $a == 0;
1853 _div($c,$trial,$base); _dec($c, $x);
1854 while (($a = _acmp($c,$trial,$x_org)) > 0)
1856 # print STDERR "# big _log_int at ", _str($c,$x), "\n";
1857 _div($c,$trial,$base); _dec($c, $x);
1859 # result is now exact (a == 0), or too small (a < 0)
1860 return ($x, $a == 0 ? 1 : 0);
1863 # else: result was to small
1864 _mul($c,$trial,$base);
1866 # did we now get the right result?
1867 $a = _acmp($c,$trial,$x_org);
1869 if ($a == 0) # yes, exactly
1874 return ($x,0) if $a > 0;
1876 # Result still too small (we should come here only if the estimate above
1877 # was very off base):
1879 # Now let the normal trial run obtain the real result
1880 # Simple loop that increments $x by 2 in each step, possible overstepping
1883 my $base_mul = _mul($c, _copy($c,$base), $base); # $base * $base
1885 while (($a = _acmp($c,$trial,$x_org)) < 0)
1887 # print STDERR "# small _log_int at ", _str($c,$x), "\n";
1888 _mul($c,$trial,$base_mul); _add($c, $x, [2]);
1894 # overstepped the result
1896 _div($c,$trial,$base);
1897 $a = _acmp($c,$trial,$x_org);
1902 $exact = 0 if $a != 0; # a = -1 => not exact result, a = 0 => exact
1905 ($x,$exact); # return result
1909 use constant DEBUG => 0;
1911 sub steps { $steps };
1915 # square-root of $x in place
1916 # Compute a guess of the result (by rule of thumb), then improve it via
1920 if (scalar @$x == 1)
1922 # fits into one Perl scalar, so result can be computed directly
1923 $x->[0] = int(sqrt($x->[0]));
1926 my $y = _copy($c,$x);
1927 # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
1928 # since our guess will "grow"
1929 my $l = int((_len($c,$x)-1) / 2);
1931 my $lastelem = $x->[-1]; # for guess
1932 my $elems = scalar @$x - 1;
1933 # not enough digits, but could have more?
1934 if ((length($lastelem) <= 3) && ($elems > 1))
1936 # right-align with zero pad
1937 my $len = length($lastelem) & 1;
1938 print "$lastelem => " if DEBUG;
1939 $lastelem .= substr($x->[-2] . '0' x $BASE_LEN,0,$BASE_LEN);
1940 # former odd => make odd again, or former even to even again
1941 $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
1942 print "$lastelem\n" if DEBUG;
1945 # construct $x (instead of _lsft($c,$x,$l,10)
1946 my $r = $l % $BASE_LEN; # 10000 00000 00000 00000 ($BASE_LEN=5)
1947 $l = int($l / $BASE_LEN);
1948 print "l = $l " if DEBUG;
1950 splice @$x,$l; # keep ref($x), but modify it
1952 # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
1954 # 14400 00000 => sqrt(14400) => guess first digits to be 120
1955 # 144000 000000 => sqrt(144000) => guess 379
1957 print "$lastelem (elems $elems) => " if DEBUG;
1958 $lastelem = $lastelem / 10 if ($elems & 1 == 1); # odd or even?
1959 my $g = sqrt($lastelem); $g =~ s/\.//; # 2.345 => 2345
1960 $r -= 1 if $elems & 1 == 0; # 70 => 7
1962 # padd with zeros if result is too short
1963 $x->[$l--] = int(substr($g . '0' x $r,0,$r+1));
1964 print "now ",$x->[-1] if DEBUG;
1965 print " would have been ", int('1' . '0' x $r),"\n" if DEBUG;
1967 # If @$x > 1, we could compute the second elem of the guess, too, to create
1968 # an even better guess. Not implemented yet. Does it improve performance?
1969 $x->[$l--] = 0 while ($l >= 0); # all other digits of guess are zero
1971 print "start x= ",_str($c,$x),"\n" if DEBUG;
1974 my $lastlast = _zero();
1975 $steps = 0 if DEBUG;
1976 while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
1979 $lastlast = _copy($c,$last);
1980 $last = _copy($c,$x);
1981 _add($c,$x, _div($c,_copy($c,$y),$x));
1983 print " x= ",_str($c,$x),"\n" if DEBUG;
1985 print "\nsteps in sqrt: $steps, " if DEBUG;
1986 _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0; # overshot?
1987 print " final ",$x->[-1],"\n" if DEBUG;
1993 # take n'th root of $x in place (n >= 3)
1996 if (scalar @$x == 1)
2000 # result will always be smaller than 2 so trunc to 1 at once
2005 # fits into one Perl scalar, so result can be computed directly
2006 # cannot use int() here, because it rounds wrongly (try
2007 # (81 ** 3) ** (1/3) to see what I mean)
2008 #$x->[0] = int( $x->[0] ** (1 / $n->[0]) );
2009 # round to 8 digits, then truncate result to integer
2010 $x->[0] = int ( sprintf ("%.8f", $x->[0] ** (1 / $n->[0]) ) );
2015 # we know now that X is more than one element long
2017 # if $n is a power of two, we can repeatedly take sqrt($X) and find the
2018 # proper result, because sqrt(sqrt($x)) == root($x,4)
2019 my $b = _as_bin($c,$n);
2020 if ($b =~ /0b1(0+)$/)
2022 my $count = CORE::length($1); # 0b100 => len('00') => 2
2023 my $cnt = $count; # counter for loop
2024 unshift (@$x, 0); # add one element, together with one
2025 # more below in the loop this makes 2
2028 # 'inflate' $X by adding one element, basically computing
2029 # $x * $BASE * $BASE. This gives us more $BASE_LEN digits for result
2030 # since len(sqrt($X)) approx == len($x) / 2.
2032 # calculate sqrt($x), $x is now one element to big, again. In the next
2033 # round we make that two, again.
2036 # $x is now one element to big, so truncate result by removing it
2041 # trial computation by starting with 2,4,8,16 etc until we overstep
2045 # while still to do more than X steps
2049 while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
2051 _mul ($c, $step, [2]);
2052 _add ($c, $trial, $step);
2056 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) == 0)
2058 @$x = @$trial; # make copy while preserving ref to $x
2061 # overstepped, so go back on step
2062 _sub($c, $trial, $step);
2063 } while (scalar @$step > 1 || $step->[0] > 128);
2067 # add two, because $trial cannot be exactly the result (otherwise we would
2068 # already have found it)
2069 _add($c, $trial, $step);
2071 # and now add more and more (2,4,6,8,10 etc)
2072 while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
2074 _add ($c, $trial, $step);
2077 # hit not exactly? (overstepped)
2078 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
2083 # hit not exactly? (overstepped)
2084 # 80 too small, 81 slightly too big, 82 too big
2085 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
2090 @$x = @$trial; # make copy while preserving ref to $x
2096 ##############################################################################
2103 # the shortcut makes equal, large numbers _really_ fast, and makes only a
2104 # very small performance drop for small numbers (e.g. something with less
2105 # than 32 bit) Since we optimize for large numbers, this is enabled.
2106 return $x if _acmp($c,$x,$y) == 0; # shortcut
2108 my $m = _one(); my ($xr,$yr);
2109 my $mask = $AND_MASK;
2112 my $y1 = _copy($c,$y); # make copy
2116 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2118 ($x1, $xr) = _div($c,$x1,$mask);
2119 ($y1, $yr) = _div($c,$y1,$mask);
2121 # make ints() from $xr, $yr
2122 # this is when the AND_BITS are greater than $BASE and is slower for
2123 # small (<256 bits) numbers, but faster for large numbers. Disabled
2124 # due to KISS principle
2126 # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2127 # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2128 # _add($c,$x, _mul($c, _new( $c, ($xrr & $yrr) ), $m) );
2130 # 0+ due to '&' doesn't work in strings
2131 _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
2141 return _zero() if _acmp($c,$x,$y) == 0; # shortcut (see -and)
2143 my $m = _one(); my ($xr,$yr);
2144 my $mask = $XOR_MASK;
2147 my $y1 = _copy($c,$y); # make copy
2151 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2153 ($x1, $xr) = _div($c,$x1,$mask);
2154 ($y1, $yr) = _div($c,$y1,$mask);
2155 # make ints() from $xr, $yr (see _and())
2156 #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2157 #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2158 #_add($c,$x, _mul($c, _new( $c, ($xrr ^ $yrr) ), $m) );
2160 # 0+ due to '^' doesn't work in strings
2161 _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
2164 # the loop stops when the shorter of the two numbers is exhausted
2165 # the remainder of the longer one will survive bit-by-bit, so we simple
2166 # multiply-add it in
2167 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
2168 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
2177 return $x if _acmp($c,$x,$y) == 0; # shortcut (see _and)
2179 my $m = _one(); my ($xr,$yr);
2180 my $mask = $OR_MASK;
2183 my $y1 = _copy($c,$y); # make copy
2187 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2189 ($x1, $xr) = _div($c,$x1,$mask);
2190 ($y1, $yr) = _div($c,$y1,$mask);
2191 # make ints() from $xr, $yr (see _and())
2192 # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2193 # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2194 # _add($c,$x, _mul($c, _new( $c, ($xrr | $yrr) ), $m) );
2196 # 0+ due to '|' doesn't work in strings
2197 _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
2200 # the loop stops when the shorter of the two numbers is exhausted
2201 # the remainder of the longer one will survive bit-by-bit, so we simple
2202 # multiply-add it in
2203 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
2204 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
2211 # convert a decimal number to hex (ref to array, return ref to string)
2214 # fits into one element (handle also 0x0 case)
2215 return sprintf("0x%x",$x->[0]) if @$x == 1;
2217 my $x1 = _copy($c,$x);
2220 my ($xr, $h, $x10000);
2223 $x10000 = [ 0x10000 ]; $h = 'h4';
2227 $x10000 = [ 0x1000 ]; $h = 'h3';
2229 while (@$x1 != 1 || $x1->[0] != 0) # _is_zero()
2231 ($x1, $xr) = _div($c,$x1,$x10000);
2232 $es .= unpack($h,pack('V',$xr->[0]));
2235 $es =~ s/^[0]+//; # strip leading zeros
2236 '0x' . $es; # return result prepended with 0x
2241 # convert a decimal number to bin (ref to array, return ref to string)
2244 # fits into one element (and Perl recent enough), handle also 0b0 case
2245 # handle zero case for older Perls
2246 if ($] <= 5.005 && @$x == 1 && $x->[0] == 0)
2248 my $t = '0b0'; return $t;
2250 if (@$x == 1 && $] >= 5.006)
2252 my $t = sprintf("0b%b",$x->[0]);
2255 my $x1 = _copy($c,$x);
2258 my ($xr, $b, $x10000);
2261 $x10000 = [ 0x10000 ]; $b = 'b16';
2265 $x10000 = [ 0x1000 ]; $b = 'b12';
2267 while (!(@$x1 == 1 && $x1->[0] == 0)) # _is_zero()
2269 ($x1, $xr) = _div($c,$x1,$x10000);
2270 $es .= unpack($b,pack('v',$xr->[0]));
2273 $es =~ s/^[0]+//; # strip leading zeros
2274 '0b' . $es; # return result prepended with 0b
2279 # convert a decimal number to octal (ref to array, return ref to string)
2282 # fits into one element (handle also 0 case)
2283 return sprintf("0%o",$x->[0]) if @$x == 1;
2285 my $x1 = _copy($c,$x);
2289 my $x1000 = [ 0100000 ];
2290 while (@$x1 != 1 || $x1->[0] != 0) # _is_zero()
2292 ($x1, $xr) = _div($c,$x1,$x1000);
2293 $es .= reverse sprintf("%05o", $xr->[0]);
2296 $es =~ s/^[0]+//; # strip leading zeros
2297 '0' . $es; # return result prepended with 0
2302 # convert a octal number to decimal (string, return ref to array)
2305 # for older Perls, play safe
2306 my $m = [ 0100000 ];
2307 my $d = 5; # 5 digits at a time
2312 my $len = int( (length($os)-1)/$d ); # $d digit parts, w/o the '0'
2313 my $val; my $i = -$d;
2316 $val = substr($os,$i,$d); # get oct digits
2317 $val = CORE::oct($val);
2319 my $adder = [ $val ];
2320 _add ($c, $x, _mul ($c, $adder, $mul ) ) if $val != 0;
2321 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
2328 # convert a hex number to decimal (string, return ref to array)
2331 my $m = _new($c, 0x10000000); # 28 bit at a time (<32 bit!)
2332 my $d = 7; # 7 digits at a time
2335 # for older Perls, play safe
2336 $m = [ 0x10000 ]; # 16 bit at a time (<32 bit!)
2337 $d = 4; # 4 digits at a time
2343 my $len = int( (length($hs)-2)/$d ); # $d digit parts, w/o the '0x'
2344 my $val; my $i = -$d;
2347 $val = substr($hs,$i,$d); # get hex digits
2348 $val =~ s/^0x// if $len == 0; # for last part only because
2349 $val = CORE::hex($val); # hex does not like wrong chars
2351 my $adder = [ $val ];
2352 # if the resulting number was to big to fit into one element, create a
2353 # two-element version (bug found by Mark Lakata - Thanx!)
2354 if (CORE::length($val) > $BASE_LEN)
2356 $adder = _new($c,$val);
2358 _add ($c, $x, _mul ($c, $adder, $mul ) ) if $val != 0;
2359 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
2366 # convert a hex number to decimal (string, return ref to array)
2369 # instead of converting X (8) bit at a time, it is faster to "convert" the
2370 # number to hex, and then call _from_hex.
2373 $hs =~ s/^[+-]?0b//; # remove sign and 0b
2374 my $l = length($hs); # bits
2375 $hs = '0' x (8-($l % 8)) . $hs if ($l % 8) != 0; # padd left side w/ 0
2376 my $h = '0x' . unpack('H*', pack ('B*', $hs)); # repack as hex
2381 ##############################################################################
2382 # special modulus functions
2386 # modular multiplicative inverse
2390 if (_is_zero($c, $y)) {
2391 return (undef, undef);
2395 if (_is_one($c, $y)) {
2396 return (_zero($c), '+');
2401 my $a = _copy($c,$y);
2402 my $b = _copy($c,$x);
2404 # Euclid's Algorithm for bgcd(), only that we calc bgcd() ($a) and the result
2405 # ($u) at the same time. See comments in BigInt for why this works.
2409 ($a, $q, $b) = ($b, _div($c, $a, $b)); # step 1
2410 last if _is_zero($c, $b);
2412 my $t = _add($c, # step 2:
2413 _mul($c, _copy($c, $v), $q) , # t = v * q
2421 # if the gcd is not 1, then return NaN
2422 return (undef, undef) unless _is_one($c, $a);
2424 ($v, $sign == 1 ? '+' : '-');
2429 # modulus of power ($x ** $y) % $z
2430 my ($c,$num,$exp,$mod) = @_;
2432 # a^b (mod 1) = 0 for all a and b
2433 if (_is_one($c,$mod))
2439 # 0^a (mod m) = 0 if m != 0, a != 0
2440 # 0^0 (mod m) = 1 if m != 0
2441 if (_is_zero($c, $num)) {
2442 if (_is_zero($c, $exp)) {
2450 # $num = _mod($c,$num,$mod); # this does not make it faster
2452 my $acc = _copy($c,$num); my $t = _one();
2454 my $expbin = _as_bin($c,$exp); $expbin =~ s/^0b//;
2455 my $len = length($expbin);
2458 if ( substr($expbin,$len,1) eq '1') # is_odd
2461 $t = _mod($c,$t,$mod);
2464 $acc = _mod($c,$acc,$mod);
2471 # Greatest common divisor.
2473 my ($c, $x, $y) = @_;
2476 # gcd(0,a) = a, if a != 0
2478 if (@$x == 1 && $x->[0] == 0) {
2479 if (@$y == 1 && $y->[0] == 0) {
2487 # Until $y is zero ...
2489 until (@$y == 1 && $y->[0] == 0) {
2491 # Compute remainder.
2499 $y = $tmp; # no deref here; that would modify input $y
2505 ##############################################################################
2506 ##############################################################################
2515 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
2519 This library provides support for big integer calculations. It is not
2520 intended to be used by other modules. Other modules which support the same
2521 API (see below) can also be used to support Math::BigInt, like
2522 Math::BigInt::GMP and Math::BigInt::Pari.
2526 In this library, the numbers are represented in base B = 10**N, where N is
2527 the largest possible value that does not cause overflow in the intermediate
2528 computations. The base B elements are stored in an array, with the least
2529 significant element stored in array element zero. There are no leading zero
2530 elements, except a single zero element when the number is zero.
2532 For instance, if B = 10000, the number 1234567890 is represented internally
2533 as [3456, 7890, 12].
2535 =head1 THE Math::BigInt API
2537 In order to allow for multiple big integer libraries, Math::BigInt was
2538 rewritten to use a plug-in library for core math routines. Any module which
2539 conforms to the API can be used by Math::BigInt by using this in your program:
2541 use Math::BigInt lib => 'libname';
2543 'libname' is either the long name, like 'Math::BigInt::Pari', or only the short
2544 version, like 'Pari'.
2546 =head2 General Notes
2548 A library only needs to deal with unsigned big integers. Testing of input
2549 parameter validity is done by the caller, so there is no need to worry about
2550 underflow (e.g., in C<_sub()> and C<_dec()>) nor about division by zero (e.g.,
2551 in C<_div()>) or similar cases.
2553 For some methods, the first parameter can be modified. That includes the
2554 possibility that you return a reference to a completely different object
2555 instead. Although keeping the reference and just changing its contents is
2556 preferred over creating and returning a different reference.
2558 Return values are always objects, strings, Perl scalars, or true/false for
2559 comparison routines.
2561 =head2 API version 1
2563 The following methods must be defined in order to support the use by
2564 Math::BigInt v1.70 or later.
2570 =item I<api_version()>
2572 Return API version as a Perl scalar, 1 for Math::BigInt v1.70, 2 for
2583 Convert a string representing an unsigned decimal number to an object
2584 representing the same number. The input is normalize, i.e., it matches
2589 Return an object representing the number zero.
2593 Return an object representing the number one.
2597 Return an object representing the number two.
2601 Return an object representing the number ten.
2603 =item I<_from_bin(STR)>
2605 Return an object given a string representing a binary number. The input has a
2606 '0b' prefix and matches the regular expression C<^0[bB](0|1[01]*)$>.
2608 =item I<_from_oct(STR)>
2610 Return an object given a string representing an octal number. The input has a
2611 '0' prefix and matches the regular expression C<^0[1-7]*$>.
2613 =item I<_from_hex(STR)>
2615 Return an object given a string representing a hexadecimal number. The input
2616 has a '0x' prefix and matches the regular expression
2617 C<^0x(0|[1-9a-fA-F][\da-fA-F]*)$>.
2621 =head3 Mathematical functions
2623 Each of these methods may modify the first input argument, except I<_bgcd()>,
2624 which shall not modify any input argument, and I<_sub()> which may modify the
2625 second input argument.
2629 =item I<_add(OBJ1, OBJ2)>
2631 Returns the result of adding OBJ2 to OBJ1.
2633 =item I<_mul(OBJ1, OBJ2)>
2635 Returns the result of multiplying OBJ2 and OBJ1.
2637 =item I<_div(OBJ1, OBJ2)>
2639 Returns the result of dividing OBJ1 by OBJ2 and truncating the result to an
2642 =item I<_sub(OBJ1, OBJ2, FLAG)>
2644 =item I<_sub(OBJ1, OBJ2)>
2646 Returns the result of subtracting OBJ2 by OBJ1. If C<flag> is false or omitted,
2647 OBJ1 might be modified. If C<flag> is true, OBJ2 might be modified.
2651 Decrement OBJ by one.
2655 Increment OBJ by one.
2657 =item I<_mod(OBJ1, OBJ2)>
2659 Return OBJ1 modulo OBJ2, i.e., the remainder after dividing OBJ1 by OBJ2.
2663 Return the square root of the object, truncated to integer.
2665 =item I<_root(OBJ, N)>
2667 Return Nth root of the object, truncated to int. N is E<gt>= 3.
2671 Return factorial of object (1*2*3*4*...).
2673 =item I<_pow(OBJ1, OBJ2)>
2675 Return OBJ1 to the power of OBJ2. By convention, 0**0 = 1.
2677 =item I<_modinv(OBJ1, OBJ2)>
2679 Return modular multiplicative inverse, i.e., return OBJ3 so that
2681 (OBJ3 * OBJ1) % OBJ2 = 1 % OBJ2
2683 The result is returned as two arguments. If the modular multiplicative
2684 inverse does not exist, both arguments are undefined. Otherwise, the
2685 arguments are a number (object) and its sign ("+" or "-").
2687 The output value, with its sign, must either be a positive value in the
2688 range 1,2,...,OBJ2-1 or the same value subtracted OBJ2. For instance, if the
2689 input arguments are objects representing the numbers 7 and 5, the method
2690 must either return an object representing the number 3 and a "+" sign, since
2691 (3*7) % 5 = 1 % 5, or an object representing the number 2 and "-" sign,
2692 since (-2*7) % 5 = 1 % 5.
2694 =item I<_modpow(OBJ1, OBJ2, OBJ3)>
2696 Return modular exponentiation, (OBJ1 ** OBJ2) % OBJ3.
2698 =item I<_rsft(OBJ, N, B)>
2700 Shift object N digits right in base B and return the resulting object. This is
2701 equivalent to performing integer division by B**N and discarding the remainder,
2702 except that it might be much faster, depending on how the number is represented
2705 For instance, if the object $obj represents the hexadecimal number 0xabcde,
2706 then C<_rsft($obj, 2, 16)> returns an object representing the number 0xabc. The
2707 "remainer", 0xde, is discarded and not returned.
2709 =item I<_lsft(OBJ, N, B)>
2711 Shift the object N digits left in base B. This is equivalent to multiplying by
2712 B**N, except that it might be much faster, depending on how the number is
2713 represented internally.
2715 =item I<_log_int(OBJ, B)>
2717 Return integer log of OBJ to base BASE. This method has two output arguments,
2718 the OBJECT and a STATUS. The STATUS is Perl scalar; it is 1 if OBJ is the exact
2719 result, 0 if the result was truncted to give OBJ, and undef if it is unknown
2720 whether OBJ is the exact result.
2722 =item I<_gcd(OBJ1, OBJ2)>
2724 Return the greatest common divisor of OBJ1 and OBJ2.
2728 =head3 Bitwise operators
2730 Each of these methods may modify the first input argument.
2734 =item I<_and(OBJ1, OBJ2)>
2736 Return bitwise and. If necessary, the smallest number is padded with leading
2739 =item I<_or(OBJ1, OBJ2)>
2741 Return bitwise or. If necessary, the smallest number is padded with leading
2744 =item I<_xor(OBJ1, OBJ2)>
2746 Return bitwise exclusive or. If necessary, the smallest number is padded
2751 =head3 Boolean operators
2755 =item I<_is_zero(OBJ)>
2757 Returns a true value if OBJ is zero, and false value otherwise.
2759 =item I<_is_one(OBJ)>
2761 Returns a true value if OBJ is one, and false value otherwise.
2763 =item I<_is_two(OBJ)>
2765 Returns a true value if OBJ is two, and false value otherwise.
2767 =item I<_is_ten(OBJ)>
2769 Returns a true value if OBJ is ten, and false value otherwise.
2771 =item I<_is_even(OBJ)>
2773 Return a true value if OBJ is an even integer, and a false value otherwise.
2775 =item I<_is_odd(OBJ)>
2777 Return a true value if OBJ is an even integer, and a false value otherwise.
2779 =item I<_acmp(OBJ1, OBJ2)>
2781 Compare OBJ1 and OBJ2 and return -1, 0, or 1, if OBJ1 is less than, equal
2782 to, or larger than OBJ2, respectively.
2786 =head3 String conversion
2792 Return a string representing the object. The returned string should have no
2793 leading zeros, i.e., it should match C<^(0|[1-9]\d*)$>.
2795 =item I<_as_bin(OBJ)>
2797 Return the binary string representation of the number. The string must have a
2800 =item I<_as_oct(OBJ)>
2802 Return the octal string representation of the number. The string must have
2805 Note: This method was required from Math::BigInt version 1.78, but the required
2806 API version number was not incremented, so there are older libraries that
2807 support API version 1, but do not support C<_as_oct()>.
2809 =item I<_as_hex(OBJ)>
2811 Return the hexadecimal string representation of the number. The string must
2816 =head3 Numeric conversion
2822 Given an object, return a Perl scalar number (int/float) representing this
2827 =head3 Miscellaneous
2833 Return a true copy of the object.
2837 Returns the number of the decimal digits in the number. The output is a
2840 =item I<_zeros(OBJ)>
2842 Return the number of trailing decimal zeros. The output is a Perl scalar.
2844 =item I<_digit(OBJ, N)>
2846 Return the Nth digit as a Perl scalar. N is a Perl scalar, where zero refers to
2847 the rightmost (least significant) digit, and negative values count from the
2848 left (most significant digit). If $obj represents the number 123, then
2849 I<_digit($obj, 0)> is 3 and I<_digit(123, -1)> is 1.
2851 =item I<_check(OBJ)>
2853 Return a true value if the object is OK, and a false value otherwise. This is a
2854 check routine to test the internal state of the object for corruption.
2858 =head2 API version 2
2860 The following methods are required for an API version of 2 or greater.
2868 Return an object representing the number 10**N where N E<gt>= 0 is a Perl
2873 =head3 Mathematical functions
2877 =item I<_nok(OBJ1, OBJ2)>
2879 Return the binomial coefficient OBJ1 over OBJ1.
2883 =head3 Miscellaneous
2889 Return the approximate number of decimal digits of the object. The
2890 output is one Perl scalar. This estimate must be greater than or equal
2891 to what C<_len()> returns.
2895 =head2 API optional methods
2897 The following methods are optional, and can be defined if the underlying lib
2898 has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
2899 slow) fallback routines to emulate these:
2901 =head3 Signed bitwise operators.
2903 Each of these methods may modify the first input argument.
2907 =item I<_signed_or(OBJ1, OBJ2, SIGN1, SIGN2)>
2909 Return the signed bitwise or.
2911 =item I<_signed_and(OBJ1, OBJ2, SIGN1, SIGN2)>
2913 Return the signed bitwise and.
2915 =item I<_signed_xor(OBJ1, OBJ2, SIGN1, SIGN2)>
2917 Return the signed bitwise exclusive or.
2921 =head1 WRAP YOUR OWN
2923 If you want to port your own favourite c-lib for big numbers to the
2924 Math::BigInt interface, you can take any of the already existing modules as
2925 a rough guideline. You should really wrap up the latest BigInt and BigFloat
2926 testsuites with your module, and replace in them any of the following:
2932 use Math::BigInt lib => 'yourlib';
2934 This way you ensure that your library really works 100% within Math::BigInt.
2938 Please report any bugs or feature requests to
2939 C<bug-math-bigint at rt.cpan.org>, or through the web interface at
2940 L<https://rt.cpan.org/Ticket/Create.html?Queue=Math-BigInt>
2942 We will be notified, and then you'll automatically be notified of progress on
2943 your bug as I make changes.
2947 You can find documentation for this module with the perldoc command.
2949 perldoc Math::BigInt::Calc
2951 You can also look for information at:
2955 =item * RT: CPAN's request tracker
2957 L<https://rt.cpan.org/Public/Dist/Display.html?Name=Math-BigInt>
2959 =item * AnnoCPAN: Annotated CPAN documentation
2961 L<http://annocpan.org/dist/Math-BigInt>
2963 =item * CPAN Ratings
2965 L<http://cpanratings.perl.org/dist/Math-BigInt>
2969 L<http://search.cpan.org/dist/Math-BigInt/>
2971 =item * CPAN Testers Matrix
2973 L<http://matrix.cpantesters.org/?dist=Math-BigInt>
2975 =item * The Bignum mailing list
2979 =item * Post to mailing list
2981 C<bignum at lists.scsys.co.uk>
2983 =item * View mailing list
2985 L<http://lists.scsys.co.uk/pipermail/bignum/>
2987 =item * Subscribe/Unsubscribe
2989 L<http://lists.scsys.co.uk/cgi-bin/mailman/listinfo/bignum>
2997 This program is free software; you may redistribute it and/or modify it under
2998 the same terms as Perl itself.
3006 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
3011 Separated from BigInt and shaped API with the help of John Peacock.
3015 Fixed, speed-up, streamlined and enhanced by Tels 2001 - 2007.
3019 API documentation corrected and extended by Peter John Acklam,
3020 E<lt>pjacklam@online.noE<gt>
3026 L<Math::BigInt>, L<Math::BigFloat>,
3027 L<Math::BigInt::GMP>, L<Math::BigInt::FastCalc> and L<Math::BigInt::Pari>.