#$@%#$%!

Ngiler SH3LL 360
: / usr / share / perl5 / vendor_perl / Math / BigInt /

Information Server
MySQL : OFF Perl : OFF CURL : ON WGET : OFF PKEXEC : OFF

Directive	Local Value
IP Address	162.0.209.101
System	Linux premium169.web-hosting.com 4.18.0-553.50.1.lve.el8.x86_64 #1 SMP Thu Apr 17 19:10:24 UTC 2025 x86_64
User	luxeelil
PHP Version	8.0.30
Software	LiteSpeed
Doc root	Writable

Edit File : Calc.pm \| Size : 75.06 KB
package Math::BigInt::Calc; use 5.006001; use strict; use warnings; use Carp; use Math::BigInt::Lib; our $VERSION = '1.999811'; our @ISA = ('Math::BigInt::Lib'); # Package to store unsigned big integers in decimal and do math with them # Internally the numbers are stored in an array with at least 1 element, no # leading zero parts (except the first) and in base 1eX where X is determined # automatically at loading time to be the maximum possible value # todo: # - fully remove funky $# stuff in div() (maybe - that code scares me...) # USE_MUL: due to problems on certain os (os390, posix-bc) "* 1e-5" is used # instead of "/ 1e5" at some places, (marked with USE_MUL). Other platforms # BS2000, some Crays need USE_DIV instead. # The BEGIN block is used to determine which of the two variants gives the # correct result. # Beware of things like: # $i = $i * $y + $car; $car = int($i / $BASE); $i = $i % $BASE; # This works on x86, but fails on ARM (SA1100, iPAQ) due to who knows what # reasons. So, use this instead (slower, but correct): # $i = $i * $y + $car; $car = int($i / $BASE); $i -= $BASE * $car; ############################################################################## # global constants, flags and accessory # announce that we are compatible with MBI v1.83 and up sub api_version () { 2; } # constants for easier life my ($BASE, $BASE_LEN, $RBASE, $MAX_VAL); my ($AND_BITS, $XOR_BITS, $OR_BITS); my ($AND_MASK, $XOR_MASK, $OR_MASK); sub _base_len { # Set/get the BASE_LEN and assorted other, related values. # Used only by the testsuite, the set variant is used only by the BEGIN # block below: my ($class, $b, $int) = @_; if (defined $b) { # avoid redefinitions undef &_mul; undef &_div; if ($] >= 5.008 && $int && $b > 7) { $BASE_LEN = $b; _mul = \&_mul_use_div_64; _div = \&_div_use_div_64; $BASE = int("1e" . $BASE_LEN); $MAX_VAL = $BASE-1; return $BASE_LEN unless wantarray; return ($BASE_LEN, $BASE, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN, $MAX_VAL); } # find whether we can use mul or div in mul()/div() $BASE_LEN = $b + 1; my $caught = 0; while (--$BASE_LEN > 5) { $BASE = int("1e" . $BASE_LEN); $RBASE = abs('1e-' . $BASE_LEN); # see USE_MUL $caught = 0; $caught += 1 if (int($BASE * $RBASE) != 1); # should be 1 $caught += 2 if (int($BASE / $BASE) != 1); # should be 1 last if $caught != 3; } $BASE = int("1e" . $BASE_LEN); $RBASE = abs('1e-' . $BASE_LEN); # see USE_MUL $MAX_VAL = $BASE-1; # ($caught & 1) != 0 => cannot use MUL # ($caught & 2) != 0 => cannot use DIV if ($caught == 2) # 2 { # must USE_MUL since we cannot use DIV _mul = \&_mul_use_mul; _div = \&_div_use_mul; } else # 0 or 1 { # can USE_DIV instead _mul = \&_mul_use_div; _div = \&_div_use_div; } } return $BASE_LEN unless wantarray; return ($BASE_LEN, $BASE, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN, $MAX_VAL); } sub _new { # Given a string representing an integer, returns a reference to an array # of integers, where each integer represents a chunk of the original input # integer. my ($class, $str) = @_; #unless ($str =~ /^([1-9]\d\|0)\z/) { # require Carp; # Carp::croak("Invalid input string '$str'"); #} my $input_len = length($str) - 1; # Shortcut for small numbers. return bless [ $str ], $class if $input_len < $BASE_LEN; my $format = "a" . (($input_len % $BASE_LEN) + 1); $format .= $] < 5.008 ? "a$BASE_LEN" x int($input_len / $BASE_LEN) : "(a$BASE_LEN)"; my $self = [ reverse(map { 0 + $_ } unpack($format, $str)) ]; return bless $self, $class; } BEGIN { # from Daniel Pfeiffer: determine largest group of digits that is precisely # multipliable with itself plus carry # Test now changed to expect the proper pattern, not a result off by 1 or 2 my ($e, $num) = 3; # lowest value we will use is 3+1-1 = 3 do { $num = '9' x ++$e; $num = $num + 1; } while $num =~ /9{$e}0{$e}/; # must be a certain pattern $e--; # last test failed, so retract one step # the limits below brush the problems with the test above under the rug: # the test should be able to find the proper $e automatically $e = 5 if $^O =~ /^uts/; # UTS get's some special treatment $e = 5 if $^O =~ /^unicos/; # unicos is also problematic (6 seems to work # there, but we play safe) my $int = 0; if ($e > 7) { use integer; my $e1 = 7; $num = 7; do { $num = ('9' x ++$e1) + 0; $num = $num + 1; } while ("$num" =~ /9{$e1}0{$e1}/); # must be a certain pattern $e1--; # last test failed, so retract one step if ($e1 > 7) { $int = 1; $e = $e1; } } __PACKAGE__ -> _base_len($e, $int); # set and store use integer; # find out how many bits _and, _or and _xor can take (old default = 16) # I don't think anybody has yet 128 bit scalars, so let's play safe. local $^W = 0; # don't warn about 'nonportable number' $AND_BITS = 15; $XOR_BITS = 15; $OR_BITS = 15; # find max bits, we will not go higher than numberofbits that fit into $BASE # to make _and etc simpler (and faster for smaller, slower for large numbers) my $max = 16; while (2 $max < $BASE) { $max++; } { no integer; $max = 16 if $] < 5.006; # older Perls might not take >16 too well } my ($x, $y, $z); do { $AND_BITS++; $x = CORE::oct('0b' . '1' x $AND_BITS); $y = $x & $x; $z = (2 $AND_BITS) - 1; } while ($AND_BITS < $max && $x == $z && $y == $x); $AND_BITS --; # retreat one step do { $XOR_BITS++; $x = CORE::oct('0b' . '1' x $XOR_BITS); $y = $x ^ 0; $z = (2 $XOR_BITS) - 1; } while ($XOR_BITS < $max && $x == $z && $y == $x); $XOR_BITS --; # retreat one step do { $OR_BITS++; $x = CORE::oct('0b' . '1' x $OR_BITS); $y = $x \| $x; $z = (2 $OR_BITS) - 1; } while ($OR_BITS < $max && $x == $z && $y == $x); $OR_BITS--; # retreat one step $AND_MASK = __PACKAGE__->_new(( 2 $AND_BITS )); $XOR_MASK = __PACKAGE__->_new(( 2 $XOR_BITS )); $OR_MASK = __PACKAGE__->_new(( 2 ** $OR_BITS )); # We can compute the approximate length no faster than the real length: _alen = \&_len; } ############################################################################### sub _zero { # create a zero my $class = shift; return bless [ 0 ], $class; } sub _one { # create a one my $class = shift; return bless [ 1 ], $class; } sub _two { # create a two my $class = shift; return bless [ 2 ], $class; } sub _ten { # create a 10 my $class = shift; bless [ 10 ], $class; } sub _1ex { # create a 1Ex my $class = shift; my $rem = $_[0] % $BASE_LEN; # remainder my $parts = $_[0] / $BASE_LEN; # parts # 000000, 000000, 100 bless [ (0) x $parts, '1' . ('0' x $rem) ], $class; } sub _copy { # make a true copy my $class = shift; return bless [ @{ $_[0] } ], $class; } # catch and throw away sub import { } ############################################################################## # convert back to string and number sub _str { # Convert number from internal base 1eN format to string format. Internal # format is always normalized, i.e., no leading zeros. my $ary = $_[1]; my $idx = $#$ary; # index of last element if ($idx < 0) { # should not happen require Carp; Carp::croak("$_[1] has no elements"); } # Handle first one differently, since it should not have any leading zeros. my $ret = int($ary->[$idx]); if ($idx > 0) { # Interestingly, the pre-padd method uses more time. # The old grep variant takes longer (14 vs. 10 sec). my $z = '0' x ($BASE_LEN - 1); while (--$idx >= 0) { $ret .= substr($z . $ary->[$idx], -$BASE_LEN); } } $ret; } sub _num { # Make a Perl scalar number (int/float) from a BigInt object. my $x = $_[1]; return $x->[0] if @$x == 1; # below $BASE # Start with the most significant element and work towards the least # significant element. Avoid multiplying "inf" (which happens if the number # overflows) with "0" (if there are zero elements in $x) since this gives # "nan" which propagates to the output. my $num = 0; for (my $i = $#$x ; $i >= 0 ; --$i) { $num = $BASE; $num += $x -> [$i]; } return $num; } ############################################################################## # actual math code sub _add { # (ref to int_num_array, ref to int_num_array) # # Routine to add two base 1eX numbers stolen from Knuth Vol 2 Algorithm A # pg 231. There are separate routines to add and sub as per Knuth pg 233. # This routine modifies array x, but not y. my ($c, $x, $y) = @_; # $x + 0 => $x return $x if @$y == 1 && $y->[0] == 0; # 0 + $y => $y->copy if (@$x == 1 && $x->[0] == 0) { @$x = @$y; return $x; } # For each in Y, add Y to X and carry. If after that, something is left in # X, foreach in X add carry to X and then return X, carry. Trades one # "$j++" for having to shift arrays. my $i; my $car = 0; my $j = 0; for $i (@$y) { $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0; $j++; } while ($car != 0) { $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++; } $x; } sub _inc { # (ref to int_num_array, ref to int_num_array) # Add 1 to $x, modify $x in place my ($c, $x) = @_; for my $i (@$x) { return $x if ($i += 1) < $BASE; # early out $i = 0; # overflow, next } push @$x, 1 if $x->[-1] == 0; # last overflowed, so extend $x; } sub _dec { # (ref to int_num_array, ref to int_num_array) # Sub 1 from $x, modify $x in place my ($c, $x) = @_; my $MAX = $BASE - 1; # since MAX_VAL based on BASE for my $i (@$x) { last if ($i -= 1) >= 0; # early out $i = $MAX; # underflow, next } pop @$x if $x->[-1] == 0 && @$x > 1; # last underflowed (but leave 0) $x; } sub _sub { # (ref to int_num_array, ref to int_num_array, swap) # # Subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y # subtract Y from X by modifying x in place my ($c, $sx, $sy, $s) = @_; my $car = 0; my $i; my $j = 0; if (!$s) { for $i (@$sx) { last unless defined $sy->[$j] \|\| $car; $i += $BASE if $car = (($i -= ($sy->[$j] \|\| 0) + $car) < 0); $j++; } # might leave leading zeros, so fix that return __strip_zeros($sx); } for $i (@$sx) { # We can't do an early out if $x < $y, since we need to copy the high # chunks from $y. Found by Bob Mathews. #last unless defined $sy->[$j] \|\| $car; $sy->[$j] += $BASE if $car = ($sy->[$j] = $i - ($sy->[$j] \|\| 0) - $car) < 0; $j++; } # might leave leading zeros, so fix that __strip_zeros($sy); } sub _mul_use_mul { # (ref to int_num_array, ref to int_num_array) # multiply two numbers in internal representation # modifies first arg, second need not be different from first my ($c, $xv, $yv) = @_; if (@$yv == 1) { # shortcut for two very short numbers (improved by Nathan Zook) # works also if xv and yv are the same reference, and handles also $x == 0 if (@$xv == 1) { if (($xv->[0] = $yv->[0]) >= $BASE) { $xv->[0] = $xv->[0] - ($xv->[1] = int($xv->[0] $RBASE)) * $BASE; } ; return $xv; } # $x * 0 => 0 if ($yv->[0] == 0) { @$xv = (0); return $xv; } # multiply a large number a by a single element one, so speed up my $y = $yv->[0]; my $car = 0; foreach my $i (@$xv) { $i = $i * $y + $car; $car = int($i * $RBASE); $i -= $car * $BASE; } push @$xv, $car if $car != 0; return $xv; } # shortcut for result $x == 0 => result = 0 return $xv if @$xv == 1 && $xv->[0] == 0; # since multiplying $x with $x fails, make copy in this case $yv = [ @$xv ] if $xv == $yv; # same references? my @prod = (); my ($prod, $car, $cty, $xi, $yi); for $xi (@$xv) { $car = 0; $cty = 0; # slow variant # for $yi (@$yv) # { # $prod = $xi * $yi + ($prod[$cty] \|\| 0) + $car; # $prod[$cty++] = # $prod - ($car = int($prod * RBASE)) * $BASE; # see USE_MUL # } # $prod[$cty] += $car if $car; # need really to check for 0? # $xi = shift @prod; # faster variant # looping through this if $xi == 0 is silly - so optimize it away! $xi = (shift @prod \|\| 0), next if $xi == 0; for $yi (@$yv) { $prod = $xi * $yi + ($prod[$cty] \|\| 0) + $car; ## this is actually a tad slower ## $prod = $prod[$cty]; $prod += ($car + $xi * $yi); # no \|\|0 here $prod[$cty++] = $prod - ($car = int($prod * $RBASE)) * $BASE; # see USE_MUL } $prod[$cty] += $car if $car; # need really to check for 0? $xi = shift @prod \|\| 0; # \|\| 0 makes v5.005_3 happy } push @$xv, @prod; # can't have leading zeros # __strip_zeros($xv); $xv; } sub _mul_use_div_64 { # (ref to int_num_array, ref to int_num_array) # multiply two numbers in internal representation # modifies first arg, second need not be different from first # works for 64 bit integer with "use integer" my ($c, $xv, $yv) = @_; use integer; if (@$yv == 1) { # shortcut for two small numbers, also handles $x == 0 if (@$xv == 1) { # shortcut for two very short numbers (improved by Nathan Zook) # works also if xv and yv are the same reference, and handles also $x == 0 if (($xv->[0] = $yv->[0]) >= $BASE) { $xv->[0] = $xv->[0] - ($xv->[1] = $xv->[0] / $BASE) $BASE; } return $xv; } # $x * 0 => 0 if ($yv->[0] == 0) { @$xv = (0); return $xv; } # multiply a large number a by a single element one, so speed up my $y = $yv->[0]; my $car = 0; foreach my $i (@$xv) { #$i = $i * $y + $car; $car = $i / $BASE; $i -= $car * $BASE; $i = $i * $y + $car; $i -= ($car = $i / $BASE) * $BASE; } push @$xv, $car if $car != 0; return $xv; } # shortcut for result $x == 0 => result = 0 return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) ); # since multiplying $x with $x fails, make copy in this case $yv = $c->_copy($xv) if $xv == $yv; # same references? my @prod = (); my ($prod, $car, $cty, $xi, $yi); for $xi (@$xv) { $car = 0; $cty = 0; # looping through this if $xi == 0 is silly - so optimize it away! $xi = (shift @prod \|\| 0), next if $xi == 0; for $yi (@$yv) { $prod = $xi * $yi + ($prod[$cty] \|\| 0) + $car; $prod[$cty++] = $prod - ($car = $prod / $BASE) * $BASE; } $prod[$cty] += $car if $car; # need really to check for 0? $xi = shift @prod \|\| 0; # \|\| 0 makes v5.005_3 happy } push @$xv, @prod; $xv; } sub _mul_use_div { # (ref to int_num_array, ref to int_num_array) # multiply two numbers in internal representation # modifies first arg, second need not be different from first my ($c, $xv, $yv) = @_; if (@$yv == 1) { # shortcut for two small numbers, also handles $x == 0 if (@$xv == 1) { # shortcut for two very short numbers (improved by Nathan Zook) # works also if xv and yv are the same reference, and handles also $x == 0 if (($xv->[0] = $yv->[0]) >= $BASE) { $xv->[0] = $xv->[0] - ($xv->[1] = int($xv->[0] / $BASE)) $BASE; } ; return $xv; } # $x * 0 => 0 if ($yv->[0] == 0) { @$xv = (0); return $xv; } # multiply a large number a by a single element one, so speed up my $y = $yv->[0]; my $car = 0; foreach my $i (@$xv) { $i = $i * $y + $car; $car = int($i / $BASE); $i -= $car * $BASE; # This (together with use integer;) does not work on 32-bit Perls #$i = $i * $y + $car; $i -= ($car = $i / $BASE) * $BASE; } push @$xv, $car if $car != 0; return $xv; } # shortcut for result $x == 0 => result = 0 return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) ); # since multiplying $x with $x fails, make copy in this case $yv = $c->_copy($xv) if $xv == $yv; # same references? my @prod = (); my ($prod, $car, $cty, $xi, $yi); for $xi (@$xv) { $car = 0; $cty = 0; # looping through this if $xi == 0 is silly - so optimize it away! $xi = (shift @prod \|\| 0), next if $xi == 0; for $yi (@$yv) { $prod = $xi * $yi + ($prod[$cty] \|\| 0) + $car; $prod[$cty++] = $prod - ($car = int($prod / $BASE)) * $BASE; } $prod[$cty] += $car if $car; # need really to check for 0? $xi = shift @prod \|\| 0; # \|\| 0 makes v5.005_3 happy } push @$xv, @prod; # can't have leading zeros # __strip_zeros($xv); $xv; } sub _div_use_mul { # ref to array, ref to array, modify first array and return remainder if # in list context # see comments in _div_use_div() for more explanations my ($c, $x, $yorg) = @_; # the general div algorithm here is about O(NN) and thus quite slow, so # we first check for some special cases and use shortcuts to handle them. # This works, because we store the numbers in a chunked format where each # element contains 5..7 digits (depending on system). # if both numbers have only one element: if (@$x == 1 && @$yorg == 1) { # shortcut, $yorg and $x are two small numbers if (wantarray) { my $rem = [ $x->[0] % $yorg->[0] ]; bless $rem, $c; $x->[0] = int($x->[0] / $yorg->[0]); return ($x, $rem); } else { $x->[0] = int($x->[0] / $yorg->[0]); return $x; } } # if x has more than one, but y has only one element: if (@$yorg == 1) { my $rem; $rem = $c->_mod($c->_copy($x), $yorg) if wantarray; # shortcut, $y is < $BASE my $j = @$x; my $r = 0; my $y = $yorg->[0]; my $b; while ($j-- > 0) { $b = $r $BASE + $x->[$j]; $x->[$j] = int($b/$y); $r = $b % $y; } pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero return ($x, $rem) if wantarray; return $x; } # now x and y have more than one element # check whether y has more elements than x, if yet, the result will be 0 if (@$yorg > @$x) { my $rem; $rem = $c->_copy($x) if wantarray; # make copy @$x = 0; # set to 0 return ($x, $rem) if wantarray; # including remainder? return $x; # only x, which is [0] now } # check whether the numbers have the same number of elements, in that case # the result will fit into one element and can be computed efficiently if (@$yorg == @$x) { # if $yorg has more digits than $x (it's leading element is longer than # the one from $x), the result will also be 0: if (length(int($yorg->[-1])) > length(int($x->[-1]))) { my $rem = $c->_copy($x) if wantarray; # make copy @$x = 0; # set to 0 return ($x, $rem) if wantarray; # including remainder? return $x; } # now calculate $x / $yorg if (length(int($yorg->[-1])) == length(int($x->[-1]))) { # same length, so make full compare my $a = 0; my $j = @$x - 1; # manual way (abort if unequal, good for early ne) while ($j >= 0) { last if ($a = $x->[$j] - $yorg->[$j]); $j--; } # $a contains the result of the compare between X and Y # a < 0: x < y, a == 0: x == y, a > 0: x > y if ($a <= 0) { # a = 0 => x == y => rem 0 # a < 0 => x < y => rem = x my $rem = $a == 0 ? $c->_zero() : $c->_copy($x); @$x = 0; # if $a < 0 $x->[0] = 1 if $a == 0; # $x == $y return ($x, $rem) if wantarray; return $x; } # $x >= $y, so proceed normally } } # all other cases: my $y = $c->_copy($yorg); # always make copy to preserve my ($car, $bar, $prd, $dd, $xi, $yi, @q, $v2, $v1, @d, $tmp, $q, $u2, $u1, $u0); $car = $bar = $prd = 0; if (($dd = int($BASE / ($y->[-1] + 1))) != 1) { for $xi (@$x) { $xi = $xi * $dd + $car; $xi -= ($car = int($xi * $RBASE)) * $BASE; # see USE_MUL } push(@$x, $car); $car = 0; for $yi (@$y) { $yi = $yi * $dd + $car; $yi -= ($car = int($yi * $RBASE)) * $BASE; # see USE_MUL } } else { push(@$x, 0); } @q = (); ($v2, $v1) = @$y[-2, -1]; $v2 = 0 unless $v2; while ($#$x > $#$y) { ($u2, $u1, $u0) = @$x[-3 .. -1]; $u2 = 0 unless $u2; #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n" # if $v1 == 0; $q = (($u0 == $v1) ? $MAX_VAL : int(($u0 * $BASE + $u1) / $v1)); --$q while ($v2 * $q > ($u0 * $BASE + $u1 - $q * $v1) * $BASE + $u2); if ($q) { ($car, $bar) = (0, 0); for ($yi = 0, $xi = $#$x - $#$y-1; $yi <= $#$y; ++$yi, ++$xi) { $prd = $q * $y->[$yi] + $car; $prd -= ($car = int($prd * $RBASE)) * $BASE; # see USE_MUL $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0)); } if ($x->[-1] < $car + $bar) { $car = 0; --$q; for ($yi = 0, $xi = $#$x - $#$y-1; $yi <= $#$y; ++$yi, ++$xi) { $x->[$xi] -= $BASE if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE)); } } } pop(@$x); unshift(@q, $q); } if (wantarray) { my $d = bless [], $c; if ($dd != 1) { $car = 0; for $xi (reverse @$x) { $prd = $car * $BASE + $xi; $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL unshift(@$d, $tmp); } } else { @$d = @$x; } @$x = @q; __strip_zeros($x); __strip_zeros($d); return ($x, $d); } @$x = @q; __strip_zeros($x); $x; } sub _div_use_div_64 { # ref to array, ref to array, modify first array and return remainder if # in list context # This version works on 64 bit integers my ($c, $x, $yorg) = @_; use integer; # the general div algorithm here is about O(NN) and thus quite slow, so # we first check for some special cases and use shortcuts to handle them. # This works, because we store the numbers in a chunked format where each # element contains 5..7 digits (depending on system). # if both numbers have only one element: if (@$x == 1 && @$yorg == 1) { # shortcut, $yorg and $x are two small numbers if (wantarray) { my $rem = [ $x->[0] % $yorg->[0] ]; bless $rem, $c; $x->[0] = int($x->[0] / $yorg->[0]); return ($x, $rem); } else { $x->[0] = int($x->[0] / $yorg->[0]); return $x; } } # if x has more than one, but y has only one element: if (@$yorg == 1) { my $rem; $rem = $c->_mod($c->_copy($x), $yorg) if wantarray; # shortcut, $y is < $BASE my $j = @$x; my $r = 0; my $y = $yorg->[0]; my $b; while ($j-- > 0) { $b = $r $BASE + $x->[$j]; $x->[$j] = int($b/$y); $r = $b % $y; } pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero return ($x, $rem) if wantarray; return $x; } # now x and y have more than one element # check whether y has more elements than x, if yet, the result will be 0 if (@$yorg > @$x) { my $rem; $rem = $c->_copy($x) if wantarray; # make copy @$x = 0; # set to 0 return ($x, $rem) if wantarray; # including remainder? return $x; # only x, which is [0] now } # check whether the numbers have the same number of elements, in that case # the result will fit into one element and can be computed efficiently if (@$yorg == @$x) { my $rem; # if $yorg has more digits than $x (it's leading element is longer than # the one from $x), the result will also be 0: if (length(int($yorg->[-1])) > length(int($x->[-1]))) { $rem = $c->_copy($x) if wantarray; # make copy @$x = 0; # set to 0 return ($x, $rem) if wantarray; # including remainder? return $x; } # now calculate $x / $yorg if (length(int($yorg->[-1])) == length(int($x->[-1]))) { # same length, so make full compare my $a = 0; my $j = @$x - 1; # manual way (abort if unequal, good for early ne) while ($j >= 0) { last if ($a = $x->[$j] - $yorg->[$j]); $j--; } # $a contains the result of the compare between X and Y # a < 0: x < y, a == 0: x == y, a > 0: x > y if ($a <= 0) { $rem = $c->_zero(); # a = 0 => x == y => rem 0 $rem = $c->_copy($x) if $a != 0; # a < 0 => x < y => rem = x @$x = 0; # if $a < 0 $x->[0] = 1 if $a == 0; # $x == $y return ($x, $rem) if wantarray; # including remainder? return $x; } # $x >= $y, so proceed normally } } # all other cases: my $y = $c->_copy($yorg); # always make copy to preserve my ($car, $bar, $prd, $dd, $xi, $yi, @q, $v2, $v1, @d, $tmp, $q, $u2, $u1, $u0); $car = $bar = $prd = 0; if (($dd = int($BASE / ($y->[-1] + 1))) != 1) { for $xi (@$x) { $xi = $xi * $dd + $car; $xi -= ($car = int($xi / $BASE)) * $BASE; } push(@$x, $car); $car = 0; for $yi (@$y) { $yi = $yi * $dd + $car; $yi -= ($car = int($yi / $BASE)) * $BASE; } } else { push(@$x, 0); } # @q will accumulate the final result, $q contains the current computed # part of the final result @q = (); ($v2, $v1) = @$y[-2, -1]; $v2 = 0 unless $v2; while ($#$x > $#$y) { ($u2, $u1, $u0) = @$x[-3..-1]; $u2 = 0 unless $u2; #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n" # if $v1 == 0; $q = (($u0 == $v1) ? $MAX_VAL : int(($u0 * $BASE + $u1) / $v1)); --$q while ($v2 * $q > ($u0 * $BASE +$ u1- $q$v1) $BASE + $u2); if ($q) { ($car, $bar) = (0, 0); for ($yi = 0, $xi = $#$x - $#$y - 1; $yi <= $#$y; ++$yi, ++$xi) { $prd = $q * $y->[$yi] + $car; $prd -= ($car = int($prd / $BASE)) * $BASE; $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0)); } if ($x->[-1] < $car + $bar) { $car = 0; --$q; for ($yi = 0, $xi = $#$x - $#$y - 1; $yi <= $#$y; ++$yi, ++$xi) { $x->[$xi] -= $BASE if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE)); } } } pop(@$x); unshift(@q, $q); } if (wantarray) { my $d = bless [], $c; if ($dd != 1) { $car = 0; for $xi (reverse @$x) { $prd = $car * $BASE + $xi; $car = $prd - ($tmp = int($prd / $dd)) * $dd; unshift(@$d, $tmp); } } else { @$d = @$x; } @$x = @q; __strip_zeros($x); __strip_zeros($d); return ($x, $d); } @$x = @q; __strip_zeros($x); $x; } sub _div_use_div { # ref to array, ref to array, modify first array and return remainder if # in list context my ($c, $x, $yorg) = @_; # the general div algorithm here is about O(NN) and thus quite slow, so # we first check for some special cases and use shortcuts to handle them. # This works, because we store the numbers in a chunked format where each # element contains 5..7 digits (depending on system). # if both numbers have only one element: if (@$x == 1 && @$yorg == 1) { # shortcut, $yorg and $x are two small numbers if (wantarray) { my $rem = [ $x->[0] % $yorg->[0] ]; bless $rem, $c; $x->[0] = int($x->[0] / $yorg->[0]); return ($x, $rem); } else { $x->[0] = int($x->[0] / $yorg->[0]); return $x; } } # if x has more than one, but y has only one element: if (@$yorg == 1) { my $rem; $rem = $c->_mod($c->_copy($x), $yorg) if wantarray; # shortcut, $y is < $BASE my $j = @$x; my $r = 0; my $y = $yorg->[0]; my $b; while ($j-- > 0) { $b = $r $BASE + $x->[$j]; $x->[$j] = int($b/$y); $r = $b % $y; } pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero return ($x, $rem) if wantarray; return $x; } # now x and y have more than one element # check whether y has more elements than x, if yet, the result will be 0 if (@$yorg > @$x) { my $rem; $rem = $c->_copy($x) if wantarray; # make copy @$x = 0; # set to 0 return ($x, $rem) if wantarray; # including remainder? return $x; # only x, which is [0] now } # check whether the numbers have the same number of elements, in that case # the result will fit into one element and can be computed efficiently if (@$yorg == @$x) { my $rem; # if $yorg has more digits than $x (it's leading element is longer than # the one from $x), the result will also be 0: if (length(int($yorg->[-1])) > length(int($x->[-1]))) { $rem = $c->_copy($x) if wantarray; # make copy @$x = 0; # set to 0 return ($x, $rem) if wantarray; # including remainder? return $x; } # now calculate $x / $yorg if (length(int($yorg->[-1])) == length(int($x->[-1]))) { # same length, so make full compare my $a = 0; my $j = @$x - 1; # manual way (abort if unequal, good for early ne) while ($j >= 0) { last if ($a = $x->[$j] - $yorg->[$j]); $j--; } # $a contains the result of the compare between X and Y # a < 0: x < y, a == 0: x == y, a > 0: x > y if ($a <= 0) { $rem = $c->_zero(); # a = 0 => x == y => rem 0 $rem = $c->_copy($x) if $a != 0; # a < 0 => x < y => rem = x @$x = 0; $x->[0] = 0; # if $a < 0 $x->[0] = 1 if $a == 0; # $x == $y return ($x, $rem) if wantarray; # including remainder? return $x; } # $x >= $y, so proceed normally } } # all other cases: my $y = $c->_copy($yorg); # always make copy to preserve my ($car, $bar, $prd, $dd, $xi, $yi, @q, $v2, $v1, @d, $tmp, $q, $u2, $u1, $u0); $car = $bar = $prd = 0; if (($dd = int($BASE / ($y->[-1] + 1))) != 1) { for $xi (@$x) { $xi = $xi * $dd + $car; $xi -= ($car = int($xi / $BASE)) * $BASE; } push(@$x, $car); $car = 0; for $yi (@$y) { $yi = $yi * $dd + $car; $yi -= ($car = int($yi / $BASE)) * $BASE; } } else { push(@$x, 0); } # @q will accumulate the final result, $q contains the current computed # part of the final result @q = (); ($v2, $v1) = @$y[-2, -1]; $v2 = 0 unless $v2; while ($#$x > $#$y) { ($u2, $u1, $u0) = @$x[-3..-1]; $u2 = 0 unless $u2; #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n" # if $v1 == 0; $q = (($u0 == $v1) ? $MAX_VAL : int(($u0 * $BASE + $u1) / $v1)); --$q while ($v2 * $q > ($u0 * $BASE + $u1 - $q * $v1) * $BASE + $u2); if ($q) { ($car, $bar) = (0, 0); for ($yi = 0, $xi = $#$x - $#$y - 1; $yi <= $#$y; ++$yi, ++$xi) { $prd = $q * $y->[$yi] + $car; $prd -= ($car = int($prd / $BASE)) * $BASE; $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0)); } if ($x->[-1] < $car + $bar) { $car = 0; --$q; for ($yi = 0, $xi = $#$x - $#$y - 1; $yi <= $#$y; ++$yi, ++$xi) { $x->[$xi] -= $BASE if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE)); } } } pop(@$x); unshift(@q, $q); } if (wantarray) { my $d = bless [], $c; if ($dd != 1) { $car = 0; for $xi (reverse @$x) { $prd = $car * $BASE + $xi; $car = $prd - ($tmp = int($prd / $dd)) * $dd; unshift(@$d, $tmp); } } else { @$d = @$x; } @$x = @q; __strip_zeros($x); __strip_zeros($d); return ($x, $d); } @$x = @q; __strip_zeros($x); $x; } ############################################################################## # testing sub _acmp { # Internal absolute post-normalized compare (ignore signs) # ref to array, ref to array, return <0, 0, >0 # Arrays must have at least one entry; this is not checked for. my ($c, $cx, $cy) = @_; # shortcut for short numbers return (($cx->[0] <=> $cy->[0]) <=> 0) if @$cx == 1 && @$cy == 1; # fast comp based on number of array elements (aka pseudo-length) my $lxy = (@$cx - @$cy) # or length of first element if same number of elements (aka difference 0) \|\| # need int() here because sometimes the last element is '00018' vs '18' (length(int($cx->[-1])) - length(int($cy->[-1]))); return -1 if $lxy < 0; # already differs, ret return 1 if $lxy > 0; # ditto # manual way (abort if unequal, good for early ne) my $a; my $j = @$cx; while (--$j >= 0) { last if $a = $cx->[$j] - $cy->[$j]; } $a <=> 0; } sub _len { # compute number of digits in base 10 # int() because add/sub sometimes leaves strings (like '00005') instead of # '5' in this place, thus causing length() to report wrong length my $cx = $_[1]; (@$cx - 1) * $BASE_LEN + length(int($cx->[-1])); } sub _digit { # Return the nth digit. Zero is rightmost, so _digit(123, 0) gives 3. # Negative values count from the left, so _digit(123, -1) gives 1. my ($c, $x, $n) = @_; my $len = _len('', $x); $n += $len if $n < 0; # -1 last, -2 second-to-last # Math::BigInt::Calc returns 0 if N is out of range, but this is not done # by the other backend libraries. return "0" if $n < 0 \|\| $n >= $len; # return 0 for digits out of range my $elem = int($n / $BASE_LEN); # index of array element my $digit = $n % $BASE_LEN; # index of digit within the element substr("0" x $BASE_LEN . "$x->[$elem]", -1 - $digit, 1); } sub _zeros { # Return number of trailing zeros in decimal. # Check each array element for having 0 at end as long as elem == 0 # Upon finding a elem != 0, stop. my $x = $_[1]; return 0 if @$x == 1 && $x->[0] == 0; my $zeros = 0; foreach my $elem (@$x) { if ($elem != 0) { $elem =~ /[^0](0)\z/; $zeros += length($1); # count trailing zeros last; # early out } $zeros += $BASE_LEN; } $zeros; } ############################################################################## # _is_ routines sub _is_zero { # return true if arg is zero @{$_[1]} == 1 && $_[1]->[0] == 0 ? 1 : 0; } sub _is_even { # return true if arg is even $_[1]->[0] & 1 ? 0 : 1; } sub _is_odd { # return true if arg is odd $_[1]->[0] & 1 ? 1 : 0; } sub _is_one { # return true if arg is one @{$_[1]} == 1 && $_[1]->[0] == 1 ? 1 : 0; } sub _is_two { # return true if arg is two @{$_[1]} == 1 && $_[1]->[0] == 2 ? 1 : 0; } sub _is_ten { # return true if arg is ten @{$_[1]} == 1 && $_[1]->[0] == 10 ? 1 : 0; } sub __strip_zeros { # Internal normalization function that strips leading zeros from the array. # Args: ref to array my $x = shift; push @$x, 0 if @$x == 0; # div might return empty results, so fix it return $x if @$x == 1; # early out #print "strip: cnt $cnt i $i\n"; # '0', '3', '4', '0', '0', # 0 1 2 3 4 # cnt = 5, i = 4 # i = 4 # i = 3 # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos) # >= 1: skip first part (this can be zero) my $i = $#$x; while ($i > 0) { last if $x->[$i] != 0; $i--; } $i++; splice(@$x, $i) if $i < @$x; $x; } ############################################################################### # check routine to test internal state for corruptions sub _check { # used by the test suite my ($class, $x) = @_; my $msg = $class -> SUPER::_check($x); return $msg if $msg; my $n; eval { $n = @$x }; return "Not an array reference" unless $@ eq ''; return "Reference to an empty array" unless $n > 0; # The following fails with Math::BigInt::FastCalc because a # Math::BigInt::FastCalc "object" is an unblessed array ref. # #return 0 unless ref($x) eq $class; for (my $i = 0 ; $i <= $#$x ; ++ $i) { my $e = $x -> [$i]; return "Element at index $i is undefined" unless defined $e; return "Element at index $i is a '" . ref($e) . "', which is not a scalar" unless ref($e) eq ""; # It would be better to use the regex /^([1-9]\d\|0)\z/, but that fails # in Math::BigInt::FastCalc, because it sometimes creates array # elements like "000000". return "Element at index $i is '$e', which does not look like an" . " normal integer" unless $e =~ /^\d+\z/; return "Element at index $i is '$e', which is not smaller than" . " the base '$BASE'" if $e >= $BASE; return "Element at index $i (last element) is zero" if $#$x > 0 && $i == $#$x && $e == 0; } return 0; } ############################################################################### sub _mod { # if possible, use mod shortcut my ($c, $x, $yo) = @_; # slow way since $y too big if (@$yo > 1) { my ($xo, $rem) = $c->_div($x, $yo); @$x = @$rem; return $x; } my $y = $yo->[0]; # if both are single element arrays if (@$x == 1) { $x->[0] %= $y; return $x; } # if @$x has more than one element, but @$y is a single element my $b = $BASE % $y; if ($b == 0) { # when BASE % Y == 0 then (B BASE) % Y == 0 # (B * BASE) % $y + A % Y => A % Y # so need to consider only last element: O(1) $x->[0] %= $y; } elsif ($b == 1) { # else need to go through all elements in @$x: O(N), but loop is a bit # simplified my $r = 0; foreach (@$x) { $r = ($r + $_) % $y; # not much faster, but heh... #$r += $_ % $y; $r %= $y; } $r = 0 if $r == $y; $x->[0] = $r; } else { # else need to go through all elements in @$x: O(N) my $r = 0; my $bm = 1; foreach (@$x) { $r = ($_ * $bm + $r) % $y; $bm = ($bm * $b) % $y; #$r += ($_ % $y) * $bm; #$bm = $b; #$bm %= $y; #$r %= $y; } $r = 0 if $r == $y; $x->[0] = $r; } @$x = $x->[0]; # keep one element of @$x return $x; } ############################################################################## # shifts sub _rsft { my ($c, $x, $y, $n) = @_; if ($n != 10) { $n = $c->_new($n); return scalar $c->_div($x, $c->_pow($n, $y)); } # shortcut (faster) for shifting by 10) # multiples of $BASE_LEN my $dst = 0; # destination my $src = $c->_num($y); # as normal int my $xlen = (@$x - 1) $BASE_LEN + length(int($x->[-1])); if ($src >= $xlen or ($src == $xlen and !defined $x->[1])) { # 12345 67890 shifted right by more than 10 digits => 0 splice(@$x, 1); # leave only one element $x->[0] = 0; # set to zero return $x; } my $rem = $src % $BASE_LEN; # remainder to shift $src = int($src / $BASE_LEN); # source if ($rem == 0) { splice(@$x, 0, $src); # even faster, 38.4 => 39.3 } else { my $len = @$x - $src; # elems to go my $vd; my $z = '0' x $BASE_LEN; $x->[ @$x ] = 0; # avoid \|\| 0 test inside loop while ($dst < $len) { $vd = $z . $x->[$src]; $vd = substr($vd, -$BASE_LEN, $BASE_LEN - $rem); $src++; $vd = substr($z . $x->[$src], -$rem, $rem) . $vd; $vd = substr($vd, -$BASE_LEN, $BASE_LEN) if length($vd) > $BASE_LEN; $x->[$dst] = int($vd); $dst++; } splice(@$x, $dst) if $dst > 0; # kill left-over array elems pop @$x if $x->[-1] == 0 && @$x > 1; # kill last element if 0 } # else rem == 0 $x; } sub _lsft { my ($c, $x, $n, $b) = @_; return $x if $c->_is_zero($x); # Handle the special case when the base is a power of 10. Don't check # whether log($b)/log(10) is an integer, because log(1000)/log(10) is not # exactly 3. my $log10 = sprintf "%.0f", log($b) / log(10); if ($b == 10 $log10) { $b = 10; $n = $c->_mul($n, $c->_new($log10)); # shortcut (faster) for shifting by 10) since we are in base 10eX # multiples of $BASE_LEN: my $src = @$x; # source my $len = $c->_num($n); # shift-len as normal int my $rem = $len % $BASE_LEN; # remainder to shift my $dst = $src + int($len / $BASE_LEN); # destination my $vd; # further speedup $x->[$src] = 0; # avoid first \|\|0 for speed my $z = '0' x $BASE_LEN; while ($src >= 0) { $vd = $x->[$src]; $vd = $z . $vd; $vd = substr($vd, -$BASE_LEN + $rem, $BASE_LEN - $rem); $vd .= $src > 0 ? substr($z . $x->[$src - 1], -$BASE_LEN, $rem) : '0' x $rem; $vd = substr($vd, -$BASE_LEN, $BASE_LEN) if length($vd) > $BASE_LEN; $x->[$dst] = int($vd); $dst--; $src--; } # set lowest parts to 0 while ($dst >= 0) { $x->[$dst--] = 0; } # fix spurious last zero element splice @$x, -1 if $x->[-1] == 0; return $x; } else { $b = $c->_new($b); #print $c->_str($b); return $c->_mul($x, $c->_pow($b, $n)); } } sub _pow { # power of $x to $y # ref to array, ref to array, return ref to array my ($c, $cx, $cy) = @_; if (@$cy == 1 && $cy->[0] == 0) { splice(@$cx, 1); $cx->[0] = 1; # y == 0 => x => 1 return $cx; } if ((@$cx == 1 && $cx->[0] == 1) \|\| # x == 1 (@$cy == 1 && $cy->[0] == 1)) # or y == 1 { return $cx; } if (@$cx == 1 && $cx->[0] == 0) { splice (@$cx, 1); $cx->[0] = 0; # 0 y => 0 (if not y <= 0) return $cx; } my $pow2 = $c->_one(); my $y_bin = $c->_as_bin($cy); $y_bin =~ s/^0b//; my $len = length($y_bin); while (--$len > 0) { $c->_mul($pow2, $cx) if substr($y_bin, $len, 1) eq '1'; # is odd? $c->_mul($cx, $cx); } $c->_mul($cx, $pow2); $cx; } sub _nok { # Return binomial coefficient (n over k). # Given refs to arrays, return ref to array. # First input argument is modified. my ($c, $n, $k) = @_; # If k > n/2, or, equivalently, 2k > n, compute nok(n, k) as # nok(n, n-k), to minimize the number if iterations in the loop. { my $twok = $c->_mul($c->_two(), $c->_copy($k)); # 2 k if ($c->_acmp($twok, $n) > 0) { # if 2k > n $k = $c->_sub($c->_copy($n), $k); # k = n - k } } # Example: # # / 7 \ 7! 1234 * 567 5 * 6 * 7 6 7 # \| \| = --------- = --------------- = --------- = 5 * - * - # \ 3 / (7-3)! 3! 1234 123 1 * 2 * 3 2 3 if ($c->_is_zero($k)) { @$n = 1; } else { # Make a copy of the original n, since we'll be modifying n in-place. my $n_orig = $c->_copy($n); # n = 5, f = 6, d = 2 (cf. example above) $c->_sub($n, $k); $c->_inc($n); my $f = $c->_copy($n); $c->_inc($f); my $d = $c->_two(); # while f <= n (the original n, that is) ... while ($c->_acmp($f, $n_orig) <= 0) { # n = (n * f / d) == 5 * 6 / 2 (cf. example above) $c->_mul($n, $f); $c->_div($n, $d); # f = 7, d = 3 (cf. example above) $c->_inc($f); $c->_inc($d); } } return $n; } my @factorials = ( 1, 1, 2, 23, 234, 2345, 23456, 234567, ); sub _fac { # factorial of $x # ref to array, return ref to array my ($c, $cx) = @_; if ((@$cx == 1) && ($cx->[0] <= 7)) { $cx->[0] = $factorials[$cx->[0]]; # 0 => 1, 1 => 1, 2 => 2 etc. return $cx; } if ((@$cx == 1) && # we do this only if $x >= 12 and $x <= 7000 ($cx->[0] >= 12 && $cx->[0] < 7000)) { # Calculate (k-j) (k-j+1) ... k .. (k+j-1) * (k + j) # See http://blogten.blogspot.com/2007/01/calculating-n.html # The above series can be expressed as factors: # k * k - (j - i) * 2 # We cache kk, and calculate (j j) as the sum of the first j odd integers # This will not work when N exceeds the storage of a Perl scalar, however, # in this case the algorithm would be way too slow to terminate, anyway. # As soon as the last element of $cx is 0, we split it up and remember # how many zeors we got so far. The reason is that n! will accumulate # zeros at the end rather fast. my $zero_elements = 0; # If n is even, set n = n -1 my $k = $c->_num($cx); my $even = 1; if (($k & 1) == 0) { $even = $k; $k --; } # set k to the center point $k = ($k + 1) / 2; # print "k $k even: $even\n"; # now calculate k * k my $k2 = $k * $k; my $odd = 1; my $sum = 1; my $i = $k - 1; # keep reference to x my $new_x = $c->_new($k * $even); @$cx = @$new_x; if ($cx->[0] == 0) { $zero_elements ++; shift @$cx; } # print STDERR "x = ", $c->_str($cx), "\n"; my $BASE2 = int(sqrt($BASE))-1; my $j = 1; while ($j <= $i) { my $m = ($k2 - $sum); $odd += 2; $sum += $odd; $j++; while ($j <= $i && ($m < $BASE2) && (($k2 - $sum) < $BASE2)) { $m = ($k2 - $sum); $odd += 2; $sum += $odd; $j++; # print STDERR "\n k2 $k2 m $m sum $sum odd $odd\n"; sleep(1); } if ($m < $BASE) { $c->_mul($cx, [$m]); } else { $c->_mul($cx, $c->_new($m)); } if ($cx->[0] == 0) { $zero_elements ++; shift @$cx; } # print STDERR "Calculate $k2 - $sum = $m (x = ", $c->_str($cx), ")\n"; } # multiply in the zeros again unshift @$cx, (0) x $zero_elements; return $cx; } # go forward until $base is exceeded limit is either $x steps (steps == 100 # means a result always too high) or $base. my $steps = 100; $steps = $cx->[0] if @$cx == 1; my $r = 2; my $cf = 3; my $step = 2; my $last = $r; while ($r $cf < $BASE && $step < $steps) { $last = $r; $r = $cf++; $step++; } if ((@$cx == 1) && $step == $cx->[0]) { # completely done, so keep reference to $x and return $cx->[0] = $r; return $cx; } # now we must do the left over steps my $n; # steps still to do if (@$cx == 1) { $n = $cx->[0]; } else { $n = $c->_copy($cx); } # Set $cx to the last result below $BASE (but keep ref to $x) $cx->[0] = $last; splice (@$cx, 1); # As soon as the last element of $cx is 0, we split it up and remember # how many zeors we got so far. The reason is that n! will accumulate # zeros at the end rather fast. my $zero_elements = 0; # do left-over steps fit into a scalar? if (ref $n eq 'ARRAY') { # No, so use slower inc() & cmp() # ($n is at least $BASE here) my $base_2 = int(sqrt($BASE)) - 1; #print STDERR "base_2: $base_2\n"; while ($step < $base_2) { if ($cx->[0] == 0) { $zero_elements ++; shift @$cx; } my $b = $step ($step + 1); $step += 2; $c->_mul($cx, [$b]); } $step = [$step]; while ($c->_acmp($step, $n) <= 0) { if ($cx->[0] == 0) { $zero_elements ++; shift @$cx; } $c->_mul($cx, $step); $c->_inc($step); } } else { # Yes, so we can speed it up slightly # print "# left over steps $n\n"; my $base_4 = int(sqrt(sqrt($BASE))) - 2; #print STDERR "base_4: $base_4\n"; my $n4 = $n - 4; while ($step < $n4 && $step < $base_4) { if ($cx->[0] == 0) { $zero_elements ++; shift @$cx; } my $b = $step * ($step + 1); $step += 2; $b = $step ($step + 1); $step += 2; $c->_mul($cx, [$b]); } my $base_2 = int(sqrt($BASE)) - 1; my $n2 = $n - 2; #print STDERR "base_2: $base_2\n"; while ($step < $n2 && $step < $base_2) { if ($cx->[0] == 0) { $zero_elements ++; shift @$cx; } my $b = $step * ($step + 1); $step += 2; $c->_mul($cx, [$b]); } # do what's left over while ($step <= $n) { $c->_mul($cx, [$step]); $step++; if ($cx->[0] == 0) { $zero_elements ++; shift @$cx; } } } # multiply in the zeros again unshift @$cx, (0) x $zero_elements; $cx; # return result } sub _log_int { # calculate integer log of $x to base $base # ref to array, ref to array - return ref to array my ($c, $x, $base) = @_; # X == 0 => NaN return if @$x == 1 && $x->[0] == 0; # BASE 0 or 1 => NaN return if @$base == 1 && $base->[0] < 2; # X == 1 => 0 (is exact) if (@$x == 1 && $x->[0] == 1) { @$x = 0; return $x, 1; } my $cmp = $c->_acmp($x, $base); # X == BASE => 1 (is exact) if ($cmp == 0) { @$x = 1; return $x, 1; } # 1 < X < BASE => 0 (is truncated) if ($cmp < 0) { @$x = 0; return $x, 0; } my $x_org = $c->_copy($x); # preserve x # Compute a guess for the result based on: # $guess = int ( length_in_base_10(X) / ( log(base) / log(10) ) ) my $len = $c->_len($x_org); my $log = log($base->[-1]) / log(10); # for each additional element in $base, we add $BASE_LEN to the result, # based on the observation that log($BASE, 10) is BASE_LEN and # log(xy) == log(x) + log(y): $log += (@$base - 1) $BASE_LEN; # calculate now a guess based on the values obtained above: my $res = int($len / $log); @$x = $res; my $trial = $c->_pow($c->_copy($base), $x); my $acmp = $c->_acmp($trial, $x_org); # Did we get the exact result? return $x, 1 if $acmp == 0; # Too small? while ($acmp < 0) { $c->_mul($trial, $base); $c->_inc($x); $acmp = $c->_acmp($trial, $x_org); } # Too big? while ($acmp > 0) { $c->_div($trial, $base); $c->_dec($x); $acmp = $c->_acmp($trial, $x_org); } return $x, 1 if $acmp == 0; # result is exact return $x, 0; # result is too small } # for debugging: use constant DEBUG => 0; my $steps = 0; sub steps { $steps }; sub _sqrt { # square-root of $x in place # Compute a guess of the result (by rule of thumb), then improve it via # Newton's method. my ($c, $x) = @_; if (@$x == 1) { # fits into one Perl scalar, so result can be computed directly $x->[0] = int(sqrt($x->[0])); return $x; } my $y = $c->_copy($x); # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess # since our guess will "grow" my $l = int(($c->_len($x)-1) / 2); my $lastelem = $x->[-1]; # for guess my $elems = @$x - 1; # not enough digits, but could have more? if ((length($lastelem) <= 3) && ($elems > 1)) { # right-align with zero pad my $len = length($lastelem) & 1; print "$lastelem => " if DEBUG; $lastelem .= substr($x->[-2] . '0' x $BASE_LEN, 0, $BASE_LEN); # former odd => make odd again, or former even to even again $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len; print "$lastelem\n" if DEBUG; } # construct $x (instead of $c->_lsft($x, $l, 10) my $r = $l % $BASE_LEN; # 10000 00000 00000 00000 ($BASE_LEN=5) $l = int($l / $BASE_LEN); print "l = $l " if DEBUG; splice @$x, $l; # keep ref($x), but modify it # we make the first part of the guess not '1000...0' but int(sqrt($lastelem)) # that gives us: # 14400 00000 => sqrt(14400) => guess first digits to be 120 # 144000 000000 => sqrt(144000) => guess 379 print "$lastelem (elems $elems) => " if DEBUG; $lastelem = $lastelem / 10 if ($elems & 1 == 1); # odd or even? my $g = sqrt($lastelem); $g =~ s/\.//; # 2.345 => 2345 $r -= 1 if $elems & 1 == 0; # 70 => 7 # padd with zeros if result is too short $x->[$l--] = int(substr($g . '0' x $r, 0, $r+1)); print "now ", $x->[-1] if DEBUG; print " would have been ", int('1' . '0' x $r), "\n" if DEBUG; # If @$x > 1, we could compute the second elem of the guess, too, to create # an even better guess. Not implemented yet. Does it improve performance? $x->[$l--] = 0 while ($l >= 0); # all other digits of guess are zero print "start x= ", $c->_str($x), "\n" if DEBUG; my $two = $c->_two(); my $last = $c->_zero(); my $lastlast = $c->_zero(); $steps = 0 if DEBUG; while ($c->_acmp($last, $x) != 0 && $c->_acmp($lastlast, $x) != 0) { $steps++ if DEBUG; $lastlast = $c->_copy($last); $last = $c->_copy($x); $c->_add($x, $c->_div($c->_copy($y), $x)); $c->_div($x, $two ); print " x= ", $c->_str($x), "\n" if DEBUG; } print "\nsteps in sqrt: $steps, " if DEBUG; $c->_dec($x) if $c->_acmp($y, $c->_mul($c->_copy($x), $x)) < 0; # overshot? print " final ", $x->[-1], "\n" if DEBUG; $x; } sub _root { # Take n'th root of $x in place. my ($c, $x, $n) = @_; # Small numbers. if (@$x == 1 && @$n == 1) { # Result can be computed directly. Adjust initial result for numerical # errors, e.g., int(1000(1/3)) is 2, not 3. my $y = int($x->[0] (1 / $n->[0])); my $yp1 = $y + 1; $y = $yp1 if $yp1 ** $n->[0] == $x->[0]; $x->[0] = $y; return $x; } # If x <= n, the result is always (truncated to) 1. if ((@$x > 1 \|\| $x -> [0] > 0) && # if x is non-zero ... $c -> _acmp($x, $n) <= 0) # ... and x <= n { my $one = $x -> _one(); @$x = @$one; return $x; } # If $n is a power of two, take sqrt($x) repeatedly, e.g., root($x, 4) = # sqrt(sqrt($x)), root($x, 8) = sqrt(sqrt(sqrt($x))). my $b = $c -> _as_bin($n); if ($b =~ /0b1(0+)$/) { my $count = length($1); # 0b100 => len('00') => 2 my $cnt = $count; # counter for loop unshift @$x, 0; # add one element, together with one # more below in the loop this makes 2 while ($cnt-- > 0) { # 'Inflate' $x by adding one element, basically computing # $x * $BASE * $BASE. This gives us more $BASE_LEN digits for # result since len(sqrt($X)) approx == len($x) / 2. unshift @$x, 0; # Calculate sqrt($x), $x is now one element to big, again. In the # next round we make that two, again. $c -> _sqrt($x); } # $x is now one element too big, so truncate result by removing it. shift @$x; return $x; } my $DEBUG = 0; # Now the general case. This works by finding an initial guess. If this # guess is incorrect, a relatively small delta is chosen. This delta is # used to find a lower and upper limit for the correct value. The delta is # doubled in each iteration. When a lower and upper limit is found, # bisection is applied to narrow down the region until we have the correct # value. # Split x into mantissa and exponent in base 10, so that # # x = xm * 10^xe, where 0 < xm < 1 and xe is an integer my $x_str = $c -> _str($x); my $xm = "." . $x_str; my $xe = length($x_str); # From this we compute the base 10 logarithm of x # # log_10(x) = log_10(xm) + log_10(xe^10) # = log(xm)/log(10) + xe # # and then the base 10 logarithm of y, where y = x^(1/n) # # log_10(y) = log_10(x)/n my $log10x = log($xm) / log(10) + $xe; my $log10y = $log10x / $c -> _num($n); # And from this we compute ym and ye, the mantissa and exponent (in # base 10) of y, where 1 < ym <= 10 and ye is an integer. my $ye = int $log10y; my $ym = 10 ** ($log10y - $ye); # Finally, we scale the mantissa and exponent to incraese the integer # part of ym, before building the string representing our guess of y. if ($DEBUG) { print "\n"; print "xm = $xm\n"; print "xe = $xe\n"; print "log10x = $log10x\n"; print "log10y = $log10y\n"; print "ym = $ym\n"; print "ye = $ye\n"; print "\n"; } my $d = $ye < 15 ? $ye : 15; $ym = 10 * $d; $ye -= $d; my $y_str = sprintf('%.0f', $ym) . "0" x $ye; my $y = $c -> _new($y_str); if ($DEBUG) { print "ym = $ym\n"; print "ye = $ye\n"; print "\n"; print "y_str = $y_str (initial guess)\n"; print "\n"; } # See if our guess y is correct. my $trial = $c -> _pow($c -> _copy($y), $n); my $acmp = $c -> _acmp($trial, $x); if ($acmp == 0) { @$x = @$y; return $x; } # Find a lower and upper limit for the correct value of y. Start off with a # delta value that is approximately the size of the accuracy of the guess. my $lower; my $upper; my $delta = $c -> _new("1" . ("0" x $ye)); my $two = $c -> _two(); if ($acmp < 0) { $lower = $y; while ($acmp < 0) { $upper = $c -> _add($c -> _copy($lower), $delta); if ($DEBUG) { print "lower = $lower\n"; print "upper = $upper\n"; print "delta = $delta\n"; print "\n"; } $acmp = $c -> _acmp($c -> _pow($c -> _copy($upper), $n), $x); if ($acmp == 0) { @$x = @$upper; return $x; } $delta = $c -> _mul($delta, $two); } } elsif ($acmp > 0) { $upper = $y; my $zero = $c -> _zero(); while ($acmp > 0) { if ($c -> _acmp($upper, $delta) <= 0) { $lower = $c -> _zero(); last; } $lower = $c -> _sub($c -> _copy($upper), $delta); if ($DEBUG) { print "lower = $lower\n"; print "upper = $upper\n"; print "delta = $delta\n"; print "\n"; } $acmp = $c -> _acmp($c -> _pow($c -> _copy($lower), $n), $x); if ($acmp == 0) { @$x = @$lower; return $x; } $delta = $c -> _mul($delta, $two); } } # Use bisection to narrow down the interval. my $one = $c -> _one(); { $delta = $c -> _sub($c -> _copy($upper), $lower); if ($c -> _acmp($delta, $one) <= 0) { @$x = @$lower; return $x; } if ($DEBUG) { print "lower = $lower\n"; print "upper = $upper\n"; print "delta = $delta\n"; print "\n"; } $delta = $c -> _div($delta, $two); my $middle = $c -> _add($c -> _copy($lower), $delta); $acmp = $c -> _acmp($c -> _pow($c -> _copy($middle), $n), $x); if ($acmp < 0) { $lower = $middle; } elsif ($acmp > 0) { $upper = $middle; } else { @$x = @$middle; return $x; } redo; } $x; } ############################################################################## # binary stuff sub _and { my ($c, $x, $y) = @_; # the shortcut makes equal, large numbers _really_ fast, and makes only a # very small performance drop for small numbers (e.g. something with less # than 32 bit) Since we optimize for large numbers, this is enabled. return $x if $c->_acmp($x, $y) == 0; # shortcut my $m = $c->_one(); my ($xr, $yr); my $mask = $AND_MASK; my $x1 = $c->_copy($x); my $y1 = $c->_copy($y); my $z = $c->_zero(); use integer; until ($c->_is_zero($x1) \|\| $c->_is_zero($y1)) { ($x1, $xr) = $c->_div($x1, $mask); ($y1, $yr) = $c->_div($y1, $mask); $c->_add($z, $c->_mul([ 0 + $xr->[0] & 0 + $yr->[0] ], $m)); $c->_mul($m, $mask); } @$x = @$z; return $x; } sub _xor { my ($c, $x, $y) = @_; return $c->_zero() if $c->_acmp($x, $y) == 0; # shortcut (see -and) my $m = $c->_one(); my ($xr, $yr); my $mask = $XOR_MASK; my $x1 = $c->_copy($x); my $y1 = $c->_copy($y); # make copy my $z = $c->_zero(); use integer; until ($c->_is_zero($x1) \|\| $c->_is_zero($y1)) { ($x1, $xr) = $c->_div($x1, $mask); ($y1, $yr) = $c->_div($y1, $mask); # make ints() from $xr, $yr (see _and()) #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b = $BASE; } #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ $b; $b = $BASE; } #$c->_add($x, $c->_mul($c->_new($xrr ^ $yrr)), $m) ); $c->_add($z, $c->_mul([ 0 + $xr->[0] ^ 0 + $yr->[0] ], $m)); $c->_mul($m, $mask); } # the loop stops when the shorter of the two numbers is exhausted # the remainder of the longer one will survive bit-by-bit, so we simple # multiply-add it in $c->_add($z, $c->_mul($x1, $m) ) if !$c->_is_zero($x1); $c->_add($z, $c->_mul($y1, $m) ) if !$c->_is_zero($y1); @$x = @$z; return $x; } sub _or { my ($c, $x, $y) = @_; return $x if $c->_acmp($x, $y) == 0; # shortcut (see _and) my $m = $c->_one(); my ($xr, $yr); my $mask = $OR_MASK; my $x1 = $c->_copy($x); my $y1 = $c->_copy($y); # make copy my $z = $c->_zero(); use integer; until ($c->_is_zero($x1) \|\| $c->_is_zero($y1)) { ($x1, $xr) = $c->_div($x1, $mask); ($y1, $yr) = $c->_div($y1, $mask); # make ints() from $xr, $yr (see _and()) # $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ $b; $b = $BASE; } # $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ $b; $b = $BASE; } # $c->_add($x, $c->_mul(_new( $c, ($xrr \| $yrr) ), $m) ); $c->_add($z, $c->_mul([ 0 + $xr->[0] \| 0 + $yr->[0] ], $m)); $c->_mul($m, $mask); } # the loop stops when the shorter of the two numbers is exhausted # the remainder of the longer one will survive bit-by-bit, so we simple # multiply-add it in $c->_add($z, $c->_mul($x1, $m) ) if !$c->_is_zero($x1); $c->_add($z, $c->_mul($y1, $m) ) if !$c->_is_zero($y1); @$x = @$z; return $x; } sub _as_hex { # convert a decimal number to hex (ref to array, return ref to string) my ($c, $x) = @_; # fits into one element (handle also 0x0 case) return sprintf("0x%x", $x->[0]) if @$x == 1; my $x1 = $c->_copy($x); my $es = ''; my ($xr, $h, $x10000); if ($] >= 5.006) { $x10000 = [ 0x10000 ]; $h = 'h4'; } else { $x10000 = [ 0x1000 ]; $h = 'h3'; } while (@$x1 != 1 \|\| $x1->[0] != 0) # _is_zero() { ($x1, $xr) = $c->_div($x1, $x10000); $es .= unpack($h, pack('V', $xr->[0])); } $es = reverse $es; $es =~ s/^[0]+//; # strip leading zeros '0x' . $es; # return result prepended with 0x } sub _as_bin { # convert a decimal number to bin (ref to array, return ref to string) my ($c, $x) = @_; # fits into one element (and Perl recent enough), handle also 0b0 case # handle zero case for older Perls if ($] <= 5.005 && @$x == 1 && $x->[0] == 0) { my $t = '0b0'; return $t; } if (@$x == 1 && $] >= 5.006) { my $t = sprintf("0b%b", $x->[0]); return $t; } my $x1 = $c->_copy($x); my $es = ''; my ($xr, $b, $x10000); if ($] >= 5.006) { $x10000 = [ 0x10000 ]; $b = 'b16'; } else { $x10000 = [ 0x1000 ]; $b = 'b12'; } while (!(@$x1 == 1 && $x1->[0] == 0)) # _is_zero() { ($x1, $xr) = $c->_div($x1, $x10000); $es .= unpack($b, pack('v', $xr->[0])); } $es = reverse $es; $es =~ s/^[0]+//; # strip leading zeros '0b' . $es; # return result prepended with 0b } sub _as_oct { # convert a decimal number to octal (ref to array, return ref to string) my ($c, $x) = @_; # fits into one element (handle also 0 case) return sprintf("0%o", $x->[0]) if @$x == 1; my $x1 = $c->_copy($x); my $es = ''; my $xr; my $x1000 = [ 0100000 ]; while (@$x1 != 1 \|\| $x1->[0] != 0) # _is_zero() { ($x1, $xr) = $c->_div($x1, $x1000); $es .= reverse sprintf("%05o", $xr->[0]); } $es = reverse $es; $es =~ s/^0+//; # strip leading zeros '0' . $es; # return result prepended with 0 } sub _from_oct { # convert a octal number to decimal (string, return ref to array) my ($c, $os) = @_; # for older Perls, play safe my $m = [ 0100000 ]; my $d = 5; # 5 digits at a time my $mul = $c->_one(); my $x = $c->_zero(); my $len = int((length($os) - 1) / $d); # $d digit parts, w/o the '0' my $val; my $i = -$d; while ($len >= 0) { $val = substr($os, $i, $d); # get oct digits $val = CORE::oct($val); $i -= $d; $len --; my $adder = [ $val ]; $c->_add($x, $c->_mul($adder, $mul)) if $val != 0; $c->_mul($mul, $m) if $len >= 0; # skip last mul } $x; } sub _from_hex { # convert a hex number to decimal (string, return ref to array) my ($c, $hs) = @_; my $m = $c->_new(0x10000000); # 28 bit at a time (<32 bit!) my $d = 7; # 7 digits at a time my $mul = $c->_one(); my $x = $c->_zero(); my $len = int((length($hs) - 2) / $d); # $d digit parts, w/o the '0x' my $val; my $i = -$d; while ($len >= 0) { $val = substr($hs, $i, $d); # get hex digits $val =~ s/^0x// if $len == 0; # for last part only because $val = CORE::hex($val); # hex does not like wrong chars $i -= $d; $len --; my $adder = [ $val ]; # if the resulting number was to big to fit into one element, create a # two-element version (bug found by Mark Lakata - Thanx!) if (CORE::length($val) > $BASE_LEN) { $adder = $c->_new($val); } $c->_add($x, $c->_mul($adder, $mul)) if $val != 0; $c->_mul($mul, $m) if $len >= 0; # skip last mul } $x; } sub _from_bin { # convert a hex number to decimal (string, return ref to array) my ($c, $bs) = @_; # instead of converting X (8) bit at a time, it is faster to "convert" the # number to hex, and then call _from_hex. my $hs = $bs; $hs =~ s/^[+-]?0b//; # remove sign and 0b my $l = length($hs); # bits $hs = '0' x (8 - ($l % 8)) . $hs if ($l % 8) != 0; # padd left side w/ 0 my $h = '0x' . unpack('H', pack ('B', $hs)); # repack as hex $c->_from_hex($h); } ############################################################################## # special modulus functions sub _modinv { # modular multiplicative inverse my ($c, $x, $y) = @_; # modulo zero if ($c->_is_zero($y)) { return undef, undef; } # modulo one if ($c->_is_one($y)) { return $c->_zero(), '+'; } my $u = $c->_zero(); my $v = $c->_one(); my $a = $c->_copy($y); my $b = $c->_copy($x); # Euclid's Algorithm for bgcd(), only that we calc bgcd() ($a) and the result # ($u) at the same time. See comments in BigInt for why this works. my $q; my $sign = 1; { ($a, $q, $b) = ($b, $c->_div($a, $b)); # step 1 last if $c->_is_zero($b); my $t = $c->_add( # step 2: $c->_mul($c->_copy($v), $q), # t = v q $u); # + u $u = $v; # u = v $v = $t; # v = t $sign = -$sign; redo; } # if the gcd is not 1, then return NaN return (undef, undef) unless $c->_is_one($a); ($v, $sign == 1 ? '+' : '-'); } sub _modpow { # modulus of power ($x $y) % $z my ($c, $num, $exp, $mod) = @_; # a^b (mod 1) = 0 for all a and b if ($c->_is_one($mod)) { @$num = 0; return $num; } # 0^a (mod m) = 0 if m != 0, a != 0 # 0^0 (mod m) = 1 if m != 0 if ($c->_is_zero($num)) { if ($c->_is_zero($exp)) { @$num = 1; } else { @$num = 0; } return $num; } # $num = $c->_mod($num, $mod); # this does not make it faster my $acc = $c->_copy($num); my $t = $c->_one(); my $expbin = $c->_as_bin($exp); $expbin =~ s/^0b//; my $len = length($expbin); while (--$len >= 0) { if (substr($expbin, $len, 1) eq '1') { # is_odd $t = $c->_mul($t, $acc); $t = $c->_mod($t, $mod); } $acc = $c->_mul($acc, $acc); $acc = $c->_mod($acc, $mod); } @$num = @$t; $num; } sub _gcd { # Greatest common divisor. my ($c, $x, $y) = @_; # gcd(0, 0) = 0 # gcd(0, a) = a, if a != 0 if (@$x == 1 && $x->[0] == 0) { if (@$y == 1 && $y->[0] == 0) { @$x = 0; } else { @$x = @$y; } return $x; } # Until $y is zero ... until (@$y == 1 && $y->[0] == 0) { # Compute remainder. $c->_mod($x, $y); # Swap $x and $y. my $tmp = $c->_copy($x); @$x = @$y; $y = $tmp; # no deref here; that would modify input $y } return $x; } 1; =pod =head1 NAME Math::BigInt::Calc - Pure Perl module to support Math::BigInt =head1 SYNOPSIS # to use it with Math::BigInt use Math::BigInt lib => 'Calc'; # to use it with Math::BigFloat use Math::BigFloat lib => 'Calc'; # to use it with Math::BigRat use Math::BigRat lib => 'Calc'; =head1 DESCRIPTION Math::BigInt::Calc inherits from Math::BigInt::Lib. In this library, the numbers are represented in base B = 10N, where N is the largest possible value that does not cause overflow in the intermediate computations. The base B elements are stored in an array, with the least significant element stored in array element zero. There are no leading zero elements, except a single zero element when the number is zero. For instance, if B = 10000, the number 1234567890 is represented internally as [7890, 3456, 12]. =head1 SEE ALSO L<Math::BigInt::Lib> for a description of the API. Alternative libraries L<Math::BigInt::FastCalc>, L<Math::BigInt::GMP>, and L<Math::BigInt::Pari>. Some of the modules that use these libraries L<Math::BigInt>, L<Math::BigFloat>, and L<Math::BigRat>. =cut Back

