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 in div() (maybe - that code scares me...)
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 # announce that we are compatible with MBI v1.70 and up
36 sub api_version () { 1; }
38 # constants for easier life
39 my ($BASE,$BASE_LEN,$MBASE,$RBASE,$MAX_VAL,$BASE_LEN_SMALL);
40 my ($AND_BITS,$XOR_BITS,$OR_BITS);
41 my ($AND_MASK,$XOR_MASK,$OR_MASK);
45 # set/get the BASE_LEN and assorted other, connected values
46 # used only be the testsuite, set is used only by the BEGIN block below
52 # find whether we can use mul or div or none in mul()/div()
53 # (in last case reduce BASE_LEN_SMALL)
54 $BASE_LEN_SMALL = $b+1;
56 while (--$BASE_LEN_SMALL > 5)
58 $MBASE = int("1e".$BASE_LEN_SMALL);
59 $RBASE = abs('1e-'.$BASE_LEN_SMALL); # see USE_MUL
61 $caught += 1 if (int($MBASE * $RBASE) != 1); # should be 1
62 $caught += 2 if (int($MBASE / $MBASE) != 1); # should be 1
65 # BASE_LEN is used for anything else than mul()/div()
66 $BASE_LEN = $BASE_LEN_SMALL;
67 $BASE_LEN = shift if (defined $_[0]); # one more arg?
68 $BASE = int("1e".$BASE_LEN);
70 $MBASE = int("1e".$BASE_LEN_SMALL);
71 $RBASE = abs('1e-'.$BASE_LEN_SMALL); # see USE_MUL
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 # must USE_MUL since we cannot use DIV
86 *{_mul} = \&_mul_use_mul;
87 *{_div} = \&_div_use_mul;
92 *{_mul} = \&_mul_use_div;
93 *{_div} = \&_div_use_div;
96 return $BASE_LEN unless wantarray;
97 return ($BASE_LEN, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN_SMALL, $MAX_VAL, $BASE);
102 # (ref to string) return ref to num_array
103 # Convert a number from string format (without sign) to internal base
104 # 1ex format. Assumes normalized value as input.
105 my $il = length($_[1])-1;
107 # < BASE_LEN due len-1 above
108 return [ int($_[1]) ] if $il < $BASE_LEN; # shortcut for short numbers
110 # this leaves '00000' instead of int 0 and will be corrected after any op
111 [ reverse(unpack("a" . ($il % $BASE_LEN+1)
112 . ("a$BASE_LEN" x ($il / $BASE_LEN)), $_[1])) ];
117 # from Daniel Pfeiffer: determine largest group of digits that is precisely
118 # multipliable with itself plus carry
119 # Test now changed to expect the proper pattern, not a result off by 1 or 2
120 my ($e, $num) = 3; # lowest value we will use is 3+1-1 = 3
123 $num = ('9' x ++$e) + 0;
125 } while ("$num" =~ /9{$e}0{$e}/); # must be a certain pattern
126 $e--; # last test failed, so retract one step
127 # the limits below brush the problems with the test above under the rug:
128 # the test should be able to find the proper $e automatically
129 $e = 5 if $^O =~ /^uts/; # UTS get's some special treatment
130 $e = 5 if $^O =~ /^unicos/; # unicos is also problematic (6 seems to work
131 # there, but we play safe)
132 $e = 5 if $] < 5.006; # cap, for older Perls
133 $e = 7 if $e > 7; # cap, for VMS, OS/390 and other 64 bit systems
134 # 8 fails inside random testsuite, so take 7
136 __PACKAGE__->_base_len($e); # set and store
139 # find out how many bits _and, _or and _xor can take (old default = 16)
140 # I don't think anybody has yet 128 bit scalars, so let's play safe.
141 local $^W = 0; # don't warn about 'nonportable number'
142 $AND_BITS = 15; $XOR_BITS = 15; $OR_BITS = 15;
144 # find max bits, we will not go higher than numberofbits that fit into $BASE
145 # to make _and etc simpler (and faster for smaller, slower for large numbers)
147 while (2 ** $max < $BASE) { $max++; }
150 $max = 16 if $] < 5.006; # older Perls might not take >16 too well
155 $x = oct('0b' . '1' x $AND_BITS); $y = $x & $x;
156 $z = (2 ** $AND_BITS) - 1;
157 } while ($AND_BITS < $max && $x == $z && $y == $x);
158 $AND_BITS --; # retreat one step
161 $x = oct('0b' . '1' x $XOR_BITS); $y = $x ^ 0;
162 $z = (2 ** $XOR_BITS) - 1;
163 } while ($XOR_BITS < $max && $x == $z && $y == $x);
164 $XOR_BITS --; # retreat one step
167 $x = oct('0b' . '1' x $OR_BITS); $y = $x | $x;
168 $z = (2 ** $OR_BITS) - 1;
169 } while ($OR_BITS < $max && $x == $z && $y == $x);
170 $OR_BITS --; # retreat one step
172 $AND_MASK = __PACKAGE__->_new( ( 2 ** $AND_BITS ));
173 $XOR_MASK = __PACKAGE__->_new( ( 2 ** $XOR_BITS ));
174 $OR_MASK = __PACKAGE__->_new( ( 2 ** $OR_BITS ));
177 ###############################################################################
193 # create a two (used internally for shifting)
199 # create a 10 (used internally for shifting)
209 # catch and throw away
212 ##############################################################################
213 # convert back to string and number
217 # (ref to BINT) return num_str
218 # Convert number from internal base 100000 format to string format.
219 # internal format is always normalized (no leading zeros, "-0" => "+0")
222 my $l = scalar @$ar; # number of parts
223 if ($l < 1) # should not happen
226 Carp::croak("$_[1] has no elements");
230 # handle first one different to strip leading zeros from it (there are no
231 # leading zero parts in internal representation)
232 $l --; $ret .= int($ar->[$l]); $l--;
233 # Interestingly, the pre-padd method uses more time
234 # the old grep variant takes longer (14 vs. 10 sec)
235 my $z = '0' x ($BASE_LEN-1);
238 $ret .= substr($z.$ar->[$l],-$BASE_LEN); # fastest way I could think of
246 # Make a number (scalar int/float) from a BigInt object
249 return 0+$x->[0] if scalar @$x == 1; # below $BASE
254 $num += $fac*$_; $fac *= $BASE;
259 ##############################################################################
264 # (ref to int_num_array, ref to int_num_array)
265 # routine to add two base 1eX numbers
266 # stolen from Knuth Vol 2 Algorithm A pg 231
267 # there are separate routines to add and sub as per Knuth pg 233
268 # This routine clobbers up array x, but not y.
272 return $x if (@$y == 1) && $y->[0] == 0; # $x + 0 => $x
273 if ((@$x == 1) && $x->[0] == 0) # 0 + $y => $y->copy
275 # twice as slow as $x = [ @$y ], but necc. to retain $x as ref :(
276 @$x = @$y; return $x;
279 # for each in Y, add Y to X and carry. If after that, something is left in
280 # X, foreach in X add carry to X and then return X, carry
281 # Trades one "$j++" for having to shift arrays
282 my $i; my $car = 0; my $j = 0;
285 $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
290 $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
297 # (ref to int_num_array, ref to int_num_array)
298 # Add 1 to $x, modify $x in place
303 return $x if (($i += 1) < $BASE); # early out
304 $i = 0; # overflow, next
306 push @$x,1 if (($x->[-1] || 0) == 0); # last overflowed, so extend
312 # (ref to int_num_array, ref to int_num_array)
313 # Sub 1 from $x, modify $x in place
316 my $MAX = $BASE-1; # since MAX_VAL based on MBASE
319 last if (($i -= 1) >= 0); # early out
320 $i = $MAX; # underflow, next
322 pop @$x if $x->[-1] == 0 && @$x > 1; # last underflowed (but leave 0)
328 # (ref to int_num_array, ref to int_num_array, swap)
329 # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
330 # subtract Y from X by modifying x in place
331 my ($c,$sx,$sy,$s) = @_;
333 my $car = 0; my $i; my $j = 0;
338 last unless defined $sy->[$j] || $car;
339 $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
341 # might leave leading zeros, so fix that
342 return __strip_zeros($sx);
346 # we can't do an early out if $x is < than $y, since we
347 # need to copy the high chunks from $y. Found by Bob Mathews.
348 #last unless defined $sy->[$j] || $car;
350 if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
353 # might leave leading zeros, so fix that
359 # (ref to int_num_array, ref to int_num_array)
360 # multiply two numbers in internal representation
361 # modifies first arg, second need not be different from first
362 my ($c,$xv,$yv) = @_;
366 # shortcut for two very short numbers (improved by Nathan Zook)
367 # works also if xv and yv are the same reference, and handles also $x == 0
370 if (($xv->[0] *= $yv->[0]) >= $MBASE)
372 $xv->[0] = $xv->[0] - ($xv->[1] = int($xv->[0] * $RBASE)) * $MBASE;
382 # multiply a large number a by a single element one, so speed up
383 my $y = $yv->[0]; my $car = 0;
386 $i = $i * $y + $car; $car = int($i * $RBASE); $i -= $car * $MBASE;
388 push @$xv, $car if $car != 0;
391 # shortcut for result $x == 0 => result = 0
392 return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
394 # since multiplying $x with $x fails, make copy in this case
395 $yv = [@$xv] if $xv == $yv; # same references?
397 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
406 # $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
408 # $prod - ($car = int($prod * RBASE)) * $MBASE; # see USE_MUL
410 # $prod[$cty] += $car if $car; # need really to check for 0?
414 # looping through this if $xi == 0 is silly - so optimize it away!
415 $xi = (shift @prod || 0), next if $xi == 0;
418 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
419 ## this is actually a tad slower
420 ## $prod = $prod[$cty]; $prod += ($car + $xi * $yi); # no ||0 here
422 $prod - ($car = int($prod * $RBASE)) * $MBASE; # see USE_MUL
424 $prod[$cty] += $car if $car; # need really to check for 0?
425 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
434 # (ref to int_num_array, ref to int_num_array)
435 # multiply two numbers in internal representation
436 # modifies first arg, second need not be different from first
437 my ($c,$xv,$yv) = @_;
441 # shortcut for two small numbers, also handles $x == 0
444 # shortcut for two very short numbers (improved by Nathan Zook)
445 # works also if xv and yv are the same reference, and handles also $x == 0
446 if (($xv->[0] *= $yv->[0]) >= $MBASE)
449 $xv->[0] - ($xv->[1] = int($xv->[0] / $MBASE)) * $MBASE;
459 # multiply a large number a by a single element one, so speed up
460 my $y = $yv->[0]; my $car = 0;
463 $i = $i * $y + $car; $car = int($i / $MBASE); $i -= $car * $MBASE;
465 push @$xv, $car if $car != 0;
468 # shortcut for result $x == 0 => result = 0
469 return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
471 # since multiplying $x with $x fails, make copy in this case
472 $yv = [@$xv] if $xv == $yv; # same references?
474 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
478 # looping through this if $xi == 0 is silly - so optimize it away!
479 $xi = (shift @prod || 0), next if $xi == 0;
482 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
484 $prod - ($car = int($prod / $MBASE)) * $MBASE;
486 $prod[$cty] += $car if $car; # need really to check for 0?
487 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
496 # ref to array, ref to array, modify first array and return remainder if
499 # see comments in _div_use_div() for more explanations
501 my ($c,$x,$yorg) = @_;
503 # the general div algorithmn here is about O(N*N) and thus quite slow, so
504 # we first check for some special cases and use shortcuts to handle them.
506 # This works, because we store the numbers in a chunked format where each
507 # element contains 5..7 digits (depending on system).
509 # if both numbers have only one element:
510 if (@$x == 1 && @$yorg == 1)
512 # shortcut, $yorg and $x are two small numbers
515 my $r = [ $x->[0] % $yorg->[0] ];
516 $x->[0] = int($x->[0] / $yorg->[0]);
521 $x->[0] = int($x->[0] / $yorg->[0]);
526 # if x has more than one, but y has only one element:
530 $rem = _mod($c,[ @$x ],$yorg) if wantarray;
532 # shortcut, $y is < $BASE
533 my $j = scalar @$x; my $r = 0;
534 my $y = $yorg->[0]; my $b;
537 $b = $r * $MBASE + $x->[$j];
538 $x->[$j] = int($b/$y);
541 pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
542 return ($x,$rem) if wantarray;
546 # now x and y have more than one element
548 # check whether y has more elements than x, if yet, the result will be 0
552 $rem = [@$x] if wantarray; # make copy
553 splice (@$x,1); # keep ref to original array
554 $x->[0] = 0; # set to 0
555 return ($x,$rem) if wantarray; # including remainder?
556 return $x; # only x, which is [0] now
558 # check whether the numbers have the same number of elements, in that case
559 # the result will fit into one element and can be computed efficiently
563 # if $yorg has more digits than $x (it's leading element is longer than
564 # the one from $x), the result will also be 0:
565 if (length(int($yorg->[-1])) > length(int($x->[-1])))
567 $rem = [@$x] if wantarray; # make copy
568 splice (@$x,1); # keep ref to org array
569 $x->[0] = 0; # set to 0
570 return ($x,$rem) if wantarray; # including remainder?
573 # now calculate $x / $yorg
574 if (length(int($yorg->[-1])) == length(int($x->[-1])))
576 # same length, so make full compare
578 my $a = 0; my $j = scalar @$x - 1;
579 # manual way (abort if unequal, good for early ne)
582 last if ($a = $x->[$j] - $yorg->[$j]); $j--;
584 # $a contains the result of the compare between X and Y
585 # a < 0: x < y, a == 0: x == y, a > 0: x > y
588 $rem = [ 0 ]; # a = 0 => x == y => rem 0
589 $rem = [@$x] if $a != 0; # a < 0 => x < y => rem = x
590 splice(@$x,1); # keep single element
591 $x->[0] = 0; # if $a < 0
592 $x->[0] = 1 if $a == 0; # $x == $y
593 return ($x,$rem) if wantarray;
596 # $x >= $y, so proceed normally
602 my $y = [ @$yorg ]; # always make copy to preserve
604 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
606 $car = $bar = $prd = 0;
607 if (($dd = int($MBASE/($y->[-1]+1))) != 1)
611 $xi = $xi * $dd + $car;
612 $xi -= ($car = int($xi * $RBASE)) * $MBASE; # see USE_MUL
614 push(@$x, $car); $car = 0;
617 $yi = $yi * $dd + $car;
618 $yi -= ($car = int($yi * $RBASE)) * $MBASE; # see USE_MUL
625 @q = (); ($v2,$v1) = @$y[-2,-1];
629 ($u2,$u1,$u0) = @$x[-3..-1];
631 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
633 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
634 --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
637 ($car, $bar) = (0,0);
638 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
640 $prd = $q * $y->[$yi] + $car;
641 $prd -= ($car = int($prd * $RBASE)) * $MBASE; # see USE_MUL
642 $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
644 if ($x->[-1] < $car + $bar)
647 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
650 if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $MBASE));
663 for $xi (reverse @$x)
665 $prd = $car * $MBASE + $xi;
666 $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
687 # ref to array, ref to array, modify first array and return remainder if
689 my ($c,$x,$yorg) = @_;
691 # the general div algorithmn here is about O(N*N) and thus quite slow, so
692 # we first check for some special cases and use shortcuts to handle them.
694 # This works, because we store the numbers in a chunked format where each
695 # element contains 5..7 digits (depending on system).
697 # if both numbers have only one element:
698 if (@$x == 1 && @$yorg == 1)
700 # shortcut, $yorg and $x are two small numbers
703 my $r = [ $x->[0] % $yorg->[0] ];
704 $x->[0] = int($x->[0] / $yorg->[0]);
709 $x->[0] = int($x->[0] / $yorg->[0]);
713 # if x has more than one, but y has only one element:
717 $rem = _mod($c,[ @$x ],$yorg) if wantarray;
719 # shortcut, $y is < $BASE
720 my $j = scalar @$x; my $r = 0;
721 my $y = $yorg->[0]; my $b;
724 $b = $r * $MBASE + $x->[$j];
725 $x->[$j] = int($b/$y);
728 pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
729 return ($x,$rem) if wantarray;
732 # now x and y have more than one element
734 # check whether y has more elements than x, if yet, the result will be 0
738 $rem = [@$x] if wantarray; # make copy
739 splice (@$x,1); # keep ref to original array
740 $x->[0] = 0; # set to 0
741 return ($x,$rem) if wantarray; # including remainder?
742 return $x; # only x, which is [0] now
744 # check whether the numbers have the same number of elements, in that case
745 # the result will fit into one element and can be computed efficiently
749 # if $yorg has more digits than $x (it's leading element is longer than
750 # the one from $x), the result will also be 0:
751 if (length(int($yorg->[-1])) > length(int($x->[-1])))
753 $rem = [@$x] if wantarray; # make copy
754 splice (@$x,1); # keep ref to org array
755 $x->[0] = 0; # set to 0
756 return ($x,$rem) if wantarray; # including remainder?
759 # now calculate $x / $yorg
761 if (length(int($yorg->[-1])) == length(int($x->[-1])))
763 # same length, so make full compare
765 my $a = 0; my $j = scalar @$x - 1;
766 # manual way (abort if unequal, good for early ne)
769 last if ($a = $x->[$j] - $yorg->[$j]); $j--;
771 # $a contains the result of the compare between X and Y
772 # a < 0: x < y, a == 0: x == y, a > 0: x > y
775 $rem = [ 0 ]; # a = 0 => x == y => rem 0
776 $rem = [@$x] if $a != 0; # a < 0 => x < y => rem = x
777 splice(@$x,1); # keep single element
778 $x->[0] = 0; # if $a < 0
779 $x->[0] = 1 if $a == 0; # $x == $y
780 return ($x,$rem) if wantarray; # including remainder?
783 # $x >= $y, so proceed normally
790 my $y = [ @$yorg ]; # always make copy to preserve
792 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
794 $car = $bar = $prd = 0;
795 if (($dd = int($MBASE/($y->[-1]+1))) != 1)
799 $xi = $xi * $dd + $car;
800 $xi -= ($car = int($xi / $MBASE)) * $MBASE;
802 push(@$x, $car); $car = 0;
805 $yi = $yi * $dd + $car;
806 $yi -= ($car = int($yi / $MBASE)) * $MBASE;
814 # @q will accumulate the final result, $q contains the current computed
815 # part of the final result
817 @q = (); ($v2,$v1) = @$y[-2,-1];
821 ($u2,$u1,$u0) = @$x[-3..-1];
823 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
825 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
826 --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
829 ($car, $bar) = (0,0);
830 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
832 $prd = $q * $y->[$yi] + $car;
833 $prd -= ($car = int($prd / $MBASE)) * $MBASE;
834 $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
836 if ($x->[-1] < $car + $bar)
839 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
842 if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $MBASE));
846 pop(@$x); unshift(@q, $q);
854 for $xi (reverse @$x)
856 $prd = $car * $MBASE + $xi;
857 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
876 ##############################################################################
881 # internal absolute post-normalized compare (ignore signs)
882 # ref to array, ref to array, return <0, 0, >0
883 # arrays must have at least one entry; this is not checked for
884 my ($c,$cx,$cy) = @_;
886 # shortcut for short numbers
887 return (($cx->[0] <=> $cy->[0]) <=> 0)
888 if scalar @$cx == scalar @$cy && scalar @$cx == 1;
890 # fast comp based on number of array elements (aka pseudo-length)
891 my $lxy = (scalar @$cx - scalar @$cy)
892 # or length of first element if same number of elements (aka difference 0)
894 # need int() here because sometimes the last element is '00018' vs '18'
895 (length(int($cx->[-1])) - length(int($cy->[-1])));
896 return -1 if $lxy < 0; # already differs, ret
897 return 1 if $lxy > 0; # ditto
899 # manual way (abort if unequal, good for early ne)
900 my $a; my $j = scalar @$cx;
903 last if ($a = $cx->[$j] - $cy->[$j]);
910 # compute number of digits
912 # int() because add/sub sometimes leaves strings (like '00005') instead of
913 # '5' in this place, thus causing length() to report wrong length
916 (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
921 # return the nth digit, negative values count backward
922 # zero is rightmost, so _digit(123,0) will give 3
925 my $len = _len('',$x);
927 $n = $len+$n if $n < 0; # -1 last, -2 second-to-last
928 $n = abs($n); # if negative was too big
929 $len--; $n = $len if $n > $len; # n to big?
931 my $elem = int($n / $BASE_LEN); # which array element
932 my $digit = $n % $BASE_LEN; # which digit in this element
933 $elem = '0' x $BASE_LEN . @$x[$elem]; # get element padded with 0's
934 substr($elem,-$digit-1,1);
939 # return amount of trailing zeros in decimal
940 # check each array elem in _m for having 0 at end as long as elem == 0
941 # Upon finding a elem != 0, stop
944 return 0 if scalar @$x == 1 && $x->[0] == 0;
946 my $zeros = 0; my $elem;
951 $elem = "$e"; # preserve x
952 $elem =~ s/.*?(0*$)/$1/; # strip anything not zero
953 $zeros *= $BASE_LEN; # elems * 5
954 $zeros += length($elem); # count trailing zeros
957 $zeros ++; # real else branch: 50% slower!
962 ##############################################################################
967 # return true if arg is zero
968 (((scalar @{$_[1]} == 1) && ($_[1]->[0] == 0))) <=> 0;
973 # return true if arg is even
974 (!($_[1]->[0] & 1)) <=> 0;
979 # return true if arg is even
980 (($_[1]->[0] & 1)) <=> 0;
985 # return true if arg is one
986 (scalar @{$_[1]} == 1) && ($_[1]->[0] == 1) <=> 0;
991 # return true if arg is two
992 (scalar @{$_[1]} == 1) && ($_[1]->[0] == 2) <=> 0;
997 # return true if arg is ten
998 (scalar @{$_[1]} == 1) && ($_[1]->[0] == 10) <=> 0;
1003 # internal normalization function that strips leading zeros from the array
1004 # args: ref to array
1007 my $cnt = scalar @$s; # get count of parts
1009 push @$s,0 if $i < 0; # div might return empty results, so fix it
1011 return $s if @$s == 1; # early out
1013 #print "strip: cnt $cnt i $i\n";
1014 # '0', '3', '4', '0', '0',
1019 # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
1020 # >= 1: skip first part (this can be zero)
1021 while ($i > 0) { last if $s->[$i] != 0; $i--; }
1022 $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
1026 ###############################################################################
1027 # check routine to test internal state for corruptions
1031 # used by the test suite
1034 return "$x is not a reference" if !ref($x);
1036 # are all parts are valid?
1037 my $i = 0; my $j = scalar @$x; my ($e,$try);
1040 $e = $x->[$i]; $e = 'undef' unless defined $e;
1041 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
1042 last if $e !~ /^[+]?[0-9]+$/;
1043 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
1044 last if "$e" !~ /^[+]?[0-9]+$/;
1045 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
1046 last if '' . "$e" !~ /^[+]?[0-9]+$/;
1047 $try = ' < 0 || >= $BASE; '."($x, $e)";
1048 last if $e <0 || $e >= $BASE;
1049 # this test is disabled, since new/bnorm and certain ops (like early out
1050 # in add/sub) are allowed/expected to leave '00000' in some elements
1051 #$try = '=~ /^00+/; '."($x, $e)";
1052 #last if $e =~ /^00+/;
1055 return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
1060 ###############################################################################
1064 # if possible, use mod shortcut
1065 my ($c,$x,$yo) = @_;
1067 # slow way since $y to big
1068 if (scalar @$yo > 1)
1070 my ($xo,$rem) = _div($c,$x,$yo);
1075 # both are single element arrays
1076 if (scalar @$x == 1)
1082 # @y is a single element, but @x has more than one element
1086 # when BASE % Y == 0 then (B * BASE) % Y == 0
1087 # (B * BASE) % $y + A % Y => A % Y
1088 # so need to consider only last element: O(1)
1093 # else need to go through all elements: O(N), but loop is a bit simplified
1097 $r = ($r + $_) % $y; # not much faster, but heh...
1098 #$r += $_ % $y; $r %= $y;
1105 # else need to go through all elements: O(N)
1106 my $r = 0; my $bm = 1;
1109 $r = ($_ * $bm + $r) % $y;
1110 $bm = ($bm * $b) % $y;
1112 #$r += ($_ % $y) * $bm;
1120 splice (@$x,1); # keep one element of $x
1124 ##############################################################################
1129 my ($c,$x,$y,$n) = @_;
1133 $n = _new($c,$n); return _div($c,$x, _pow($c,$n,$y));
1136 # shortcut (faster) for shifting by 10)
1137 # multiples of $BASE_LEN
1138 my $dst = 0; # destination
1139 my $src = _num($c,$y); # as normal int
1140 my $xlen = (@$x-1)*$BASE_LEN+length(int($x->[-1])); # len of x in digits
1141 if ($src >= $xlen or ($src == $xlen and ! defined $x->[1]))
1143 # 12345 67890 shifted right by more than 10 digits => 0
1144 splice (@$x,1); # leave only one element
1145 $x->[0] = 0; # set to zero
1148 my $rem = $src % $BASE_LEN; # remainder to shift
1149 $src = int($src / $BASE_LEN); # source
1152 splice (@$x,0,$src); # even faster, 38.4 => 39.3
1156 my $len = scalar @$x - $src; # elems to go
1157 my $vd; my $z = '0'x $BASE_LEN;
1158 $x->[scalar @$x] = 0; # avoid || 0 test inside loop
1161 $vd = $z.$x->[$src];
1162 $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
1164 $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
1165 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1166 $x->[$dst] = int($vd);
1169 splice (@$x,$dst) if $dst > 0; # kill left-over array elems
1170 pop @$x if $x->[-1] == 0 && @$x > 1; # kill last element if 0
1177 my ($c,$x,$y,$n) = @_;
1181 $n = _new($c,$n); return _mul($c,$x, _pow($c,$n,$y));
1184 # shortcut (faster) for shifting by 10) since we are in base 10eX
1185 # multiples of $BASE_LEN:
1186 my $src = scalar @$x; # source
1187 my $len = _num($c,$y); # shift-len as normal int
1188 my $rem = $len % $BASE_LEN; # remainder to shift
1189 my $dst = $src + int($len/$BASE_LEN); # destination
1190 my $vd; # further speedup
1191 $x->[$src] = 0; # avoid first ||0 for speed
1192 my $z = '0' x $BASE_LEN;
1195 $vd = $x->[$src]; $vd = $z.$vd;
1196 $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
1197 $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
1198 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1199 $x->[$dst] = int($vd);
1202 # set lowest parts to 0
1203 while ($dst >= 0) { $x->[$dst--] = 0; }
1204 # fix spurios last zero element
1205 splice @$x,-1 if $x->[-1] == 0;
1212 # ref to array, ref to array, return ref to array
1213 my ($c,$cx,$cy) = @_;
1215 if (scalar @$cy == 1 && $cy->[0] == 0)
1217 splice (@$cx,1); $cx->[0] = 1; # y == 0 => x => 1
1220 if ((scalar @$cx == 1 && $cx->[0] == 1) || # x == 1
1221 (scalar @$cy == 1 && $cy->[0] == 1)) # or y == 1
1225 if (scalar @$cx == 1 && $cx->[0] == 0)
1227 splice (@$cx,1); $cx->[0] = 0; # 0 ** y => 0 (if not y <= 0)
1233 my $y_bin = _as_bin($c,$cy); $y_bin =~ s/^0b//;
1234 my $len = length($y_bin);
1237 _mul($c,$pow2,$cx) if substr($y_bin,$len,1) eq '1'; # is odd?
1248 # ref to array, return ref to array
1251 if ((@$cx == 1) && ($cx->[0] <= 2))
1253 $cx->[0] ||= 1; # 0 => 1, 1 => 1, 2 => 2
1257 # go forward until $base is exceeded
1258 # limit is either $x steps (steps == 100 means a result always too high) or
1260 my $steps = 100; $steps = $cx->[0] if @$cx == 1;
1261 my $r = 2; my $cf = 3; my $step = 2; my $last = $r;
1262 while ($r*$cf < $BASE && $step < $steps)
1264 $last = $r; $r *= $cf++; $step++;
1266 if ((@$cx == 1) && $step == $cx->[0])
1268 # completely done, so keep reference to $x and return
1273 # now we must do the left over steps
1274 my $n; # steps still to do
1275 if (scalar @$cx == 1)
1284 $cx->[0] = $last; splice (@$cx,1); # keep ref to $x
1285 my $zero_elements = 0;
1287 # do left-over steps fit into a scalar?
1288 if (ref $n eq 'ARRAY')
1290 # No, so use slower inc() & cmp()
1292 while (_acmp($step,$n) <= 0)
1294 # as soon as the last element of $cx is 0, we split it up and remember
1295 # how many zeors we got so far. The reason is that n! will accumulate
1296 # zeros at the end rather fast.
1299 $zero_elements ++; shift @$cx;
1301 _mul($c,$cx,$step); _inc($c,$step);
1306 # Yes, so we can speed it up slightly
1309 # When the last element of $cx is 0, we split it up and remember
1310 # how many we got so far. The reason is that n! will accumulate
1311 # zeros at the end rather fast.
1314 $zero_elements ++; shift @$cx;
1316 _mul($c,$cx,[$step]); $step++;
1319 # multiply in the zeros again
1320 while ($zero_elements-- > 0)
1324 $cx; # return result
1327 #############################################################################
1331 # calculate integer log of $x to base $base
1332 # ref to array, ref to array - return ref to array
1333 my ($c,$x,$base) = @_;
1336 return if (scalar @$x == 1 && $x->[0] == 0);
1337 # BASE 0 or 1 => NaN
1338 return if (scalar @$base == 1 && $base->[0] < 2);
1339 my $cmp = _acmp($c,$x,$base); # X == BASE => 1
1342 splice (@$x,1); $x->[0] = 1;
1348 splice (@$x,1); $x->[0] = 0;
1352 # this trial multiplication is very fast, even for large counts (like for
1353 # 2 ** 1024, since this still requires only 1024 very fast steps
1354 # (multiplication of a large number by a very small number is very fast))
1355 my $x_org = _copy($c,$x); # preserve x
1356 splice(@$x,1); $x->[0] = 1; # keep ref to $x
1358 my $trial = _copy($c,$base);
1360 # XXX TODO this only works if $base has only one element
1361 if (scalar @$base == 1)
1363 # compute int ( length_in_base_10(X) / ( log(base) / log(10) ) )
1364 my $len = _len($c,$x_org);
1365 my $res = int($len / (log($base->[0]) / log(10))) || 1; # avoid $res == 0
1368 $trial = _pow ($c, _copy($c, $base), $x);
1369 my $a = _acmp($x,$trial,$x_org);
1370 return ($x,1) if $a == 0;
1371 # we now know that $res is too small
1374 _mul($c,$trial,$base); _add($c, $x, [1]);
1379 _div($c,$trial,$base); _sub($c, $x, [1]);
1381 # did we now get the right result?
1382 $a = _acmp($x,$trial,$x_org);
1383 return ($x,1) if $a == 0; # yes, exactly
1387 _div($c,$trial,$base); _sub($c, $x, [1]);
1391 # simple loop that increments $x by two in each step, possible overstepping
1392 # the real result by one
1395 my $base_mul = _mul($c, _copy($c,$base), $base);
1397 while (($a = _acmp($c,$trial,$x_org)) < 0)
1399 _mul($c,$trial,$base_mul); _add($c, $x, [2]);
1405 # overstepped the result
1407 _div($c,$trial,$base);
1408 $a = _acmp($c,$trial,$x_org);
1413 $exact = 0 if $a != 0;
1416 ($x,$exact); # return result
1420 use constant DEBUG => 0;
1422 sub steps { $steps };
1426 # square-root of $x in place
1427 # Compute a guess of the result (by rule of thumb), then improve it via
1431 if (scalar @$x == 1)
1433 # fit's into one Perl scalar, so result can be computed directly
1434 $x->[0] = int(sqrt($x->[0]));
1437 my $y = _copy($c,$x);
1438 # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
1439 # since our guess will "grow"
1440 my $l = int((_len($c,$x)-1) / 2);
1442 my $lastelem = $x->[-1]; # for guess
1443 my $elems = scalar @$x - 1;
1444 # not enough digits, but could have more?
1445 if ((length($lastelem) <= 3) && ($elems > 1))
1447 # right-align with zero pad
1448 my $len = length($lastelem) & 1;
1449 print "$lastelem => " if DEBUG;
1450 $lastelem .= substr($x->[-2] . '0' x $BASE_LEN,0,$BASE_LEN);
1451 # former odd => make odd again, or former even to even again
1452 $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
1453 print "$lastelem\n" if DEBUG;
1456 # construct $x (instead of _lsft($c,$x,$l,10)
1457 my $r = $l % $BASE_LEN; # 10000 00000 00000 00000 ($BASE_LEN=5)
1458 $l = int($l / $BASE_LEN);
1459 print "l = $l " if DEBUG;
1461 splice @$x,$l; # keep ref($x), but modify it
1463 # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
1465 # 14400 00000 => sqrt(14400) => guess first digits to be 120
1466 # 144000 000000 => sqrt(144000) => guess 379
1468 print "$lastelem (elems $elems) => " if DEBUG;
1469 $lastelem = $lastelem / 10 if ($elems & 1 == 1); # odd or even?
1470 my $g = sqrt($lastelem); $g =~ s/\.//; # 2.345 => 2345
1471 $r -= 1 if $elems & 1 == 0; # 70 => 7
1473 # padd with zeros if result is too short
1474 $x->[$l--] = int(substr($g . '0' x $r,0,$r+1));
1475 print "now ",$x->[-1] if DEBUG;
1476 print " would have been ", int('1' . '0' x $r),"\n" if DEBUG;
1478 # If @$x > 1, we could compute the second elem of the guess, too, to create
1479 # an even better guess. Not implemented yet. Does it improve performance?
1480 $x->[$l--] = 0 while ($l >= 0); # all other digits of guess are zero
1482 print "start x= ",_str($c,$x),"\n" if DEBUG;
1485 my $lastlast = _zero();
1486 $steps = 0 if DEBUG;
1487 while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
1490 $lastlast = _copy($c,$last);
1491 $last = _copy($c,$x);
1492 _add($c,$x, _div($c,_copy($c,$y),$x));
1494 print " x= ",_str($c,$x),"\n" if DEBUG;
1496 print "\nsteps in sqrt: $steps, " if DEBUG;
1497 _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0; # overshot?
1498 print " final ",$x->[-1],"\n" if DEBUG;
1504 # take n'th root of $x in place (n >= 3)
1507 if (scalar @$x == 1)
1511 # result will always be smaller than 2 so trunc to 1 at once
1516 # fit's into one Perl scalar, so result can be computed directly
1517 # cannot use int() here, because it rounds wrongly (try
1518 # (81 ** 3) ** (1/3) to see what I mean)
1519 #$x->[0] = int( $x->[0] ** (1 / $n->[0]) );
1520 # round to 8 digits, then truncate result to integer
1521 $x->[0] = int ( sprintf ("%.8f", $x->[0] ** (1 / $n->[0]) ) );
1526 # we know now that X is more than one element long
1528 # if $n is a power of two, we can repeatedly take sqrt($X) and find the
1529 # proper result, because sqrt(sqrt($x)) == root($x,4)
1530 my $b = _as_bin($c,$n);
1531 if ($b =~ /0b1(0+)$/)
1533 my $count = CORE::length($1); # 0b100 => len('00') => 2
1534 my $cnt = $count; # counter for loop
1535 unshift (@$x, 0); # add one element, together with one
1536 # more below in the loop this makes 2
1539 # 'inflate' $X by adding one element, basically computing
1540 # $x * $BASE * $BASE. This gives us more $BASE_LEN digits for result
1541 # since len(sqrt($X)) approx == len($x) / 2.
1543 # calculate sqrt($x), $x is now one element to big, again. In the next
1544 # round we make that two, again.
1547 # $x is now one element to big, so truncate result by removing it
1552 # trial computation by starting with 2,4,8,16 etc until we overstep
1556 # while still to do more than X steps
1560 while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
1562 _mul ($c, $step, [2]);
1563 _add ($c, $trial, $step);
1567 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) == 0)
1569 @$x = @$trial; # make copy while preserving ref to $x
1572 # overstepped, so go back on step
1573 _sub($c, $trial, $step);
1574 } while (scalar @$step > 1 || $step->[0] > 128);
1578 # add two, because $trial cannot be exactly the result (otherwise we would
1579 # alrady have found it)
1580 _add($c, $trial, $step);
1582 # and now add more and more (2,4,6,8,10 etc)
1583 while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
1585 _add ($c, $trial, $step);
1588 # hit not exactly? (overstepped)
1589 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
1594 # hit not exactly? (overstepped)
1595 # 80 too small, 81 slightly too big, 82 too big
1596 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
1601 @$x = @$trial; # make copy while preserving ref to $x
1607 ##############################################################################
1614 # the shortcut makes equal, large numbers _really_ fast, and makes only a
1615 # very small performance drop for small numbers (e.g. something with less
1616 # than 32 bit) Since we optimize for large numbers, this is enabled.
1617 return $x if _acmp($c,$x,$y) == 0; # shortcut
1619 my $m = _one(); my ($xr,$yr);
1620 my $mask = $AND_MASK;
1623 my $y1 = _copy($c,$y); # make copy
1627 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1629 ($x1, $xr) = _div($c,$x1,$mask);
1630 ($y1, $yr) = _div($c,$y1,$mask);
1632 # make ints() from $xr, $yr
1633 # this is when the AND_BITS are greater than $BASE and is slower for
1634 # small (<256 bits) numbers, but faster for large numbers. Disabled
1635 # due to KISS principle
1637 # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1638 # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1639 # _add($c,$x, _mul($c, _new( $c, ($xrr & $yrr) ), $m) );
1641 # 0+ due to '&' doesn't work in strings
1642 _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
1652 return _zero() if _acmp($c,$x,$y) == 0; # shortcut (see -and)
1654 my $m = _one(); my ($xr,$yr);
1655 my $mask = $XOR_MASK;
1658 my $y1 = _copy($c,$y); # make copy
1662 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1664 ($x1, $xr) = _div($c,$x1,$mask);
1665 ($y1, $yr) = _div($c,$y1,$mask);
1666 # make ints() from $xr, $yr (see _and())
1667 #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1668 #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1669 #_add($c,$x, _mul($c, _new( $c, ($xrr ^ $yrr) ), $m) );
1671 # 0+ due to '^' doesn't work in strings
1672 _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
1675 # the loop stops when the shorter of the two numbers is exhausted
1676 # the remainder of the longer one will survive bit-by-bit, so we simple
1677 # multiply-add it in
1678 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1679 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1688 return $x if _acmp($c,$x,$y) == 0; # shortcut (see _and)
1690 my $m = _one(); my ($xr,$yr);
1691 my $mask = $OR_MASK;
1694 my $y1 = _copy($c,$y); # make copy
1698 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1700 ($x1, $xr) = _div($c,$x1,$mask);
1701 ($y1, $yr) = _div($c,$y1,$mask);
1702 # make ints() from $xr, $yr (see _and())
1703 # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1704 # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1705 # _add($c,$x, _mul($c, _new( $c, ($xrr | $yrr) ), $m) );
1707 # 0+ due to '|' doesn't work in strings
1708 _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
1711 # the loop stops when the shorter of the two numbers is exhausted
1712 # the remainder of the longer one will survive bit-by-bit, so we simple
1713 # multiply-add it in
1714 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1715 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1722 # convert a decimal number to hex (ref to array, return ref to string)
1725 # fit's into one element (handle also 0x0 case)
1726 return sprintf("0x%x",$x->[0]) if @$x == 1;
1728 my $x1 = _copy($c,$x);
1731 my ($xr, $h, $x10000);
1734 $x10000 = [ 0x10000 ]; $h = 'h4';
1738 $x10000 = [ 0x1000 ]; $h = 'h3';
1740 while (@$x1 != 1 || $x1->[0] != 0) # _is_zero()
1742 ($x1, $xr) = _div($c,$x1,$x10000);
1743 $es .= unpack($h,pack('v',$xr->[0])); # XXX TODO: why pack('v',...)?
1746 $es =~ s/^[0]+//; # strip leading zeros
1747 '0x' . $es; # return result prepended with 0x
1752 # convert a decimal number to bin (ref to array, return ref to string)
1755 # fit's into one element (and Perl recent enough), handle also 0b0 case
1756 # handle zero case for older Perls
1757 if ($] <= 5.005 && @$x == 1 && $x->[0] == 0)
1759 my $t = '0b0'; return $t;
1761 if (@$x == 1 && $] >= 5.006)
1763 my $t = sprintf("0b%b",$x->[0]);
1766 my $x1 = _copy($c,$x);
1769 my ($xr, $b, $x10000);
1772 $x10000 = [ 0x10000 ]; $b = 'b16';
1776 $x10000 = [ 0x1000 ]; $b = 'b12';
1778 while (!(@$x1 == 1 && $x1->[0] == 0)) # _is_zero()
1780 ($x1, $xr) = _div($c,$x1,$x10000);
1781 $es .= unpack($b,pack('v',$xr->[0])); # XXX TODO: why pack('v',...)?
1782 # $es .= unpack($b,$xr->[0]);
1785 $es =~ s/^[0]+//; # strip leading zeros
1786 '0b' . $es; # return result prepended with 0b
1791 # convert a hex number to decimal (ref to string, return ref to array)
1794 my $m = _new($c, 0x10000000); # 28 bit at a time (<32 bit!)
1795 my $d = 7; # 7 digits at a time
1798 # for older Perls, play safe
1799 $m = [ 0x10000 ]; # 16 bit at a time (<32 bit!)
1800 $d = 4; # 4 digits at a time
1806 my $len = int( (length($hs)-2)/$d ); # $d digit parts, w/o the '0x'
1807 my $val; my $i = -$d;
1810 $val = substr($hs,$i,$d); # get hex digits
1811 $val =~ s/^[+-]?0x// if $len == 0; # for last part only because
1812 $val = hex($val); # hex does not like wrong chars
1814 my $adder = [ $val ];
1815 # if the resulting number was to big to fit into one element, create a
1816 # two-element version (bug found by Mark Lakata - Thanx!)
1817 if (CORE::length($val) > $BASE_LEN)
1819 $adder = _new($c,$val);
1821 _add ($c, $x, _mul ($c, $adder, $mul ) ) if $val != 0;
1822 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
1829 # convert a hex number to decimal (ref to string, return ref to array)
1832 # instead of converting X (8) bit at a time, it is faster to "convert" the
1833 # number to hex, and then call _from_hex.
1836 $hs =~ s/^[+-]?0b//; # remove sign and 0b
1837 my $l = length($hs); # bits
1838 $hs = '0' x (8-($l % 8)) . $hs if ($l % 8) != 0; # padd left side w/ 0
1839 my $h = '0x' . unpack('H*', pack ('B*', $hs)); # repack as hex
1844 ##############################################################################
1845 # special modulus functions
1852 my $u = _zero($c); my $u1 = _one($c);
1853 my $a = _copy($c,$y); my $b = _copy($c,$x);
1855 # Euclid's Algorithm for bgcd(), only that we calc bgcd() ($a) and the
1856 # result ($u) at the same time. See comments in BigInt for why this works.
1858 ($a, $q, $b) = ($b, _div($c,$a,$b)); # step 1
1860 while (!_is_zero($c,$b))
1862 my $t = _add($c, # step 2:
1863 _mul($c,_copy($c,$u1), $q) , # t = u1 * q
1865 $u = $u1; # u = u1, u1 = t
1868 ($a, $q, $b) = ($b, _div($c,$a,$b)); # step 1
1871 # if the gcd is not 1, then return NaN
1872 return (undef,undef) unless _is_one($c,$a);
1874 ($u1, $sign == 1 ? '+' : '-');
1879 # modulus of power ($x ** $y) % $z
1880 my ($c,$num,$exp,$mod) = @_;
1882 # in the trivial case,
1883 if (_is_one($c,$mod))
1885 splice @$num,0,1; $num->[0] = 0;
1888 if ((scalar @$num == 1) && (($num->[0] == 0) || ($num->[0] == 1)))
1894 # $num = _mod($c,$num,$mod); # this does not make it faster
1896 my $acc = _copy($c,$num); my $t = _one();
1898 my $expbin = _as_bin($c,$exp); $expbin =~ s/^0b//;
1899 my $len = length($expbin);
1902 if ( substr($expbin,$len,1) eq '1') # is_odd
1905 $t = _mod($c,$t,$mod);
1908 $acc = _mod($c,$acc,$mod);
1916 # greatest common divisor
1919 while ( (scalar @$y != 1) || ($y->[0] != 0) ) # while ($y != 0)
1921 my $t = _copy($c,$y);
1922 $y = _mod($c, $x, $y);
1928 ##############################################################################
1929 ##############################################################################
1936 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
1940 Provides support for big integer calculations. Not intended to be used by other
1941 modules. Other modules which sport the same functions can also be used to support
1942 Math::BigInt, like Math::BigInt::GMP or Math::BigInt::Pari.
1946 In order to allow for multiple big integer libraries, Math::BigInt was
1947 rewritten to use library modules for core math routines. Any module which
1948 follows the same API as this can be used instead by using the following:
1950 use Math::BigInt lib => 'libname';
1952 'libname' is either the long name ('Math::BigInt::Pari'), or only the short
1953 version like 'Pari'.
1959 The following functions MUST be defined in order to support the use by
1960 Math::BigInt v1.70 or later:
1962 api_version() return API version, minimum 1 for v1.70
1963 _new(string) return ref to new object from ref to decimal string
1964 _zero() return a new object with value 0
1965 _one() return a new object with value 1
1966 _two() return a new object with value 2
1967 _ten() return a new object with value 10
1969 _str(obj) return ref to a string representing the object
1970 _num(obj) returns a Perl integer/floating point number
1971 NOTE: because of Perl numeric notation defaults,
1972 the _num'ified obj may lose accuracy due to
1973 machine-dependent floating point size limitations
1975 _add(obj,obj) Simple addition of two objects
1976 _mul(obj,obj) Multiplication of two objects
1977 _div(obj,obj) Division of the 1st object by the 2nd
1978 In list context, returns (result,remainder).
1979 NOTE: this is integer math, so no
1980 fractional part will be returned.
1981 The second operand will be not be 0, so no need to
1983 _sub(obj,obj) Simple subtraction of 1 object from another
1984 a third, optional parameter indicates that the params
1985 are swapped. In this case, the first param needs to
1986 be preserved, while you can destroy the second.
1987 sub (x,y,1) => return x - y and keep x intact!
1988 _dec(obj) decrement object by one (input is guaranteed to be > 0)
1989 _inc(obj) increment object by one
1992 _acmp(obj,obj) <=> operator for objects (return -1, 0 or 1)
1994 _len(obj) returns count of the decimal digits of the object
1995 _digit(obj,n) returns the n'th decimal digit of object
1997 _is_one(obj) return true if argument is 1
1998 _is_two(obj) return true if argument is 2
1999 _is_ten(obj) return true if argument is 10
2000 _is_zero(obj) return true if argument is 0
2001 _is_even(obj) return true if argument is even (0,2,4,6..)
2002 _is_odd(obj) return true if argument is odd (1,3,5,7..)
2004 _copy return a ref to a true copy of the object
2006 _check(obj) check whether internal representation is still intact
2007 return 0 for ok, otherwise error message as string
2009 _from_hex(str) return ref to new object from ref to hexadecimal string
2010 _from_bin(str) return ref to new object from ref to binary string
2012 _as_hex(str) return string containing the value as
2013 unsigned hex string, with the '0x' prepended.
2014 Leading zeros must be stripped.
2015 _as_bin(str) Like as_hex, only as binary string containing only
2016 zeros and ones. Leading zeros must be stripped and a
2017 '0b' must be prepended.
2019 _rsft(obj,N,B) shift object in base B by N 'digits' right
2020 _lsft(obj,N,B) shift object in base B by N 'digits' left
2022 _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
2023 Note: XOR, AND and OR pad with zeros if size mismatches
2024 _and(obj1,obj2) AND (bit-wise) object 1 with object 2
2025 _or(obj1,obj2) OR (bit-wise) object 1 with object 2
2027 _mod(obj,obj) Return remainder of div of the 1st by the 2nd object
2028 _sqrt(obj) return the square root of object (truncated to int)
2029 _root(obj) return the n'th (n >= 3) root of obj (truncated to int)
2030 _fac(obj) return factorial of object 1 (1*2*3*4..)
2031 _pow(obj,obj) return object 1 to the power of object 2
2032 return undef for NaN
2033 _zeros(obj) return number of trailing decimal zeros
2034 _modinv return inverse modulus
2035 _modpow return modulus of power ($x ** $y) % $z
2036 _log_int(X,N) calculate integer log() of X in base N
2037 X >= 0, N >= 0 (return undef for NaN)
2038 returns (RESULT, EXACT) where EXACT is:
2039 1 : result is exactly RESULT
2040 0 : result was truncated to RESULT
2041 undef : unknown whether result is exactly RESULT
2042 _gcd(obj,obj) return Greatest Common Divisor of two objects
2044 The following functions are optional, and can be defined if the underlying lib
2045 has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
2046 slow) fallback routines to emulate these:
2053 Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
2056 So the library needs only to deal with unsigned big integers. Testing of input
2057 parameter validity is done by the caller, so you need not worry about
2058 underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by zero or similar
2061 The first parameter can be modified, that includes the possibility that you
2062 return a reference to a completely different object instead. Although keeping
2063 the reference and just changing it's contents is preferred over creating and
2064 returning a different reference.
2066 Return values are always references to objects, strings, or true/false for
2067 comparison routines.
2069 =head1 WRAP YOUR OWN
2071 If you want to port your own favourite c-lib for big numbers to the
2072 Math::BigInt interface, you can take any of the already existing modules as
2073 a rough guideline. You should really wrap up the latest BigInt and BigFloat
2074 testsuites with your module, and replace in them any of the following:
2080 use Math::BigInt lib => 'yourlib';
2082 This way you ensure that your library really works 100% within Math::BigInt.
2086 This program is free software; you may redistribute it and/or modify it under
2087 the same terms as Perl itself.
2091 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
2093 Seperated from BigInt and shaped API with the help of John Peacock.
2095 Fixed, speed-up, streamlined and enhanced by Tels 2001 - 2005.
2099 L<Math::BigInt>, L<Math::BigFloat>, L<Math::BigInt::BitVect>,
2100 L<Math::BigInt::GMP>, L<Math::BigInt::FastCalc> and L<Math::BigInt::Pari>.