1 package Math::BigInt::Calc;
5 # use warnings; # dont use warnings for older Perls
11 # Package to store unsigned big integers in decimal and do math with them
13 # Internally the numbers are stored in an array with at least 1 element, no
14 # leading zero parts (except the first) and in base 1eX where X is determined
15 # automatically at loading time to be the maximum possible value
18 # - fully remove funky $# stuff (maybe)
20 # USE_MUL: due to problems on certain os (os390, posix-bc) "* 1e-5" is used
21 # instead of "/ 1e5" at some places, (marked with USE_MUL). Other platforms
22 # BS2000, some Crays need USE_DIV instead.
23 # The BEGIN block is used to determine which of the two variants gives the
26 # Beware of things like:
27 # $i = $i * $y + $car; $car = int($i / $MBASE); $i = $i % $MBASE;
28 # This works on x86, but fails on ARM (SA1100, iPAQ) due to whoknows what
29 # reasons. So, use this instead (slower, but correct):
30 # $i = $i * $y + $car; $car = int($i / $MBASE); $i -= $MBASE * $car;
32 ##############################################################################
33 # global constants, flags and accessory
35 # announce that we are compatible with MBI v1.70 and up
36 sub api_version () { 1; }
38 # constants for easier life
40 my ($MBASE,$BASE,$RBASE,$BASE_LEN,$MAX_VAL,$BASE_LEN2,$BASE_LEN_SMALL);
41 my ($AND_BITS,$XOR_BITS,$OR_BITS);
42 my ($AND_MASK,$XOR_MASK,$OR_MASK);
46 # set/get the BASE_LEN and assorted other, connected values
47 # used only be the testsuite, set is used only by the BEGIN block below
53 # find whether we can use mul or div or none in mul()/div()
54 # (in last case reduce BASE_LEN_SMALL)
55 $BASE_LEN_SMALL = $b+1;
57 while (--$BASE_LEN_SMALL > 5)
59 $MBASE = int("1e".$BASE_LEN_SMALL);
60 $RBASE = abs('1e-'.$BASE_LEN_SMALL); # see USE_MUL
62 $caught += 1 if (int($MBASE * $RBASE) != 1); # should be 1
63 $caught += 2 if (int($MBASE / $MBASE) != 1); # should be 1
66 # BASE_LEN is used for anything else than mul()/div()
67 $BASE_LEN = $BASE_LEN_SMALL;
68 $BASE_LEN = shift if (defined $_[0]); # one more arg?
69 $BASE = int("1e".$BASE_LEN);
71 $BASE_LEN2 = int($BASE_LEN_SMALL / 2); # for mul shortcut
72 $MBASE = int("1e".$BASE_LEN_SMALL);
73 $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);
102 # from Daniel Pfeiffer: determine largest group of digits that is precisely
103 # multipliable with itself plus carry
104 # Test now changed to expect the proper pattern, not a result off by 1 or 2
105 my ($e, $num) = 3; # lowest value we will use is 3+1-1 = 3
108 $num = ('9' x ++$e) + 0;
110 } while ("$num" =~ /9{$e}0{$e}/); # must be a certain pattern
111 $e--; # last test failed, so retract one step
112 # the limits below brush the problems with the test above under the rug:
113 # the test should be able to find the proper $e automatically
114 $e = 5 if $^O =~ /^uts/; # UTS get's some special treatment
115 $e = 5 if $^O =~ /^unicos/; # unicos is also problematic (6 seems to work
116 # there, but we play safe)
117 $e = 5 if $] < 5.006; # cap, for older Perls
118 $e = 7 if $e > 7; # cap, for VMS, OS/390 and other 64 bit systems
119 # 8 fails inside random testsuite, so take 7
121 # determine how many digits fit into an integer and can be safely added
122 # together plus carry w/o causing an overflow
126 ############################################################################
127 # the next block is no longer important
129 ## this below detects 15 on a 64 bit system, because after that it becomes
130 ## 1e16 and not 1000000 :/ I can make it detect 18, but then I get a lot of
131 ## test failures. Ugh! (Tomake detect 18: uncomment lines marked with *)
133 #my $bi = 5; # approx. 16 bit
134 #$num = int('9' x $bi);
136 ## while ( ($num+$num+1) eq '1' . '9' x $bi) # *
137 #while ( int($num+$num+1) eq '1' . '9' x $bi)
139 # $bi++; $num = int('9' x $bi);
140 # # $bi++; $num *= 10; $num += 9; # *
142 #$bi--; # back off one step
143 # by setting them equal, we ignore the findings and use the default
144 # one-size-fits-all approach from former versions
145 my $bi = $e; # XXX, this should work always
147 __PACKAGE__->_base_len($e,$bi); # set and store
149 # find out how many bits _and, _or and _xor can take (old default = 16)
150 # I don't think anybody has yet 128 bit scalars, so let's play safe.
151 local $^W = 0; # don't warn about 'nonportable number'
152 $AND_BITS = 15; $XOR_BITS = 15; $OR_BITS = 15;
154 # find max bits, we will not go higher than numberofbits that fit into $BASE
155 # to make _and etc simpler (and faster for smaller, slower for large numbers)
157 while (2 ** $max < $BASE) { $max++; }
160 $max = 16 if $] < 5.006; # older Perls might not take >16 too well
165 $x = oct('0b' . '1' x $AND_BITS); $y = $x & $x;
166 $z = (2 ** $AND_BITS) - 1;
167 } while ($AND_BITS < $max && $x == $z && $y == $x);
168 $AND_BITS --; # retreat one step
171 $x = oct('0b' . '1' x $XOR_BITS); $y = $x ^ 0;
172 $z = (2 ** $XOR_BITS) - 1;
173 } while ($XOR_BITS < $max && $x == $z && $y == $x);
174 $XOR_BITS --; # retreat one step
177 $x = oct('0b' . '1' x $OR_BITS); $y = $x | $x;
178 $z = (2 ** $OR_BITS) - 1;
179 } while ($OR_BITS < $max && $x == $z && $y == $x);
180 $OR_BITS --; # retreat one step
184 ###############################################################################
188 # (ref to string) return ref to num_array
189 # Convert a number from string format (without sign) to internal base
190 # 1ex format. Assumes normalized value as input.
191 my $il = length($_[1])-1;
193 # < BASE_LEN due len-1 above
194 return [ int($_[1]) ] if $il < $BASE_LEN; # shortcut for short numbers
196 # this leaves '00000' instead of int 0 and will be corrected after any op
197 [ reverse(unpack("a" . ($il % $BASE_LEN+1)
198 . ("a$BASE_LEN" x ($il / $BASE_LEN)), $_[1])) ];
203 $AND_MASK = __PACKAGE__->_new( ( 2 ** $AND_BITS ));
204 $XOR_MASK = __PACKAGE__->_new( ( 2 ** $XOR_BITS ));
205 $OR_MASK = __PACKAGE__->_new( ( 2 ** $OR_BITS ));
222 # create a two (used internally for shifting)
228 # create a 10 (used internally for shifting)
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")
252 my $l = scalar @$ar; # number of parts
253 return $nan if $l < 1; # should not happen
255 # handle first one different to strip leading zeros from it (there are no
256 # leading zero parts in internal representation)
257 $l --; $ret .= int($ar->[$l]); $l--;
258 # Interestingly, the pre-padd method uses more time
259 # the old grep variant takes longer (14 vs. 10 sec)
260 my $z = '0' x ($BASE_LEN-1);
263 $ret .= substr($z.$ar->[$l],-$BASE_LEN); # fastest way I could think of
271 # Make a number (scalar int/float) from a BigInt object
274 return 0+$x->[0] if scalar @$x == 1; # below $BASE
279 $num += $fac*$_; $fac *= $BASE;
284 ##############################################################################
289 # (ref to int_num_array, ref to int_num_array)
290 # routine to add two base 1eX numbers
291 # stolen from Knuth Vol 2 Algorithm A pg 231
292 # there are separate routines to add and sub as per Knuth pg 233
293 # This routine clobbers up array x, but not y.
297 return $x if (@$y == 1) && $y->[0] == 0; # $x + 0 => $x
298 if ((@$x == 1) && $x->[0] == 0) # 0 + $y => $y->copy
300 # twice as slow as $x = [ @$y ], but necc. to retain $x as ref :(
301 @$x = @$y; return $x;
304 # for each in Y, add Y to X and carry. If after that, something is left in
305 # X, foreach in X add carry to X and then return X, carry
306 # Trades one "$j++" for having to shift arrays
307 my $i; my $car = 0; my $j = 0;
310 $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
315 $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
322 # (ref to int_num_array, ref to int_num_array)
323 # Add 1 to $x, modify $x in place
328 return $x if (($i += 1) < $BASE); # early out
329 $i = 0; # overflow, next
331 push @$x,1 if ($x->[-1] == 0); # last overflowed, so extend
337 # (ref to int_num_array, ref to int_num_array)
338 # Sub 1 from $x, modify $x in place
341 my $MAX = $BASE-1; # since MAX_VAL based on MBASE
344 last if (($i -= 1) >= 0); # early out
345 $i = $MAX; # underflow, next
347 pop @$x if $x->[-1] == 0 && @$x > 1; # last underflowed (but leave 0)
353 # (ref to int_num_array, ref to int_num_array, swap)
354 # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
355 # subtract Y from X by modifying x in place
356 my ($c,$sx,$sy,$s) = @_;
358 my $car = 0; my $i; my $j = 0;
363 last unless defined $sy->[$j] || $car;
364 $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
366 # might leave leading zeros, so fix that
367 return __strip_zeros($sx);
371 # we can't do an early out if $x is < than $y, since we
372 # need to copy the high chunks from $y. Found by Bob Mathews.
373 #last unless defined $sy->[$j] || $car;
375 if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
378 # might leave leading zeros, so fix that
384 # (ref to int_num_array, ref to int_num_array)
385 # multiply two numbers in internal representation
386 # modifies first arg, second need not be different from first
387 my ($c,$xv,$yv) = @_;
391 # shortcut for two very short numbers (improved by Nathan Zook)
392 # works also if xv and yv are the same reference, and handles also $x == 0
395 if (($xv->[0] *= $yv->[0]) >= $MBASE)
397 $xv->[0] = $xv->[0] - ($xv->[1] = int($xv->[0] * $RBASE)) * $MBASE;
407 # multiply a large number a by a single element one, so speed up
408 my $y = $yv->[0]; my $car = 0;
411 $i = $i * $y + $car; $car = int($i * $RBASE); $i -= $car * $MBASE;
413 push @$xv, $car if $car != 0;
416 # shortcut for result $x == 0 => result = 0
417 return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
419 # since multiplying $x with $x fails, make copy in this case
420 $yv = [@$xv] if $xv == $yv; # same references?
422 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
431 # $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
433 # $prod - ($car = int($prod * RBASE)) * $MBASE; # see USE_MUL
435 # $prod[$cty] += $car if $car; # need really to check for 0?
439 # looping through this if $xi == 0 is silly - so optimize it away!
440 $xi = (shift @prod || 0), next if $xi == 0;
443 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
444 ## this is actually a tad slower
445 ## $prod = $prod[$cty]; $prod += ($car + $xi * $yi); # no ||0 here
447 $prod - ($car = int($prod * $RBASE)) * $MBASE; # see USE_MUL
449 $prod[$cty] += $car if $car; # need really to check for 0?
450 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
459 # (ref to int_num_array, ref to int_num_array)
460 # multiply two numbers in internal representation
461 # modifies first arg, second need not be different from first
462 my ($c,$xv,$yv) = @_;
466 # shortcut for two small numbers, also handles $x == 0
469 # shortcut for two very short numbers (improved by Nathan Zook)
470 # works also if xv and yv are the same reference, and handles also $x == 0
471 if (($xv->[0] *= $yv->[0]) >= $MBASE)
474 $xv->[0] - ($xv->[1] = int($xv->[0] / $MBASE)) * $MBASE;
484 # multiply a large number a by a single element one, so speed up
485 my $y = $yv->[0]; my $car = 0;
488 $i = $i * $y + $car; $car = int($i / $MBASE); $i -= $car * $MBASE;
490 push @$xv, $car if $car != 0;
493 # shortcut for result $x == 0 => result = 0
494 return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
496 # since multiplying $x with $x fails, make copy in this case
497 $yv = [@$xv] if $xv == $yv; # same references?
499 my @prod = (); my ($prod,$car,$cty,$xi,$yi);
503 # looping through this if $xi == 0 is silly - so optimize it away!
504 $xi = (shift @prod || 0), next if $xi == 0;
507 $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
509 $prod - ($car = int($prod / $MBASE)) * $MBASE;
511 $prod[$cty] += $car if $car; # need really to check for 0?
512 $xi = shift @prod || 0; # || 0 makes v5.005_3 happy
521 # ref to array, ref to array, modify first array and return remainder if
524 # see comments in _div_use_div() for more explanations
526 my ($c,$x,$yorg) = @_;
528 # the general div algorithmn here is about O(N*N) and thus quite slow, so
529 # we first check for some special cases and use shortcuts to handle them.
531 # This works, because we store the numbers in a chunked format where each
532 # element contains 5..7 digits (depending on system).
534 # if both numbers have only one element:
535 if (@$x == 1 && @$yorg == 1)
537 # shortcut, $yorg and $x are two small numbers
540 my $r = [ $x->[0] % $yorg->[0] ];
541 $x->[0] = int($x->[0] / $yorg->[0]);
546 $x->[0] = int($x->[0] / $yorg->[0]);
551 # if x has more than one, but y has only one element:
555 $rem = _mod($c,[ @$x ],$yorg) if wantarray;
557 # shortcut, $y is < $BASE
558 my $j = scalar @$x; my $r = 0;
559 my $y = $yorg->[0]; my $b;
562 $b = $r * $MBASE + $x->[$j];
563 $x->[$j] = int($b/$y);
566 pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
567 return ($x,$rem) if wantarray;
571 # now x and y have more than one element
573 # check whether y has more elements than x, if yet, the result will be 0
577 $rem = [@$x] if wantarray; # make copy
578 splice (@$x,1); # keep ref to original array
579 $x->[0] = 0; # set to 0
580 return ($x,$rem) if wantarray; # including remainder?
581 return $x; # only x, which is [0] now
583 # check whether the numbers have the same number of elements, in that case
584 # the result will fit into one element and can be computed efficiently
588 # if $yorg has more digits than $x (it's leading element is longer than
589 # the one from $x), the result will also be 0:
590 if (length(int($yorg->[-1])) > length(int($x->[-1])))
592 $rem = [@$x] if wantarray; # make copy
593 splice (@$x,1); # keep ref to org array
594 $x->[0] = 0; # set to 0
595 return ($x,$rem) if wantarray; # including remainder?
598 # now calculate $x / $yorg
599 if (length(int($yorg->[-1])) == length(int($x->[-1])))
601 # same length, so make full compare, and if equal, return 1
602 # hm, same lengths, but same contents? So we need to check all parts:
603 my $a = 0; my $j = scalar @$x - 1;
604 # manual way (abort if unequal, good for early ne)
607 last if ($a = $x->[$j] - $yorg->[$j]); $j--;
609 # $a contains the result of the compare between X and Y
610 # a < 0: x < y, a == 0 => x == y, a > 0: x > y
615 $rem = [ 0 ]; # a = 0 => x == y => rem 1
616 $rem = [@$x] if $a != 0; # a < 0 => x < y => rem = x
618 splice(@$x,1); # keep single element
619 $x->[0] = 0; # if $a < 0
625 return ($x,$rem) if wantarray;
628 # $x >= $y, proceed normally
634 my $y = [ @$yorg ]; # always make copy to preserve
636 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
638 $car = $bar = $prd = 0;
639 if (($dd = int($MBASE/($y->[-1]+1))) != 1)
643 $xi = $xi * $dd + $car;
644 $xi -= ($car = int($xi * $RBASE)) * $MBASE; # see USE_MUL
646 push(@$x, $car); $car = 0;
649 $yi = $yi * $dd + $car;
650 $yi -= ($car = int($yi * $RBASE)) * $MBASE; # see USE_MUL
657 @q = (); ($v2,$v1) = @$y[-2,-1];
661 ($u2,$u1,$u0) = @$x[-3..-1];
663 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
665 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
666 --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
669 ($car, $bar) = (0,0);
670 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
672 $prd = $q * $y->[$yi] + $car;
673 $prd -= ($car = int($prd * $RBASE)) * $MBASE; # see USE_MUL
674 $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
676 if ($x->[-1] < $car + $bar)
679 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
682 if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $MBASE));
695 for $xi (reverse @$x)
697 $prd = $car * $MBASE + $xi;
698 $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
719 # ref to array, ref to array, modify first array and return remainder if
721 my ($c,$x,$yorg) = @_;
723 # the general div algorithmn here is about O(N*N) and thus quite slow, so
724 # we first check for some special cases and use shortcuts to handle them.
726 # This works, because we store the numbers in a chunked format where each
727 # element contains 5..7 digits (depending on system).
729 # if both numbers have only one element:
730 if (@$x == 1 && @$yorg == 1)
732 # shortcut, $yorg and $x are two small numbers
735 my $r = [ $x->[0] % $yorg->[0] ];
736 $x->[0] = int($x->[0] / $yorg->[0]);
741 $x->[0] = int($x->[0] / $yorg->[0]);
745 # if x has more than one, but y has only one element:
749 $rem = _mod($c,[ @$x ],$yorg) if wantarray;
751 # shortcut, $y is < $BASE
752 my $j = scalar @$x; my $r = 0;
753 my $y = $yorg->[0]; my $b;
756 $b = $r * $MBASE + $x->[$j];
757 $x->[$j] = int($b/$y);
760 pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
761 return ($x,$rem) if wantarray;
764 # now x and y have more than one element
766 # check whether y has more elements than x, if yet, the result will be 0
770 $rem = [@$x] if wantarray; # make copy
771 splice (@$x,1); # keep ref to original array
772 $x->[0] = 0; # set to 0
773 return ($x,$rem) if wantarray; # including remainder?
774 return $x; # only x, which is [0] now
776 # check whether the numbers have the same number of elements, in that case
777 # the result will fit into one element and can be computed efficiently
781 # if $yorg has more digits than $x (it's leading element is longer than
782 # the one from $x), the result will also be 0:
783 if (length(int($yorg->[-1])) > length(int($x->[-1])))
785 $rem = [@$x] if wantarray; # make copy
786 splice (@$x,1); # keep ref to org array
787 $x->[0] = 0; # set to 0
788 return ($x,$rem) if wantarray; # including remainder?
791 # now calculate $x / $yorg
793 if (length(int($yorg->[-1])) == length(int($x->[-1])))
795 # same length, so make full compare, and if equal, return 1
796 # hm, same lengths, but same contents? So we need to check all parts:
797 my $a = 0; my $j = scalar @$x - 1;
798 # manual way (abort if unequal, good for early ne)
801 last if ($a = $x->[$j] - $yorg->[$j]); $j--;
803 # $a contains the result of the compare between X and Y
804 # a < 0: x < y, a == 0 => x == y, a > 0: x > y
809 $rem = [ 0 ]; # a = 0 => x == y => rem 1
810 $rem = [@$x] if $a != 0; # a < 0 => x < y => rem = x
812 splice(@$x,1); # keep single element
813 $x->[0] = 0; # if $a < 0
819 return ($x,$rem) if wantarray;
822 # $x >= $y, so proceed normally
828 my $y = [ @$yorg ]; # always make copy to preserve
830 my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
832 $car = $bar = $prd = 0;
833 if (($dd = int($MBASE/($y->[-1]+1))) != 1)
837 $xi = $xi * $dd + $car;
838 $xi -= ($car = int($xi / $MBASE)) * $MBASE;
840 push(@$x, $car); $car = 0;
843 $yi = $yi * $dd + $car;
844 $yi -= ($car = int($yi / $MBASE)) * $MBASE;
852 # @q will accumulate the final result, $q contains the current computed
853 # part of the final result
855 @q = (); ($v2,$v1) = @$y[-2,-1];
859 ($u2,$u1,$u0) = @$x[-3..-1];
861 #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
863 $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
864 --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
867 ($car, $bar) = (0,0);
868 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
870 $prd = $q * $y->[$yi] + $car;
871 $prd -= ($car = int($prd / $MBASE)) * $MBASE;
872 $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
874 if ($x->[-1] < $car + $bar)
877 for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
880 if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $MBASE));
884 pop(@$x); unshift(@q, $q);
892 for $xi (reverse @$x)
894 $prd = $car * $MBASE + $xi;
895 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
914 ##############################################################################
919 # internal absolute post-normalized compare (ignore signs)
920 # ref to array, ref to array, return <0, 0, >0
921 # arrays must have at least one entry; this is not checked for
922 my ($c,$cx,$cy) = @_;
924 # shortcut for short numbers
925 return (($cx->[0] <=> $cy->[0]) <=> 0)
926 if scalar @$cx == scalar @$cy && scalar @$cx == 1;
928 # fast comp based on number of array elements (aka pseudo-length)
929 my $lxy = (scalar @$cx - scalar @$cy)
930 # or length of first element if same number of elements (aka difference 0)
932 # need int() here because sometimes the last element is '00018' vs '18'
933 (length(int($cx->[-1])) - length(int($cy->[-1])));
934 return -1 if $lxy < 0; # already differs, ret
935 return 1 if $lxy > 0; # ditto
937 # manual way (abort if unequal, good for early ne)
938 my $a; my $j = scalar @$cx;
941 last if ($a = $cx->[$j] - $cy->[$j]);
948 # compute number of digits
950 # int() because add/sub sometimes leaves strings (like '00005') instead of
951 # '5' in this place, thus causing length() to report wrong length
954 (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
959 # return the nth digit, negative values count backward
960 # zero is rightmost, so _digit(123,0) will give 3
963 my $len = _len('',$x);
965 $n = $len+$n if $n < 0; # -1 last, -2 second-to-last
966 $n = abs($n); # if negative was too big
967 $len--; $n = $len if $n > $len; # n to big?
969 my $elem = int($n / $BASE_LEN); # which array element
970 my $digit = $n % $BASE_LEN; # which digit in this element
971 $elem = '0000'.@$x[$elem]; # get element padded with 0's
972 substr($elem,-$digit-1,1);
977 # return amount of trailing zeros in decimal
978 # check each array elem in _m for having 0 at end as long as elem == 0
979 # Upon finding a elem != 0, stop
982 return 0 if scalar @$x == 1 && $x->[0] == 0;
984 my $zeros = 0; my $elem;
989 $elem = "$e"; # preserve x
990 $elem =~ s/.*?(0*$)/$1/; # strip anything not zero
991 $zeros *= $BASE_LEN; # elems * 5
992 $zeros += length($elem); # count trailing zeros
995 $zeros ++; # real else branch: 50% slower!
1000 ##############################################################################
1005 # return true if arg is zero
1006 (((scalar @{$_[1]} == 1) && ($_[1]->[0] == 0))) <=> 0;
1011 # return true if arg is even
1012 (!($_[1]->[0] & 1)) <=> 0;
1017 # return true if arg is even
1018 (($_[1]->[0] & 1)) <=> 0;
1023 # return true if arg is one
1024 (scalar @{$_[1]} == 1) && ($_[1]->[0] == 1) <=> 0;
1029 # return true if arg is two
1030 (scalar @{$_[1]} == 1) && ($_[1]->[0] == 2) <=> 0;
1035 # return true if arg is ten
1036 (scalar @{$_[1]} == 1) && ($_[1]->[0] == 10) <=> 0;
1041 # internal normalization function that strips leading zeros from the array
1042 # args: ref to array
1045 my $cnt = scalar @$s; # get count of parts
1047 push @$s,0 if $i < 0; # div might return empty results, so fix it
1049 return $s if @$s == 1; # early out
1051 #print "strip: cnt $cnt i $i\n";
1052 # '0', '3', '4', '0', '0',
1057 # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
1058 # >= 1: skip first part (this can be zero)
1059 while ($i > 0) { last if $s->[$i] != 0; $i--; }
1060 $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
1064 ###############################################################################
1065 # check routine to test internal state for corruptions
1069 # used by the test suite
1072 return "$x is not a reference" if !ref($x);
1074 # are all parts are valid?
1075 my $i = 0; my $j = scalar @$x; my ($e,$try);
1078 $e = $x->[$i]; $e = 'undef' unless defined $e;
1079 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
1080 last if $e !~ /^[+]?[0-9]+$/;
1081 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
1082 last if "$e" !~ /^[+]?[0-9]+$/;
1083 $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
1084 last if '' . "$e" !~ /^[+]?[0-9]+$/;
1085 $try = ' < 0 || >= $BASE; '."($x, $e)";
1086 last if $e <0 || $e >= $BASE;
1087 # this test is disabled, since new/bnorm and certain ops (like early out
1088 # in add/sub) are allowed/expected to leave '00000' in some elements
1089 #$try = '=~ /^00+/; '."($x, $e)";
1090 #last if $e =~ /^00+/;
1093 return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
1098 ###############################################################################
1102 # if possible, use mod shortcut
1103 my ($c,$x,$yo) = @_;
1105 # slow way since $y to big
1106 if (scalar @$yo > 1)
1108 my ($xo,$rem) = _div($c,$x,$yo);
1113 # both are single element arrays
1114 if (scalar @$x == 1)
1120 # @y is a single element, but @x has more than one element
1124 # when BASE % Y == 0 then (B * BASE) % Y == 0
1125 # (B * BASE) % $y + A % Y => A % Y
1126 # so need to consider only last element: O(1)
1131 # else need to go through all elements: O(N), but loop is a bit simplified
1135 $r = ($r + $_) % $y; # not much faster, but heh...
1136 #$r += $_ % $y; $r %= $y;
1143 # else need to go through all elements: O(N)
1144 my $r = 0; my $bm = 1;
1147 $r = ($_ * $bm + $r) % $y;
1148 $bm = ($bm * $b) % $y;
1150 #$r += ($_ % $y) * $bm;
1158 splice (@$x,1); # keep one element of $x
1162 ##############################################################################
1167 my ($c,$x,$y,$n) = @_;
1171 $n = _new($c,$n); return _div($c,$x, _pow($c,$n,$y));
1174 # shortcut (faster) for shifting by 10)
1175 # multiples of $BASE_LEN
1176 my $dst = 0; # destination
1177 my $src = _num($c,$y); # as normal int
1178 my $xlen = (@$x-1)*$BASE_LEN+length(int($x->[-1])); # len of x in digits
1179 if ($src > $xlen or ($src == $xlen and ! defined $x->[1]))
1181 # 12345 67890 shifted right by more than 10 digits => 0
1182 splice (@$x,1); # leave only one element
1183 $x->[0] = 0; # set to zero
1186 my $rem = $src % $BASE_LEN; # remainder to shift
1187 $src = int($src / $BASE_LEN); # source
1190 splice (@$x,0,$src); # even faster, 38.4 => 39.3
1194 my $len = scalar @$x - $src; # elems to go
1195 my $vd; my $z = '0'x $BASE_LEN;
1196 $x->[scalar @$x] = 0; # avoid || 0 test inside loop
1199 $vd = $z.$x->[$src];
1200 $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
1202 $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
1203 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1204 $x->[$dst] = int($vd);
1207 splice (@$x,$dst) if $dst > 0; # kill left-over array elems
1208 pop @$x if $x->[-1] == 0 && @$x > 1; # kill last element if 0
1215 my ($c,$x,$y,$n) = @_;
1219 $n = _new($c,$n); return _mul($c,$x, _pow($c,$n,$y));
1222 # shortcut (faster) for shifting by 10) since we are in base 10eX
1223 # multiples of $BASE_LEN:
1224 my $src = scalar @$x; # source
1225 my $len = _num($c,$y); # shift-len as normal int
1226 my $rem = $len % $BASE_LEN; # remainder to shift
1227 my $dst = $src + int($len/$BASE_LEN); # destination
1228 my $vd; # further speedup
1229 $x->[$src] = 0; # avoid first ||0 for speed
1230 my $z = '0' x $BASE_LEN;
1233 $vd = $x->[$src]; $vd = $z.$vd;
1234 $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
1235 $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
1236 $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1237 $x->[$dst] = int($vd);
1240 # set lowest parts to 0
1241 while ($dst >= 0) { $x->[$dst--] = 0; }
1242 # fix spurios last zero element
1243 splice @$x,-1 if $x->[-1] == 0;
1250 # ref to array, ref to array, return ref to array
1251 my ($c,$cx,$cy) = @_;
1253 if (scalar @$cy == 1 && $cy->[0] == 0)
1255 splice (@$cx,1); $cx->[0] = 1; # y == 0 => x => 1
1258 if ((scalar @$cx == 1 && $cx->[0] == 1) || # x == 1
1259 (scalar @$cy == 1 && $cy->[0] == 1)) # or y == 1
1263 if (scalar @$cx == 1 && $cx->[0] == 0)
1265 splice (@$cx,1); $cx->[0] = 0; # 0 ** y => 0 (if not y <= 0)
1271 my $y_bin = _as_bin($c,$cy); $y_bin =~ s/^0b//;
1272 my $len = length($y_bin);
1275 _mul($c,$pow2,$cx) if substr($y_bin,$len,1) eq '1'; # is odd?
1286 # ref to array, return ref to array
1289 if ((@$cx == 1) && ($cx->[0] <= 2))
1291 $cx->[0] ||= 1; # 0 => 1, 1 => 1, 2 => 2
1295 # go forward until $base is exceeded
1296 # limit is either $x steps (steps == 100 means a result always too high) or
1298 my $steps = 100; $steps = $cx->[0] if @$cx == 1;
1299 my $r = 2; my $cf = 3; my $step = 2; my $last = $r;
1300 while ($r*$cf < $BASE && $step < $steps)
1302 $last = $r; $r *= $cf++; $step++;
1304 if ((@$cx == 1) && $step == $cx->[0])
1306 # completely done, so keep reference to $x and return
1311 # now we must do the left over steps
1312 my $n; # steps still to do
1313 if (scalar @$cx == 1)
1322 $cx->[0] = $last; splice (@$cx,1); # keep ref to $x
1323 my $zero_elements = 0;
1325 # do left-over steps fit into a scalar?
1326 if (ref $n eq 'ARRAY')
1328 # No, so use slower inc() & cmp()
1330 while (_acmp($step,$n) <= 0)
1332 # as soon as the last element of $cx is 0, we split it up and remember
1333 # how many zeors we got so far. The reason is that n! will accumulate
1334 # zeros at the end rather fast.
1337 $zero_elements ++; shift @$cx;
1339 _mul($c,$cx,$step); _inc($c,$step);
1344 # Yes, so we can speed it up slightly
1347 # When the last element of $cx is 0, we split it up and remember
1348 # how many we got so far. The reason is that n! will accumulate
1349 # zeros at the end rather fast.
1352 $zero_elements ++; shift @$cx;
1354 _mul($c,$cx,[$step]); $step++;
1357 # multiply in the zeros again
1358 while ($zero_elements-- > 0)
1362 $cx; # return result
1365 #############################################################################
1369 # calculate integer log of $x to base $base
1370 # ref to array, ref to array - return ref to array
1371 my ($c,$x,$base) = @_;
1374 return if (scalar @$x == 1 && $x->[0] == 0);
1375 # BASE 0 or 1 => NaN
1376 return if (scalar @$base == 1 && $base->[0] < 2);
1377 my $cmp = _acmp($c,$x,$base); # X == BASE => 1
1380 splice (@$x,1); $x->[0] = 1;
1386 splice (@$x,1); $x->[0] = 0;
1390 # this trial multiplication is very fast, even for large counts (like for
1391 # 2 ** 1024, since this still requires only 1024 very fast steps
1392 # (multiplication of a large number by a very small number is very fast))
1393 my $x_org = _copy($c,$x); # preserve x
1394 splice(@$x,1); $x->[0] = 1; # keep ref to $x
1396 my $trial = _copy($c,$base);
1398 # XXX TODO this only works if $base has only one element
1399 if (scalar @$base == 1)
1401 # compute int ( length_in_base_10(X) / ( log(base) / log(10) ) )
1402 my $len = _len($c,$x_org);
1403 my $res = int($len / (log($base->[0]) / log(10))) || 1; # avoid $res == 0
1406 $trial = _pow ($c, _copy($c, $base), $x);
1407 my $a = _acmp($x,$trial,$x_org);
1408 return ($x,1) if $a == 0;
1409 # we now know that $res is too small
1412 _mul($c,$trial,$base); _add($c, $x, [1]);
1417 _div($c,$trial,$base); _sub($c, $x, [1]);
1419 # did we now get the right result?
1420 $a = _acmp($x,$trial,$x_org);
1421 return ($x,1) if $a == 0; # yes, exactly
1425 _div($c,$trial,$base); _sub($c, $x, [1]);
1429 # simple loop that increments $x by two in each step, possible overstepping
1430 # the real result by one
1433 my $base_mul = _mul($c, _copy($c,$base), $base);
1435 while (($a = _acmp($c,$trial,$x_org)) < 0)
1437 _mul($c,$trial,$base_mul); _add($c, $x, [2]);
1443 # overstepped the result
1445 _div($c,$trial,$base);
1446 $a = _acmp($c,$trial,$x_org);
1451 $exact = 0 if $a != 0;
1454 ($x,$exact); # return result
1458 use constant DEBUG => 0;
1460 sub steps { $steps };
1464 # square-root of $x in place
1465 # Compute a guess of the result (by rule of thumb), then improve it via
1469 if (scalar @$x == 1)
1471 # fit's into one Perl scalar, so result can be computed directly
1472 $x->[0] = int(sqrt($x->[0]));
1475 my $y = _copy($c,$x);
1476 # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
1477 # since our guess will "grow"
1478 my $l = int((_len($c,$x)-1) / 2);
1480 my $lastelem = $x->[-1]; # for guess
1481 my $elems = scalar @$x - 1;
1482 # not enough digits, but could have more?
1483 if ((length($lastelem) <= 3) && ($elems > 1))
1485 # right-align with zero pad
1486 my $len = length($lastelem) & 1;
1487 print "$lastelem => " if DEBUG;
1488 $lastelem .= substr($x->[-2] . '0' x $BASE_LEN,0,$BASE_LEN);
1489 # former odd => make odd again, or former even to even again
1490 $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
1491 print "$lastelem\n" if DEBUG;
1494 # construct $x (instead of _lsft($c,$x,$l,10)
1495 my $r = $l % $BASE_LEN; # 10000 00000 00000 00000 ($BASE_LEN=5)
1496 $l = int($l / $BASE_LEN);
1497 print "l = $l " if DEBUG;
1499 splice @$x,$l; # keep ref($x), but modify it
1501 # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
1503 # 14400 00000 => sqrt(14400) => guess first digits to be 120
1504 # 144000 000000 => sqrt(144000) => guess 379
1506 print "$lastelem (elems $elems) => " if DEBUG;
1507 $lastelem = $lastelem / 10 if ($elems & 1 == 1); # odd or even?
1508 my $g = sqrt($lastelem); $g =~ s/\.//; # 2.345 => 2345
1509 $r -= 1 if $elems & 1 == 0; # 70 => 7
1511 # padd with zeros if result is too short
1512 $x->[$l--] = int(substr($g . '0' x $r,0,$r+1));
1513 print "now ",$x->[-1] if DEBUG;
1514 print " would have been ", int('1' . '0' x $r),"\n" if DEBUG;
1516 # If @$x > 1, we could compute the second elem of the guess, too, to create
1517 # an even better guess. Not implemented yet. Does it improve performance?
1518 $x->[$l--] = 0 while ($l >= 0); # all other digits of guess are zero
1520 print "start x= ",_str($c,$x),"\n" if DEBUG;
1523 my $lastlast = _zero();
1524 $steps = 0 if DEBUG;
1525 while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
1528 $lastlast = _copy($c,$last);
1529 $last = _copy($c,$x);
1530 _add($c,$x, _div($c,_copy($c,$y),$x));
1532 print " x= ",_str($c,$x),"\n" if DEBUG;
1534 print "\nsteps in sqrt: $steps, " if DEBUG;
1535 _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0; # overshot?
1536 print " final ",$x->[-1],"\n" if DEBUG;
1542 # take n'th root of $x in place (n >= 3)
1545 if (scalar @$x == 1)
1549 # result will always be smaller than 2 so trunc to 1 at once
1554 # fit's into one Perl scalar, so result can be computed directly
1555 # cannot use int() here, because it rounds wrongly (try
1556 # (81 ** 3) ** (1/3) to see what I mean)
1557 #$x->[0] = int( $x->[0] ** (1 / $n->[0]) );
1558 # round to 8 digits, then truncate result to integer
1559 $x->[0] = int ( sprintf ("%.8f", $x->[0] ** (1 / $n->[0]) ) );
1564 # we know now that X is more than one element long
1566 # if $n is a power of two, we can repeatedly take sqrt($X) and find the
1567 # proper result, because sqrt(sqrt($x)) == root($x,4)
1568 my $b = _as_bin($c,$n);
1569 if ($b =~ /0b1(0+)$/)
1571 my $count = CORE::length($1); # 0b100 => len('00') => 2
1572 my $cnt = $count; # counter for loop
1573 unshift (@$x, 0); # add one element, together with one
1574 # more below in the loop this makes 2
1577 # 'inflate' $X by adding one element, basically computing
1578 # $x * $BASE * $BASE. This gives us more $BASE_LEN digits for result
1579 # since len(sqrt($X)) approx == len($x) / 2.
1581 # calculate sqrt($x), $x is now one element to big, again. In the next
1582 # round we make that two, again.
1585 # $x is now one element to big, so truncate result by removing it
1590 # trial computation by starting with 2,4,8,16 etc until we overstep
1594 # while still to do more than X steps
1598 while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
1600 _mul ($c, $step, [2]);
1601 _add ($c, $trial, $step);
1605 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) == 0)
1607 @$x = @$trial; # make copy while preserving ref to $x
1610 # overstepped, so go back on step
1611 _sub($c, $trial, $step);
1612 } while (scalar @$step > 1 || $step->[0] > 128);
1616 # add two, because $trial cannot be exactly the result (otherwise we would
1617 # alrady have found it)
1618 _add($c, $trial, $step);
1620 # and now add more and more (2,4,6,8,10 etc)
1621 while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
1623 _add ($c, $trial, $step);
1626 # hit not exactly? (overstepped)
1627 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
1632 # hit not exactly? (overstepped)
1633 # 80 too small, 81 slightly too big, 82 too big
1634 if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
1639 @$x = @$trial; # make copy while preserving ref to $x
1645 ##############################################################################
1652 # the shortcut makes equal, large numbers _really_ fast, and makes only a
1653 # very small performance drop for small numbers (e.g. something with less
1654 # than 32 bit) Since we optimize for large numbers, this is enabled.
1655 return $x if _acmp($c,$x,$y) == 0; # shortcut
1657 my $m = _one(); my ($xr,$yr);
1658 my $mask = $AND_MASK;
1661 my $y1 = _copy($c,$y); # make copy
1665 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1667 ($x1, $xr) = _div($c,$x1,$mask);
1668 ($y1, $yr) = _div($c,$y1,$mask);
1670 # make ints() from $xr, $yr
1671 # this is when the AND_BITS are greater than $BASE and is slower for
1672 # small (<256 bits) numbers, but faster for large numbers. Disabled
1673 # due to KISS principle
1675 # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1676 # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1677 # _add($c,$x, _mul($c, _new( $c, ($xrr & $yrr) ), $m) );
1679 # 0+ due to '&' doesn't work in strings
1680 _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
1690 return _zero() if _acmp($c,$x,$y) == 0; # shortcut (see -and)
1692 my $m = _one(); my ($xr,$yr);
1693 my $mask = $XOR_MASK;
1696 my $y1 = _copy($c,$y); # make copy
1700 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1702 ($x1, $xr) = _div($c,$x1,$mask);
1703 ($y1, $yr) = _div($c,$y1,$mask);
1704 # make ints() from $xr, $yr (see _and())
1705 #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1706 #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1707 #_add($c,$x, _mul($c, _new( $c, ($xrr ^ $yrr) ), $m) );
1709 # 0+ due to '^' doesn't work in strings
1710 _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
1713 # the loop stops when the shorter of the two numbers is exhausted
1714 # the remainder of the longer one will survive bit-by-bit, so we simple
1715 # multiply-add it in
1716 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1717 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1726 return $x if _acmp($c,$x,$y) == 0; # shortcut (see _and)
1728 my $m = _one(); my ($xr,$yr);
1729 my $mask = $OR_MASK;
1732 my $y1 = _copy($c,$y); # make copy
1736 while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1738 ($x1, $xr) = _div($c,$x1,$mask);
1739 ($y1, $yr) = _div($c,$y1,$mask);
1740 # make ints() from $xr, $yr (see _and())
1741 # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1742 # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1743 # _add($c,$x, _mul($c, _new( $c, ($xrr | $yrr) ), $m) );
1745 # 0+ due to '|' doesn't work in strings
1746 _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
1749 # the loop stops when the shorter of the two numbers is exhausted
1750 # the remainder of the longer one will survive bit-by-bit, so we simple
1751 # multiply-add it in
1752 _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1753 _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1760 # convert a decimal number to hex (ref to array, return ref to string)
1763 # fit's into one element (handle also 0x0 case)
1766 my $t = sprintf("0x%x",$x->[0]);
1770 my $x1 = _copy($c,$x);
1773 my ($xr, $h, $x10000);
1776 $x10000 = [ 0x10000 ]; $h = 'h4';
1780 $x10000 = [ 0x1000 ]; $h = 'h3';
1782 # while (! _is_zero($c,$x1))
1783 while (@$x1 != 1 || $x1->[0] != 0) # _is_zero()
1785 ($x1, $xr) = _div($c,$x1,$x10000);
1786 $es .= unpack($h,pack('v',$xr->[0])); # XXX TODO: why pack('v',...)?
1789 $es =~ s/^[0]+//; # strip leading zeros
1796 # convert a decimal number to bin (ref to array, return ref to string)
1799 # fit's into one element (and Perl recent enough), handle also 0b0 case
1800 # handle zero case for older Perls
1801 if ($] <= 5.005 && @$x == 1 && $x->[0] == 0)
1803 my $t = '0b0'; return $t;
1805 if (@$x == 1 && $] >= 5.006)
1807 my $t = sprintf("0b%b",$x->[0]);
1810 my $x1 = _copy($c,$x);
1813 my ($xr, $b, $x10000);
1816 $x10000 = [ 0x10000 ]; $b = 'b16';
1820 $x10000 = [ 0x1000 ]; $b = 'b12';
1822 # while (! _is_zero($c,$x1))
1823 while (!(@$x1 == 1 && $x1->[0] == 0)) # _is_zero()
1825 ($x1, $xr) = _div($c,$x1,$x10000);
1826 $es .= unpack($b,pack('v',$xr->[0])); # XXX TODO: why pack('v',...)?
1827 # $es .= unpack($b,$xr->[0]);
1830 $es =~ s/^[0]+//; # strip leading zeros
1837 # convert a hex number to decimal (ref to string, return ref to array)
1841 my $m = [ 0x10000 ]; # 16 bit at a time
1844 my $len = length($hs)-2;
1845 $len = int($len/4); # 4-digit parts, w/o '0x'
1846 my $val; my $i = -4;
1849 $val = substr($hs,$i,4);
1850 $val =~ s/^[+-]?0x// if $len == 0; # for last part only because
1851 $val = hex($val); # hex does not like wrong chars
1853 _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1854 _mul ($c, $mul, $m ) if $len >= 0; # skip last mul
1861 # convert a hex number to decimal (ref to string, return ref to array)
1864 # instead of converting X (8) bit at a time, it is faster to "convert" the
1865 # number to hex, and then call _from_hex.
1868 $hs =~ s/^[+-]?0b//; # remove sign and 0b
1869 my $l = length($hs); # bits
1870 $hs = '0' x (8-($l % 8)) . $hs if ($l % 8) != 0; # padd left side w/ 0
1871 my $h = unpack('H*', pack ('B*', $hs)); # repack as hex
1873 $c->_from_hex('0x'.$h);
1876 ##############################################################################
1877 # special modulus functions
1884 my $u = _zero($c); my $u1 = _one($c);
1885 my $a = _copy($c,$y); my $b = _copy($c,$x);
1887 # Euclid's Algorithm for bgcd(), only that we calc bgcd() ($a) and the
1888 # result ($u) at the same time. See comments in BigInt for why this works.
1890 ($a, $q, $b) = ($b, _div($c,$a,$b)); # step 1
1892 while (!_is_zero($c,$b))
1894 my $t = _add($c, # step 2:
1895 _mul($c,_copy($c,$u1), $q) , # t = u1 * q
1897 $u = $u1; # u = u1, u1 = t
1900 ($a, $q, $b) = ($b, _div($c,$a,$b)); # step 1
1903 # if the gcd is not 1, then return NaN
1904 return (undef,undef) unless _is_one($c,$a);
1906 $sign = $sign == 1 ? '+' : '-';
1912 # modulus of power ($x ** $y) % $z
1913 my ($c,$num,$exp,$mod) = @_;
1915 # in the trivial case,
1916 if (_is_one($c,$mod))
1918 splice @$num,0,1; $num->[0] = 0;
1921 if ((scalar @$num == 1) && (($num->[0] == 0) || ($num->[0] == 1)))
1927 # $num = _mod($c,$num,$mod); # this does not make it faster
1929 my $acc = _copy($c,$num); my $t = _one();
1931 my $expbin = _as_bin($c,$exp); $expbin =~ s/^0b//;
1932 my $len = length($expbin);
1935 if ( substr($expbin,$len,1) eq '1') # is_odd
1938 $t = _mod($c,$t,$mod);
1941 $acc = _mod($c,$acc,$mod);
1949 # greatest common divisor
1952 while (! _is_zero($c,$y))
1954 my $t = _copy($c,$y);
1955 $y = _mod($c, $x, $y);
1961 ##############################################################################
1962 ##############################################################################
1969 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
1973 Provides support for big integer calculations. Not intended to be used by other
1974 modules. Other modules which sport the same functions can also be used to support
1975 Math::BigInt, like Math::BigInt::GMP or Math::BigInt::Pari.
1979 In order to allow for multiple big integer libraries, Math::BigInt was
1980 rewritten to use library modules for core math routines. Any module which
1981 follows the same API as this can be used instead by using the following:
1983 use Math::BigInt lib => 'libname';
1985 'libname' is either the long name ('Math::BigInt::Pari'), or only the short
1986 version like 'Pari'.
1992 The following functions MUST be defined in order to support the use by
1993 Math::BigInt v1.70 or later:
1995 api_version() return API version, minimum 1 for v1.70
1996 _new(string) return ref to new object from ref to decimal string
1997 _zero() return a new object with value 0
1998 _one() return a new object with value 1
1999 _two() return a new object with value 2
2000 _ten() return a new object with value 10
2002 _str(obj) return ref to a string representing the object
2003 _num(obj) returns a Perl integer/floating point number
2004 NOTE: because of Perl numeric notation defaults,
2005 the _num'ified obj may lose accuracy due to
2006 machine-dependend floating point size limitations
2008 _add(obj,obj) Simple addition of two objects
2009 _mul(obj,obj) Multiplication of two objects
2010 _div(obj,obj) Division of the 1st object by the 2nd
2011 In list context, returns (result,remainder).
2012 NOTE: this is integer math, so no
2013 fractional part will be returned.
2014 The second operand will be not be 0, so no need to
2016 _sub(obj,obj) Simple subtraction of 1 object from another
2017 a third, optional parameter indicates that the params
2018 are swapped. In this case, the first param needs to
2019 be preserved, while you can destroy the second.
2020 sub (x,y,1) => return x - y and keep x intact!
2021 _dec(obj) decrement object by one (input is garant. to be > 0)
2022 _inc(obj) increment object by one
2025 _acmp(obj,obj) <=> operator for objects (return -1, 0 or 1)
2027 _len(obj) returns count of the decimal digits of the object
2028 _digit(obj,n) returns the n'th decimal digit of object
2030 _is_one(obj) return true if argument is 1
2031 _is_two(obj) return true if argument is 2
2032 _is_ten(obj) return true if argument is 10
2033 _is_zero(obj) return true if argument is 0
2034 _is_even(obj) return true if argument is even (0,2,4,6..)
2035 _is_odd(obj) return true if argument is odd (1,3,5,7..)
2037 _copy return a ref to a true copy of the object
2039 _check(obj) check whether internal representation is still intact
2040 return 0 for ok, otherwise error message as string
2042 _from_hex(str) return ref to new object from ref to hexadecimal string
2043 _from_bin(str) return ref to new object from ref to binary string
2045 _as_hex(str) return string containing the value as
2046 unsigned hex string, with the '0x' prepended.
2047 Leading zeros must be stripped.
2048 _as_bin(str) Like as_hex, only as binary string containing only
2049 zeros and ones. Leading zeros must be stripped and a
2050 '0b' must be prepended.
2052 _rsft(obj,N,B) shift object in base B by N 'digits' right
2053 _lsft(obj,N,B) shift object in base B by N 'digits' left
2055 _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
2056 Note: XOR, AND and OR pad with zeros if size mismatches
2057 _and(obj1,obj2) AND (bit-wise) object 1 with object 2
2058 _or(obj1,obj2) OR (bit-wise) object 1 with object 2
2060 _mod(obj,obj) Return remainder of div of the 1st by the 2nd object
2061 _sqrt(obj) return the square root of object (truncated to int)
2062 _root(obj) return the n'th (n >= 3) root of obj (truncated to int)
2063 _fac(obj) return factorial of object 1 (1*2*3*4..)
2064 _pow(obj,obj) return object 1 to the power of object 2
2065 return undef for NaN
2066 _zeros(obj) return number of trailing decimal zeros
2067 _modinv return inverse modulus
2068 _modpow return modulus of power ($x ** $y) % $z
2069 _log_int(X,N) calculate integer log() of X in base N
2070 X >= 0, N >= 0 (return undef for NaN)
2071 returns (RESULT, EXACT) where EXACT is:
2072 1 : result is exactly RESULT
2073 0 : result was truncated to RESULT
2074 undef : unknown whether result is exactly RESULT
2075 _gcd(obj,obj) return Greatest Common Divisor of two objects
2077 The following functions are optional, and can be defined if the underlying lib
2078 has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
2079 slow) fallback routines to emulate these:
2086 Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
2089 So the library needs only to deal with unsigned big integers. Testing of input
2090 parameter validity is done by the caller, so you need not worry about
2091 underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by zero or similar
2094 The first parameter can be modified, that includes the possibility that you
2095 return a reference to a completely different object instead. Although keeping
2096 the reference and just changing it's contents is prefered over creating and
2097 returning a different reference.
2099 Return values are always references to objects, strings, or true/false for
2100 comparisation routines.
2102 =head1 WRAP YOUR OWN
2104 If you want to port your own favourite c-lib for big numbers to the
2105 Math::BigInt interface, you can take any of the already existing modules as
2106 a rough guideline. You should really wrap up the latest BigInt and BigFloat
2107 testsuites with your module, and replace in them any of the following:
2113 use Math::BigInt lib => 'yourlib';
2115 This way you ensure that your library really works 100% within Math::BigInt.
2119 This program is free software; you may redistribute it and/or modify it under
2120 the same terms as Perl itself.
2124 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
2126 Seperated from BigInt and shaped API with the help of John Peacock.
2127 Fixed, sped-up and enhanced by Tels http://bloodgate.com 2001-2003.
2128 Further streamlining (api_version 1) by Tels 2004.
2132 L<Math::BigInt>, L<Math::BigFloat>, L<Math::BigInt::BitVect>,
2133 L<Math::BigInt::GMP>, L<Math::BigInt::FastCalc> and L<Math::BigInt::Pari>.