Edit File : Calc.pm | Size : 75.06 KB

package Math::BigInt::Calc;

use 5.006001;
use strict;
use warnings;

use Carp;
use Math::BigInt::Lib;

our $VERSION = '1.999811';

our @ISA = ('Math::BigInt::Lib');

# Package to store unsigned big integers in decimal and do math with them

# Internally the numbers are stored in an array with at least 1 element, no
# leading zero parts (except the first) and in base 1eX where X is determined
# automatically at loading time to be the maximum possible value

# todo:
# - fully remove funky $# stuff in div() (maybe - that code scares me...)

# USE_MUL: due to problems on certain os (os390, posix-bc) "* 1e-5" is used
# instead of "/ 1e5" at some places, (marked with USE_MUL). Other platforms
# BS2000, some Crays need USE_DIV instead.
# The BEGIN block is used to determine which of the two variants gives the
# correct result.

# Beware of things like:
# $i = $i * $y + $car; $car = int($i / $BASE); $i = $i % $BASE;
# This works on x86, but fails on ARM (SA1100, iPAQ) due to who knows what
# reasons. So, use this instead (slower, but correct):
# $i = $i * $y + $car; $car = int($i / $BASE); $i -= $BASE * $car;

##############################################################################
# global constants, flags and accessory

# announce that we are compatible with MBI v1.83 and up
sub api_version () { 2; }

