dist/Math-BigInt/lib/Math/BigInt/Calc.pm

   1 package Math::BigInt::Calc;
   2
   3 use 5.006002;
   4 use strict;
   5 # use warnings; # dont use warnings for older Perls
   6
   7 our $VERSION = '0.58';
   8
   9 # Package to store unsigned big integers in decimal and do math with them
  10
  11 # Internally the numbers are stored in an array with at least 1 element, no
  12 # leading zero parts (except the first) and in base 1eX where X is determined
  13 # automatically at loading time to be the maximum possible value
  14
  15 # todo:
  16 # - fully remove funky $# stuff in div() (maybe - that code scares me...)
  17
  18 # USE_MUL: due to problems on certain os (os390, posix-bc) "* 1e-5" is used
  19 # instead of "/ 1e5" at some places, (marked with USE_MUL). Other platforms
  20 # BS2000, some Crays need USE_DIV instead.
  21 # The BEGIN block is used to determine which of the two variants gives the
  22 # correct result.
  23
  24 # Beware of things like:
  25 # $i = $i * $y + $car; $car = int($i / $BASE); $i = $i % $BASE;
  26 # This works on x86, but fails on ARM (SA1100, iPAQ) due to whoknows what
  27 # reasons. So, use this instead (slower, but correct):
  28 # $i = $i * $y + $car; $car = int($i / $BASE); $i -= $BASE * $car;
  29
  30 ##############################################################################
  31 # global constants, flags and accessory
  32
  33 # announce that we are compatible with MBI v1.83 and up
  34 sub api_version () { 2; }
  35
  36 # constants for easier life
  37 my ($BASE,$BASE_LEN,$RBASE,$MAX_VAL);
  38 my ($AND_BITS,$XOR_BITS,$OR_BITS);
  39 my ($AND_MASK,$XOR_MASK,$OR_MASK);
  40
  41 sub _base_len
  42   {
  43   # Set/get the BASE_LEN and assorted other, connected values.
  44   # Used only by the testsuite, the set variant is used only by the BEGIN
  45   # block below:
  46   shift;
  47
  48   my ($b, $int) = @_;
  49   if (defined $b)
  50     {
  51     # avoid redefinitions
  52     undef &_mul;
  53     undef &_div;
  54
  55     if ($] >= 5.008 && $int && $b > 7)
  56       {
  57       $BASE_LEN = $b;
  58       *_mul = \&_mul_use_div_64;
  59       *_div = \&_div_use_div_64;
  60       $BASE = int("1e".$BASE_LEN);
  61       $MAX_VAL = $BASE-1;
  62       return $BASE_LEN unless wantarray;
  63       return ($BASE_LEN, $BASE, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN, $MAX_VAL,);
  64       }
  65
  66     # find whether we can use mul or div in mul()/div()
  67     $BASE_LEN = $b+1;
  68     my $caught = 0;
  69     while (--$BASE_LEN > 5)
  70       {
  71       $BASE = int("1e".$BASE_LEN);
  72       $RBASE = abs('1e-'.$BASE_LEN);                    # see USE_MUL
  73       $caught = 0;
  74       $caught += 1 if (int($BASE * $RBASE) != 1);       # should be 1
  75       $caught += 2 if (int($BASE / $BASE) != 1);        # should be 1
  76       last if $caught != 3;
  77       }
  78     $BASE = int("1e".$BASE_LEN);
  79     $RBASE = abs('1e-'.$BASE_LEN);                      # see USE_MUL
  80     $MAX_VAL = $BASE-1;
  81
  82     # ($caught & 1) != 0 => cannot use MUL
  83     # ($caught & 2) != 0 => cannot use DIV
  84     if ($caught == 2)                           # 2
  85       {
  86       # must USE_MUL since we cannot use DIV
  87       *_mul = \&_mul_use_mul;
  88       *_div = \&_div_use_mul;
  89       }
  90     else                                        # 0 or 1
  91       {
  92       # can USE_DIV instead
  93       *_mul = \&_mul_use_div;
  94       *_div = \&_div_use_div;
  95       }
  96     }
  97   return $BASE_LEN unless wantarray;
  98   return ($BASE_LEN, $BASE, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN, $MAX_VAL);
  99   }
 100
 101 sub _new
 102   {
 103   # (ref to string) return ref to num_array
 104   # Convert a number from string format (without sign) to internal base
 105   # 1ex format. Assumes normalized value as input.
 106   my $il = length($_[1])-1;
 107
 108   # < BASE_LEN due len-1 above
 109   return [ int($_[1]) ] if $il < $BASE_LEN;     # shortcut for short numbers
 110
 111   # this leaves '00000' instead of int 0 and will be corrected after any op
 112   [ reverse(unpack("a" . ($il % $BASE_LEN+1)
 113     . ("a$BASE_LEN" x ($il / $BASE_LEN)), $_[1])) ];
 114   }
 115
 116 BEGIN
 117   {
 118   # from Daniel Pfeiffer: determine largest group of digits that is precisely
 119   # multipliable with itself plus carry
 120   # Test now changed to expect the proper pattern, not a result off by 1 or 2
 121   my ($e, $num) = 3;    # lowest value we will use is 3+1-1 = 3
 122   do
 123     {
 124     $num = ('9' x ++$e) + 0;
 125     $num *= $num + 1.0;
 126     } while ("$num" =~ /9{$e}0{$e}/);   # must be a certain pattern
 127   $e--;                                 # last test failed, so retract one step
 128   # the limits below brush the problems with the test above under the rug:
 129   # the test should be able to find the proper $e automatically
 130   $e = 5 if $^O =~ /^uts/;      # UTS get's some special treatment
 131   $e = 5 if $^O =~ /^unicos/;   # unicos is also problematic (6 seems to work
 132                                 # there, but we play safe)
 133
 134   my $int = 0;
 135   if ($e > 7)
 136     {
 137     use integer;
 138     my $e1 = 7;
 139     $num = 7;
 140     do
 141       {
 142       $num = ('9' x ++$e1) + 0;
 143       $num *= $num + 1;
 144       } while ("$num" =~ /9{$e1}0{$e1}/);       # must be a certain pattern
 145     $e1--;                                      # last test failed, so retract one step
 146     if ($e1 > 7)
 147       {
 148       $int = 1; $e = $e1;
 149       }
 150     }
 151
 152   __PACKAGE__->_base_len($e,$int);      # set and store
 153
 154   use integer;
 155   # find out how many bits _and, _or and _xor can take (old default = 16)
 156   # I don't think anybody has yet 128 bit scalars, so let's play safe.
 157   local $^W = 0;        # don't warn about 'nonportable number'
 158   $AND_BITS = 15; $XOR_BITS = 15; $OR_BITS = 15;
 159
 160   # find max bits, we will not go higher than numberofbits that fit into $BASE
 161   # to make _and etc simpler (and faster for smaller, slower for large numbers)
 162   my $max = 16;
 163   while (2 ** $max < $BASE) { $max++; }
 164   {
 165     no integer;
 166     $max = 16 if $] < 5.006;    # older Perls might not take >16 too well
 167   }
 168   my ($x,$y,$z);
 169   do {
 170     $AND_BITS++;
 171     $x = CORE::oct('0b' . '1' x $AND_BITS); $y = $x & $x;
 172     $z = (2 ** $AND_BITS) - 1;
 173     } while ($AND_BITS < $max && $x == $z && $y == $x);
 174   $AND_BITS --;                                         # retreat one step
 175   do {
 176     $XOR_BITS++;
 177     $x = CORE::oct('0b' . '1' x $XOR_BITS); $y = $x ^ 0;
 178     $z = (2 ** $XOR_BITS) - 1;
 179     } while ($XOR_BITS < $max && $x == $z && $y == $x);
 180   $XOR_BITS --;                                         # retreat one step
 181   do {
 182     $OR_BITS++;
 183     $x = CORE::oct('0b' . '1' x $OR_BITS); $y = $x | $x;
 184     $z = (2 ** $OR_BITS) - 1;
 185     } while ($OR_BITS < $max && $x == $z && $y == $x);
 186   $OR_BITS --;                                          # retreat one step
 187
 188   $AND_MASK = __PACKAGE__->_new( ( 2 ** $AND_BITS ));
 189   $XOR_MASK = __PACKAGE__->_new( ( 2 ** $XOR_BITS ));
 190   $OR_MASK = __PACKAGE__->_new( ( 2 ** $OR_BITS ));
 191
 192   # We can compute the approximate length no faster than the real length:
 193   *_alen = \&_len;
 194   }
 195
 196 ###############################################################################
 197
 198 sub _zero
 199   {
 200   # create a zero
 201   [ 0 ];
 202   }
 203
 204 sub _one
 205   {
 206   # create a one
 207   [ 1 ];
 208   }
 209
 210 sub _two
 211   {
 212   # create a two (used internally for shifting)
 213   [ 2 ];
 214   }
 215
 216 sub _ten
 217   {
 218   # create a 10 (used internally for shifting)
 219   [ 10 ];
 220   }
 221
 222 sub _1ex
 223   {
 224   # create a 1Ex
 225   my $rem = $_[1] % $BASE_LEN;          # remainder
 226   my $parts = $_[1] / $BASE_LEN;        # parts
 227
 228   # 000000, 000000, 100
 229   [ (0) x $parts, '1' . ('0' x $rem) ];
 230   }
 231
 232 sub _copy
 233   {
 234   # make a true copy
 235   [ @{$_[1]} ];
 236   }
 237
 238 # catch and throw away
 239 sub import { }
 240
 241 ##############################################################################
 242 # convert back to string and number
 243
 244 sub _str
 245   {
 246   # (ref to BINT) return num_str
 247   # Convert number from internal base 100000 format to string format.
 248   # internal format is always normalized (no leading zeros, "-0" => "+0")
 249   my $ar = $_[1];
 250
 251   my $l = scalar @$ar;                          # number of parts
 252   if ($l < 1)                                   # should not happen
 253     {
 254     require Carp;
 255     Carp::croak("$_[1] has no elements");
 256     }
 257
 258   my $ret = "";
 259   # handle first one different to strip leading zeros from it (there are no
 260   # leading zero parts in internal representation)
 261   $l --; $ret .= int($ar->[$l]); $l--;
 262   # Interestingly, the pre-padd method uses more time
 263   # the old grep variant takes longer (14 vs. 10 sec)
 264   my $z = '0' x ($BASE_LEN-1);
 265   while ($l >= 0)
 266     {
 267     $ret .= substr($z.$ar->[$l],-$BASE_LEN); # fastest way I could think of
 268     $l--;
 269     }
 270   $ret;
 271   }
 272
 273 sub _num
 274   {
 275   # Make a number (scalar int/float) from a BigInt object
 276   my $x = $_[1];
 277
 278   return 0+$x->[0] if scalar @$x == 1;  # below $BASE
 279   my $fac = 1;
 280   my $num = 0;
 281   foreach (@$x)
 282     {
 283     $num += $fac*$_; $fac *= $BASE;
 284     }
 285   $num;
 286   }
 287
 288 ##############################################################################
 289 # actual math code
 290
 291 sub _add
 292   {
 293   # (ref to int_num_array, ref to int_num_array)
 294   # routine to add two base 1eX numbers
 295   # stolen from Knuth Vol 2 Algorithm A pg 231
 296   # there are separate routines to add and sub as per Knuth pg 233
 297   # This routine clobbers up array x, but not y.
 298
 299   my ($c,$x,$y) = @_;
 300
 301   return $x if (@$y == 1) && $y->[0] == 0;              # $x + 0 => $x
 302   if ((@$x == 1) && $x->[0] == 0)                       # 0 + $y => $y->copy
 303     {
 304     # twice as slow as $x = [ @$y ], but nec. to retain $x as ref :(
 305     @$x = @$y; return $x;
 306     }
 307
 308   # for each in Y, add Y to X and carry. If after that, something is left in
 309   # X, foreach in X add carry to X and then return X, carry
 310   # Trades one "$j++" for having to shift arrays
 311   my $i; my $car = 0; my $j = 0;
 312   for $i (@$y)
 313     {
 314     $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
 315     $j++;
 316     }
 317   while ($car != 0)
 318     {
 319     $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
 320     }
 321   $x;
 322   }
 323
 324 sub _inc
 325   {
 326   # (ref to int_num_array, ref to int_num_array)
 327   # Add 1 to $x, modify $x in place
 328   my ($c,$x) = @_;
 329
 330   for my $i (@$x)
 331     {
 332     return $x if (($i += 1) < $BASE);           # early out
 333     $i = 0;                                     # overflow, next
 334     }
 335   push @$x,1 if (($x->[-1] || 0) == 0);         # last overflowed, so extend
 336   $x;
 337   }
 338
 339 sub _dec
 340   {
 341   # (ref to int_num_array, ref to int_num_array)
 342   # Sub 1 from $x, modify $x in place
 343   my ($c,$x) = @_;
 344
 345   my $MAX = $BASE-1;                            # since MAX_VAL based on BASE
 346   for my $i (@$x)
 347     {
 348     last if (($i -= 1) >= 0);                   # early out
 349     $i = $MAX;                                  # underflow, next
 350     }
 351   pop @$x if $x->[-1] == 0 && @$x > 1;          # last underflowed (but leave 0)
 352   $x;
 353   }
 354
 355 sub _sub
 356   {
 357   # (ref to int_num_array, ref to int_num_array, swap)
 358   # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
 359   # subtract Y from X by modifying x in place
 360   my ($c,$sx,$sy,$s) = @_;
 361
 362   my $car = 0; my $i; my $j = 0;
 363   if (!$s)
 364     {
 365     for $i (@$sx)
 366       {
 367       last unless defined $sy->[$j] || $car;
 368       $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
 369       }
 370     # might leave leading zeros, so fix that
 371     return __strip_zeros($sx);
 372     }
 373   for $i (@$sx)
 374     {
 375     # we can't do an early out if $x is < than $y, since we
 376     # need to copy the high chunks from $y. Found by Bob Mathews.
 377     #last unless defined $sy->[$j] || $car;
 378     $sy->[$j] += $BASE
 379      if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
 380     $j++;
 381     }
 382   # might leave leading zeros, so fix that
 383   __strip_zeros($sy);
 384   }
 385
 386 sub _mul_use_mul
 387   {
 388   # (ref to int_num_array, ref to int_num_array)
 389   # multiply two numbers in internal representation
 390   # modifies first arg, second need not be different from first
 391   my ($c,$xv,$yv) = @_;
 392
 393   if (@$yv == 1)
 394     {
 395     # shortcut for two very short numbers (improved by Nathan Zook)
 396     # works also if xv and yv are the same reference, and handles also $x == 0
 397     if (@$xv == 1)
 398       {
 399       if (($xv->[0] *= $yv->[0]) >= $BASE)
 400          {
 401          $xv->[0] = $xv->[0] - ($xv->[1] = int($xv->[0] * $RBASE)) * $BASE;
 402          };
 403       return $xv;
 404       }
 405     # $x * 0 => 0
 406     if ($yv->[0] == 0)
 407       {
 408       @$xv = (0);
 409       return $xv;
 410       }
 411     # multiply a large number a by a single element one, so speed up
 412     my $y = $yv->[0]; my $car = 0;
 413     foreach my $i (@$xv)
 414       {
 415       $i = $i * $y + $car; $car = int($i * $RBASE); $i -= $car * $BASE;
 416       }
 417     push @$xv, $car if $car != 0;
 418     return $xv;
 419     }
 420   # shortcut for result $x == 0 => result = 0
 421   return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
 422
 423   # since multiplying $x with $x fails, make copy in this case
 424   $yv = [@$xv] if $xv == $yv;   # same references?
 425
 426   my @prod = (); my ($prod,$car,$cty,$xi,$yi);
 427
 428   for $xi (@$xv)
 429     {
 430     $car = 0; $cty = 0;
 431
 432     # slow variant
 433 #    for $yi (@$yv)
 434 #      {
 435 #      $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
 436 #      $prod[$cty++] =
 437 #       $prod - ($car = int($prod * RBASE)) * $BASE;  # see USE_MUL
 438 #      }
 439 #    $prod[$cty] += $car if $car; # need really to check for 0?
 440 #    $xi = shift @prod;
 441
 442     # faster variant
 443     # looping through this if $xi == 0 is silly - so optimize it away!
 444     $xi = (shift @prod || 0), next if $xi == 0;
 445     for $yi (@$yv)
 446       {
 447       $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
 448 ##     this is actually a tad slower
 449 ##        $prod = $prod[$cty]; $prod += ($car + $xi * $yi);     # no ||0 here
 450       $prod[$cty++] =
 451        $prod - ($car = int($prod * $RBASE)) * $BASE;  # see USE_MUL
 452       }
 453     $prod[$cty] += $car if $car; # need really to check for 0?
 454     $xi = shift @prod || 0;     # || 0 makes v5.005_3 happy
 455     }
 456   push @$xv, @prod;
 457   # can't have leading zeros
 458 #  __strip_zeros($xv);
 459   $xv;
 460   }
 461
 462 sub _mul_use_div_64
 463   {
 464   # (ref to int_num_array, ref to int_num_array)
 465   # multiply two numbers in internal representation
 466   # modifies first arg, second need not be different from first
 467   # works for 64 bit integer with "use integer"
 468   my ($c,$xv,$yv) = @_;
 469
 470   use integer;
 471   if (@$yv == 1)
 472     {
 473     # shortcut for two small numbers, also handles $x == 0
 474     if (@$xv == 1)
 475       {
 476       # shortcut for two very short numbers (improved by Nathan Zook)
 477       # works also if xv and yv are the same reference, and handles also $x == 0
 478       if (($xv->[0] *= $yv->[0]) >= $BASE)
 479           {
 480           $xv->[0] =
 481               $xv->[0] - ($xv->[1] = $xv->[0] / $BASE) * $BASE;
 482           };
 483       return $xv;
 484       }
 485     # $x * 0 => 0
 486     if ($yv->[0] == 0)
 487       {
 488       @$xv = (0);
 489       return $xv;
 490       }
 491     # multiply a large number a by a single element one, so speed up
 492     my $y = $yv->[0]; my $car = 0;
 493     foreach my $i (@$xv)
 494       {
 495       #$i = $i * $y + $car; $car = $i / $BASE; $i -= $car * $BASE;
 496       $i = $i * $y + $car; $i -= ($car = $i / $BASE) * $BASE;
 497       }
 498     push @$xv, $car if $car != 0;
 499     return $xv;
 500     }
 501   # shortcut for result $x == 0 => result = 0
 502   return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
 503
 504   # since multiplying $x with $x fails, make copy in this case
 505   $yv = [@$xv] if $xv == $yv;   # same references?
 506
 507   my @prod = (); my ($prod,$car,$cty,$xi,$yi);
 508   for $xi (@$xv)
 509     {
 510     $car = 0; $cty = 0;
 511     # looping through this if $xi == 0 is silly - so optimize it away!
 512     $xi = (shift @prod || 0), next if $xi == 0;
 513     for $yi (@$yv)
 514       {
 515       $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
 516       $prod[$cty++] = $prod - ($car = $prod / $BASE) * $BASE;
 517       }
 518     $prod[$cty] += $car if $car; # need really to check for 0?
 519     $xi = shift @prod || 0;     # || 0 makes v5.005_3 happy
 520     }
 521   push @$xv, @prod;
 522   $xv;
 523   }
 524
 525 sub _mul_use_div
 526   {
 527   # (ref to int_num_array, ref to int_num_array)
 528   # multiply two numbers in internal representation
 529   # modifies first arg, second need not be different from first
 530   my ($c,$xv,$yv) = @_;
 531
 532   if (@$yv == 1)
 533     {
 534     # shortcut for two small numbers, also handles $x == 0
 535     if (@$xv == 1)
 536       {
 537       # shortcut for two very short numbers (improved by Nathan Zook)
 538       # works also if xv and yv are the same reference, and handles also $x == 0
 539       if (($xv->[0] *= $yv->[0]) >= $BASE)
 540           {
 541           $xv->[0] =
 542               $xv->[0] - ($xv->[1] = int($xv->[0] / $BASE)) * $BASE;
 543           };
 544       return $xv;
 545       }
 546     # $x * 0 => 0
 547     if ($yv->[0] == 0)
 548       {
 549       @$xv = (0);
 550       return $xv;
 551       }
 552     # multiply a large number a by a single element one, so speed up
 553     my $y = $yv->[0]; my $car = 0;
 554     foreach my $i (@$xv)
 555       {
 556       $i = $i * $y + $car; $car = int($i / $BASE); $i -= $car * $BASE;
 557       # This (together with use integer;) does not work on 32-bit Perls
 558       #$i = $i * $y + $car; $i -= ($car = $i / $BASE) * $BASE;
 559       }
 560     push @$xv, $car if $car != 0;
 561     return $xv;
 562     }
 563   # shortcut for result $x == 0 => result = 0
 564   return $xv if ( ((@$xv == 1) && ($xv->[0] == 0)) );
 565
 566   # since multiplying $x with $x fails, make copy in this case
 567   $yv = [@$xv] if $xv == $yv;   # same references?
 568
 569   my @prod = (); my ($prod,$car,$cty,$xi,$yi);
 570   for $xi (@$xv)
 571     {
 572     $car = 0; $cty = 0;
 573     # looping through this if $xi == 0 is silly - so optimize it away!
 574     $xi = (shift @prod || 0), next if $xi == 0;
 575     for $yi (@$yv)
 576       {
 577       $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
 578       $prod[$cty++] = $prod - ($car = int($prod / $BASE)) * $BASE;
 579       }
 580     $prod[$cty] += $car if $car; # need really to check for 0?
 581     $xi = shift @prod || 0;     # || 0 makes v5.005_3 happy
 582     }
 583   push @$xv, @prod;
 584   # can't have leading zeros
 585 #  __strip_zeros($xv);
 586   $xv;
 587   }
 588
 589 sub _div_use_mul
 590   {
 591   # ref to array, ref to array, modify first array and return remainder if
 592   # in list context
 593
 594   # see comments in _div_use_div() for more explanations
 595
 596   my ($c,$x,$yorg) = @_;
 597
 598   # the general div algorithm here is about O(N*N) and thus quite slow, so
 599   # we first check for some special cases and use shortcuts to handle them.
 600
 601   # This works, because we store the numbers in a chunked format where each
 602   # element contains 5..7 digits (depending on system).
 603
 604   # if both numbers have only one element:
 605   if (@$x == 1 && @$yorg == 1)
 606     {
 607     # shortcut, $yorg and $x are two small numbers
 608     if (wantarray)
 609       {
 610       my $r = [ $x->[0] % $yorg->[0] ];
 611       $x->[0] = int($x->[0] / $yorg->[0]);
 612       return ($x,$r);
 613       }
 614     else
 615       {
 616       $x->[0] = int($x->[0] / $yorg->[0]);
 617       return $x;
 618       }
 619     }
 620
 621   # if x has more than one, but y has only one element:
 622   if (@$yorg == 1)
 623     {
 624     my $rem;
 625     $rem = _mod($c,[ @$x ],$yorg) if wantarray;
 626
 627     # shortcut, $y is < $BASE
 628     my $j = scalar @$x; my $r = 0;
 629     my $y = $yorg->[0]; my $b;
 630     while ($j-- > 0)
 631       {
 632       $b = $r * $BASE + $x->[$j];
 633       $x->[$j] = int($b/$y);
 634       $r = $b % $y;
 635       }
 636     pop @$x if @$x > 1 && $x->[-1] == 0;        # splice up a leading zero
 637     return ($x,$rem) if wantarray;
 638     return $x;
 639     }
 640
 641   # now x and y have more than one element
 642
 643   # check whether y has more elements than x, if yet, the result will be 0
 644   if (@$yorg > @$x)
 645     {
 646     my $rem;
 647     $rem = [@$x] if wantarray;                  # make copy
 648     splice (@$x,1);                             # keep ref to original array
 649     $x->[0] = 0;                                # set to 0
 650     return ($x,$rem) if wantarray;              # including remainder?
 651     return $x;                                  # only x, which is [0] now
 652     }
 653   # check whether the numbers have the same number of elements, in that case
 654   # the result will fit into one element and can be computed efficiently
 655   if (@$yorg == @$x)
 656     {
 657     my $rem;
 658     # if $yorg has more digits than $x (it's leading element is longer than
 659     # the one from $x), the result will also be 0:
 660     if (length(int($yorg->[-1])) > length(int($x->[-1])))
 661       {
 662       $rem = [@$x] if wantarray;                # make copy
 663       splice (@$x,1);                           # keep ref to org array
 664       $x->[0] = 0;                              # set to 0
 665       return ($x,$rem) if wantarray;            # including remainder?
 666       return $x;
 667       }
 668     # now calculate $x / $yorg
 669     if (length(int($yorg->[-1])) == length(int($x->[-1])))
 670       {
 671       # same length, so make full compare
 672
 673       my $a = 0; my $j = scalar @$x - 1;
 674       # manual way (abort if unequal, good for early ne)
 675       while ($j >= 0)
 676         {
 677         last if ($a = $x->[$j] - $yorg->[$j]); $j--;
 678         }
 679       # $a contains the result of the compare between X and Y
 680       # a < 0: x < y, a == 0: x == y, a > 0: x > y
 681       if ($a <= 0)
 682         {
 683         $rem = [ 0 ];                   # a = 0 => x == y => rem 0
 684         $rem = [@$x] if $a != 0;        # a < 0 => x < y => rem = x
 685         splice(@$x,1);                  # keep single element
 686         $x->[0] = 0;                    # if $a < 0
 687         $x->[0] = 1 if $a == 0;         # $x == $y
 688         return ($x,$rem) if wantarray;
 689         return $x;
 690         }
 691       # $x >= $y, so proceed normally
 692       }
 693     }
 694
 695   # all other cases:
 696
 697   my $y = [ @$yorg ];                           # always make copy to preserve
 698
 699   my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
 700
 701   $car = $bar = $prd = 0;
 702   if (($dd = int($BASE/($y->[-1]+1))) != 1)
 703     {
 704     for $xi (@$x)
 705       {
 706       $xi = $xi * $dd + $car;
 707       $xi -= ($car = int($xi * $RBASE)) * $BASE;        # see USE_MUL
 708       }
 709     push(@$x, $car); $car = 0;
 710     for $yi (@$y)
 711       {
 712       $yi = $yi * $dd + $car;
 713       $yi -= ($car = int($yi * $RBASE)) * $BASE;        # see USE_MUL
 714       }
 715     }
 716   else
 717     {
 718     push(@$x, 0);
 719     }
 720   @q = (); ($v2,$v1) = @$y[-2,-1];
 721   $v2 = 0 unless $v2;
 722   while ($#$x > $#$y)
 723     {
 724     ($u2,$u1,$u0) = @$x[-3..-1];
 725     $u2 = 0 unless $u2;
 726     #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
 727     # if $v1 == 0;
 728     $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
 729     --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
 730     if ($q)
 731       {
 732       ($car, $bar) = (0,0);
 733       for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
 734         {
 735         $prd = $q * $y->[$yi] + $car;
 736         $prd -= ($car = int($prd * $RBASE)) * $BASE;    # see USE_MUL
 737         $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
 738         }
 739       if ($x->[-1] < $car + $bar)
 740         {
 741         $car = 0; --$q;
 742         for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
 743           {
 744           $x->[$xi] -= $BASE
 745            if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
 746           }
 747         }
 748       }
 749     pop(@$x);
 750     unshift(@q, $q);
 751     }
 752   if (wantarray)
 753     {
 754     @d = ();
 755     if ($dd != 1)
 756       {
 757       $car = 0;
 758       for $xi (reverse @$x)
 759         {
 760         $prd = $car * $BASE + $xi;
 761         $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
 762         unshift(@d, $tmp);
 763         }
 764       }
 765     else
 766       {
 767       @d = @$x;
 768       }
 769     @$x = @q;
 770     my $d = \@d;
 771     __strip_zeros($x);
 772     __strip_zeros($d);
 773     return ($x,$d);
 774     }
 775   @$x = @q;
 776   __strip_zeros($x);
 777   $x;
 778   }
 779
 780 sub _div_use_div_64
 781   {
 782   # ref to array, ref to array, modify first array and return remainder if
 783   # in list context
 784   # This version works on 64 bit integers
 785   my ($c,$x,$yorg) = @_;
 786
 787   use integer;
 788   # the general div algorithm here is about O(N*N) and thus quite slow, so
 789   # we first check for some special cases and use shortcuts to handle them.
 790
 791   # This works, because we store the numbers in a chunked format where each
 792   # element contains 5..7 digits (depending on system).
 793
 794   # if both numbers have only one element:
 795   if (@$x == 1 && @$yorg == 1)
 796     {
 797     # shortcut, $yorg and $x are two small numbers
 798     if (wantarray)
 799       {
 800       my $r = [ $x->[0] % $yorg->[0] ];
 801       $x->[0] = int($x->[0] / $yorg->[0]);
 802       return ($x,$r);
 803       }
 804     else
 805       {
 806       $x->[0] = int($x->[0] / $yorg->[0]);
 807       return $x;
 808       }
 809     }
 810   # if x has more than one, but y has only one element:
 811   if (@$yorg == 1)
 812     {
 813     my $rem;
 814     $rem = _mod($c,[ @$x ],$yorg) if wantarray;
 815
 816     # shortcut, $y is < $BASE
 817     my $j = scalar @$x; my $r = 0;
 818     my $y = $yorg->[0]; my $b;
 819     while ($j-- > 0)
 820       {
 821       $b = $r * $BASE + $x->[$j];
 822       $x->[$j] = int($b/$y);
 823       $r = $b % $y;
 824       }
 825     pop @$x if @$x > 1 && $x->[-1] == 0;        # splice up a leading zero
 826     return ($x,$rem) if wantarray;
 827     return $x;
 828     }
 829   # now x and y have more than one element
 830
 831   # check whether y has more elements than x, if yet, the result will be 0
 832   if (@$yorg > @$x)
 833     {
 834     my $rem;
 835     $rem = [@$x] if wantarray;                  # make copy
 836     splice (@$x,1);                             # keep ref to original array
 837     $x->[0] = 0;                                # set to 0
 838     return ($x,$rem) if wantarray;              # including remainder?
 839     return $x;                                  # only x, which is [0] now
 840     }
 841   # check whether the numbers have the same number of elements, in that case
 842   # the result will fit into one element and can be computed efficiently
 843   if (@$yorg == @$x)
 844     {
 845     my $rem;
 846     # if $yorg has more digits than $x (it's leading element is longer than
 847     # the one from $x), the result will also be 0:
 848     if (length(int($yorg->[-1])) > length(int($x->[-1])))
 849       {
 850       $rem = [@$x] if wantarray;                # make copy
 851       splice (@$x,1);                           # keep ref to org array
 852       $x->[0] = 0;                              # set to 0
 853       return ($x,$rem) if wantarray;            # including remainder?
 854       return $x;
 855       }
 856     # now calculate $x / $yorg
 857
 858     if (length(int($yorg->[-1])) == length(int($x->[-1])))
 859       {
 860       # same length, so make full compare
 861
 862       my $a = 0; my $j = scalar @$x - 1;
 863       # manual way (abort if unequal, good for early ne)
 864       while ($j >= 0)
 865         {
 866         last if ($a = $x->[$j] - $yorg->[$j]); $j--;
 867         }
 868       # $a contains the result of the compare between X and Y
 869       # a < 0: x < y, a == 0: x == y, a > 0: x > y
 870       if ($a <= 0)
 871         {
 872         $rem = [ 0 ];                   # a = 0 => x == y => rem 0
 873         $rem = [@$x] if $a != 0;        # a < 0 => x < y => rem = x
 874         splice(@$x,1);                  # keep single element
 875         $x->[0] = 0;                    # if $a < 0
 876         $x->[0] = 1 if $a == 0;         # $x == $y
 877         return ($x,$rem) if wantarray;  # including remainder?
 878         return $x;
 879         }
 880       # $x >= $y, so proceed normally
 881
 882       }
 883     }
 884
 885   # all other cases:
 886
 887   my $y = [ @$yorg ];                           # always make copy to preserve
 888
 889   my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
 890
 891   $car = $bar = $prd = 0;
 892   if (($dd = int($BASE/($y->[-1]+1))) != 1)
 893     {
 894     for $xi (@$x)
 895       {
 896       $xi = $xi * $dd + $car;
 897       $xi -= ($car = int($xi / $BASE)) * $BASE;
 898       }
 899     push(@$x, $car); $car = 0;
 900     for $yi (@$y)
 901       {
 902       $yi = $yi * $dd + $car;
 903       $yi -= ($car = int($yi / $BASE)) * $BASE;
 904       }
 905     }
 906   else
 907     {
 908     push(@$x, 0);
 909     }
 910
 911   # @q will accumulate the final result, $q contains the current computed
 912   # part of the final result
 913
 914   @q = (); ($v2,$v1) = @$y[-2,-1];
 915   $v2 = 0 unless $v2;
 916   while ($#$x > $#$y)
 917     {
 918     ($u2,$u1,$u0) = @$x[-3..-1];
 919     $u2 = 0 unless $u2;
 920     #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
 921     # if $v1 == 0;
 922     $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
 923     --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
 924     if ($q)
 925       {
 926       ($car, $bar) = (0,0);
 927       for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
 928         {
 929         $prd = $q * $y->[$yi] + $car;
 930         $prd -= ($car = int($prd / $BASE)) * $BASE;
 931         $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
 932         }
 933       if ($x->[-1] < $car + $bar)
 934         {
 935         $car = 0; --$q;
 936         for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
 937           {
 938           $x->[$xi] -= $BASE
 939            if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
 940           }
 941         }
 942       }
 943     pop(@$x); unshift(@q, $q);
 944     }
 945   if (wantarray)
 946     {
 947     @d = ();
 948     if ($dd != 1)
 949       {
 950       $car = 0;
 951       for $xi (reverse @$x)
 952         {
 953         $prd = $car * $BASE + $xi;
 954         $car = $prd - ($tmp = int($prd / $dd)) * $dd;
 955         unshift(@d, $tmp);
 956         }
 957       }
 958     else
 959       {
 960       @d = @$x;
 961       }
 962     @$x = @q;
 963     my $d = \@d;
 964     __strip_zeros($x);
 965     __strip_zeros($d);
 966     return ($x,$d);
 967     }
 968   @$x = @q;
 969   __strip_zeros($x);
 970   $x;
 971   }
 972
 973 sub _div_use_div
 974   {
 975   # ref to array, ref to array, modify first array and return remainder if
 976   # in list context
 977   my ($c,$x,$yorg) = @_;
 978
 979   # the general div algorithm here is about O(N*N) and thus quite slow, so
 980   # we first check for some special cases and use shortcuts to handle them.
 981
 982   # This works, because we store the numbers in a chunked format where each
 983   # element contains 5..7 digits (depending on system).
 984
 985   # if both numbers have only one element:
 986   if (@$x == 1 && @$yorg == 1)
 987     {
 988     # shortcut, $yorg and $x are two small numbers
 989     if (wantarray)
 990       {
 991       my $r = [ $x->[0] % $yorg->[0] ];
 992       $x->[0] = int($x->[0] / $yorg->[0]);
 993       return ($x,$r);
 994       }
 995     else
 996       {
 997       $x->[0] = int($x->[0] / $yorg->[0]);
 998       return $x;
 999       }
1000     }
1001   # if x has more than one, but y has only one element:
1002   if (@$yorg == 1)
1003     {
1004     my $rem;
1005     $rem = _mod($c,[ @$x ],$yorg) if wantarray;
1006
1007     # shortcut, $y is < $BASE
1008     my $j = scalar @$x; my $r = 0;
1009     my $y = $yorg->[0]; my $b;
1010     while ($j-- > 0)
1011       {
1012       $b = $r * $BASE + $x->[$j];
1013       $x->[$j] = int($b/$y);
1014       $r = $b % $y;
1015       }
1016     pop @$x if @$x > 1 && $x->[-1] == 0;        # splice up a leading zero
1017     return ($x,$rem) if wantarray;
1018     return $x;
1019     }
1020   # now x and y have more than one element
1021
1022   # check whether y has more elements than x, if yet, the result will be 0
1023   if (@$yorg > @$x)
1024     {
1025     my $rem;
1026     $rem = [@$x] if wantarray;                  # make copy
1027     splice (@$x,1);                             # keep ref to original array
1028     $x->[0] = 0;                                # set to 0
1029     return ($x,$rem) if wantarray;              # including remainder?
1030     return $x;                                  # only x, which is [0] now
1031     }
1032   # check whether the numbers have the same number of elements, in that case
1033   # the result will fit into one element and can be computed efficiently
1034   if (@$yorg == @$x)
1035     {
1036     my $rem;
1037     # if $yorg has more digits than $x (it's leading element is longer than
1038     # the one from $x), the result will also be 0:
1039     if (length(int($yorg->[-1])) > length(int($x->[-1])))
1040       {
1041       $rem = [@$x] if wantarray;                # make copy
1042       splice (@$x,1);                           # keep ref to org array
1043       $x->[0] = 0;                              # set to 0
1044       return ($x,$rem) if wantarray;            # including remainder?
1045       return $x;
1046       }
1047     # now calculate $x / $yorg
1048
1049     if (length(int($yorg->[-1])) == length(int($x->[-1])))
1050       {
1051       # same length, so make full compare
1052
1053       my $a = 0; my $j = scalar @$x - 1;
1054       # manual way (abort if unequal, good for early ne)
1055       while ($j >= 0)
1056         {
1057         last if ($a = $x->[$j] - $yorg->[$j]); $j--;
1058         }
1059       # $a contains the result of the compare between X and Y
1060       # a < 0: x < y, a == 0: x == y, a > 0: x > y
1061       if ($a <= 0)
1062         {
1063         $rem = [ 0 ];                   # a = 0 => x == y => rem 0
1064         $rem = [@$x] if $a != 0;        # a < 0 => x < y => rem = x
1065         splice(@$x,1);                  # keep single element
1066         $x->[0] = 0;                    # if $a < 0
1067         $x->[0] = 1 if $a == 0;         # $x == $y
1068         return ($x,$rem) if wantarray;  # including remainder?
1069         return $x;
1070         }
1071       # $x >= $y, so proceed normally
1072
1073       }
1074     }
1075
1076   # all other cases:
1077
1078   my $y = [ @$yorg ];                           # always make copy to preserve
1079
1080   my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
1081
1082   $car = $bar = $prd = 0;
1083   if (($dd = int($BASE/($y->[-1]+1))) != 1)
1084     {
1085     for $xi (@$x)
1086       {
1087       $xi = $xi * $dd + $car;
1088       $xi -= ($car = int($xi / $BASE)) * $BASE;
1089       }
1090     push(@$x, $car); $car = 0;
1091     for $yi (@$y)
1092       {
1093       $yi = $yi * $dd + $car;
1094       $yi -= ($car = int($yi / $BASE)) * $BASE;
1095       }
1096     }
1097   else
1098     {
1099     push(@$x, 0);
1100     }
1101
1102   # @q will accumulate the final result, $q contains the current computed
1103   # part of the final result
1104
1105   @q = (); ($v2,$v1) = @$y[-2,-1];
1106   $v2 = 0 unless $v2;
1107   while ($#$x > $#$y)
1108     {
1109     ($u2,$u1,$u0) = @$x[-3..-1];
1110     $u2 = 0 unless $u2;
1111     #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
1112     # if $v1 == 0;
1113     $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$BASE+$u1)/$v1));
1114     --$q while ($v2*$q > ($u0*$BASE+$u1-$q*$v1)*$BASE+$u2);
1115     if ($q)
1116       {
1117       ($car, $bar) = (0,0);
1118       for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
1119         {
1120         $prd = $q * $y->[$yi] + $car;
1121         $prd -= ($car = int($prd / $BASE)) * $BASE;
1122         $x->[$xi] += $BASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
1123         }
1124       if ($x->[-1] < $car + $bar)
1125         {
1126         $car = 0; --$q;
1127         for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
1128           {
1129           $x->[$xi] -= $BASE
1130            if ($car = (($x->[$xi] += $y->[$yi] + $car) >= $BASE));
1131           }
1132         }
1133       }
1134     pop(@$x); unshift(@q, $q);
1135     }
1136   if (wantarray)
1137     {
1138     @d = ();
1139     if ($dd != 1)
1140       {
1141       $car = 0;
1142       for $xi (reverse @$x)
1143         {
1144         $prd = $car * $BASE + $xi;
1145         $car = $prd - ($tmp = int($prd / $dd)) * $dd;
1146         unshift(@d, $tmp);
1147         }
1148       }
1149     else
1150       {
1151       @d = @$x;
1152       }
1153     @$x = @q;
1154     my $d = \@d;
1155     __strip_zeros($x);
1156     __strip_zeros($d);
1157     return ($x,$d);
1158     }
1159   @$x = @q;
1160   __strip_zeros($x);
1161   $x;
1162   }
1163
1164 ##############################################################################
1165 # testing
1166
1167 sub _acmp
1168   {
1169   # internal absolute post-normalized compare (ignore signs)
1170   # ref to array, ref to array, return <0, 0, >0
1171   # arrays must have at least one entry; this is not checked for
1172   my ($c,$cx,$cy) = @_;
1173
1174   # shortcut for short numbers
1175   return (($cx->[0] <=> $cy->[0]) <=> 0)
1176    if scalar @$cx == scalar @$cy && scalar @$cx == 1;
1177
1178   # fast comp based on number of array elements (aka pseudo-length)
1179   my $lxy = (scalar @$cx - scalar @$cy)
1180   # or length of first element if same number of elements (aka difference 0)
1181     ||
1182   # need int() here because sometimes the last element is '00018' vs '18'
1183    (length(int($cx->[-1])) - length(int($cy->[-1])));
1184   return -1 if $lxy < 0;                                # already differs, ret
1185   return 1 if $lxy > 0;                                 # ditto
1186
1187   # manual way (abort if unequal, good for early ne)
1188   my $a; my $j = scalar @$cx;
1189   while (--$j >= 0)
1190     {
1191     last if ($a = $cx->[$j] - $cy->[$j]);
1192     }
1193   $a <=> 0;
1194   }
1195
1196 sub _len
1197   {
1198   # compute number of digits in base 10
1199
1200   # int() because add/sub sometimes leaves strings (like '00005') instead of
1201   # '5' in this place, thus causing length() to report wrong length
1202   my $cx = $_[1];
1203
1204   (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
1205   }
1206
1207 sub _digit
1208   {
1209   # Return the nth digit. Zero is rightmost, so _digit(123,0) gives 3.
1210   # Negative values count from the left, so _digit(123, -1) gives 1.
1211   my ($c,$x,$n) = @_;
1212
1213   my $len = _len('',$x);
1214
1215   $n += $len if $n < 0;                 # -1 last, -2 second-to-last
1216   return "0" if $n < 0 || $n >= $len;   # return 0 for digits out of range
1217
1218   my $elem = int($n / $BASE_LEN);       # which array element
1219   my $digit = $n % $BASE_LEN;           # which digit in this element
1220   substr("$x->[$elem]", -$digit-1, 1);
1221   }
1222
1223 sub _zeros
1224   {
1225   # return amount of trailing zeros in decimal
1226   # check each array elem in _m for having 0 at end as long as elem == 0
1227   # Upon finding a elem != 0, stop
1228   my $x = $_[1];
1229
1230   return 0 if scalar @$x == 1 && $x->[0] == 0;
1231
1232   my $zeros = 0; my $elem;
1233   foreach my $e (@$x)
1234     {
1235     if ($e != 0)
1236       {
1237       $elem = "$e";                             # preserve x
1238       $elem =~ s/.*?(0*$)/$1/;                  # strip anything not zero
1239       $zeros *= $BASE_LEN;                      # elems * 5
1240       $zeros += length($elem);                  # count trailing zeros
1241       last;                                     # early out
1242       }
1243     $zeros ++;                                  # real else branch: 50% slower!
1244     }
1245   $zeros;
1246   }
1247
1248 ##############################################################################
1249 # _is_* routines
1250
1251 sub _is_zero
1252   {
1253   # return true if arg is zero
1254   (((scalar @{$_[1]} == 1) && ($_[1]->[0] == 0))) <=> 0;
1255   }
1256
1257 sub _is_even
1258   {
1259   # return true if arg is even
1260   (!($_[1]->[0] & 1)) <=> 0;
1261   }
1262
1263 sub _is_odd
1264   {
1265   # return true if arg is odd
1266   (($_[1]->[0] & 1)) <=> 0;
1267   }
1268
1269 sub _is_one
1270   {
1271   # return true if arg is one
1272   (scalar @{$_[1]} == 1) && ($_[1]->[0] == 1) <=> 0;
1273   }
1274
1275 sub _is_two
1276   {
1277   # return true if arg is two
1278   (scalar @{$_[1]} == 1) && ($_[1]->[0] == 2) <=> 0;
1279   }
1280
1281 sub _is_ten
1282   {
1283   # return true if arg is ten
1284   (scalar @{$_[1]} == 1) && ($_[1]->[0] == 10) <=> 0;
1285   }
1286
1287 sub __strip_zeros
1288   {
1289   # internal normalization function that strips leading zeros from the array
1290   # args: ref to array
1291   my $s = shift;
1292
1293   my $cnt = scalar @$s; # get count of parts
1294   my $i = $cnt-1;
1295   push @$s,0 if $i < 0;         # div might return empty results, so fix it
1296
1297   return $s if @$s == 1;                # early out
1298
1299   #print "strip: cnt $cnt i $i\n";
1300   # '0', '3', '4', '0', '0',
1301   #  0    1    2    3    4
1302   # cnt = 5, i = 4
1303   # i = 4
1304   # i = 3
1305   # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
1306   # >= 1: skip first part (this can be zero)
1307   while ($i > 0) { last if $s->[$i] != 0; $i--; }
1308   $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
1309   $s;
1310   }
1311
1312 ###############################################################################
1313 # check routine to test internal state for corruptions
1314
1315 sub _check
1316   {
1317   # used by the test suite
1318   my $x = $_[1];
1319
1320   return "$x is not a reference" if !ref($x);
1321
1322   # are all parts are valid?
1323   my $i = 0; my $j = scalar @$x; my ($e,$try);
1324   while ($i < $j)
1325     {
1326     $e = $x->[$i]; $e = 'undef' unless defined $e;
1327     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
1328     last if $e !~ /^[+]?[0-9]+$/;
1329     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
1330     last if "$e" !~ /^[+]?[0-9]+$/;
1331     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
1332     last if '' . "$e" !~ /^[+]?[0-9]+$/;
1333     $try = ' < 0 || >= $BASE; '."($x, $e)";
1334     last if $e <0 || $e >= $BASE;
1335     # this test is disabled, since new/bnorm and certain ops (like early out
1336     # in add/sub) are allowed/expected to leave '00000' in some elements
1337     #$try = '=~ /^00+/; '."($x, $e)";
1338     #last if $e =~ /^00+/;
1339     $i++;
1340     }
1341   return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
1342   0;
1343   }
1344
1345
1346 ###############################################################################
1347
1348 sub _mod
1349   {
1350   # if possible, use mod shortcut
1351   my ($c,$x,$yo) = @_;
1352
1353   # slow way since $y to big
1354   if (scalar @$yo > 1)
1355     {
1356     my ($xo,$rem) = _div($c,$x,$yo);
1357     return $rem;
1358     }
1359
1360   my $y = $yo->[0];
1361   # both are single element arrays
1362   if (scalar @$x == 1)
1363     {
1364     $x->[0] %= $y;
1365     return $x;
1366     }
1367
1368   # @y is a single element, but @x has more than one element
1369   my $b = $BASE % $y;
1370   if ($b == 0)
1371     {
1372     # when BASE % Y == 0 then (B * BASE) % Y == 0
1373     # (B * BASE) % $y + A % Y => A % Y
1374     # so need to consider only last element: O(1)
1375     $x->[0] %= $y;
1376     }
1377   elsif ($b == 1)
1378     {
1379     # else need to go through all elements: O(N), but loop is a bit simplified
1380     my $r = 0;
1381     foreach (@$x)
1382       {
1383       $r = ($r + $_) % $y;              # not much faster, but heh...
1384       #$r += $_ % $y; $r %= $y;
1385       }
1386     $r = 0 if $r == $y;
1387     $x->[0] = $r;
1388     }
1389   else
1390     {
1391     # else need to go through all elements: O(N)
1392     my $r = 0; my $bm = 1;
1393     foreach (@$x)
1394       {
1395       $r = ($_ * $bm + $r) % $y;
1396       $bm = ($bm * $b) % $y;
1397
1398       #$r += ($_ % $y) * $bm;
1399       #$bm *= $b;
1400       #$bm %= $y;
1401       #$r %= $y;
1402       }
1403     $r = 0 if $r == $y;
1404     $x->[0] = $r;
1405     }
1406   splice (@$x,1);               # keep one element of $x
1407   $x;
1408   }
1409
1410 ##############################################################################
1411 # shifts
1412
1413 sub _rsft
1414   {
1415   my ($c,$x,$y,$n) = @_;
1416
1417   if ($n != 10)
1418     {
1419     $n = _new($c,$n); return _div($c,$x, _pow($c,$n,$y));
1420     }
1421
1422   # shortcut (faster) for shifting by 10)
1423   # multiples of $BASE_LEN
1424   my $dst = 0;                          # destination
1425   my $src = _num($c,$y);                # as normal int
1426   my $xlen = (@$x-1)*$BASE_LEN+length(int($x->[-1]));  # len of x in digits
1427   if ($src >= $xlen or ($src == $xlen and ! defined $x->[1]))
1428     {
1429     # 12345 67890 shifted right by more than 10 digits => 0
1430     splice (@$x,1);                    # leave only one element
1431     $x->[0] = 0;                       # set to zero
1432     return $x;
1433     }
1434   my $rem = $src % $BASE_LEN;           # remainder to shift
1435   $src = int($src / $BASE_LEN);         # source
1436   if ($rem == 0)
1437     {
1438     splice (@$x,0,$src);                # even faster, 38.4 => 39.3
1439     }
1440   else
1441     {
1442     my $len = scalar @$x - $src;        # elems to go
1443     my $vd; my $z = '0'x $BASE_LEN;
1444     $x->[scalar @$x] = 0;               # avoid || 0 test inside loop
1445     while ($dst < $len)
1446       {
1447       $vd = $z.$x->[$src];
1448       $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
1449       $src++;
1450       $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
1451       $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1452       $x->[$dst] = int($vd);
1453       $dst++;
1454       }
1455     splice (@$x,$dst) if $dst > 0;              # kill left-over array elems
1456     pop @$x if $x->[-1] == 0 && @$x > 1;        # kill last element if 0
1457     } # else rem == 0
1458   $x;
1459   }
1460
1461 sub _lsft
1462   {
1463   my ($c,$x,$y,$n) = @_;
1464
1465   if ($n != 10)
1466     {
1467     $n = _new($c,$n); return _mul($c,$x, _pow($c,$n,$y));
1468     }
1469
1470   # shortcut (faster) for shifting by 10) since we are in base 10eX
1471   # multiples of $BASE_LEN:
1472   my $src = scalar @$x;                 # source
1473   my $len = _num($c,$y);                # shift-len as normal int
1474   my $rem = $len % $BASE_LEN;           # remainder to shift
1475   my $dst = $src + int($len/$BASE_LEN); # destination
1476   my $vd;                               # further speedup
1477   $x->[$src] = 0;                       # avoid first ||0 for speed
1478   my $z = '0' x $BASE_LEN;
1479   while ($src >= 0)
1480     {
1481     $vd = $x->[$src]; $vd = $z.$vd;
1482     $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
1483     $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
1484     $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1485     $x->[$dst] = int($vd);
1486     $dst--; $src--;
1487     }
1488   # set lowest parts to 0
1489   while ($dst >= 0) { $x->[$dst--] = 0; }
1490   # fix spurios last zero element
1491   splice @$x,-1 if $x->[-1] == 0;
1492   $x;
1493   }
1494
1495 sub _pow
1496   {
1497   # power of $x to $y
1498   # ref to array, ref to array, return ref to array
1499   my ($c,$cx,$cy) = @_;
1500
1501   if (scalar @$cy == 1 && $cy->[0] == 0)
1502     {
1503     splice (@$cx,1); $cx->[0] = 1;              # y == 0 => x => 1
1504     return $cx;
1505     }
1506   if ((scalar @$cx == 1 && $cx->[0] == 1) ||    #    x == 1
1507       (scalar @$cy == 1 && $cy->[0] == 1))      # or y == 1
1508     {
1509     return $cx;
1510     }
1511   if (scalar @$cx == 1 && $cx->[0] == 0)
1512     {
1513     splice (@$cx,1); $cx->[0] = 0;              # 0 ** y => 0 (if not y <= 0)
1514     return $cx;
1515     }
1516
1517   my $pow2 = _one();
1518
1519   my $y_bin = _as_bin($c,$cy); $y_bin =~ s/^0b//;
1520   my $len = length($y_bin);
1521   while (--$len > 0)
1522     {
1523     _mul($c,$pow2,$cx) if substr($y_bin,$len,1) eq '1';         # is odd?
1524     _mul($c,$cx,$cx);
1525     }
1526
1527   _mul($c,$cx,$pow2);
1528   $cx;
1529   }
1530
1531 sub _nok
1532   {
1533   # n over k
1534   # ref to array, return ref to array
1535   my ($c,$n,$k) = @_;
1536
1537   # ( 7 )       7!       1*2*3*4 * 5*6*7   5 * 6 * 7       6   7
1538   # ( - ) = --------- =  --------------- = --------- = 5 * - * -
1539   # ( 3 )   (7-3)! 3!    1*2*3*4 * 1*2*3   1 * 2 * 3       2   3
1540
1541   if (!_is_zero($c,$k))
1542     {
1543     my $x = _copy($c,$n);
1544     _sub($c,$n,$k);
1545     _inc($c,$n);
1546     my $f = _copy($c,$n); _inc($c,$f);          # n = 5, f = 6, d = 2
1547     my $d = _two($c);
1548     while (_acmp($c,$f,$x) <= 0)                # f <= n ?
1549       {
1550       # n = (n * f / d) == 5 * 6 / 2
1551       $n = _mul($c,$n,$f); $n = _div($c,$n,$d);
1552       # f = 7, d = 3
1553       _inc($c,$f); _inc($c,$d);
1554       }
1555     }
1556   else
1557     {
1558     # keep ref to $n and set it to 1
1559     splice (@$n,1); $n->[0] = 1;
1560     }
1561   $n;
1562   }
1563
1564 my @factorials = (
1565   1,
1566   1,
1567   2,
1568   2*3,
1569   2*3*4,
1570   2*3*4*5,
1571   2*3*4*5*6,
1572   2*3*4*5*6*7,
1573 );
1574
1575 sub _fac
1576   {
1577   # factorial of $x
1578   # ref to array, return ref to array
1579   my ($c,$cx) = @_;
1580
1581   if ((@$cx == 1) && ($cx->[0] <= 7))
1582     {
1583     $cx->[0] = $factorials[$cx->[0]];           # 0 => 1, 1 => 1, 2 => 2 etc.
1584     return $cx;
1585     }
1586
1587   if ((@$cx == 1) &&            # we do this only if $x >= 12 and $x <= 7000
1588       ($cx->[0] >= 12 && $cx->[0] < 7000))
1589     {
1590
1591   # Calculate (k-j) * (k-j+1) ... k .. (k+j-1) * (k + j)
1592   # See http://blogten.blogspot.com/2007/01/calculating-n.html
1593   # The above series can be expressed as factors:
1594   #   k * k - (j - i) * 2
1595   # We cache k*k, and calculate (j * j) as the sum of the first j odd integers
1596
1597   # This will not work when N exceeds the storage of a Perl scalar, however,
1598   # in this case the algorithm would be way to slow to terminate, anyway.
1599
1600   # As soon as the last element of $cx is 0, we split it up and remember
1601   # how many zeors we got so far. The reason is that n! will accumulate
1602   # zeros at the end rather fast.
1603   my $zero_elements = 0;
1604
1605   # If n is even, set n = n -1
1606   my $k = _num($c,$cx); my $even = 1;
1607   if (($k & 1) == 0)
1608     {
1609     $even = $k; $k --;
1610     }
1611   # set k to the center point
1612   $k = ($k + 1) / 2;
1613 #  print "k $k even: $even\n";
1614   # now calculate k * k
1615   my $k2 = $k * $k;
1616   my $odd = 1; my $sum = 1;
1617   my $i = $k - 1;
1618   # keep reference to x
1619   my $new_x = _new($c, $k * $even);
1620   @$cx = @$new_x;
1621   if ($cx->[0] == 0)
1622     {
1623     $zero_elements ++; shift @$cx;
1624     }
1625 #  print STDERR "x = ", _str($c,$cx),"\n";
1626   my $BASE2 = int(sqrt($BASE))-1;
1627   my $j = 1;
1628   while ($j <= $i)
1629     {
1630     my $m = ($k2 - $sum); $odd += 2; $sum += $odd; $j++;
1631     while ($j <= $i && ($m < $BASE2) && (($k2 - $sum) < $BASE2))
1632       {
1633       $m *= ($k2 - $sum);
1634       $odd += 2; $sum += $odd; $j++;
1635 #      print STDERR "\n k2 $k2 m $m sum $sum odd $odd\n"; sleep(1);
1636       }
1637     if ($m < $BASE)
1638       {
1639       _mul($c,$cx,[$m]);
1640       }
1641     else
1642       {
1643       _mul($c,$cx,$c->_new($m));
1644       }
1645     if ($cx->[0] == 0)
1646       {
1647       $zero_elements ++; shift @$cx;
1648       }
1649 #    print STDERR "Calculate $k2 - $sum = $m (x = ", _str($c,$cx),")\n";
1650     }
1651   # multiply in the zeros again
1652   unshift @$cx, (0) x $zero_elements;
1653   return $cx;
1654   }
1655
1656   # go forward until $base is exceeded
1657   # limit is either $x steps (steps == 100 means a result always too high) or
1658   # $base.
1659   my $steps = 100; $steps = $cx->[0] if @$cx == 1;
1660   my $r = 2; my $cf = 3; my $step = 2; my $last = $r;
1661   while ($r*$cf < $BASE && $step < $steps)
1662     {
1663     $last = $r; $r *= $cf++; $step++;
1664     }
1665   if ((@$cx == 1) && $step == $cx->[0])
1666     {
1667     # completely done, so keep reference to $x and return
1668     $cx->[0] = $r;
1669     return $cx;
1670     }
1671
1672   # now we must do the left over steps
1673   my $n;                                        # steps still to do
1674   if (scalar @$cx == 1)
1675     {
1676     $n = $cx->[0];
1677     }
1678   else
1679     {
1680     $n = _copy($c,$cx);
1681     }
1682
1683   # Set $cx to the last result below $BASE (but keep ref to $x)
1684   $cx->[0] = $last; splice (@$cx,1);
1685   # As soon as the last element of $cx is 0, we split it up and remember
1686   # how many zeors we got so far. The reason is that n! will accumulate
1687   # zeros at the end rather fast.
1688   my $zero_elements = 0;
1689
1690   # do left-over steps fit into a scalar?
1691   if (ref $n eq 'ARRAY')
1692     {
1693     # No, so use slower inc() & cmp()
1694     # ($n is at least $BASE here)
1695     my $base_2 = int(sqrt($BASE)) - 1;
1696     #print STDERR "base_2: $base_2\n";
1697     while ($step < $base_2)
1698       {
1699       if ($cx->[0] == 0)
1700         {
1701         $zero_elements ++; shift @$cx;
1702         }
1703       my $b = $step * ($step + 1); $step += 2;
1704       _mul($c,$cx,[$b]);
1705       }
1706     $step = [$step];
1707     while (_acmp($c,$step,$n) <= 0)
1708       {
1709       if ($cx->[0] == 0)
1710         {
1711         $zero_elements ++; shift @$cx;
1712         }
1713       _mul($c,$cx,$step); _inc($c,$step);
1714       }
1715     }
1716   else
1717     {
1718     # Yes, so we can speed it up slightly
1719
1720 #    print "# left over steps $n\n";
1721
1722     my $base_4 = int(sqrt(sqrt($BASE))) - 2;
1723     #print STDERR "base_4: $base_4\n";
1724     my $n4 = $n - 4;
1725     while ($step < $n4 && $step < $base_4)
1726       {
1727       if ($cx->[0] == 0)
1728         {
1729         $zero_elements ++; shift @$cx;
1730         }
1731       my $b = $step * ($step + 1); $step += 2; $b *= $step * ($step + 1); $step += 2;
1732       _mul($c,$cx,[$b]);
1733       }
1734     my $base_2 = int(sqrt($BASE)) - 1;
1735     my $n2 = $n - 2;
1736     #print STDERR "base_2: $base_2\n";
1737     while ($step < $n2 && $step < $base_2)
1738       {
1739       if ($cx->[0] == 0)
1740         {
1741         $zero_elements ++; shift @$cx;
1742         }
1743       my $b = $step * ($step + 1); $step += 2;
1744       _mul($c,$cx,[$b]);
1745       }
1746     # do what's left over
1747     while ($step <= $n)
1748       {
1749       _mul($c,$cx,[$step]); $step++;
1750       if ($cx->[0] == 0)
1751         {
1752         $zero_elements ++; shift @$cx;
1753         }
1754       }
1755     }
1756   # multiply in the zeros again
1757   unshift @$cx, (0) x $zero_elements;
1758   $cx;                  # return result
1759   }
1760
1761 #############################################################################
1762
1763 sub _log_int
1764   {
1765   # calculate integer log of $x to base $base
1766   # ref to array, ref to array - return ref to array
1767   my ($c,$x,$base) = @_;
1768
1769   # X == 0 => NaN
1770   return if (scalar @$x == 1 && $x->[0] == 0);
1771   # BASE 0 or 1 => NaN
1772   return if (scalar @$base == 1 && $base->[0] < 2);
1773   my $cmp = _acmp($c,$x,$base); # X == BASE => 1
1774   if ($cmp == 0)
1775     {
1776     splice (@$x,1); $x->[0] = 1;
1777     return ($x,1)
1778     }
1779   # X < BASE
1780   if ($cmp < 0)
1781     {
1782     splice (@$x,1); $x->[0] = 0;
1783     return ($x,undef);
1784     }
1785
1786   my $x_org = _copy($c,$x);             # preserve x
1787   splice(@$x,1); $x->[0] = 1;           # keep ref to $x
1788
1789   # Compute a guess for the result based on:
1790   # $guess = int ( length_in_base_10(X) / ( log(base) / log(10) ) )
1791   my $len = _len($c,$x_org);
1792   my $log = log($base->[-1]) / log(10);
1793
1794   # for each additional element in $base, we add $BASE_LEN to the result,
1795   # based on the observation that log($BASE,10) is BASE_LEN and
1796   # log(x*y) == log(x) + log(y):
1797   $log += ((scalar @$base)-1) * $BASE_LEN;
1798
1799   # calculate now a guess based on the values obtained above:
1800   my $res = int($len / $log);
1801
1802   $x->[0] = $res;
1803   my $trial = _pow ($c, _copy($c, $base), $x);
1804   my $a = _acmp($c,$trial,$x_org);
1805
1806 #  print STDERR "# trial ", _str($c,$x)," was: $a (0 = exact, -1 too small, +1 too big)\n";
1807
1808   # found an exact result?
1809   return ($x,1) if $a == 0;
1810
1811   if ($a > 0)
1812     {
1813     # or too big
1814     _div($c,$trial,$base); _dec($c, $x);
1815     while (($a = _acmp($c,$trial,$x_org)) > 0)
1816       {
1817 #      print STDERR "# big _log_int at ", _str($c,$x), "\n";
1818       _div($c,$trial,$base); _dec($c, $x);
1819       }
1820     # result is now exact (a == 0), or too small (a < 0)
1821     return ($x, $a == 0 ? 1 : 0);
1822     }
1823
1824   # else: result was to small
1825   _mul($c,$trial,$base);
1826
1827   # did we now get the right result?
1828   $a = _acmp($c,$trial,$x_org);
1829
1830   if ($a == 0)                          # yes, exactly
1831     {
1832     _inc($c, $x);
1833     return ($x,1);
1834     }
1835   return ($x,0) if $a > 0;
1836
1837   # Result still too small (we should come here only if the estimate above
1838   # was very off base):
1839
1840   # Now let the normal trial run obtain the real result
1841   # Simple loop that increments $x by 2 in each step, possible overstepping
1842   # the real result
1843
1844   my $base_mul = _mul($c, _copy($c,$base), $base);      # $base * $base
1845
1846   while (($a = _acmp($c,$trial,$x_org)) < 0)
1847     {
1848 #    print STDERR "# small _log_int at ", _str($c,$x), "\n";
1849     _mul($c,$trial,$base_mul); _add($c, $x, [2]);
1850     }
1851
1852   my $exact = 1;
1853   if ($a > 0)
1854     {
1855     # overstepped the result
1856     _dec($c, $x);
1857     _div($c,$trial,$base);
1858     $a = _acmp($c,$trial,$x_org);
1859     if ($a > 0)
1860       {
1861       _dec($c, $x);
1862       }
1863     $exact = 0 if $a != 0;              # a = -1 => not exact result, a = 0 => exact
1864     }
1865
1866   ($x,$exact);                          # return result
1867   }
1868
1869 # for debugging:
1870   use constant DEBUG => 0;
1871   my $steps = 0;
1872   sub steps { $steps };
1873
1874 sub _sqrt
1875   {
1876   # square-root of $x in place
1877   # Compute a guess of the result (by rule of thumb), then improve it via
1878   # Newton's method.
1879   my ($c,$x) = @_;
1880
1881   if (scalar @$x == 1)
1882     {
1883     # fits into one Perl scalar, so result can be computed directly
1884     $x->[0] = int(sqrt($x->[0]));
1885     return $x;
1886     }
1887   my $y = _copy($c,$x);
1888   # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
1889   # since our guess will "grow"
1890   my $l = int((_len($c,$x)-1) / 2);
1891
1892   my $lastelem = $x->[-1];                                      # for guess
1893   my $elems = scalar @$x - 1;
1894   # not enough digits, but could have more?
1895   if ((length($lastelem) <= 3) && ($elems > 1))
1896     {
1897     # right-align with zero pad
1898     my $len = length($lastelem) & 1;
1899     print "$lastelem => " if DEBUG;
1900     $lastelem .= substr($x->[-2] . '0' x $BASE_LEN,0,$BASE_LEN);
1901     # former odd => make odd again, or former even to even again
1902     $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
1903     print "$lastelem\n" if DEBUG;
1904     }
1905
1906   # construct $x (instead of _lsft($c,$x,$l,10)
1907   my $r = $l % $BASE_LEN;       # 10000 00000 00000 00000 ($BASE_LEN=5)
1908   $l = int($l / $BASE_LEN);
1909   print "l =  $l " if DEBUG;
1910
1911   splice @$x,$l;                # keep ref($x), but modify it
1912
1913   # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
1914   # that gives us:
1915   # 14400 00000 => sqrt(14400) => guess first digits to be 120
1916   # 144000 000000 => sqrt(144000) => guess 379
1917
1918   print "$lastelem (elems $elems) => " if DEBUG;
1919   $lastelem = $lastelem / 10 if ($elems & 1 == 1);              # odd or even?
1920   my $g = sqrt($lastelem); $g =~ s/\.//;                        # 2.345 => 2345
1921   $r -= 1 if $elems & 1 == 0;                                   # 70 => 7
1922
1923   # padd with zeros if result is too short
1924   $x->[$l--] = int(substr($g . '0' x $r,0,$r+1));
1925   print "now ",$x->[-1] if DEBUG;
1926   print " would have been ", int('1' . '0' x $r),"\n" if DEBUG;
1927
1928   # If @$x > 1, we could compute the second elem of the guess, too, to create
1929   # an even better guess. Not implemented yet. Does it improve performance?
1930   $x->[$l--] = 0 while ($l >= 0);       # all other digits of guess are zero
1931
1932   print "start x= ",_str($c,$x),"\n" if DEBUG;
1933   my $two = _two();
1934   my $last = _zero();
1935   my $lastlast = _zero();
1936   $steps = 0 if DEBUG;
1937   while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
1938     {
1939     $steps++ if DEBUG;
1940     $lastlast = _copy($c,$last);
1941     $last = _copy($c,$x);
1942     _add($c,$x, _div($c,_copy($c,$y),$x));
1943     _div($c,$x, $two );
1944     print " x= ",_str($c,$x),"\n" if DEBUG;
1945     }
1946   print "\nsteps in sqrt: $steps, " if DEBUG;
1947   _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0;     # overshot?
1948   print " final ",$x->[-1],"\n" if DEBUG;
1949   $x;
1950   }
1951
1952 sub _root
1953   {
1954   # take n'th root of $x in place (n >= 3)
1955   my ($c,$x,$n) = @_;
1956
1957   if (scalar @$x == 1)
1958     {
1959     if (scalar @$n > 1)
1960       {
1961       # result will always be smaller than 2 so trunc to 1 at once
1962       $x->[0] = 1;
1963       }
1964     else
1965       {
1966       # fits into one Perl scalar, so result can be computed directly
1967       # cannot use int() here, because it rounds wrongly (try
1968       # (81 ** 3) ** (1/3) to see what I mean)
1969       #$x->[0] = int( $x->[0] ** (1 / $n->[0]) );
1970       # round to 8 digits, then truncate result to integer
1971       $x->[0] = int ( sprintf ("%.8f", $x->[0] ** (1 / $n->[0]) ) );
1972       }
1973     return $x;
1974     }
1975
1976   # we know now that X is more than one element long
1977
1978   # if $n is a power of two, we can repeatedly take sqrt($X) and find the
1979   # proper result, because sqrt(sqrt($x)) == root($x,4)
1980   my $b = _as_bin($c,$n);
1981   if ($b =~ /0b1(0+)$/)
1982     {
1983     my $count = CORE::length($1);       # 0b100 => len('00') => 2
1984     my $cnt = $count;                   # counter for loop
1985     unshift (@$x, 0);                   # add one element, together with one
1986                                         # more below in the loop this makes 2
1987     while ($cnt-- > 0)
1988       {
1989       # 'inflate' $X by adding one element, basically computing
1990       # $x * $BASE * $BASE. This gives us more $BASE_LEN digits for result
1991       # since len(sqrt($X)) approx == len($x) / 2.
1992       unshift (@$x, 0);
1993       # calculate sqrt($x), $x is now one element to big, again. In the next
1994       # round we make that two, again.
1995       _sqrt($c,$x);
1996       }
1997     # $x is now one element to big, so truncate result by removing it
1998     splice (@$x,0,1);
1999     }
2000   else
2001     {
2002     # trial computation by starting with 2,4,8,16 etc until we overstep
2003     my $step;
2004     my $trial = _two();
2005
2006     # while still to do more than X steps
2007     do
2008       {
2009       $step = _two();
2010       while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
2011         {
2012         _mul ($c, $step, [2]);
2013         _add ($c, $trial, $step);
2014         }
2015
2016       # hit exactly?
2017       if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) == 0)
2018         {
2019         @$x = @$trial;                  # make copy while preserving ref to $x
2020         return $x;
2021         }
2022       # overstepped, so go back on step
2023       _sub($c, $trial, $step);
2024       } while (scalar @$step > 1 || $step->[0] > 128);
2025
2026     # reset step to 2
2027     $step = _two();
2028     # add two, because $trial cannot be exactly the result (otherwise we would
2029     # already have found it)
2030     _add($c, $trial, $step);
2031
2032     # and now add more and more (2,4,6,8,10 etc)
2033     while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
2034       {
2035       _add ($c, $trial, $step);
2036       }
2037
2038     # hit not exactly? (overstepped)
2039     if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
2040       {
2041       _dec($c,$trial);
2042       }
2043
2044     # hit not exactly? (overstepped)
2045     # 80 too small, 81 slightly too big, 82 too big
2046     if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
2047       {
2048       _dec ($c, $trial);
2049       }
2050
2051     @$x = @$trial;                      # make copy while preserving ref to $x
2052     return $x;
2053     }
2054   $x;
2055   }
2056
2057 ##############################################################################
2058 # binary stuff
2059
2060 sub _and
2061   {
2062   my ($c,$x,$y) = @_;
2063
2064   # the shortcut makes equal, large numbers _really_ fast, and makes only a
2065   # very small performance drop for small numbers (e.g. something with less
2066   # than 32 bit) Since we optimize for large numbers, this is enabled.
2067   return $x if _acmp($c,$x,$y) == 0;            # shortcut
2068
2069   my $m = _one(); my ($xr,$yr);
2070   my $mask = $AND_MASK;
2071
2072   my $x1 = $x;
2073   my $y1 = _copy($c,$y);                        # make copy
2074   $x = _zero();
2075   my ($b,$xrr,$yrr);
2076   use integer;
2077   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2078     {
2079     ($x1, $xr) = _div($c,$x1,$mask);
2080     ($y1, $yr) = _div($c,$y1,$mask);
2081
2082     # make ints() from $xr, $yr
2083     # this is when the AND_BITS are greater than $BASE and is slower for
2084     # small (<256 bits) numbers, but faster for large numbers. Disabled
2085     # due to KISS principle
2086
2087 #    $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2088 #    $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2089 #    _add($c,$x, _mul($c, _new( $c, ($xrr & $yrr) ), $m) );
2090
2091     # 0+ due to '&' doesn't work in strings
2092     _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
2093     _mul($c,$m,$mask);
2094     }
2095   $x;
2096   }
2097
2098 sub _xor
2099   {
2100   my ($c,$x,$y) = @_;
2101
2102   return _zero() if _acmp($c,$x,$y) == 0;       # shortcut (see -and)
2103
2104   my $m = _one(); my ($xr,$yr);
2105   my $mask = $XOR_MASK;
2106
2107   my $x1 = $x;
2108   my $y1 = _copy($c,$y);                        # make copy
2109   $x = _zero();
2110   my ($b,$xrr,$yrr);
2111   use integer;
2112   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2113     {
2114     ($x1, $xr) = _div($c,$x1,$mask);
2115     ($y1, $yr) = _div($c,$y1,$mask);
2116     # make ints() from $xr, $yr (see _and())
2117     #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2118     #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2119     #_add($c,$x, _mul($c, _new( $c, ($xrr ^ $yrr) ), $m) );
2120
2121     # 0+ due to '^' doesn't work in strings
2122     _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
2123     _mul($c,$m,$mask);
2124     }
2125   # the loop stops when the shorter of the two numbers is exhausted
2126   # the remainder of the longer one will survive bit-by-bit, so we simple
2127   # multiply-add it in
2128   _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
2129   _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
2130
2131   $x;
2132   }
2133
2134 sub _or
2135   {
2136   my ($c,$x,$y) = @_;
2137
2138   return $x if _acmp($c,$x,$y) == 0;            # shortcut (see _and)
2139
2140   my $m = _one(); my ($xr,$yr);
2141   my $mask = $OR_MASK;
2142
2143   my $x1 = $x;
2144   my $y1 = _copy($c,$y);                        # make copy
2145   $x = _zero();
2146   my ($b,$xrr,$yrr);
2147   use integer;
2148   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2149     {
2150     ($x1, $xr) = _div($c,$x1,$mask);
2151     ($y1, $yr) = _div($c,$y1,$mask);
2152     # make ints() from $xr, $yr (see _and())
2153 #    $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2154 #    $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2155 #    _add($c,$x, _mul($c, _new( $c, ($xrr | $yrr) ), $m) );
2156
2157     # 0+ due to '|' doesn't work in strings
2158     _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
2159     _mul($c,$m,$mask);
2160     }
2161   # the loop stops when the shorter of the two numbers is exhausted
2162   # the remainder of the longer one will survive bit-by-bit, so we simple
2163   # multiply-add it in
2164   _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
2165   _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
2166
2167   $x;
2168   }
2169
2170 sub _as_hex
2171   {
2172   # convert a decimal number to hex (ref to array, return ref to string)
2173   my ($c,$x) = @_;
2174
2175   # fits into one element (handle also 0x0 case)
2176   return sprintf("0x%x",$x->[0]) if @$x == 1;
2177
2178   my $x1 = _copy($c,$x);
2179
2180   my $es = '';
2181   my ($xr, $h, $x10000);
2182   if ($] >= 5.006)
2183     {
2184     $x10000 = [ 0x10000 ]; $h = 'h4';
2185     }
2186   else
2187     {
2188     $x10000 = [ 0x1000 ]; $h = 'h3';
2189     }
2190   while (@$x1 != 1 || $x1->[0] != 0)            # _is_zero()
2191     {
2192     ($x1, $xr) = _div($c,$x1,$x10000);
2193     $es .= unpack($h,pack('V',$xr->[0]));
2194     }
2195   $es = reverse $es;
2196   $es =~ s/^[0]+//;   # strip leading zeros
2197   '0x' . $es;                                   # return result prepended with 0x
2198   }
2199
2200 sub _as_bin
2201   {
2202   # convert a decimal number to bin (ref to array, return ref to string)
2203   my ($c,$x) = @_;
2204
2205   # fits into one element (and Perl recent enough), handle also 0b0 case
2206   # handle zero case for older Perls
2207   if ($] <= 5.005 && @$x == 1 && $x->[0] == 0)
2208     {
2209     my $t = '0b0'; return $t;
2210     }
2211   if (@$x == 1 && $] >= 5.006)
2212     {
2213     my $t = sprintf("0b%b",$x->[0]);
2214     return $t;
2215     }
2216   my $x1 = _copy($c,$x);
2217
2218   my $es = '';
2219   my ($xr, $b, $x10000);
2220   if ($] >= 5.006)
2221     {
2222     $x10000 = [ 0x10000 ]; $b = 'b16';
2223     }
2224   else
2225     {
2226     $x10000 = [ 0x1000 ]; $b = 'b12';
2227     }
2228   while (!(@$x1 == 1 && $x1->[0] == 0))         # _is_zero()
2229     {
2230     ($x1, $xr) = _div($c,$x1,$x10000);
2231     $es .= unpack($b,pack('v',$xr->[0]));
2232     }
2233   $es = reverse $es;
2234   $es =~ s/^[0]+//;   # strip leading zeros
2235   '0b' . $es;                                   # return result prepended with 0b
2236   }
2237
2238 sub _as_oct
2239   {
2240   # convert a decimal number to octal (ref to array, return ref to string)
2241   my ($c,$x) = @_;
2242
2243   # fits into one element (handle also 0 case)
2244   return sprintf("0%o",$x->[0]) if @$x == 1;
2245
2246   my $x1 = _copy($c,$x);
2247
2248   my $es = '';
2249   my $xr;
2250   my $x1000 = [ 0100000 ];
2251   while (@$x1 != 1 || $x1->[0] != 0)            # _is_zero()
2252     {
2253     ($x1, $xr) = _div($c,$x1,$x1000);
2254     $es .= reverse sprintf("%05o", $xr->[0]);
2255     }
2256   $es = reverse $es;
2257   $es =~ s/^[0]+//;   # strip leading zeros
2258   '0' . $es;                                    # return result prepended with 0
2259   }
2260
2261 sub _from_oct
2262   {
2263   # convert a octal number to decimal (string, return ref to array)
2264   my ($c,$os) = @_;
2265
2266   # for older Perls, play safe
2267   my $m = [ 0100000 ];
2268   my $d = 5;                                    # 5 digits at a time
2269
2270   my $mul = _one();
2271   my $x = _zero();
2272
2273   my $len = int( (length($os)-1)/$d );          # $d digit parts, w/o the '0'
2274   my $val; my $i = -$d;
2275   while ($len >= 0)
2276     {
2277     $val = substr($os,$i,$d);                   # get oct digits
2278     $val = CORE::oct($val);
2279     $i -= $d; $len --;
2280     my $adder = [ $val ];
2281     _add ($c, $x, _mul ($c, $adder, $mul ) ) if $val != 0;
2282     _mul ($c, $mul, $m ) if $len >= 0;          # skip last mul
2283     }
2284   $x;
2285   }
2286
2287 sub _from_hex
2288   {
2289   # convert a hex number to decimal (string, return ref to array)
2290   my ($c,$hs) = @_;
2291
2292   my $m = _new($c, 0x10000000);                 # 28 bit at a time (<32 bit!)
2293   my $d = 7;                                    # 7 digits at a time
2294   if ($] <= 5.006)
2295     {
2296     # for older Perls, play safe
2297     $m = [ 0x10000 ];                           # 16 bit at a time (<32 bit!)
2298     $d = 4;                                     # 4 digits at a time
2299     }
2300
2301   my $mul = _one();
2302   my $x = _zero();
2303
2304   my $len = int( (length($hs)-2)/$d );          # $d digit parts, w/o the '0x'
2305   my $val; my $i = -$d;
2306   while ($len >= 0)
2307     {
2308     $val = substr($hs,$i,$d);                   # get hex digits
2309     $val =~ s/^0x// if $len == 0;               # for last part only because
2310     $val = CORE::hex($val);                     # hex does not like wrong chars
2311     $i -= $d; $len --;
2312     my $adder = [ $val ];
2313     # if the resulting number was to big to fit into one element, create a
2314     # two-element version (bug found by Mark Lakata - Thanx!)
2315     if (CORE::length($val) > $BASE_LEN)
2316       {
2317       $adder = _new($c,$val);
2318       }
2319     _add ($c, $x, _mul ($c, $adder, $mul ) ) if $val != 0;
2320     _mul ($c, $mul, $m ) if $len >= 0;          # skip last mul
2321     }
2322   $x;
2323   }
2324
2325 sub _from_bin
2326   {
2327   # convert a hex number to decimal (string, return ref to array)
2328   my ($c,$bs) = @_;
2329
2330   # instead of converting X (8) bit at a time, it is faster to "convert" the
2331   # number to hex, and then call _from_hex.
2332
2333   my $hs = $bs;
2334   $hs =~ s/^[+-]?0b//;                                  # remove sign and 0b
2335   my $l = length($hs);                                  # bits
2336   $hs = '0' x (8-($l % 8)) . $hs if ($l % 8) != 0;      # padd left side w/ 0
2337   my $h = '0x' . unpack('H*', pack ('B*', $hs));        # repack as hex
2338
2339   $c->_from_hex($h);
2340   }
2341
2342 ##############################################################################
2343 # special modulus functions
2344
2345 sub _modinv
2346   {
2347   # modular inverse
2348   my ($c,$x,$y) = @_;
2349
2350   my $u = _zero($c); my $u1 = _one($c);
2351   my $a = _copy($c,$y); my $b = _copy($c,$x);
2352
2353   # Euclid's Algorithm for bgcd(), only that we calc bgcd() ($a) and the
2354   # result ($u) at the same time. See comments in BigInt for why this works.
2355   my $q;
2356   ($a, $q, $b) = ($b, _div($c,$a,$b));          # step 1
2357   my $sign = 1;
2358   while (!_is_zero($c,$b))
2359     {
2360     my $t = _add($c,                            # step 2:
2361        _mul($c,_copy($c,$u1), $q) ,             #  t =  u1 * q
2362        $u );                                    #     + u
2363     $u = $u1;                                   #  u = u1, u1 = t
2364     $u1 = $t;
2365     $sign = -$sign;
2366     ($a, $q, $b) = ($b, _div($c,$a,$b));        # step 1
2367     }
2368
2369   # if the gcd is not 1, then return NaN
2370   return (undef,undef) unless _is_one($c,$a);
2371
2372   ($u1, $sign == 1 ? '+' : '-');
2373   }
2374
2375 sub _modpow
2376   {
2377   # modulus of power ($x ** $y) % $z
2378   my ($c,$num,$exp,$mod) = @_;
2379
2380   # a^b (mod 1) = 0 for all a and b
2381   if (_is_one($c,$mod))
2382     {
2383         @$num = 0;
2384         return $num;
2385     }
2386
2387   # 0^a (mod m) = 0 if m != 0, a != 0
2388   # 0^0 (mod m) = 1 if m != 0
2389   if (_is_one($c, $num)) {
2390       if (_is_zero($c, $exp)) {
2391           @$num = 1;
2392       } else {
2393           @$num = 0;
2394       }
2395       return $num;
2396   }
2397
2398 #  $num = _mod($c,$num,$mod);   # this does not make it faster
2399
2400   my $acc = _copy($c,$num); my $t = _one();
2401
2402   my $expbin = _as_bin($c,$exp); $expbin =~ s/^0b//;
2403   my $len = length($expbin);
2404   while (--$len >= 0)
2405     {
2406     if ( substr($expbin,$len,1) eq '1')                 # is_odd
2407       {
2408       _mul($c,$t,$acc);
2409       $t = _mod($c,$t,$mod);
2410       }
2411     _mul($c,$acc,$acc);
2412     $acc = _mod($c,$acc,$mod);
2413     }
2414   @$num = @$t;
2415   $num;
2416   }
2417
2418 sub _gcd
2419   {
2420   # greatest common divisor
2421   my ($c,$x,$y) = @_;
2422
2423   while ( (scalar @$y != 1) || ($y->[0] != 0) )         # while ($y != 0)
2424     {
2425     my $t = _copy($c,$y);
2426     $y = _mod($c, $x, $y);
2427     $x = $t;
2428     }
2429   $x;
2430   }
2431
2432 ##############################################################################
2433 ##############################################################################
2434
2435 1;
2436 __END__
2437
2438 =head1 NAME
2439
2440 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
2441
2442 =head1 SYNOPSIS
2443
2444 Provides support for big integer calculations. Not intended to be used by other
2445 modules. Other modules which sport the same functions can also be used to support
2446 Math::BigInt, like Math::BigInt::GMP or Math::BigInt::Pari.
2447
2448 =head1 DESCRIPTION
2449
2450 In order to allow for multiple big integer libraries, Math::BigInt was
2451 rewritten to use library modules for core math routines. Any module which
2452 follows the same API as this can be used instead by using the following:
2453
2454         use Math::BigInt lib => 'libname';
2455
2456 'libname' is either the long name ('Math::BigInt::Pari'), or only the short
2457 version like 'Pari'.
2458
2459 =head1 STORAGE
2460
2461 =head1 METHODS
2462
2463 The following functions MUST be defined in order to support the use by
2464 Math::BigInt v1.70 or later:
2465
2466         api_version()   return API version, 1 for v1.70, 2 for v1.83
2467         _new(string)    return ref to new object from ref to decimal string
2468         _zero()         return a new object with value 0
2469         _one()          return a new object with value 1
2470         _two()          return a new object with value 2
2471         _ten()          return a new object with value 10
2472
2473         _str(obj)       return ref to a string representing the object
2474         _num(obj)       returns a Perl integer/floating point number
2475                         NOTE: because of Perl numeric notation defaults,
2476                         the _num'ified obj may lose accuracy due to
2477                         machine-dependent floating point size limitations
2478
2479         _add(obj,obj)   Simple addition of two objects
2480         _mul(obj,obj)   Multiplication of two objects
2481         _div(obj,obj)   Division of the 1st object by the 2nd
2482                         In list context, returns (result,remainder).
2483                         NOTE: this is integer math, so no
2484                         fractional part will be returned.
2485                         The second operand will be not be 0, so no need to
2486                         check for that.
2487         _sub(obj,obj)   Simple subtraction of 1 object from another
2488                         a third, optional parameter indicates that the params
2489                         are swapped. In this case, the first param needs to
2490                         be preserved, while you can destroy the second.
2491                         sub (x,y,1) => return x - y and keep x intact!
2492         _dec(obj)       decrement object by one (input is guaranteed to be > 0)
2493         _inc(obj)       increment object by one
2494
2495
2496         _acmp(obj,obj)  <=> operator for objects (return -1, 0 or 1)
2497
2498         _len(obj)       returns count of the decimal digits of the object
2499         _digit(obj,n)   returns the n'th decimal digit of object
2500
2501         _is_one(obj)    return true if argument is 1
2502         _is_two(obj)    return true if argument is 2
2503         _is_ten(obj)    return true if argument is 10
2504         _is_zero(obj)   return true if argument is 0
2505         _is_even(obj)   return true if argument is even (0,2,4,6..)
2506         _is_odd(obj)    return true if argument is odd (1,3,5,7..)
2507
2508         _copy           return a ref to a true copy of the object
2509
2510         _check(obj)     check whether internal representation is still intact
2511                         return 0 for ok, otherwise error message as string
2512
2513         _from_hex(str)  return new object from a hexadecimal string
2514         _from_bin(str)  return new object from a binary string
2515         _from_oct(str)  return new object from an octal string
2516
2517         _as_hex(str)    return string containing the value as
2518                         unsigned hex string, with the '0x' prepended.
2519                         Leading zeros must be stripped.
2520         _as_bin(str)    Like as_hex, only as binary string containing only
2521                         zeros and ones. Leading zeros must be stripped and a
2522                         '0b' must be prepended.
2523
2524         _rsft(obj,N,B)  shift object in base B by N 'digits' right
2525         _lsft(obj,N,B)  shift object in base B by N 'digits' left
2526
2527         _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
2528                         Note: XOR, AND and OR pad with zeros if size mismatches
2529         _and(obj1,obj2) AND (bit-wise) object 1 with object 2
2530         _or(obj1,obj2)  OR (bit-wise) object 1 with object 2
2531
2532         _mod(obj1,obj2) Return remainder of div of the 1st by the 2nd object
2533         _sqrt(obj)      return the square root of object (truncated to int)
2534         _root(obj)      return the n'th (n >= 3) root of obj (truncated to int)
2535         _fac(obj)       return factorial of object 1 (1*2*3*4..)
2536         _pow(obj1,obj2) return object 1 to the power of object 2
2537                         return undef for NaN
2538         _zeros(obj)     return number of trailing decimal zeros
2539         _modinv         return inverse modulus
2540         _modpow         return modulus of power ($x ** $y) % $z
2541         _log_int(X,N)   calculate integer log() of X in base N
2542                         X >= 0, N >= 0 (return undef for NaN)
2543                         returns (RESULT, EXACT) where EXACT is:
2544                          1     : result is exactly RESULT
2545                          0     : result was truncated to RESULT
2546                          undef : unknown whether result is exactly RESULT
2547         _gcd(obj,obj)   return Greatest Common Divisor of two objects
2548
2549 The following functions are REQUIRED for an api_version of 2 or greater:
2550
2551         _1ex($x)        create the number 1Ex where x >= 0
2552         _alen(obj)      returns approximate count of the decimal digits of the
2553                         object. This estimate MUST always be greater or equal
2554                         to what _len() returns.
2555         _nok(n,k)       calculate n over k (binomial coefficient)
2556
2557 The following functions are optional, and can be defined if the underlying lib
2558 has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
2559 slow) fallback routines to emulate these:
2560
2561         _signed_or
2562         _signed_and
2563         _signed_xor
2564
2565 Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
2566 or '0b1101').
2567
2568 So the library needs only to deal with unsigned big integers. Testing of input
2569 parameter validity is done by the caller, so you need not worry about
2570 underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by zero or similar
2571 cases.
2572
2573 The first parameter can be modified, that includes the possibility that you
2574 return a reference to a completely different object instead. Although keeping
2575 the reference and just changing its contents is preferred over creating and
2576 returning a different reference.
2577
2578 Return values are always references to objects, strings, or true/false for
2579 comparison routines.
2580
2581 =head1 WRAP YOUR OWN
2582
2583 If you want to port your own favourite c-lib for big numbers to the
2584 Math::BigInt interface, you can take any of the already existing modules as
2585 a rough guideline. You should really wrap up the latest BigInt and BigFloat
2586 testsuites with your module, and replace in them any of the following:
2587
2588         use Math::BigInt;
2589
2590 by this:
2591
2592         use Math::BigInt lib => 'yourlib';
2593
2594 This way you ensure that your library really works 100% within Math::BigInt.
2595
2596 =head1 LICENSE
2597
2598 This program is free software; you may redistribute it and/or modify it under
2599 the same terms as Perl itself.
2600
2601 =head1 AUTHORS
2602
2603 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
2604 in late 2000.
2605 Separated from BigInt and shaped API with the help of John Peacock.
2606
2607 Fixed, speed-up, streamlined and enhanced by Tels 2001 - 2007.
2608
2609 =head1 SEE ALSO
2610
2611 L<Math::BigInt>, L<Math::BigFloat>,
2612 L<Math::BigInt::GMP>, L<Math::BigInt::FastCalc> and L<Math::BigInt::Pari>.
2613
2614 =cut