# constants for easier life
my ($BASE, $BASE_LEN, $RBASE, $MAX_VAL);
my ($AND_BITS, $XOR_BITS, $OR_BITS);
my ($AND_MASK, $XOR_MASK, $OR_MASK);

sub _base_len {
    # Set/get the BASE_LEN and assorted other, related values.
    # Used only by the testsuite, the set variant is used only by the BEGIN
    # block below:

my ($class, $b, $int) = @_;
    if (defined $b) {
        # avoid redefinitions
        undef &_mul;
        undef &_div;

if ($] >= 5.008 && $int && $b > 7) {
            $BASE_LEN = $b;
            *_mul = \&_mul_use_div_64;
            *_div = \&_div_use_div_64;
            $BASE = int("1e" . $BASE_LEN);
            $MAX_VAL = $BASE-1;
            return $BASE_LEN unless wantarray;
            return ($BASE_LEN, $BASE, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN, $MAX_VAL);
        }

# find whether we can use mul or div in mul()/div()
        $BASE_LEN = $b + 1;
        my $caught = 0;
        while (--$BASE_LEN > 5) {
            $BASE = int("1e" . $BASE_LEN);
            $RBASE = abs('1e-' . $BASE_LEN); # see USE_MUL
            $caught = 0;
            $caught += 1 if (int($BASE * $RBASE) != 1); # should be 1
            $caught += 2 if (int($BASE / $BASE) != 1);  # should be 1
            last if $caught != 3;
        }
        $BASE = int("1e" . $BASE_LEN);
        $RBASE = abs('1e-' . $BASE_LEN); # see USE_MUL
        $MAX_VAL = $BASE-1;

# ($caught & 1) != 0 => cannot use MUL
        # ($caught & 2) != 0 => cannot use DIV
        if ($caught == 2)       # 2
        {
            # must USE_MUL since we cannot use DIV
            *_mul = \&_mul_use_mul;
            *_div = \&_div_use_mul;
        } else                  # 0 or 1
        {
            # can USE_DIV instead
            *_mul = \&_mul_use_div;
            *_div = \&_div_use_div;
        }
    }
    return $BASE_LEN unless wantarray;
    return ($BASE_LEN, $BASE, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN, $MAX_VAL);
}

sub _new {
    # Given a string representing an integer, returns a reference to an array
    # of integers, where each integer represents a chunk of the original input
    # integer.

my ($class, $str) = @_;
    #unless ($str =~ /^([1-9]\d*|0)\z/) {
    #    require Carp;
    #    Carp::croak("Invalid input string '$str'");
    #}

my $input_len = length($str) - 1;

# Shortcut for small numbers.
    return bless [ $str ], $class if $input_len < $BASE_LEN;

my $format = "a" . (($input_len % $BASE_LEN) + 1);
    $format .= $] < 5.008 ? "a$BASE_LEN" x int($input_len / $BASE_LEN)
                          : "(a$BASE_LEN)*";

my $self = [ reverse(map { 0 + $_ } unpack($format, $str)) ];
    return bless $self, $class;
}

BEGIN {
    # from Daniel Pfeiffer: determine largest group of digits that is precisely
    # multipliable with itself plus carry
    # Test now changed to expect the proper pattern, not a result off by 1 or 2
    my ($e, $num) = 3;          # lowest value we will use is 3+1-1 = 3
    do {
        $num = '9' x ++$e;
        $num *= $num + 1;
    } while $num =~ /9{$e}0{$e}/; # must be a certain pattern
    $e--;                         # last test failed, so retract one step
    # the limits below brush the problems with the test above under the rug:
    # the test should be able to find the proper $e automatically
    $e = 5 if $^O =~ /^uts/;    # UTS get's some special treatment
    $e = 5 if $^O =~ /^unicos/; # unicos is also problematic (6 seems to work
                                # there, but we play safe)

my $int = 0;
    if ($e > 7) {
        use integer;
        my $e1 = 7;
        $num = 7;
        do {
            $num = ('9' x ++$e1) + 0;
            $num *= $num + 1;
        } while ("$num" =~ /9{$e1}0{$e1}/); # must be a certain pattern
        $e1--;                  # last test failed, so retract one step
        if ($e1 > 7) {
            $int = 1;
            $e = $e1;
        }
    }

__PACKAGE__ -> _base_len($e, $int);        # set and store

use integer;
    # find out how many bits _and, _or and _xor can take (old default = 16)
    # I don't think anybody has yet 128 bit scalars, so let's play safe.
    local $^W = 0;              # don't warn about 'nonportable number'
    $AND_BITS = 15;
    $XOR_BITS = 15;
    $OR_BITS  = 15;

# find max bits, we will not go higher than numberofbits that fit into $BASE
    # to make _and etc simpler (and faster for smaller, slower for large numbers)
    my $max = 16;
    while (2 ** $max < $BASE) {
        $max++;
    }
    {
        no integer;
        $max = 16 if $] < 5.006; # older Perls might not take >16 too well
    }
    my ($x, $y, $z);

do {
        $AND_BITS++;
        $x = CORE::oct('0b' . '1' x $AND_BITS);
        $y = $x & $x;
        $z = (2 ** $AND_BITS) - 1;
    } while ($AND_BITS < $max && $x == $z && $y == $x);
    $AND_BITS --;               # retreat one step

do {
        $XOR_BITS++;
        $x = CORE::oct('0b' . '1' x $XOR_BITS);
        $y = $x ^ 0;
        $z = (2 ** $XOR_BITS) - 1;
    } while ($XOR_BITS < $max && $x == $z && $y == $x);
    $XOR_BITS --;               # retreat one step

do {
        $OR_BITS++;
        $x = CORE::oct('0b' . '1' x $OR_BITS);
        $y = $x | $x;
        $z = (2 ** $OR_BITS) -  1;
    } while ($OR_BITS < $max && $x == $z && $y == $x);
    $OR_BITS--;                # retreat one step

$AND_MASK = __PACKAGE__->_new(( 2 ** $AND_BITS ));
    $XOR_MASK = __PACKAGE__->_new(( 2 ** $XOR_BITS ));
    $OR_MASK  = __PACKAGE__->_new(( 2 ** $OR_BITS  ));

# We can compute the approximate length no faster than the real length:
    *_alen = \&_len;
}

###############################################################################

sub _zero {
    # create a zero
    my $class = shift;
    return bless [ 0 ], $class;
}

sub _one {
    # create a one
    my $class = shift;
    return bless [ 1 ], $class;
}

sub _two {
    # create a two
    my $class = shift;
    return bless [ 2 ], $class;
}

sub _ten {
    # create a 10
    my $class = shift;
    bless [ 10 ], $class;
}

sub _1ex {
    # create a 1Ex
    my $class = shift;

my $rem   = $_[0] % $BASE_LEN;      # remainder
    my $parts = $_[0] / $BASE_LEN;      # parts

# 000000, 000000, 100
    bless [ (0) x $parts, '1' . ('0' x $rem) ], $class;
}

sub _copy {
    # make a true copy
    my $class = shift;
    return bless [ @{ $_[0] } ], $class;
}

# catch and throw away
sub import { }

##############################################################################
# convert back to string and number

sub _str {
    # Convert number from internal base 1eN format to string format. Internal
    # format is always normalized, i.e., no leading zeros.

my $ary = $_[1];
    my $idx = $#$ary;           # index of last element

if ($idx < 0) {             # should not happen
        require Carp;
        Carp::croak("$_[1] has no elements");
    }

# Handle first one differently, since it should not have any leading zeros.
    my $ret = int($ary->[$idx]);
    if ($idx > 0) {
        # Interestingly, the pre-padd method uses more time.
        # The old grep variant takes longer (14 vs. 10 sec).
        my $z = '0' x ($BASE_LEN - 1);
        while (--$idx >= 0) {
            $ret .= substr($z . $ary->[$idx], -$BASE_LEN);
        }
    }
    $ret;
}

sub _num {
    # Make a Perl scalar number (int/float) from a BigInt object.
    my $x = $_[1];

return $x->[0] if @$x == 1;         # below $BASE

# Start with the most significant element and work towards the least
    # significant element. Avoid multiplying "inf" (which happens if the number
    # overflows) with "0" (if there are zero elements in $x) since this gives
    # "nan" which propagates to the output.

my $num = 0;
    for (my $i = $#$x ; $i >= 0 ; --$i) {
        $num *= $BASE;
        $num += $x -> [$i];
    }
    return $num;
}

##############################################################################
# actual math code

sub _add {
    # (ref to int_num_array, ref to int_num_array)
    #
    # Routine to add two base 1eX numbers stolen from Knuth Vol 2 Algorithm A
    # pg 231. There are separate routines to add and sub as per Knuth pg 233.
    # This routine modifies array x, but not y.

my ($c, $x, $y) = @_;

# $x + 0 => $x

return $x if @$y == 1 && $y->[0] == 0;

# 0 + $y => $y->copy

if (@$x == 1 && $x->[0] == 0) {
        @$x = @$y;
        return $x;
    }

# For each in Y, add Y to X and carry. If after that, something is left in
    # X, foreach in X add carry to X and then return X, carry. Trades one
    # "$j++" for having to shift arrays.
    my $i;
    my $car = 0;
    my $j = 0;
    for $i (@$y) {
        $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
        $j++;
    }
    while ($car != 0) {
        $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0;
        $j++;
    }
    $x;
}

sub _inc {
    # (ref to int_num_array, ref to int_num_array)
    # Add 1 to $x, modify $x in place
    my ($c, $x) = @_;

for my $i (@$x) {
        return $x if ($i += 1) < $BASE; # early out
        $i = 0;                         # overflow, next
    }
    push @$x, 1 if $x->[-1] == 0;       # last overflowed, so extend
    $x;
}

sub _dec {
    # (ref to int_num_array, ref to int_num_array)
    # Sub 1 from $x, modify $x in place
    my ($c, $x) = @_;

my $MAX = $BASE - 1;                # since MAX_VAL based on BASE
    for my $i (@$x) {
        last if ($i -= 1) >= 0;         # early out
        $i = $MAX;                      # underflow, next
    }
    pop @$x if $x->[-1] == 0 && @$x > 1; # last underflowed (but leave 0)
    $x;
}

sub _sub {
    # (ref to int_num_array, ref to int_num_array, swap)
    #
    # Subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
    # subtract Y from X by modifying x in place
    my ($c, $sx, $sy, $s) = @_;

my $car = 0;
    my $i;
    my $j = 0;
    if (!$s) {
        for $i (@$sx) {
            last unless defined $sy->[$j] || $car;
            $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0);
            $j++;
        }
        # might leave leading zeros, so fix that
        return __strip_zeros($sx);
    }
    for $i (@$sx) {
        # We can't do an early out if $x < $y, since we need to copy the high
        # chunks from $y. Found by Bob Mathews.
        #last unless defined $sy->[$j] || $car;
        $sy->[$j] += $BASE
          if $car = ($sy->[$j] = $i - ($sy->[$j] || 0) - $car) < 0;
        $j++;
    }
    # might leave leading zeros, so fix that
    __strip_zeros($sy);
}

sub _mul_use_mul {
    # (ref to int_num_array, ref to int_num_array)
    # multiply two numbers in internal representation
    # modifies first arg, second need not be different from first
    my ($c, $xv, $yv) = @_;

if (@$yv == 1) {
        # shortcut for two very short numbers (improved by Nathan Zook)
        # works also if xv and yv are the same reference, and handles also $x == 0
        if (@$xv == 1) {
            if (($xv->[0] *= $yv->[0]) >= $BASE) {
                $xv->[0] = $xv->[0] - ($xv->[1] = int($xv->[0] * $RBASE)) * $BASE;
            }
            ;
            return $xv;
        }
        # $x * 0 => 0
        if ($yv->[0] == 0) {
            @$xv = (0);
            return $xv;
        }
        # multiply a large number a by a single element one, so speed up
        my $y = $yv->[0];
        my $car = 0;
        foreach my $i (@$xv) {
            $i = $i * $y + $car;
            $car = int($i * $RBASE);
            $i -= $car * $BASE;
        }
        push @$xv, $car if $car != 0;
        return $xv;
    }
    # shortcut for result $x == 0 => result = 0
    return $xv if @$xv == 1 && $xv->[0] == 0;

# since multiplying $x with $x fails, make copy in this case
    $yv = [ @$xv ] if $xv == $yv; # same references?

my @prod = ();
    my ($prod, $car, $cty, $xi, $yi);

for $xi (@$xv) {
        $car = 0;
        $cty = 0;

# slow variant
        #    for $yi (@$yv)
        #      {
        #      $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
        #      $prod[$cty++] =
        #       $prod - ($car = int($prod * RBASE)) * $BASE;  # see USE_MUL
        #      }
        #    $prod[$cty] += $car if $car; # need really to check for 0?
        #    $xi = shift @prod;

# faster variant
        # looping through this if $xi == 0 is silly - so optimize it away!
        $xi = (shift @prod || 0), next if $xi == 0;
        for $yi (@$yv) {
            $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
            ##     this is actually a tad slower
            ##        $prod = $prod[$cty]; $prod += ($car + $xi * $yi);     # no ||0 here
            $prod[$cty++] =
              $prod - ($car = int($prod * $RBASE)) * $BASE; # see USE_MUL
        }
        $prod[$cty] += $car if $car; # need really to check for 0?
        $xi = shift @prod || 0;      # || 0 makes v5.005_3 happy
    }
    push @$xv, @prod;
    # can't have leading zeros
    #  __strip_zeros($xv);
    $xv;
}

sub _mul_use_div_64 {
    # (ref to int_num_array, ref to int_num_array)
    # multiply two numbers in internal representation
    # modifies first arg, second need not be different from first
    # works for 64 bit integer with "use integer"
    my ($c, $xv, $yv) = @_;

use integer;
    if (@$yv == 1) {
        # shortcut for two small numbers, also handles $x == 0
        if (@$xv == 1) {
            # shortcut for two very short numbers (improved by Nathan Zook)
            # works also if xv and yv are the same reference, and handles also $x == 0
            if (($xv->[0] *= $yv->[0]) >= $BASE) {
                $xv->[0] =
                  $xv->[0] - ($xv->[1] = $xv->[0] / $BASE) * $BASE;
            }
            return $xv;
        }
        # $x * 0 => 0
        if ($yv->[0] == 0) {
            @$xv = (0);
            return $xv;
        }
        # multiply a large number a by a single element one, so speed up
        my $y = $yv->[0];
        my $car = 0;
        foreach my $i (@$xv) {
            #$i = $i * $y + $car; $car = $i / $BASE; $i -= $car * $BASE;
            $i = $i * $y + $car;
            $i -= ($car = $i / $BASE) * $BASE;
        }
        push @$xv, $car if $car != 0;
        return $xv;
    }
    # shortcut for result $x == 0 => result = 0
    return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );

# since multiplying $x with $x fails, make copy in this case
    $yv = $c->_copy($xv) if $xv == $yv; # same references?

my @prod = ();
    my ($prod, $car, $cty, $xi, $yi);
    for $xi (@$xv) {
        $car = 0;
        $cty = 0;
        # looping through this if $xi == 0 is silly - so optimize it away!
        $xi = (shift @prod || 0), next if $xi == 0;
        for $yi (@$yv) {
            $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
            $prod[$cty++] = $prod - ($car = $prod / $BASE) * $BASE;
        }
        $prod[$cty] += $car if $car; # need really to check for 0?
        $xi = shift @prod || 0;      # || 0 makes v5.005_3 happy
    }
    push @$xv, @prod;
    $xv;
}

sub _mul_use_div {
    # (ref to int_num_array, ref to int_num_array)
    # multiply two numbers in internal representation
    # modifies first arg, second need not be different from first
    my ($c, $xv, $yv) = @_;

if (@$yv == 1) {
        # shortcut for two small numbers, also handles $x == 0
        if (@$xv == 1) {
            # shortcut for two very short numbers (improved by Nathan Zook)
            # works also if xv and yv are the same reference, and handles also $x == 0
            if (($xv->[0] *= $yv->[0]) >= $BASE) {
                $xv->[0] =
                  $xv->[0] - ($xv->[1] = int($xv->[0] / $BASE)) * $BASE;
            }
            ;
            return $xv;
        }
        # $x * 0 => 0
        if ($yv->[0] == 0) {
            @$xv = (0);
            return $xv;
        }
        # multiply a large number a by a single element one, so speed up
        my $y = $yv->[0];
        my $car = 0;
        foreach my $i (@$xv) {
            $i = $i * $y + $car;
            $car = int($i / $BASE);
            $i -= $car * $BASE;
            # This (together with use integer;) does not work on 32-bit Perls
            #$i = $i * $y + $car; $i -= ($car = $i / $BASE) * $BASE;
        }
        push @$xv, $car if $car != 0;
        return $xv;
    }
    # shortcut for result $x == 0 => result = 0
    return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );

# since multiplying $x with $x fails, make copy in this case
    $yv = $c->_copy($xv) if $xv == $yv; # same references?

my @prod = ();
    my ($prod, $car, $cty, $xi, $yi);
    for $xi (@$xv) {
        $car = 0;
        $cty = 0;
        # looping through this if $xi == 0 is silly - so optimize it away!
        $xi = (shift @prod || 0), next if $xi == 0;
        for $yi (@$yv) {
            $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
            $prod[$cty++] = $prod - ($car = int($prod / $BASE)) * $BASE;
        }
        $prod[$cty] += $car if $car; # need really to check for 0?
        $xi = shift @prod || 0;      # || 0 makes v5.005_3 happy
    }
    push @$xv, @prod;
    # can't have leading zeros
    #  __strip_zeros($xv);
    $xv;
}

sub _div_use_mul {
    # ref to array, ref to array, modify first array and return remainder if
    # in list context

# see comments in _div_use_div() for more explanations

my ($c, $x, $yorg) = @_;

# the general div algorithm here is about O(N*N) and thus quite slow, so
    # we first check for some special cases and use shortcuts to handle them.

# This works, because we store the numbers in a chunked format where each
    # element contains 5..7 digits (depending on system).

# if both numbers have only one element:
    if (@$x == 1 && @$yorg == 1) {
        # shortcut, $yorg and $x are two small numbers
        if (wantarray) {
            my $rem = [ $x->[0] % $yorg->[0] ];
            bless $rem, $c;
            $x->[0] = int($x->[0] / $yorg->[0]);
            return ($x, $rem);
        } else {
            $x->[0] = int($x->[0] / $yorg->[0]);
            return $x;
        }
    }

# if x has more than one, but y has only one element:
    if (@$yorg == 1) {
        my $rem;
        $rem = $c->_mod($c->_copy($x), $yorg) if wantarray;

# shortcut, $y is < $BASE
        my $j = @$x;
        my $r = 0;
        my $y = $yorg->[0];
        my $b;
        while ($j-- > 0) {
            $b = $r * $BASE + $x->[$j];
            $x->[$j] = int($b/$y);
            $r = $b % $y;
        }
        pop @$x if @$x > 1 && $x->[-1] == 0; # splice up a leading zero
        return ($x, $rem) if wantarray;
        return $x;
    }

# now x and y have more than one element

# if $yorg has more digits than $x (it's leading element is longer than
        # the one from $x), the result will also be 0:
        if (length(int($yorg->[-1])) > length(int($x->[-1]))) {
            my $rem = $c->_copy($x) if wantarray;        # make copy
            @$x = 0;                            # set to 0
            return ($x, $rem) if wantarray;     # including remainder?
            return $x;
        }
        # now calculate $x / $yorg
        if (length(int($yorg->[-1])) == length(int($x->[-1]))) {
            # same length, so make full compare

my $a = 0;
            my $j = @$x - 1;
            # manual way (abort if unequal, good for early ne)
            while ($j >= 0) {
                last if ($a = $x->[$j] - $yorg->[$j]);
                $j--;
            }
            # $a contains the result of the compare between X and Y
            # a < 0: x < y, a == 0: x == y, a > 0: x > y
            if ($a <= 0) {
                # a = 0 => x == y => rem 0
                # a < 0 => x < y => rem = x
                my $rem = $a == 0 ? $c->_zero() : $c->_copy($x);
                @$x = 0;             # if $a < 0
                $x->[0] = 1 if $a == 0;  # $x == $y
                return ($x, $rem) if wantarray;
                return $x;
            }
            # $x >= $y, so proceed normally
        }
    }

# all other cases:

my $y = $c->_copy($yorg);         # always make copy to preserve

my ($car, $bar, $prd, $dd, $xi, $yi, @q, $v2, $v1, @d, $tmp, $q, $u2, $u1, $u0);

$car = $bar = $prd = 0;
    if (($dd = int($BASE / ($y->[-1] + 1))) != 1) {
        for $xi (@$x) {
            $xi = $xi * $dd + $car;
            $xi -= ($car = int($xi * $RBASE)) * $BASE; # see USE_MUL
        }
        push(@$x, $car);
        $car = 0;
        for $yi (@$y) {
            $yi = $yi * $dd + $car;
            $yi -= ($car = int($yi * $RBASE)) * $BASE; # see USE_MUL
        }
    } else {
        push(@$x, 0);
    }
    @q = ();
    ($v2, $v1) = @$y[-2, -1];
    $v2 = 0 unless $v2;
    while ($#$x > $#$y) {
        ($u2, $u1, $u0) = @$x[-3 .. -1];
        $u2 = 0 unless $u2;
        #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
        # if $v1 == 0;
        $q = (($u0 == $v1) ? $MAX_VAL : int(($u0 * $BASE + $u1) / $v1));
        --$q while ($v2 * $q > ($u0 * $BASE + $u1 - $q * $v1) * $BASE + $u2);
        if ($q) {
            ($car, $bar) = (0, 0);
            for ($yi = 0, $xi = $#$x - $#$y-1; $yi <= $#$y; ++$yi, ++$xi) {
                $prd = $q * $y->[$yi] + $car;
                $prd -= ($car = int($prd * $RBASE)) * $BASE; # see USE_MUL
                $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
            }
            if ($x->[-1] < $car + $bar) {
                $car = 0;
                --$q;
                for ($yi = 0, $xi = $#$x - $#$y-1; $yi <= $#$y; ++$yi, ++$xi) {
                    $x->[$xi] -= $BASE
                      if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
                }
            }
        }
        pop(@$x);
        unshift(@q, $q);
    }
    if (wantarray) {
        my $d = bless [], $c;
        if ($dd != 1) {
            $car = 0;
            for $xi (reverse @$x) {
                $prd = $car * $BASE + $xi;
                $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
                unshift(@$d, $tmp);
            }
        } else {
            @$d = @$x;
        }
        @$x = @q;
        __strip_zeros($x);
        __strip_zeros($d);
        return ($x, $d);
    }
    @$x = @q;
    __strip_zeros($x);
    $x;
}

sub _div_use_div_64 {
    # ref to array, ref to array, modify first array and return remainder if
    # in list context
    # This version works on 64 bit integers
    my ($c, $x, $yorg) = @_;

use integer;
    # the general div algorithm here is about O(N*N) and thus quite slow, so
    # we first check for some special cases and use shortcuts to handle them.

# This works, because we store the numbers in a chunked format where each
    # element contains 5..7 digits (depending on system).

# check whether y has more elements than x, if yet, the result will be 0
    if (@$yorg > @$x) {
        my $rem;
        $rem = $c->_copy($x) if wantarray;    # make copy
        @$x = 0;                        # set to 0
        return ($x, $rem) if wantarray; # including remainder?
        return $x;                      # only x, which is [0] now
    }
    # check whether the numbers have the same number of elements, in that case
    # the result will fit into one element and can be computed efficiently
    if (@$yorg == @$x) {
        my $rem;
        # if $yorg has more digits than $x (it's leading element is longer than
        # the one from $x), the result will also be 0:
        if (length(int($yorg->[-1])) > length(int($x->[-1]))) {
            $rem = $c->_copy($x) if wantarray;     # make copy
            @$x = 0;                          # set to 0
            return ($x, $rem) if wantarray; # including remainder?
            return $x;
        }
        # now calculate $x / $yorg

if (length(int($yorg->[-1])) == length(int($x->[-1]))) {
            # same length, so make full compare

my $a = 0;
            my $j = @$x - 1;
            # manual way (abort if unequal, good for early ne)
            while ($j >= 0) {
                last if ($a = $x->[$j] - $yorg->[$j]);
                $j--;
            }
            # $a contains the result of the compare between X and Y
            # a < 0: x < y, a == 0: x == y, a > 0: x > y
            if ($a <= 0) {
                $rem = $c->_zero();                  # a = 0 => x == y => rem 0
                $rem = $c->_copy($x) if $a != 0;       # a < 0 => x < y => rem = x
                @$x = 0;                       # if $a < 0
                $x->[0] = 1 if $a == 0;        # $x == $y
                return ($x, $rem) if wantarray; # including remainder?
                return $x;
            }
            # $x >= $y, so proceed normally
        }
    }

# all other cases:

my $y = $c->_copy($yorg);         # always make copy to preserve

my ($car, $bar, $prd, $dd, $xi, $yi, @q, $v2, $v1, @d, $tmp, $q, $u2, $u1, $u0);

$car = $bar = $prd = 0;
    if (($dd = int($BASE / ($y->[-1] + 1))) != 1) {
        for $xi (@$x) {
            $xi = $xi * $dd + $car;
            $xi -= ($car = int($xi / $BASE)) * $BASE;
        }
        push(@$x, $car);
        $car = 0;
        for $yi (@$y) {
            $yi = $yi * $dd + $car;
            $yi -= ($car = int($yi / $BASE)) * $BASE;
        }
    } else {
        push(@$x, 0);
    }

# @q will accumulate the final result, $q contains the current computed
    # part of the final result

@q = ();
    ($v2, $v1) = @$y[-2, -1];
    $v2 = 0 unless $v2;
    while ($#$x > $#$y) {
        ($u2, $u1, $u0) = @$x[-3..-1];
        $u2 = 0 unless $u2;
        #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
        # if $v1 == 0;
        $q = (($u0 == $v1) ? $MAX_VAL : int(($u0 * $BASE + $u1) / $v1));
        --$q while ($v2 * $q > ($u0 * $BASE +$ u1- $q*$v1) * $BASE + $u2);
        if ($q) {
            ($car, $bar) = (0, 0);
            for ($yi = 0, $xi = $#$x - $#$y - 1; $yi <= $#$y; ++$yi, ++$xi) {
                $prd = $q * $y->[$yi] + $car;
                $prd -= ($car = int($prd / $BASE)) * $BASE;
                $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
            }
            if ($x->[-1] < $car + $bar) {
                $car = 0;
                --$q;
                for ($yi = 0, $xi = $#$x - $#$y - 1; $yi <= $#$y; ++$yi, ++$xi) {
                    $x->[$xi] -= $BASE
                      if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
                }
            }
        }
        pop(@$x);
        unshift(@q, $q);
    }
    if (wantarray) {
        my $d = bless [], $c;
        if ($dd != 1) {
            $car = 0;
            for $xi (reverse @$x) {
                $prd = $car * $BASE + $xi;
                $car = $prd - ($tmp = int($prd / $dd)) * $dd;
                unshift(@$d, $tmp);
            }
        } else {
            @$d = @$x;
        }
        @$x = @q;
        __strip_zeros($x);
        __strip_zeros($d);
        return ($x, $d);
    }
    @$x = @q;
    __strip_zeros($x);
    $x;
}

sub _div_use_div {
    # ref to array, ref to array, modify first array and return remainder if
    # in list context
    my ($c, $x, $yorg) = @_;

# the general div algorithm here is about O(N*N) and thus quite slow, so
    # we first check for some special cases and use shortcuts to handle them.

# This works, because we store the numbers in a chunked format where each
    # element contains 5..7 digits (depending on system).

# check whether y has more elements than x, if yet, the result will be 0
    if (@$yorg > @$x) {
        my $rem;
        $rem = $c->_copy($x) if wantarray;    # make copy
        @$x = 0;                        # set to 0
        return ($x, $rem) if wantarray; # including remainder?
        return $x;                      # only x, which is [0] now
    }
    # check whether the numbers have the same number of elements, in that case
    # the result will fit into one element and can be computed efficiently
    if (@$yorg == @$x) {
        my $rem;
        # if $yorg has more digits than $x (it's leading element is longer than
        # the one from $x), the result will also be 0:
        if (length(int($yorg->[-1])) > length(int($x->[-1]))) {
            $rem = $c->_copy($x) if wantarray;        # make copy
            @$x = 0;                            # set to 0
            return ($x, $rem) if wantarray;     # including remainder?
            return $x;
        }
        # now calculate $x / $yorg

if (length(int($yorg->[-1])) == length(int($x->[-1]))) {
            # same length, so make full compare

my $a = 0;
            my $j = @$x - 1;
            # manual way (abort if unequal, good for early ne)
            while ($j >= 0) {
                last if ($a = $x->[$j] - $yorg->[$j]);
                $j--;
            }
            # $a contains the result of the compare between X and Y
            # a < 0: x < y, a == 0: x == y, a > 0: x > y
            if ($a <= 0) {
                $rem = $c->_zero();                   # a = 0 => x == y => rem 0
                $rem = $c->_copy($x) if $a != 0;      # a < 0 => x < y => rem = x
                @$x = 0;
                $x->[0] = 0;                    # if $a < 0
                $x->[0] = 1 if $a == 0;         # $x == $y
                return ($x, $rem) if wantarray; # including remainder?
                return $x;
            }
            # $x >= $y, so proceed normally

}
    }

# all other cases:

my $y = $c->_copy($yorg);         # always make copy to preserve

my ($car, $bar, $prd, $dd, $xi, $yi, @q, $v2, $v1, @d, $tmp, $q, $u2, $u1, $u0);

# @q will accumulate the final result, $q contains the current computed
    # part of the final result

@q = ();
    ($v2, $v1) = @$y[-2, -1];
    $v2 = 0 unless $v2;
    while ($#$x > $#$y) {
        ($u2, $u1, $u0) = @$x[-3..-1];
        $u2 = 0 unless $u2;
        #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
        # if $v1 == 0;
        $q = (($u0 == $v1) ? $MAX_VAL : int(($u0 * $BASE + $u1) / $v1));
        --$q while ($v2 * $q > ($u0 * $BASE + $u1 - $q * $v1) * $BASE + $u2);
        if ($q) {
            ($car, $bar) = (0, 0);
            for ($yi = 0, $xi = $#$x - $#$y - 1; $yi <= $#$y; ++$yi, ++$xi) {
                $prd = $q * $y->[$yi] + $car;
                $prd -= ($car = int($prd / $BASE)) * $BASE;
                $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
            }
            if ($x->[-1] < $car + $bar) {
                $car = 0;
                --$q;
                for ($yi = 0, $xi = $#$x - $#$y - 1; $yi <= $#$y; ++$yi, ++$xi) {
                    $x->[$xi] -= $BASE
                      if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
                }
            }
        }
        pop(@$x);
        unshift(@q, $q);
    }
    if (wantarray) {
        my $d = bless [], $c;
        if ($dd != 1) {
            $car = 0;
            for $xi (reverse @$x) {
                $prd = $car * $BASE + $xi;
                $car = $prd - ($tmp = int($prd / $dd)) * $dd;
                unshift(@$d, $tmp);
            }
        } else {
            @$d = @$x;
        }
        @$x = @q;
        __strip_zeros($x);
        __strip_zeros($d);
        return ($x, $d);
    }
    @$x = @q;
    __strip_zeros($x);
    $x;
}

##############################################################################
# testing

sub _acmp {
    # Internal absolute post-normalized compare (ignore signs)
    # ref to array, ref to array, return <0, 0, >0
    # Arrays must have at least one entry; this is not checked for.
    my ($c, $cx, $cy) = @_;

# shortcut for short numbers
    return (($cx->[0] <=> $cy->[0]) <=> 0)
      if @$cx == 1 && @$cy == 1;

# fast comp based on number of array elements (aka pseudo-length)
    my $lxy = (@$cx - @$cy)
      # or length of first element if same number of elements (aka difference 0)
      ||
        # need int() here because sometimes the last element is '00018' vs '18'
        (length(int($cx->[-1])) - length(int($cy->[-1])));

return -1 if $lxy < 0;      # already differs, ret
    return  1 if $lxy > 0;      # ditto

# manual way (abort if unequal, good for early ne)
    my $a;
    my $j = @$cx;
    while (--$j >= 0) {
        last if $a = $cx->[$j] - $cy->[$j];
    }
    $a <=> 0;
}

sub _len {
    # compute number of digits in base 10

# int() because add/sub sometimes leaves strings (like '00005') instead of
    # '5' in this place, thus causing length() to report wrong length
    my $cx = $_[1];

(@$cx - 1) * $BASE_LEN + length(int($cx->[-1]));
}

sub _digit {
    # Return the nth digit. Zero is rightmost, so _digit(123, 0) gives 3.
    # Negative values count from the left, so _digit(123, -1) gives 1.
    my ($c, $x, $n) = @_;

my $len = _len('', $x);

$n += $len if $n < 0;               # -1 last, -2 second-to-last

# Math::BigInt::Calc returns 0 if N is out of range, but this is not done
    # by the other backend libraries.

return "0" if $n < 0 || $n >= $len; # return 0 for digits out of range

my $elem = int($n / $BASE_LEN);     # index of array element
    my $digit = $n % $BASE_LEN;         # index of digit within the element
    substr("0" x $BASE_LEN . "$x->[$elem]", -1 - $digit, 1);
}

sub _zeros {
    # Return number of trailing zeros in decimal.
    # Check each array element for having 0 at end as long as elem == 0
    # Upon finding a elem != 0, stop.

my $x = $_[1];

return 0 if @$x == 1 && $x->[0] == 0;

my $zeros = 0;
    foreach my $elem (@$x) {
        if ($elem != 0) {
            $elem =~ /[^0](0*)\z/;
            $zeros += length($1);       # count trailing zeros
            last;                       # early out
        }
        $zeros += $BASE_LEN;
    }
    $zeros;
}

##############################################################################
# _is_* routines

sub _is_zero {
    # return true if arg is zero
    @{$_[1]} == 1 && $_[1]->[0] == 0 ? 1 : 0;
}

sub _is_even {
    # return true if arg is even
    $_[1]->[0] & 1 ? 0 : 1;
}

sub _is_odd {
    # return true if arg is odd
    $_[1]->[0] & 1 ? 1 : 0;
}

sub _is_one {
    # return true if arg is one
    @{$_[1]} == 1 && $_[1]->[0] == 1 ? 1 : 0;
}

sub _is_two {
    # return true if arg is two
    @{$_[1]} == 1 && $_[1]->[0] == 2 ? 1 : 0;
}

sub _is_ten {
    # return true if arg is ten
    @{$_[1]} == 1 && $_[1]->[0] == 10 ? 1 : 0;
}

sub __strip_zeros {
    # Internal normalization function that strips leading zeros from the array.
    # Args: ref to array
    my $x = shift;

push @$x, 0 if @$x == 0;    # div might return empty results, so fix it
    return $x if @$x == 1;      # early out

#print "strip: cnt $cnt i $i\n";
    # '0', '3', '4', '0', '0',
    #  0    1    2    3    4
    # cnt = 5, i = 4
    # i = 4
    # i = 3
    # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
    # >= 1: skip first part (this can be zero)

my $i = $#$x;
    while ($i > 0) {
        last if $x->[$i] != 0;
        $i--;
    }
    $i++;
    splice(@$x, $i) if $i < @$x;
    $x;
}

###############################################################################
# check routine to test internal state for corruptions

sub _check {
    # used by the test suite
    my ($class, $x) = @_;

my $msg = $class -> SUPER::_check($x);
    return $msg if $msg;

my $n;
    eval { $n = @$x };
    return "Not an array reference" unless $@ eq '';

return "Reference to an empty array" unless $n > 0;

# The following fails with Math::BigInt::FastCalc because a
    # Math::BigInt::FastCalc "object" is an unblessed array ref.
    #
    #return 0 unless ref($x) eq $class;

for (my $i = 0 ; $i <= $#$x ; ++ $i) {
        my $e = $x -> [$i];

return "Element at index $i is undefined"
          unless defined $e;

return "Element at index $i is a '" . ref($e) .
          "', which is not a scalar"
          unless ref($e) eq "";

# It would be better to use the regex /^([1-9]\d*|0)\z/, but that fails
        # in Math::BigInt::FastCalc, because it sometimes creates array
        # elements like "000000".
        return "Element at index $i is '$e', which does not look like an" .
          " normal integer" unless $e =~ /^\d+\z/;

return "Element at index $i is '$e', which is not smaller than" .
          " the base '$BASE'" if $e >= $BASE;

return "Element at index $i (last element) is zero"
          if $#$x > 0 && $i == $#$x && $e == 0;
    }

return 0;
}

###############################################################################

sub _mod {
    # if possible, use mod shortcut
    my ($c, $x, $yo) = @_;

# slow way since $y too big
    if (@$yo > 1) {
        my ($xo, $rem) = $c->_div($x, $yo);
        @$x = @$rem;
        return $x;
    }

my $y = $yo->[0];

# if both are single element arrays
    if (@$x == 1) {
        $x->[0] %= $y;
        return $x;
    }

# if @$x has more than one element, but @$y is a single element
    my $b = $BASE % $y;
    if ($b == 0) {
        # when BASE % Y == 0 then (B * BASE) % Y == 0
        # (B * BASE) % $y + A % Y => A % Y
        # so need to consider only last element: O(1)
        $x->[0] %= $y;
    } elsif ($b == 1) {
        # else need to go through all elements in @$x: O(N), but loop is a bit
        # simplified
        my $r = 0;
        foreach (@$x) {
            $r = ($r + $_) % $y; # not much faster, but heh...
            #$r += $_ % $y; $r %= $y;
        }
        $r = 0 if $r == $y;
        $x->[0] = $r;
    } else {
        # else need to go through all elements in @$x: O(N)
        my $r = 0;
        my $bm = 1;
        foreach (@$x) {
            $r = ($_ * $bm + $r) % $y;
            $bm = ($bm * $b) % $y;

#$r += ($_ % $y) * $bm;
            #$bm *= $b;
            #$bm %= $y;
            #$r %= $y;
        }
        $r = 0 if $r == $y;
        $x->[0] = $r;
    }
    @$x = $x->[0];              # keep one element of @$x
    return $x;
}

##############################################################################
# shifts

sub _rsft {
    my ($c, $x, $y, $n) = @_;

if ($n != 10) {
        $n = $c->_new($n);
        return scalar $c->_div($x, $c->_pow($n, $y));
    }

# shortcut (faster) for shifting by 10)
    # multiples of $BASE_LEN
    my $dst = 0;                # destination
    my $src = $c->_num($y);     # as normal int
    my $xlen = (@$x - 1) * $BASE_LEN + length(int($x->[-1]));
    if ($src >= $xlen or ($src == $xlen and !defined $x->[1])) {
        # 12345 67890 shifted right by more than 10 digits => 0
        splice(@$x, 1);         # leave only one element
        $x->[0] = 0;            # set to zero
        return $x;
    }
    my $rem = $src % $BASE_LEN;   # remainder to shift
    $src = int($src / $BASE_LEN); # source
    if ($rem == 0) {
        splice(@$x, 0, $src);   # even faster, 38.4 => 39.3
    } else {
        my $len = @$x - $src;   # elems to go
        my $vd;
        my $z = '0' x $BASE_LEN;
        $x->[ @$x ] = 0;          # avoid || 0 test inside loop
        while ($dst < $len) {
            $vd = $z . $x->[$src];
            $vd = substr($vd, -$BASE_LEN, $BASE_LEN - $rem);
            $src++;
            $vd = substr($z . $x->[$src], -$rem, $rem) . $vd;
            $vd = substr($vd, -$BASE_LEN, $BASE_LEN) if length($vd) > $BASE_LEN;
            $x->[$dst] = int($vd);
            $dst++;
        }
        splice(@$x, $dst) if $dst > 0;       # kill left-over array elems
        pop @$x if $x->[-1] == 0 && @$x > 1; # kill last element if 0
    }                                        # else rem == 0
    $x;
}

sub _lsft {
    my ($c, $x, $n, $b) = @_;

return $x if $c->_is_zero($x);

# Handle the special case when the base is a power of 10. Don't check
    # whether log($b)/log(10) is an integer, because log(1000)/log(10) is not
    # exactly 3.

my $log10 = sprintf "%.0f", log($b) / log(10);
    if ($b == 10 ** $log10) {
        $b = 10;
        $n = $c->_mul($n, $c->_new($log10));

# shortcut (faster) for shifting by 10) since we are in base 10eX
        # multiples of $BASE_LEN:
        my $src = @$x;                      # source
        my $len = $c->_num($n);             # shift-len as normal int
        my $rem = $len % $BASE_LEN;         # remainder to shift
        my $dst = $src + int($len / $BASE_LEN); # destination
        my $vd;                                 # further speedup
        $x->[$src] = 0;                         # avoid first ||0 for speed
        my $z = '0' x $BASE_LEN;
        while ($src >= 0) {
            $vd = $x->[$src];
            $vd = $z . $vd;
            $vd = substr($vd, -$BASE_LEN + $rem, $BASE_LEN - $rem);
            $vd .= $src > 0 ? substr($z . $x->[$src - 1], -$BASE_LEN, $rem)
              : '0' x $rem;
            $vd = substr($vd, -$BASE_LEN, $BASE_LEN) if length($vd) > $BASE_LEN;
            $x->[$dst] = int($vd);
            $dst--;
            $src--;
        }
        # set lowest parts to 0
        while ($dst >= 0) {
            $x->[$dst--] = 0;
        }
        # fix spurious last zero element
        splice @$x, -1 if $x->[-1] == 0;
        return $x;
    } else {
        $b = $c->_new($b);
        #print $c->_str($b);
        return $c->_mul($x, $c->_pow($b, $n));
    }
}

sub _pow {
    # power of $x to $y
    # ref to array, ref to array, return ref to array
    my ($c, $cx, $cy) = @_;

if (@$cy == 1 && $cy->[0] == 0) {
        splice(@$cx, 1);
        $cx->[0] = 1;           # y == 0 => x => 1
        return $cx;
    }

if ((@$cx == 1 && $cx->[0] == 1) || #    x == 1
        (@$cy == 1 && $cy->[0] == 1))   # or y == 1
    {
        return $cx;
    }

if (@$cx == 1 && $cx->[0] == 0) {
        splice (@$cx, 1);
        $cx->[0] = 0;           # 0 ** y => 0 (if not y <= 0)
        return $cx;
    }

my $pow2 = $c->_one();

my $y_bin = $c->_as_bin($cy);
    $y_bin =~ s/^0b//;
    my $len = length($y_bin);
    while (--$len > 0) {
        $c->_mul($pow2, $cx) if substr($y_bin, $len, 1) eq '1'; # is odd?
        $c->_mul($cx, $cx);
    }

$c->_mul($cx, $pow2);
    $cx;
}

sub _nok {
    # Return binomial coefficient (n over k).
    # Given refs to arrays, return ref to array.
    # First input argument is modified.

my ($c, $n, $k) = @_;

# If k > n/2, or, equivalently, 2*k > n, compute nok(n, k) as
    # nok(n, n-k), to minimize the number if iterations in the loop.

{
        my $twok = $c->_mul($c->_two(), $c->_copy($k)); # 2 * k
        if ($c->_acmp($twok, $n) > 0) {               # if 2*k > n
            $k = $c->_sub($c->_copy($n), $k);         # k = n - k
        }
    }

# Example:
    #
    # / 7 \       7!       1*2*3*4 * 5*6*7   5 * 6 * 7       6   7
    # |   | = --------- =  --------------- = --------- = 5 * - * -
    # \ 3 /   (7-3)! 3!    1*2*3*4 * 1*2*3   1 * 2 * 3       2   3

if ($c->_is_zero($k)) {
        @$n = 1;
    } else {

# Make a copy of the original n, since we'll be modifying n in-place.

my $n_orig = $c->_copy($n);

# n = 5, f = 6, d = 2 (cf. example above)

$c->_sub($n, $k);
        $c->_inc($n);

my $f = $c->_copy($n);
        $c->_inc($f);

my $d = $c->_two();

# while f <= n (the original n, that is) ...

while ($c->_acmp($f, $n_orig) <= 0) {

# n = (n * f / d) == 5 * 6 / 2 (cf. example above)

$c->_mul($n, $f);
            $c->_div($n, $d);

# f = 7, d = 3 (cf. example above)

$c->_inc($f);
            $c->_inc($d);
        }

}

return $n;
}

my @factorials = (
                  1,
                  1,
                  2,
                  2*3,
                  2*3*4,
                  2*3*4*5,
                  2*3*4*5*6,
                  2*3*4*5*6*7,
                 );

sub _fac {
    # factorial of $x
    # ref to array, return ref to array
    my ($c, $cx) = @_;

if ((@$cx == 1) && ($cx->[0] <= 7)) {
        $cx->[0] = $factorials[$cx->[0]]; # 0 => 1, 1 => 1, 2 => 2 etc.
        return $cx;
    }

if ((@$cx == 1) &&          # we do this only if $x >= 12 and $x <= 7000
        ($cx->[0] >= 12 && $cx->[0] < 7000)) {

# Calculate (k-j) * (k-j+1) ... k .. (k+j-1) * (k + j)
        # See http://blogten.blogspot.com/2007/01/calculating-n.html
        # The above series can be expressed as factors:
        #   k * k - (j - i) * 2
        # We cache k*k, and calculate (j * j) as the sum of the first j odd integers

# This will not work when N exceeds the storage of a Perl scalar, however,
        # in this case the algorithm would be way too slow to terminate, anyway.

# As soon as the last element of $cx is 0, we split it up and remember
        # how many zeors we got so far. The reason is that n! will accumulate
        # zeros at the end rather fast.
        my $zero_elements = 0;

# If n is even, set n = n -1
        my $k = $c->_num($cx);
        my $even = 1;
        if (($k & 1) == 0) {
            $even = $k;
            $k --;
        }
        # set k to the center point
        $k = ($k + 1) / 2;
        #  print "k $k even: $even\n";
        # now calculate k * k
        my $k2 = $k * $k;
        my $odd = 1;
        my $sum = 1;
        my $i = $k - 1;
        # keep reference to x
        my $new_x = $c->_new($k * $even);
        @$cx = @$new_x;
        if ($cx->[0] == 0) {
            $zero_elements ++;
            shift @$cx;
        }
        #  print STDERR "x = ", $c->_str($cx), "\n";
        my $BASE2 = int(sqrt($BASE))-1;
        my $j = 1;
        while ($j <= $i) {
            my $m = ($k2 - $sum);
            $odd += 2;
            $sum += $odd;
            $j++;
            while ($j <= $i && ($m < $BASE2) && (($k2 - $sum) < $BASE2)) {
                $m *= ($k2 - $sum);
                $odd += 2;
                $sum += $odd;
                $j++;
                #      print STDERR "\n k2 $k2 m $m sum $sum odd $odd\n"; sleep(1);
            }
            if ($m < $BASE) {
                $c->_mul($cx, [$m]);
            } else {
                $c->_mul($cx, $c->_new($m));
            }
            if ($cx->[0] == 0) {
                $zero_elements ++;
                shift @$cx;
            }
            #    print STDERR "Calculate $k2 - $sum = $m (x = ", $c->_str($cx), ")\n";
        }
        # multiply in the zeros again
        unshift @$cx, (0) x $zero_elements;
        return $cx;
    }

# go forward until $base is exceeded limit is either $x steps (steps == 100
    # means a result always too high) or $base.
    my $steps = 100;
    $steps = $cx->[0] if @$cx == 1;
    my $r = 2;
    my $cf = 3;
    my $step = 2;
    my $last = $r;
    while ($r * $cf < $BASE && $step < $steps) {
        $last = $r;
        $r *= $cf++;
        $step++;
    }
    if ((@$cx == 1) && $step == $cx->[0]) {
        # completely done, so keep reference to $x and return
        $cx->[0] = $r;
        return $cx;
    }

# now we must do the left over steps
    my $n;                      # steps still to do
    if (@$cx == 1) {
        $n = $cx->[0];
    } else {
        $n = $c->_copy($cx);
    }

# Set $cx to the last result below $BASE (but keep ref to $x)
    $cx->[0] = $last;
    splice (@$cx, 1);
    # As soon as the last element of $cx is 0, we split it up and remember
    # how many zeors we got so far. The reason is that n! will accumulate
    # zeros at the end rather fast.
    my $zero_elements = 0;

# do left-over steps fit into a scalar?
    if (ref $n eq 'ARRAY') {
        # No, so use slower inc() & cmp()
        # ($n is at least $BASE here)
        my $base_2 = int(sqrt($BASE)) - 1;
        #print STDERR "base_2: $base_2\n";
        while ($step < $base_2) {
            if ($cx->[0] == 0) {
                $zero_elements ++;
                shift @$cx;
            }
            my $b = $step * ($step + 1);
            $step += 2;
            $c->_mul($cx, [$b]);
        }
        $step = [$step];
        while ($c->_acmp($step, $n) <= 0) {
            if ($cx->[0] == 0) {
                $zero_elements ++;
                shift @$cx;
            }
            $c->_mul($cx, $step);
            $c->_inc($step);
        }
    } else {
        # Yes, so we can speed it up slightly

#    print "# left over steps $n\n";

my $base_4 = int(sqrt(sqrt($BASE))) - 2;
        #print STDERR "base_4: $base_4\n";
        my $n4 = $n - 4;
        while ($step < $n4 && $step < $base_4) {
            if ($cx->[0] == 0) {
                $zero_elements ++;
                shift @$cx;
            }
            my $b = $step * ($step + 1);
            $step += 2;
            $b *= $step * ($step + 1);
            $step += 2;
            $c->_mul($cx, [$b]);
        }
        my $base_2 = int(sqrt($BASE)) - 1;
        my $n2 = $n - 2;
        #print STDERR "base_2: $base_2\n";
        while ($step < $n2 && $step < $base_2) {
            if ($cx->[0] == 0) {
                $zero_elements ++;
                shift @$cx;
            }
            my $b = $step * ($step + 1);
            $step += 2;
            $c->_mul($cx, [$b]);
        }
        # do what's left over
        while ($step <= $n) {
            $c->_mul($cx, [$step]);
            $step++;
            if ($cx->[0] == 0) {
                $zero_elements ++;
                shift @$cx;
            }
        }
    }
    # multiply in the zeros again
    unshift @$cx, (0) x $zero_elements;
    $cx;                        # return result
}

sub _log_int {
    # calculate integer log of $x to base $base
    # ref to array, ref to array - return ref to array
    my ($c, $x, $base) = @_;

# X == 0 => NaN
    return if @$x == 1 && $x->[0] == 0;

# BASE 0 or 1 => NaN
    return if @$base == 1 && $base->[0] < 2;

# X == 1 => 0 (is exact)
    if (@$x == 1 && $x->[0] == 1) {
        @$x = 0;
        return $x, 1;
    }

my $cmp = $c->_acmp($x, $base);

# X == BASE => 1 (is exact)
    if ($cmp == 0) {
        @$x = 1;
        return $x, 1;
    }

# 1 < X < BASE => 0 (is truncated)
    if ($cmp < 0) {
        @$x = 0;
        return $x, 0;
    }

my $x_org = $c->_copy($x);  # preserve x

# Compute a guess for the result based on:
    # $guess = int ( length_in_base_10(X) / ( log(base) / log(10) ) )
    my $len = $c->_len($x_org);
    my $log = log($base->[-1]) / log(10);

# for each additional element in $base, we add $BASE_LEN to the result,
    # based on the observation that log($BASE, 10) is BASE_LEN and
    # log(x*y) == log(x) + log(y):
    $log += (@$base - 1) * $BASE_LEN;

# calculate now a guess based on the values obtained above:
    my $res = int($len / $log);

@$x = $res;
    my $trial = $c->_pow($c->_copy($base), $x);
    my $acmp = $c->_acmp($trial, $x_org);

# Did we get the exact result?

return $x, 1 if $acmp == 0;

# Too small?

while ($acmp < 0) {
        $c->_mul($trial, $base);
        $c->_inc($x);
        $acmp = $c->_acmp($trial, $x_org);
    }

# Too big?

while ($acmp > 0) {
        $c->_div($trial, $base);
        $c->_dec($x);
        $acmp = $c->_acmp($trial, $x_org);
    }

return $x, 1 if $acmp == 0;         # result is exact
    return $x, 0;                       # result is too small
}

# for debugging:
use constant DEBUG => 0;
my $steps = 0;
sub steps { $steps };

sub _sqrt {
    # square-root of $x in place
    # Compute a guess of the result (by rule of thumb), then improve it via
    # Newton's method.
    my ($c, $x) = @_;

if (@$x == 1) {
        # fits into one Perl scalar, so result can be computed directly
        $x->[0] = int(sqrt($x->[0]));
        return $x;
    }
    my $y = $c->_copy($x);
    # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
    # since our guess will "grow"
    my $l = int(($c->_len($x)-1) / 2);

my $lastelem = $x->[-1];    # for guess
    my $elems = @$x - 1;
    # not enough digits, but could have more?
    if ((length($lastelem) <= 3) && ($elems > 1)) {
        # right-align with zero pad
        my $len = length($lastelem) & 1;
        print "$lastelem => " if DEBUG;
        $lastelem .= substr($x->[-2] . '0' x $BASE_LEN, 0, $BASE_LEN);
        # former odd => make odd again, or former even to even again
        $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
        print "$lastelem\n" if DEBUG;
    }

# construct $x (instead of $c->_lsft($x, $l, 10)
    my $r = $l % $BASE_LEN;     # 10000 00000 00000 00000 ($BASE_LEN=5)
    $l = int($l / $BASE_LEN);
    print "l =  $l " if DEBUG;

splice @$x, $l;              # keep ref($x), but modify it

# we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
    # that gives us:
    # 14400 00000 => sqrt(14400) => guess first digits to be 120
    # 144000 000000 => sqrt(144000) => guess 379

print "$lastelem (elems $elems) => " if DEBUG;
    $lastelem = $lastelem / 10 if ($elems & 1 == 1); # odd or even?
    my $g = sqrt($lastelem);
    $g =~ s/\.//;               # 2.345 => 2345
    $r -= 1 if $elems & 1 == 0; # 70 => 7

# padd with zeros if result is too short
    $x->[$l--] = int(substr($g . '0' x $r, 0, $r+1));
    print "now ", $x->[-1] if DEBUG;
    print " would have been ", int('1' . '0' x $r), "\n" if DEBUG;

# If @$x > 1, we could compute the second elem of the guess, too, to create
    # an even better guess. Not implemented yet. Does it improve performance?
    $x->[$l--] = 0 while ($l >= 0); # all other digits of guess are zero

print "start x= ", $c->_str($x), "\n" if DEBUG;
    my $two = $c->_two();
    my $last = $c->_zero();
    my $lastlast = $c->_zero();
    $steps = 0 if DEBUG;
    while ($c->_acmp($last, $x) != 0 && $c->_acmp($lastlast, $x) != 0) {
        $steps++ if DEBUG;
        $lastlast = $c->_copy($last);
        $last = $c->_copy($x);
        $c->_add($x, $c->_div($c->_copy($y), $x));
        $c->_div($x, $two );
        print " x= ", $c->_str($x), "\n" if DEBUG;
    }
    print "\nsteps in sqrt: $steps, " if DEBUG;
    $c->_dec($x) if $c->_acmp($y, $c->_mul($c->_copy($x), $x)) < 0; # overshot?
    print " final ", $x->[-1], "\n" if DEBUG;
    $x;
}

sub _root {
    # Take n'th root of $x in place.

my ($c, $x, $n) = @_;

# Small numbers.

if (@$x == 1 && @$n == 1) {
        # Result can be computed directly. Adjust initial result for numerical
        # errors, e.g., int(1000**(1/3)) is 2, not 3.
        my $y = int($x->[0] ** (1 / $n->[0]));
        my $yp1 = $y + 1;
        $y = $yp1 if $yp1 ** $n->[0] == $x->[0];
        $x->[0] = $y;
        return $x;
    }

# If x <= n, the result is always (truncated to) 1.

if ((@$x > 1 || $x -> [0] > 0) &&           # if x is non-zero ...
        $c -> _acmp($x, $n) <= 0)               # ... and x <= n
    {
        my $one = $x -> _one();
        @$x = @$one;
        return $x;
    }

# If $n is a power of two, take sqrt($x) repeatedly, e.g., root($x, 4) =
    # sqrt(sqrt($x)), root($x, 8) = sqrt(sqrt(sqrt($x))).

my $b = $c -> _as_bin($n);
    if ($b =~ /0b1(0+)$/) {
        my $count = length($1);       # 0b100 => len('00') => 2
        my $cnt = $count;             # counter for loop
        unshift @$x, 0;               # add one element, together with one
                                      #   more below in the loop this makes 2
        while ($cnt-- > 0) {
            # 'Inflate' $x by adding one element, basically computing
            # $x * $BASE * $BASE. This gives us more $BASE_LEN digits for
            # result since len(sqrt($X)) approx == len($x) / 2.
            unshift @$x, 0;
            # Calculate sqrt($x), $x is now one element to big, again. In the
            # next round we make that two, again.
            $c -> _sqrt($x);
        }

# $x is now one element too big, so truncate result by removing it.
        shift @$x;

return $x;
    }

my $DEBUG = 0;

# Now the general case. This works by finding an initial guess. If this
    # guess is incorrect, a relatively small delta is chosen. This delta is
    # used to find a lower and upper limit for the correct value. The delta is
    # doubled in each iteration. When a lower and upper limit is found,
    # bisection is applied to narrow down the region until we have the correct
    # value.

# Split x into mantissa and exponent in base 10, so that
    #
    #   x = xm * 10^xe, where 0 < xm < 1 and xe is an integer

my $x_str = $c -> _str($x);
    my $xm    = "." . $x_str;
    my $xe    = length($x_str);

# From this we compute the base 10 logarithm of x
    #
    #   log_10(x) = log_10(xm) + log_10(xe^10)
    #             = log(xm)/log(10) + xe
    #
    # and then the base 10 logarithm of y, where y = x^(1/n)
    #
    #   log_10(y) = log_10(x)/n

my $log10x = log($xm) / log(10) + $xe;
    my $log10y = $log10x / $c -> _num($n);

# And from this we compute ym and ye, the mantissa and exponent (in
    # base 10) of y, where 1 < ym <= 10 and ye is an integer.

my $ye = int $log10y;
    my $ym = 10 ** ($log10y - $ye);

# Finally, we scale the mantissa and exponent to incraese the integer
    # part of ym, before building the string representing our guess of y.

if ($DEBUG) {
        print "\n";
        print "xm     = $xm\n";
        print "xe     = $xe\n";
        print "log10x = $log10x\n";
        print "log10y = $log10y\n";
        print "ym     = $ym\n";
        print "ye     = $ye\n";
        print "\n";
    }

my $d = $ye < 15 ? $ye : 15;
    $ym *= 10 ** $d;
    $ye -= $d;

my $y_str = sprintf('%.0f', $ym) . "0" x $ye;
    my $y = $c -> _new($y_str);

if ($DEBUG) {
        print "ym     = $ym\n";
        print "ye     = $ye\n";
        print "\n";
        print "y_str  = $y_str (initial guess)\n";
        print "\n";
    }

# See if our guess y is correct.

my $trial = $c -> _pow($c -> _copy($y), $n);
    my $acmp  = $c -> _acmp($trial, $x);

if ($acmp == 0) {
        @$x = @$y;
        return $x;
    }

# Find a lower and upper limit for the correct value of y. Start off with a
    # delta value that is approximately the size of the accuracy of the guess.

my $lower;
    my $upper;

my $delta = $c -> _new("1" . ("0" x $ye));
    my $two   = $c -> _two();

if ($acmp < 0) {
        $lower = $y;
        while ($acmp < 0) {
            $upper = $c -> _add($c -> _copy($lower), $delta);

if ($DEBUG) {
                print "lower  = $lower\n";
                print "upper  = $upper\n";
                print "delta  = $delta\n";
                print "\n";
            }
            $acmp  = $c -> _acmp($c -> _pow($c -> _copy($upper), $n), $x);
            if ($acmp == 0) {
                @$x = @$upper;
                return $x;
            }
            $delta = $c -> _mul($delta, $two);
        }
    }

elsif ($acmp > 0) {
        $upper = $y;
        my $zero = $c -> _zero();
        while ($acmp > 0) {
            if ($c -> _acmp($upper, $delta) <= 0) {
                $lower = $c -> _zero();
                last;
            }
            $lower = $c -> _sub($c -> _copy($upper), $delta);

if ($DEBUG) {
                print "lower  = $lower\n";
                print "upper  = $upper\n";
                print "delta  = $delta\n";
                print "\n";
            }
            $acmp  = $c -> _acmp($c -> _pow($c -> _copy($lower), $n), $x);
            if ($acmp == 0) {
                @$x = @$lower;
                return $x;
            }
            $delta = $c -> _mul($delta, $two);
        }
    }

# Use bisection to narrow down the interval.

my $one = $c -> _one();
    {

$delta = $c -> _sub($c -> _copy($upper), $lower);
        if ($c -> _acmp($delta, $one) <= 0) {
            @$x = @$lower;
            return $x;
        }

if ($DEBUG) {
            print "lower  = $lower\n";
            print "upper  = $upper\n";
            print "delta   = $delta\n";
            print "\n";
        }

$delta = $c -> _div($delta, $two);
        my $middle = $c -> _add($c -> _copy($lower), $delta);

$acmp  = $c -> _acmp($c -> _pow($c -> _copy($middle), $n), $x);
        if ($acmp < 0) {
            $lower = $middle;
        } elsif ($acmp > 0) {
            $upper = $middle;
        } else {
            @$x = @$middle;
            return $x;
        }

redo;
    }

$x;
}

##############################################################################
# binary stuff

sub _and {
    my ($c, $x, $y) = @_;

# the shortcut makes equal, large numbers _really_ fast, and makes only a
    # very small performance drop for small numbers (e.g. something with less
    # than 32 bit) Since we optimize for large numbers, this is enabled.
    return $x if $c->_acmp($x, $y) == 0; # shortcut

my $m = $c->_one();
    my ($xr, $yr);
    my $mask = $AND_MASK;

my $x1 = $c->_copy($x);
    my $y1 = $c->_copy($y);
    my $z  = $c->_zero();

use integer;
    until ($c->_is_zero($x1) || $c->_is_zero($y1)) {
        ($x1, $xr) = $c->_div($x1, $mask);
        ($y1, $yr) = $c->_div($y1, $mask);

$c->_add($z, $c->_mul([ 0 + $xr->[0] & 0 + $yr->[0] ], $m));
        $c->_mul($m, $mask);
    }

@$x = @$z;
    return $x;
}

sub _xor {
    my ($c, $x, $y) = @_;

return $c->_zero() if $c->_acmp($x, $y) == 0; # shortcut (see -and)

my $m = $c->_one();
    my ($xr, $yr);
    my $mask = $XOR_MASK;

my $x1 = $c->_copy($x);
    my $y1 = $c->_copy($y);      # make copy
    my $z  = $c->_zero();

use integer;
    until ($c->_is_zero($x1) || $c->_is_zero($y1)) {
        ($x1, $xr) = $c->_div($x1, $mask);
        ($y1, $yr) = $c->_div($y1, $mask);
        # make ints() from $xr, $yr (see _and())
        #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
        #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
        #$c->_add($x, $c->_mul($c->_new($xrr ^ $yrr)), $m) );

$c->_add($z, $c->_mul([ 0 + $xr->[0] ^ 0 + $yr->[0] ], $m));
        $c->_mul($m, $mask);
    }
    # the loop stops when the shorter of the two numbers is exhausted
    # the remainder of the longer one will survive bit-by-bit, so we simple
    # multiply-add it in
    $c->_add($z, $c->_mul($x1, $m) ) if !$c->_is_zero($x1);
    $c->_add($z, $c->_mul($y1, $m) ) if !$c->_is_zero($y1);

@$x = @$z;
    return $x;
}

sub _or {
    my ($c, $x, $y) = @_;

return $x if $c->_acmp($x, $y) == 0; # shortcut (see _and)

my $m = $c->_one();
    my ($xr, $yr);
    my $mask = $OR_MASK;

my $x1 = $c->_copy($x);
    my $y1 = $c->_copy($y);      # make copy
    my $z  = $c->_zero();

use integer;
    until ($c->_is_zero($x1) || $c->_is_zero($y1)) {
        ($x1, $xr) = $c->_div($x1, $mask);
        ($y1, $yr) = $c->_div($y1, $mask);
        # make ints() from $xr, $yr (see _and())
        #    $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
        #    $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
        #    $c->_add($x, $c->_mul(_new( $c, ($xrr | $yrr) ), $m) );

$c->_add($z, $c->_mul([ 0 + $xr->[0] | 0 + $yr->[0] ], $m));
        $c->_mul($m, $mask);
    }
    # the loop stops when the shorter of the two numbers is exhausted
    # the remainder of the longer one will survive bit-by-bit, so we simple
    # multiply-add it in
    $c->_add($z, $c->_mul($x1, $m) ) if !$c->_is_zero($x1);
    $c->_add($z, $c->_mul($y1, $m) ) if !$c->_is_zero($y1);

@$x = @$z;
    return $x;
}

sub _as_hex {
    # convert a decimal number to hex (ref to array, return ref to string)
    my ($c, $x) = @_;

# fits into one element (handle also 0x0 case)
    return sprintf("0x%x", $x->[0]) if @$x == 1;

my $x1 = $c->_copy($x);

my $es = '';
    my ($xr, $h, $x10000);
    if ($] >= 5.006) {
        $x10000 = [ 0x10000 ];
        $h = 'h4';
    } else {
        $x10000 = [ 0x1000 ];
        $h = 'h3';
    }
    while (@$x1 != 1 || $x1->[0] != 0) # _is_zero()
    {
        ($x1, $xr) = $c->_div($x1, $x10000);
        $es .= unpack($h, pack('V', $xr->[0]));
    }
    $es = reverse $es;
    $es =~ s/^[0]+//;           # strip leading zeros
    '0x' . $es;                 # return result prepended with 0x
}

sub _as_bin {
    # convert a decimal number to bin (ref to array, return ref to string)
    my ($c, $x) = @_;

# fits into one element (and Perl recent enough), handle also 0b0 case
    # handle zero case for older Perls
    if ($] <= 5.005 && @$x == 1 && $x->[0] == 0) {
        my $t = '0b0';
        return $t;
    }
    if (@$x == 1 && $] >= 5.006) {
        my $t = sprintf("0b%b", $x->[0]);
        return $t;
    }
    my $x1 = $c->_copy($x);

my $es = '';
    my ($xr, $b, $x10000);
    if ($] >= 5.006) {
        $x10000 = [ 0x10000 ];
        $b = 'b16';
    } else {
        $x10000 = [ 0x1000 ];
        $b = 'b12';
    }
    while (!(@$x1 == 1 && $x1->[0] == 0)) # _is_zero()
    {
        ($x1, $xr) = $c->_div($x1, $x10000);
        $es .= unpack($b, pack('v', $xr->[0]));
    }
    $es = reverse $es;
    $es =~ s/^[0]+//;           # strip leading zeros
    '0b' . $es;                 # return result prepended with 0b
}

sub _as_oct {
    # convert a decimal number to octal (ref to array, return ref to string)
    my ($c, $x) = @_;

# fits into one element (handle also 0 case)
    return sprintf("0%o", $x->[0]) if @$x == 1;

my $x1 = $c->_copy($x);

my $es = '';
    my $xr;
    my $x1000 = [ 0100000 ];
    while (@$x1 != 1 || $x1->[0] != 0) # _is_zero()
    {
        ($x1, $xr) = $c->_div($x1, $x1000);
        $es .= reverse sprintf("%05o", $xr->[0]);
    }
    $es = reverse $es;
    $es =~ s/^0+//;             # strip leading zeros
    '0' . $es;                  # return result prepended with 0
}

sub _from_oct {
    # convert a octal number to decimal (string, return ref to array)
    my ($c, $os) = @_;

# for older Perls, play safe
    my $m = [ 0100000 ];
    my $d = 5;                  # 5 digits at a time

my $mul = $c->_one();
    my $x = $c->_zero();

my $len = int((length($os) - 1) / $d);      # $d digit parts, w/o the '0'
    my $val;
    my $i = -$d;
    while ($len >= 0) {
        $val = substr($os, $i, $d);             # get oct digits
        $val = CORE::oct($val);
        $i -= $d;
        $len --;
        my $adder = [ $val ];
        $c->_add($x, $c->_mul($adder, $mul)) if $val != 0;
        $c->_mul($mul, $m) if $len >= 0;        # skip last mul
    }
    $x;
}

sub _from_hex {
    # convert a hex number to decimal (string, return ref to array)
    my ($c, $hs) = @_;

my $m = $c->_new(0x10000000); # 28 bit at a time (<32 bit!)
    my $d = 7;                    # 7 digits at a time
    my $mul = $c->_one();
    my $x = $c->_zero();

my $len = int((length($hs) - 2) / $d); # $d digit parts, w/o the '0x'
    my $val;
    my $i = -$d;
    while ($len >= 0) {
        $val = substr($hs, $i, $d);     # get hex digits
        $val =~ s/^0x// if $len == 0; # for last part only because
        $val = CORE::hex($val);       # hex does not like wrong chars
        $i -= $d;
        $len --;
        my $adder = [ $val ];
        # if the resulting number was to big to fit into one element, create a
        # two-element version (bug found by Mark Lakata - Thanx!)
        if (CORE::length($val) > $BASE_LEN) {
            $adder = $c->_new($val);
        }
        $c->_add($x, $c->_mul($adder, $mul)) if $val != 0;
        $c->_mul($mul, $m) if $len >= 0; # skip last mul
    }
    $x;
}

sub _from_bin {
    # convert a hex number to decimal (string, return ref to array)
    my ($c, $bs) = @_;

# instead of converting X (8) bit at a time, it is faster to "convert" the
    # number to hex, and then call _from_hex.

my $hs = $bs;
    $hs =~ s/^[+-]?0b//;                                # remove sign and 0b
    my $l = length($hs);                                # bits
    $hs = '0' x (8 - ($l % 8)) . $hs if ($l % 8) != 0;  # padd left side w/ 0
    my $h = '0x' . unpack('H*', pack ('B*', $hs));      # repack as hex

$c->_from_hex($h);
}

##############################################################################
# special modulus functions

sub _modinv {
    # modular multiplicative inverse
    my ($c, $x, $y) = @_;

# modulo zero
    if ($c->_is_zero($y)) {
        return undef, undef;
    }

# modulo one
    if ($c->_is_one($y)) {
        return $c->_zero(), '+';
    }

my $u = $c->_zero();
    my $v = $c->_one();
    my $a = $c->_copy($y);
    my $b = $c->_copy($x);

# Euclid's Algorithm for bgcd(), only that we calc bgcd() ($a) and the result
    # ($u) at the same time. See comments in BigInt for why this works.
    my $q;
    my $sign = 1;
    {
        ($a, $q, $b) = ($b, $c->_div($a, $b));          # step 1
        last if $c->_is_zero($b);

my $t = $c->_add(                               # step 2:
                         $c->_mul($c->_copy($v), $q),   #  t =   v * q
                         $u);                           #      + u
        $u = $v;                                        #  u = v
        $v = $t;                                        #  v = t
        $sign = -$sign;
        redo;
    }

# if the gcd is not 1, then return NaN
    return (undef, undef) unless $c->_is_one($a);

($v, $sign == 1 ? '+' : '-');
}

sub _modpow {
    # modulus of power ($x ** $y) % $z
    my ($c, $num, $exp, $mod) = @_;

# a^b (mod 1) = 0 for all a and b
    if ($c->_is_one($mod)) {
        @$num = 0;
        return $num;
    }

# 0^a (mod m) = 0 if m != 0, a != 0
    # 0^0 (mod m) = 1 if m != 0
    if ($c->_is_zero($num)) {
        if ($c->_is_zero($exp)) {
            @$num = 1;
        } else {
            @$num = 0;
        }
        return $num;
    }

#  $num = $c->_mod($num, $mod);   # this does not make it faster

my $acc = $c->_copy($num);
    my $t = $c->_one();

my $expbin = $c->_as_bin($exp);
    $expbin =~ s/^0b//;
    my $len = length($expbin);
    while (--$len >= 0) {
        if (substr($expbin, $len, 1) eq '1') { # is_odd
            $t = $c->_mul($t, $acc);
            $t = $c->_mod($t, $mod);
        }
        $acc = $c->_mul($acc, $acc);
        $acc = $c->_mod($acc, $mod);
    }
    @$num = @$t;
    $num;
}

sub _gcd {
    # Greatest common divisor.

my ($c, $x, $y) = @_;

# gcd(0, 0) = 0
    # gcd(0, a) = a, if a != 0

if (@$x == 1 && $x->[0] == 0) {
        if (@$y == 1 && $y->[0] == 0) {
            @$x = 0;
        } else {
            @$x = @$y;
        }
        return $x;
    }

# Until $y is zero ...

until (@$y == 1 && $y->[0] == 0) {

# Compute remainder.

$c->_mod($x, $y);

# Swap $x and $y.

my $tmp = $c->_copy($x);
        @$x = @$y;
        $y = $tmp;              # no deref here; that would modify input $y
    }

return $x;
}

=pod

=head1 NAME

Math::BigInt::Calc - Pure Perl module to support Math::BigInt

=head1 SYNOPSIS

# to use it with Math::BigInt
    use Math::BigInt lib => 'Calc';

# to use it with Math::BigFloat
    use Math::BigFloat lib => 'Calc';

# to use it with Math::BigRat
    use Math::BigRat lib => 'Calc';

=head1 DESCRIPTION

Math::BigInt::Calc inherits from Math::BigInt::Lib.

In this library, the numbers are represented in base B = 10**N, where N is the
largest possible value that does not cause overflow in the intermediate
computations. The base B elements are stored in an array, with the least
significant element stored in array element zero. There are no leading zero
elements, except a single zero element when the number is zero.

For instance, if B = 10000, the number 1234567890 is represented internally
as [7890, 3456, 12].

=head1 SEE ALSO

L<Math::BigInt::Lib> for a description of the API.

Alternative libraries L<Math::BigInt::FastCalc>, L<Math::BigInt::GMP>, and
L<Math::BigInt::Pari>.

Some of the modules that use these libraries L<Math::BigInt>,
L<Math::BigFloat>, and L<Math::BigRat>.

=cut

Back