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

   1 package Math::BigInt::Calc;
   2
   3 use 5.006;
   4 use strict;
   5 # use warnings; # dont use warnings for older Perls
   6
   7 our $VERSION = '0.52';
   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, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN, $MAX_VAL, $BASE);
  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, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN, $MAX_VAL, $BASE);
  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 lenght 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 algorithmn 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 algorithmn 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 algorithmn 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, negative values count backward
1210   # zero is rightmost, so _digit(123,0) will give 3
1211   my ($c,$x,$n) = @_;
1212
1213   my $len = _len('',$x);
1214
1215   $n = $len+$n if $n < 0;               # -1 last, -2 second-to-last
1216   $n = abs($n);                         # if negative was too big
1217   $len--; $n = $len if $n > $len;       # n to big?
1218
1219   my $elem = int($n / $BASE_LEN);       # which array element
1220   my $digit = $n % $BASE_LEN;           # which digit in this element
1221   $elem = '0' x $BASE_LEN . @$x[$elem]; # get element padded with 0's
1222   substr($elem,-$digit-1,1);
1223   }
1224
1225 sub _zeros
1226   {
1227   # return amount of trailing zeros in decimal
1228   # check each array elem in _m for having 0 at end as long as elem == 0
1229   # Upon finding a elem != 0, stop
1230   my $x = $_[1];
1231
1232   return 0 if scalar @$x == 1 && $x->[0] == 0;
1233
1234   my $zeros = 0; my $elem;
1235   foreach my $e (@$x)
1236     {
1237     if ($e != 0)
1238       {
1239       $elem = "$e";                             # preserve x
1240       $elem =~ s/.*?(0*$)/$1/;                  # strip anything not zero
1241       $zeros *= $BASE_LEN;                      # elems * 5
1242       $zeros += length($elem);                  # count trailing zeros
1243       last;                                     # early out
1244       }
1245     $zeros ++;                                  # real else branch: 50% slower!
1246     }
1247   $zeros;
1248   }
1249
1250 ##############################################################################
1251 # _is_* routines
1252
1253 sub _is_zero
1254   {
1255   # return true if arg is zero
1256   (((scalar @{$_[1]} == 1) && ($_[1]->[0] == 0))) <=> 0;
1257   }
1258
1259 sub _is_even
1260   {
1261   # return true if arg is even
1262   (!($_[1]->[0] & 1)) <=> 0;
1263   }
1264
1265 sub _is_odd
1266   {
1267   # return true if arg is even
1268   (($_[1]->[0] & 1)) <=> 0;
1269   }
1270
1271 sub _is_one
1272   {
1273   # return true if arg is one
1274   (scalar @{$_[1]} == 1) && ($_[1]->[0] == 1) <=> 0;
1275   }
1276
1277 sub _is_two
1278   {
1279   # return true if arg is two
1280   (scalar @{$_[1]} == 1) && ($_[1]->[0] == 2) <=> 0;
1281   }
1282
1283 sub _is_ten
1284   {
1285   # return true if arg is ten
1286   (scalar @{$_[1]} == 1) && ($_[1]->[0] == 10) <=> 0;
1287   }
1288
1289 sub __strip_zeros
1290   {
1291   # internal normalization function that strips leading zeros from the array
1292   # args: ref to array
1293   my $s = shift;
1294
1295   my $cnt = scalar @$s; # get count of parts
1296   my $i = $cnt-1;
1297   push @$s,0 if $i < 0;         # div might return empty results, so fix it
1298
1299   return $s if @$s == 1;                # early out
1300
1301   #print "strip: cnt $cnt i $i\n";
1302   # '0', '3', '4', '0', '0',
1303   #  0    1    2    3    4
1304   # cnt = 5, i = 4
1305   # i = 4
1306   # i = 3
1307   # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
1308   # >= 1: skip first part (this can be zero)
1309   while ($i > 0) { last if $s->[$i] != 0; $i--; }
1310   $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
1311   $s;
1312   }
1313
1314 ###############################################################################
1315 # check routine to test internal state for corruptions
1316
1317 sub _check
1318   {
1319   # used by the test suite
1320   my $x = $_[1];
1321
1322   return "$x is not a reference" if !ref($x);
1323
1324   # are all parts are valid?
1325   my $i = 0; my $j = scalar @$x; my ($e,$try);
1326   while ($i < $j)
1327     {
1328     $e = $x->[$i]; $e = 'undef' unless defined $e;
1329     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
1330     last if $e !~ /^[+]?[0-9]+$/;
1331     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
1332     last if "$e" !~ /^[+]?[0-9]+$/;
1333     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
1334     last if '' . "$e" !~ /^[+]?[0-9]+$/;
1335     $try = ' < 0 || >= $BASE; '."($x, $e)";
1336     last if $e <0 || $e >= $BASE;
1337     # this test is disabled, since new/bnorm and certain ops (like early out
1338     # in add/sub) are allowed/expected to leave '00000' in some elements
1339     #$try = '=~ /^00+/; '."($x, $e)";
1340     #last if $e =~ /^00+/;
1341     $i++;
1342     }
1343   return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
1344   0;
1345   }
1346
1347
1348 ###############################################################################
1349
1350 sub _mod
1351   {
1352   # if possible, use mod shortcut
1353   my ($c,$x,$yo) = @_;
1354
1355   # slow way since $y to big
1356   if (scalar @$yo > 1)
1357     {
1358     my ($xo,$rem) = _div($c,$x,$yo);
1359     return $rem;
1360     }
1361
1362   my $y = $yo->[0];
1363   # both are single element arrays
1364   if (scalar @$x == 1)
1365     {
1366     $x->[0] %= $y;
1367     return $x;
1368     }
1369
1370   # @y is a single element, but @x has more than one element
1371   my $b = $BASE % $y;
1372   if ($b == 0)
1373     {
1374     # when BASE % Y == 0 then (B * BASE) % Y == 0
1375     # (B * BASE) % $y + A % Y => A % Y
1376     # so need to consider only last element: O(1)
1377     $x->[0] %= $y;
1378     }
1379   elsif ($b == 1)
1380     {
1381     # else need to go through all elements: O(N), but loop is a bit simplified
1382     my $r = 0;
1383     foreach (@$x)
1384       {
1385       $r = ($r + $_) % $y;              # not much faster, but heh...
1386       #$r += $_ % $y; $r %= $y;
1387       }
1388     $r = 0 if $r == $y;
1389     $x->[0] = $r;
1390     }
1391   else
1392     {
1393     # else need to go through all elements: O(N)
1394     my $r = 0; my $bm = 1;
1395     foreach (@$x)
1396       {
1397       $r = ($_ * $bm + $r) % $y;
1398       $bm = ($bm * $b) % $y;
1399
1400       #$r += ($_ % $y) * $bm;
1401       #$bm *= $b;
1402       #$bm %= $y;
1403       #$r %= $y;
1404       }
1405     $r = 0 if $r == $y;
1406     $x->[0] = $r;
1407     }
1408   splice (@$x,1);               # keep one element of $x
1409   $x;
1410   }
1411
1412 ##############################################################################
1413 # shifts
1414
1415 sub _rsft
1416   {
1417   my ($c,$x,$y,$n) = @_;
1418
1419   if ($n != 10)
1420     {
1421     $n = _new($c,$n); return _div($c,$x, _pow($c,$n,$y));
1422     }
1423
1424   # shortcut (faster) for shifting by 10)
1425   # multiples of $BASE_LEN
1426   my $dst = 0;                          # destination
1427   my $src = _num($c,$y);                # as normal int
1428   my $xlen = (@$x-1)*$BASE_LEN+length(int($x->[-1]));  # len of x in digits
1429   if ($src >= $xlen or ($src == $xlen and ! defined $x->[1]))
1430     {
1431     # 12345 67890 shifted right by more than 10 digits => 0
1432     splice (@$x,1);                    # leave only one element
1433     $x->[0] = 0;                       # set to zero
1434     return $x;
1435     }
1436   my $rem = $src % $BASE_LEN;           # remainder to shift
1437   $src = int($src / $BASE_LEN);         # source
1438   if ($rem == 0)
1439     {
1440     splice (@$x,0,$src);                # even faster, 38.4 => 39.3
1441     }
1442   else
1443     {
1444     my $len = scalar @$x - $src;        # elems to go
1445     my $vd; my $z = '0'x $BASE_LEN;
1446     $x->[scalar @$x] = 0;               # avoid || 0 test inside loop
1447     while ($dst < $len)
1448       {
1449       $vd = $z.$x->[$src];
1450       $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
1451       $src++;
1452       $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
1453       $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1454       $x->[$dst] = int($vd);
1455       $dst++;
1456       }
1457     splice (@$x,$dst) if $dst > 0;              # kill left-over array elems
1458     pop @$x if $x->[-1] == 0 && @$x > 1;        # kill last element if 0
1459     } # else rem == 0
1460   $x;
1461   }
1462
1463 sub _lsft
1464   {
1465   my ($c,$x,$y,$n) = @_;
1466
1467   if ($n != 10)
1468     {
1469     $n = _new($c,$n); return _mul($c,$x, _pow($c,$n,$y));
1470     }
1471
1472   # shortcut (faster) for shifting by 10) since we are in base 10eX
1473   # multiples of $BASE_LEN:
1474   my $src = scalar @$x;                 # source
1475   my $len = _num($c,$y);                # shift-len as normal int
1476   my $rem = $len % $BASE_LEN;           # remainder to shift
1477   my $dst = $src + int($len/$BASE_LEN); # destination
1478   my $vd;                               # further speedup
1479   $x->[$src] = 0;                       # avoid first ||0 for speed
1480   my $z = '0' x $BASE_LEN;
1481   while ($src >= 0)
1482     {
1483     $vd = $x->[$src]; $vd = $z.$vd;
1484     $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
1485     $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
1486     $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1487     $x->[$dst] = int($vd);
1488     $dst--; $src--;
1489     }
1490   # set lowest parts to 0
1491   while ($dst >= 0) { $x->[$dst--] = 0; }
1492   # fix spurios last zero element
1493   splice @$x,-1 if $x->[-1] == 0;
1494   $x;
1495   }
1496
1497 sub _pow
1498   {
1499   # power of $x to $y
1500   # ref to array, ref to array, return ref to array
1501   my ($c,$cx,$cy) = @_;
1502
1503   if (scalar @$cy == 1 && $cy->[0] == 0)
1504     {
1505     splice (@$cx,1); $cx->[0] = 1;              # y == 0 => x => 1
1506     return $cx;
1507     }
1508   if ((scalar @$cx == 1 && $cx->[0] == 1) ||    #    x == 1
1509       (scalar @$cy == 1 && $cy->[0] == 1))      # or y == 1
1510     {
1511     return $cx;
1512     }
1513   if (scalar @$cx == 1 && $cx->[0] == 0)
1514     {
1515     splice (@$cx,1); $cx->[0] = 0;              # 0 ** y => 0 (if not y <= 0)
1516     return $cx;
1517     }
1518
1519   my $pow2 = _one();
1520
1521   my $y_bin = _as_bin($c,$cy); $y_bin =~ s/^0b//;
1522   my $len = length($y_bin);
1523   while (--$len > 0)
1524     {
1525     _mul($c,$pow2,$cx) if substr($y_bin,$len,1) eq '1';         # is odd?
1526     _mul($c,$cx,$cx);
1527     }
1528
1529   _mul($c,$cx,$pow2);
1530   $cx;
1531   }
1532
1533 sub _nok
1534   {
1535   # n over k
1536   # ref to array, return ref to array
1537   my ($c,$n,$k) = @_;
1538
1539   # ( 7 )    7!          7*6*5 * 4*3*2*1   7 * 6 * 5
1540   # ( - ) = --------- =  --------------- = ---------
1541   # ( 3 )   3! (7-3)!    3*2*1 * 4*3*2*1   3 * 2 * 1
1542
1543   # compute n - k + 2 (so we start with 5 in the example above)
1544   my $x = _copy($c,$n);
1545
1546   _sub($c,$n,$k);
1547   if (!_is_one($c,$n))
1548     {
1549     _inc($c,$n);
1550     my $f = _copy($c,$n); _inc($c,$f);          # n = 5, f = 6, d = 2
1551     my $d = _two($c);
1552     while (_acmp($c,$f,$x) <= 0)                # f < n ?
1553       {
1554       # n = (n * f / d) == 5 * 6 / 2 => n == 3
1555       $n = _mul($c,$n,$f); $n = _div($c,$n,$d);
1556       # f = 7, d = 3
1557       _inc($c,$f); _inc($c,$d);
1558       }
1559     }
1560   else
1561     {
1562     # keep ref to $n and set it to 1
1563     splice (@$n,1); $n->[0] = 1;
1564     }
1565   $n;
1566   }
1567
1568 my @factorials = (
1569   1,
1570   1,
1571   2,
1572   2*3,
1573   2*3*4,
1574   2*3*4*5,
1575   2*3*4*5*6,
1576   2*3*4*5*6*7,
1577 );
1578
1579 sub _fac
1580   {
1581   # factorial of $x
1582   # ref to array, return ref to array
1583   my ($c,$cx) = @_;
1584
1585   if ((@$cx == 1) && ($cx->[0] <= 7))
1586     {
1587     $cx->[0] = $factorials[$cx->[0]];           # 0 => 1, 1 => 1, 2 => 2 etc.
1588     return $cx;
1589     }
1590
1591   if ((@$cx == 1) &&            # we do this only if $x >= 12 and $x <= 7000
1592       ($cx->[0] >= 12 && $cx->[0] < 7000))
1593     {
1594
1595   # Calculate (k-j) * (k-j+1) ... k .. (k+j-1) * (k + j)
1596   # See http://blogten.blogspot.com/2007/01/calculating-n.html
1597   # The above series can be expressed as factors:
1598   #   k * k - (j - i) * 2
1599   # We cache k*k, and calculate (j * j) as the sum of the first j odd integers
1600
1601   # This will not work when N exceeds the storage of a Perl scalar, however,
1602   # in this case the algorithm would be way to slow to terminate, anyway.
1603
1604   # As soon as the last element of $cx is 0, we split it up and remember
1605   # how many zeors we got so far. The reason is that n! will accumulate
1606   # zeros at the end rather fast.
1607   my $zero_elements = 0;
1608
1609   # If n is even, set n = n -1
1610   my $k = _num($c,$cx); my $even = 1;
1611   if (($k & 1) == 0)
1612     {
1613     $even = $k; $k --;
1614     }
1615   # set k to the center point
1616   $k = ($k + 1) / 2;
1617 #  print "k $k even: $even\n";
1618   # now calculate k * k
1619   my $k2 = $k * $k;
1620   my $odd = 1; my $sum = 1;
1621   my $i = $k - 1;
1622   # keep reference to x
1623   my $new_x = _new($c, $k * $even);
1624   @$cx = @$new_x;
1625   if ($cx->[0] == 0)
1626     {
1627     $zero_elements ++; shift @$cx;
1628     }
1629 #  print STDERR "x = ", _str($c,$cx),"\n";
1630   my $BASE2 = int(sqrt($BASE))-1;
1631   my $j = 1;
1632   while ($j <= $i)
1633     {
1634     my $m = ($k2 - $sum); $odd += 2; $sum += $odd; $j++;
1635     while ($j <= $i && ($m < $BASE2) && (($k2 - $sum) < $BASE2))
1636       {
1637       $m *= ($k2 - $sum);
1638       $odd += 2; $sum += $odd; $j++;
1639 #      print STDERR "\n k2 $k2 m $m sum $sum odd $odd\n"; sleep(1);
1640       }
1641     if ($m < $BASE)
1642       {
1643       _mul($c,$cx,[$m]);
1644       }
1645     else
1646       {
1647       _mul($c,$cx,$c->_new($m));
1648       }
1649     if ($cx->[0] == 0)
1650       {
1651       $zero_elements ++; shift @$cx;
1652       }
1653 #    print STDERR "Calculate $k2 - $sum = $m (x = ", _str($c,$cx),")\n";
1654     }
1655   # multiply in the zeros again
1656   unshift @$cx, (0) x $zero_elements;
1657   return $cx;
1658   }
1659
1660   # go forward until $base is exceeded
1661   # limit is either $x steps (steps == 100 means a result always too high) or
1662   # $base.
1663   my $steps = 100; $steps = $cx->[0] if @$cx == 1;
1664   my $r = 2; my $cf = 3; my $step = 2; my $last = $r;
1665   while ($r*$cf < $BASE && $step < $steps)
1666     {
1667     $last = $r; $r *= $cf++; $step++;
1668     }
1669   if ((@$cx == 1) && $step == $cx->[0])
1670     {
1671     # completely done, so keep reference to $x and return
1672     $cx->[0] = $r;
1673     return $cx;
1674     }
1675
1676   # now we must do the left over steps
1677   my $n;                                        # steps still to do
1678   if (scalar @$cx == 1)
1679     {
1680     $n = $cx->[0];
1681     }
1682   else
1683     {
1684     $n = _copy($c,$cx);
1685     }
1686
1687   # Set $cx to the last result below $BASE (but keep ref to $x)
1688   $cx->[0] = $last; splice (@$cx,1);
1689   # As soon as the last element of $cx is 0, we split it up and remember
1690   # how many zeors we got so far. The reason is that n! will accumulate
1691   # zeros at the end rather fast.
1692   my $zero_elements = 0;
1693
1694   # do left-over steps fit into a scalar?
1695   if (ref $n eq 'ARRAY')
1696     {
1697     # No, so use slower inc() & cmp()
1698     # ($n is at least $BASE here)
1699     my $base_2 = int(sqrt($BASE)) - 1;
1700     #print STDERR "base_2: $base_2\n";
1701     while ($step < $base_2)
1702       {
1703       if ($cx->[0] == 0)
1704         {
1705         $zero_elements ++; shift @$cx;
1706         }
1707       my $b = $step * ($step + 1); $step += 2;
1708       _mul($c,$cx,[$b]);
1709       }
1710     $step = [$step];
1711     while (_acmp($c,$step,$n) <= 0)
1712       {
1713       if ($cx->[0] == 0)
1714         {
1715         $zero_elements ++; shift @$cx;
1716         }
1717       _mul($c,$cx,$step); _inc($c,$step);
1718       }
1719     }
1720   else
1721     {
1722     # Yes, so we can speed it up slightly
1723
1724 #    print "# left over steps $n\n";
1725
1726     my $base_4 = int(sqrt(sqrt($BASE))) - 2;
1727     #print STDERR "base_4: $base_4\n";
1728     my $n4 = $n - 4;
1729     while ($step < $n4 && $step < $base_4)
1730       {
1731       if ($cx->[0] == 0)
1732         {
1733         $zero_elements ++; shift @$cx;
1734         }
1735       my $b = $step * ($step + 1); $step += 2; $b *= $step * ($step + 1); $step += 2;
1736       _mul($c,$cx,[$b]);
1737       }
1738     my $base_2 = int(sqrt($BASE)) - 1;
1739     my $n2 = $n - 2;
1740     #print STDERR "base_2: $base_2\n";
1741     while ($step < $n2 && $step < $base_2)
1742       {
1743       if ($cx->[0] == 0)
1744         {
1745         $zero_elements ++; shift @$cx;
1746         }
1747       my $b = $step * ($step + 1); $step += 2;
1748       _mul($c,$cx,[$b]);
1749       }
1750     # do what's left over
1751     while ($step <= $n)
1752       {
1753       _mul($c,$cx,[$step]); $step++;
1754       if ($cx->[0] == 0)
1755         {
1756         $zero_elements ++; shift @$cx;
1757         }
1758       }
1759     }
1760   # multiply in the zeros again
1761   unshift @$cx, (0) x $zero_elements;
1762   $cx;                  # return result
1763   }
1764
1765 #############################################################################
1766
1767 sub _log_int
1768   {
1769   # calculate integer log of $x to base $base
1770   # ref to array, ref to array - return ref to array
1771   my ($c,$x,$base) = @_;
1772
1773   # X == 0 => NaN
1774   return if (scalar @$x == 1 && $x->[0] == 0);
1775   # BASE 0 or 1 => NaN
1776   return if (scalar @$base == 1 && $base->[0] < 2);
1777   my $cmp = _acmp($c,$x,$base); # X == BASE => 1
1778   if ($cmp == 0)
1779     {
1780     splice (@$x,1); $x->[0] = 1;
1781     return ($x,1)
1782     }
1783   # X < BASE
1784   if ($cmp < 0)
1785     {
1786     splice (@$x,1); $x->[0] = 0;
1787     return ($x,undef);
1788     }
1789
1790   my $x_org = _copy($c,$x);             # preserve x
1791   splice(@$x,1); $x->[0] = 1;           # keep ref to $x
1792
1793   # Compute a guess for the result based on:
1794   # $guess = int ( length_in_base_10(X) / ( log(base) / log(10) ) )
1795   my $len = _len($c,$x_org);
1796   my $log = log($base->[-1]) / log(10);
1797
1798   # for each additional element in $base, we add $BASE_LEN to the result,
1799   # based on the observation that log($BASE,10) is BASE_LEN and
1800   # log(x*y) == log(x) + log(y):
1801   $log += ((scalar @$base)-1) * $BASE_LEN;
1802
1803   # calculate now a guess based on the values obtained above:
1804   my $res = int($len / $log);
1805
1806   $x->[0] = $res;
1807   my $trial = _pow ($c, _copy($c, $base), $x);
1808   my $a = _acmp($c,$trial,$x_org);
1809
1810 #  print STDERR "# trial ", _str($c,$x)," was: $a (0 = exact, -1 too small, +1 too big)\n";
1811
1812   # found an exact result?
1813   return ($x,1) if $a == 0;
1814
1815   if ($a > 0)
1816     {
1817     # or too big
1818     _div($c,$trial,$base); _dec($c, $x);
1819     while (($a = _acmp($c,$trial,$x_org)) > 0)
1820       {
1821 #      print STDERR "# big _log_int at ", _str($c,$x), "\n";
1822       _div($c,$trial,$base); _dec($c, $x);
1823       }
1824     # result is now exact (a == 0), or too small (a < 0)
1825     return ($x, $a == 0 ? 1 : 0);
1826     }
1827
1828   # else: result was to small
1829   _mul($c,$trial,$base);
1830
1831   # did we now get the right result?
1832   $a = _acmp($c,$trial,$x_org);
1833
1834   if ($a == 0)                          # yes, exactly
1835     {
1836     _inc($c, $x);
1837     return ($x,1);
1838     }
1839   return ($x,0) if $a > 0;
1840
1841   # Result still too small (we should come here only if the estimate above
1842   # was very off base):
1843
1844   # Now let the normal trial run obtain the real result
1845   # Simple loop that increments $x by 2 in each step, possible overstepping
1846   # the real result
1847
1848   my $base_mul = _mul($c, _copy($c,$base), $base);      # $base * $base
1849
1850   while (($a = _acmp($c,$trial,$x_org)) < 0)
1851     {
1852 #    print STDERR "# small _log_int at ", _str($c,$x), "\n";
1853     _mul($c,$trial,$base_mul); _add($c, $x, [2]);
1854     }
1855
1856   my $exact = 1;
1857   if ($a > 0)
1858     {
1859     # overstepped the result
1860     _dec($c, $x);
1861     _div($c,$trial,$base);
1862     $a = _acmp($c,$trial,$x_org);
1863     if ($a > 0)
1864       {
1865       _dec($c, $x);
1866       }
1867     $exact = 0 if $a != 0;              # a = -1 => not exact result, a = 0 => exact
1868     }
1869
1870   ($x,$exact);                          # return result
1871   }
1872
1873 # for debugging:
1874   use constant DEBUG => 0;
1875   my $steps = 0;
1876   sub steps { $steps };
1877
1878 sub _sqrt
1879   {
1880   # square-root of $x in place
1881   # Compute a guess of the result (by rule of thumb), then improve it via
1882   # Newton's method.
1883   my ($c,$x) = @_;
1884
1885   if (scalar @$x == 1)
1886     {
1887     # fits into one Perl scalar, so result can be computed directly
1888     $x->[0] = int(sqrt($x->[0]));
1889     return $x;
1890     }
1891   my $y = _copy($c,$x);
1892   # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
1893   # since our guess will "grow"
1894   my $l = int((_len($c,$x)-1) / 2);
1895
1896   my $lastelem = $x->[-1];                                      # for guess
1897   my $elems = scalar @$x - 1;
1898   # not enough digits, but could have more?
1899   if ((length($lastelem) <= 3) && ($elems > 1))
1900     {
1901     # right-align with zero pad
1902     my $len = length($lastelem) & 1;
1903     print "$lastelem => " if DEBUG;
1904     $lastelem .= substr($x->[-2] . '0' x $BASE_LEN,0,$BASE_LEN);
1905     # former odd => make odd again, or former even to even again
1906     $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
1907     print "$lastelem\n" if DEBUG;
1908     }
1909
1910   # construct $x (instead of _lsft($c,$x,$l,10)
1911   my $r = $l % $BASE_LEN;       # 10000 00000 00000 00000 ($BASE_LEN=5)
1912   $l = int($l / $BASE_LEN);
1913   print "l =  $l " if DEBUG;
1914
1915   splice @$x,$l;                # keep ref($x), but modify it
1916
1917   # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
1918   # that gives us:
1919   # 14400 00000 => sqrt(14400) => guess first digits to be 120
1920   # 144000 000000 => sqrt(144000) => guess 379
1921
1922   print "$lastelem (elems $elems) => " if DEBUG;
1923   $lastelem = $lastelem / 10 if ($elems & 1 == 1);              # odd or even?
1924   my $g = sqrt($lastelem); $g =~ s/\.//;                        # 2.345 => 2345
1925   $r -= 1 if $elems & 1 == 0;                                   # 70 => 7
1926
1927   # padd with zeros if result is too short
1928   $x->[$l--] = int(substr($g . '0' x $r,0,$r+1));
1929   print "now ",$x->[-1] if DEBUG;
1930   print " would have been ", int('1' . '0' x $r),"\n" if DEBUG;
1931
1932   # If @$x > 1, we could compute the second elem of the guess, too, to create
1933   # an even better guess. Not implemented yet. Does it improve performance?
1934   $x->[$l--] = 0 while ($l >= 0);       # all other digits of guess are zero
1935
1936   print "start x= ",_str($c,$x),"\n" if DEBUG;
1937   my $two = _two();
1938   my $last = _zero();
1939   my $lastlast = _zero();
1940   $steps = 0 if DEBUG;
1941   while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
1942     {
1943     $steps++ if DEBUG;
1944     $lastlast = _copy($c,$last);
1945     $last = _copy($c,$x);
1946     _add($c,$x, _div($c,_copy($c,$y),$x));
1947     _div($c,$x, $two );
1948     print " x= ",_str($c,$x),"\n" if DEBUG;
1949     }
1950   print "\nsteps in sqrt: $steps, " if DEBUG;
1951   _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0;     # overshot?
1952   print " final ",$x->[-1],"\n" if DEBUG;
1953   $x;
1954   }
1955
1956 sub _root
1957   {
1958   # take n'th root of $x in place (n >= 3)
1959   my ($c,$x,$n) = @_;
1960
1961   if (scalar @$x == 1)
1962     {
1963     if (scalar @$n > 1)
1964       {
1965       # result will always be smaller than 2 so trunc to 1 at once
1966       $x->[0] = 1;
1967       }
1968     else
1969       {
1970       # fits into one Perl scalar, so result can be computed directly
1971       # cannot use int() here, because it rounds wrongly (try
1972       # (81 ** 3) ** (1/3) to see what I mean)
1973       #$x->[0] = int( $x->[0] ** (1 / $n->[0]) );
1974       # round to 8 digits, then truncate result to integer
1975       $x->[0] = int ( sprintf ("%.8f", $x->[0] ** (1 / $n->[0]) ) );
1976       }
1977     return $x;
1978     }
1979
1980   # we know now that X is more than one element long
1981
1982   # if $n is a power of two, we can repeatedly take sqrt($X) and find the
1983   # proper result, because sqrt(sqrt($x)) == root($x,4)
1984   my $b = _as_bin($c,$n);
1985   if ($b =~ /0b1(0+)$/)
1986     {
1987     my $count = CORE::length($1);       # 0b100 => len('00') => 2
1988     my $cnt = $count;                   # counter for loop
1989     unshift (@$x, 0);                   # add one element, together with one
1990                                         # more below in the loop this makes 2
1991     while ($cnt-- > 0)
1992       {
1993       # 'inflate' $X by adding one element, basically computing
1994       # $x * $BASE * $BASE. This gives us more $BASE_LEN digits for result
1995       # since len(sqrt($X)) approx == len($x) / 2.
1996       unshift (@$x, 0);
1997       # calculate sqrt($x), $x is now one element to big, again. In the next
1998       # round we make that two, again.
1999       _sqrt($c,$x);
2000       }
2001     # $x is now one element to big, so truncate result by removing it
2002     splice (@$x,0,1);
2003     }
2004   else
2005     {
2006     # trial computation by starting with 2,4,8,16 etc until we overstep
2007     my $step;
2008     my $trial = _two();
2009
2010     # while still to do more than X steps
2011     do
2012       {
2013       $step = _two();
2014       while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
2015         {
2016         _mul ($c, $step, [2]);
2017         _add ($c, $trial, $step);
2018         }
2019
2020       # hit exactly?
2021       if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) == 0)
2022         {
2023         @$x = @$trial;                  # make copy while preserving ref to $x
2024         return $x;
2025         }
2026       # overstepped, so go back on step
2027       _sub($c, $trial, $step);
2028       } while (scalar @$step > 1 || $step->[0] > 128);
2029
2030     # reset step to 2
2031     $step = _two();
2032     # add two, because $trial cannot be exactly the result (otherwise we would
2033     # alrady have found it)
2034     _add($c, $trial, $step);
2035
2036     # and now add more and more (2,4,6,8,10 etc)
2037     while (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) < 0)
2038       {
2039       _add ($c, $trial, $step);
2040       }
2041
2042     # hit not exactly? (overstepped)
2043     if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
2044       {
2045       _dec($c,$trial);
2046       }
2047
2048     # hit not exactly? (overstepped)
2049     # 80 too small, 81 slightly too big, 82 too big
2050     if (_acmp($c, _pow($c, _copy($c, $trial), $n), $x) > 0)
2051       {
2052       _dec ($c, $trial);
2053       }
2054
2055     @$x = @$trial;                      # make copy while preserving ref to $x
2056     return $x;
2057     }
2058   $x;
2059   }
2060
2061 ##############################################################################
2062 # binary stuff
2063
2064 sub _and
2065   {
2066   my ($c,$x,$y) = @_;
2067
2068   # the shortcut makes equal, large numbers _really_ fast, and makes only a
2069   # very small performance drop for small numbers (e.g. something with less
2070   # than 32 bit) Since we optimize for large numbers, this is enabled.
2071   return $x if _acmp($c,$x,$y) == 0;            # shortcut
2072
2073   my $m = _one(); my ($xr,$yr);
2074   my $mask = $AND_MASK;
2075
2076   my $x1 = $x;
2077   my $y1 = _copy($c,$y);                        # make copy
2078   $x = _zero();
2079   my ($b,$xrr,$yrr);
2080   use integer;
2081   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2082     {
2083     ($x1, $xr) = _div($c,$x1,$mask);
2084     ($y1, $yr) = _div($c,$y1,$mask);
2085
2086     # make ints() from $xr, $yr
2087     # this is when the AND_BITS are greater than $BASE and is slower for
2088     # small (<256 bits) numbers, but faster for large numbers. Disabled
2089     # due to KISS principle
2090
2091 #    $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2092 #    $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2093 #    _add($c,$x, _mul($c, _new( $c, ($xrr & $yrr) ), $m) );
2094
2095     # 0+ due to '&' doesn't work in strings
2096     _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
2097     _mul($c,$m,$mask);
2098     }
2099   $x;
2100   }
2101
2102 sub _xor
2103   {
2104   my ($c,$x,$y) = @_;
2105
2106   return _zero() if _acmp($c,$x,$y) == 0;       # shortcut (see -and)
2107
2108   my $m = _one(); my ($xr,$yr);
2109   my $mask = $XOR_MASK;
2110
2111   my $x1 = $x;
2112   my $y1 = _copy($c,$y);                        # make copy
2113   $x = _zero();
2114   my ($b,$xrr,$yrr);
2115   use integer;
2116   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2117     {
2118     ($x1, $xr) = _div($c,$x1,$mask);
2119     ($y1, $yr) = _div($c,$y1,$mask);
2120     # make ints() from $xr, $yr (see _and())
2121     #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2122     #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2123     #_add($c,$x, _mul($c, _new( $c, ($xrr ^ $yrr) ), $m) );
2124
2125     # 0+ due to '^' doesn't work in strings
2126     _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
2127     _mul($c,$m,$mask);
2128     }
2129   # the loop stops when the shorter of the two numbers is exhausted
2130   # the remainder of the longer one will survive bit-by-bit, so we simple
2131   # multiply-add it in
2132   _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
2133   _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
2134
2135   $x;
2136   }
2137
2138 sub _or
2139   {
2140   my ($c,$x,$y) = @_;
2141
2142   return $x if _acmp($c,$x,$y) == 0;            # shortcut (see _and)
2143
2144   my $m = _one(); my ($xr,$yr);
2145   my $mask = $OR_MASK;
2146
2147   my $x1 = $x;
2148   my $y1 = _copy($c,$y);                        # make copy
2149   $x = _zero();
2150   my ($b,$xrr,$yrr);
2151   use integer;
2152   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
2153     {
2154     ($x1, $xr) = _div($c,$x1,$mask);
2155     ($y1, $yr) = _div($c,$y1,$mask);
2156     # make ints() from $xr, $yr (see _and())
2157 #    $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
2158 #    $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
2159 #    _add($c,$x, _mul($c, _new( $c, ($xrr | $yrr) ), $m) );
2160
2161     # 0+ due to '|' doesn't work in strings
2162     _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
2163     _mul($c,$m,$mask);
2164     }
2165   # the loop stops when the shorter of the two numbers is exhausted
2166   # the remainder of the longer one will survive bit-by-bit, so we simple
2167   # multiply-add it in
2168   _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
2169   _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
2170
2171   $x;
2172   }
2173
2174 sub _as_hex
2175   {
2176   # convert a decimal number to hex (ref to array, return ref to string)
2177   my ($c,$x) = @_;
2178
2179   # fits into one element (handle also 0x0 case)
2180   return sprintf("0x%x",$x->[0]) if @$x == 1;
2181
2182   my $x1 = _copy($c,$x);
2183
2184   my $es = '';
2185   my ($xr, $h, $x10000);
2186   if ($] >= 5.006)
2187     {
2188     $x10000 = [ 0x10000 ]; $h = 'h4';
2189     }
2190   else
2191     {
2192     $x10000 = [ 0x1000 ]; $h = 'h3';
2193     }
2194   while (@$x1 != 1 || $x1->[0] != 0)            # _is_zero()
2195     {
2196     ($x1, $xr) = _div($c,$x1,$x10000);
2197     $es .= unpack($h,pack('V',$xr->[0]));
2198     }
2199   $es = reverse $es;
2200   $es =~ s/^[0]+//;   # strip leading zeros
2201   '0x' . $es;                                   # return result prepended with 0x
2202   }
2203
2204 sub _as_bin
2205   {
2206   # convert a decimal number to bin (ref to array, return ref to string)
2207   my ($c,$x) = @_;
2208
2209   # fits into one element (and Perl recent enough), handle also 0b0 case
2210   # handle zero case for older Perls
2211   if ($] <= 5.005 && @$x == 1 && $x->[0] == 0)
2212     {
2213     my $t = '0b0'; return $t;
2214     }
2215   if (@$x == 1 && $] >= 5.006)
2216     {
2217     my $t = sprintf("0b%b",$x->[0]);
2218     return $t;
2219     }
2220   my $x1 = _copy($c,$x);
2221
2222   my $es = '';
2223   my ($xr, $b, $x10000);
2224   if ($] >= 5.006)
2225     {
2226     $x10000 = [ 0x10000 ]; $b = 'b16';
2227     }
2228   else
2229     {
2230     $x10000 = [ 0x1000 ]; $b = 'b12';
2231     }
2232   while (!(@$x1 == 1 && $x1->[0] == 0))         # _is_zero()
2233     {
2234     ($x1, $xr) = _div($c,$x1,$x10000);
2235     $es .= unpack($b,pack('v',$xr->[0]));
2236     }
2237   $es = reverse $es;
2238   $es =~ s/^[0]+//;   # strip leading zeros
2239   '0b' . $es;                                   # return result prepended with 0b
2240   }
2241
2242 sub _as_oct
2243   {
2244   # convert a decimal number to octal (ref to array, return ref to string)
2245   my ($c,$x) = @_;
2246
2247   # fits into one element (handle also 0 case)
2248   return sprintf("0%o",$x->[0]) if @$x == 1;
2249
2250   my $x1 = _copy($c,$x);
2251
2252   my $es = '';
2253   my $xr;
2254   my $x1000 = [ 0100000 ];
2255   while (@$x1 != 1 || $x1->[0] != 0)            # _is_zero()
2256     {
2257     ($x1, $xr) = _div($c,$x1,$x1000);
2258     $es .= reverse sprintf("%05o", $xr->[0]);
2259     }
2260   $es = reverse $es;
2261   $es =~ s/^[0]+//;   # strip leading zeros
2262   '0' . $es;                                    # return result prepended with 0
2263   }
2264
2265 sub _from_oct
2266   {
2267   # convert a octal number to decimal (string, return ref to array)
2268   my ($c,$os) = @_;
2269
2270   # for older Perls, play safe
2271   my $m = [ 0100000 ];
2272   my $d = 5;                                    # 5 digits at a time
2273
2274   my $mul = _one();
2275   my $x = _zero();
2276
2277   my $len = int( (length($os)-1)/$d );          # $d digit parts, w/o the '0'
2278   my $val; my $i = -$d;
2279   while ($len >= 0)
2280     {
2281     $val = substr($os,$i,$d);                   # get oct digits
2282     $val = CORE::oct($val);
2283     $i -= $d; $len --;
2284     my $adder = [ $val ];
2285     _add ($c, $x, _mul ($c, $adder, $mul ) ) if $val != 0;
2286     _mul ($c, $mul, $m ) if $len >= 0;          # skip last mul
2287     }
2288   $x;
2289   }
2290
2291 sub _from_hex
2292   {
2293   # convert a hex number to decimal (string, return ref to array)
2294   my ($c,$hs) = @_;
2295
2296   my $m = _new($c, 0x10000000);                 # 28 bit at a time (<32 bit!)
2297   my $d = 7;                                    # 7 digits at a time
2298   if ($] <= 5.006)
2299     {
2300     # for older Perls, play safe
2301     $m = [ 0x10000 ];                           # 16 bit at a time (<32 bit!)
2302     $d = 4;                                     # 4 digits at a time
2303     }
2304
2305   my $mul = _one();
2306   my $x = _zero();
2307
2308   my $len = int( (length($hs)-2)/$d );          # $d digit parts, w/o the '0x'
2309   my $val; my $i = -$d;
2310   while ($len >= 0)
2311     {
2312     $val = substr($hs,$i,$d);                   # get hex digits
2313     $val =~ s/^0x// if $len == 0;               # for last part only because
2314     $val = CORE::hex($val);                     # hex does not like wrong chars
2315     $i -= $d; $len --;
2316     my $adder = [ $val ];
2317     # if the resulting number was to big to fit into one element, create a
2318     # two-element version (bug found by Mark Lakata - Thanx!)
2319     if (CORE::length($val) > $BASE_LEN)
2320       {
2321       $adder = _new($c,$val);
2322       }
2323     _add ($c, $x, _mul ($c, $adder, $mul ) ) if $val != 0;
2324     _mul ($c, $mul, $m ) if $len >= 0;          # skip last mul
2325     }
2326   $x;
2327   }
2328
2329 sub _from_bin
2330   {
2331   # convert a hex number to decimal (string, return ref to array)
2332   my ($c,$bs) = @_;
2333
2334   # instead of converting X (8) bit at a time, it is faster to "convert" the
2335   # number to hex, and then call _from_hex.
2336
2337   my $hs = $bs;
2338   $hs =~ s/^[+-]?0b//;                                  # remove sign and 0b
2339   my $l = length($hs);                                  # bits
2340   $hs = '0' x (8-($l % 8)) . $hs if ($l % 8) != 0;      # padd left side w/ 0
2341   my $h = '0x' . unpack('H*', pack ('B*', $hs));        # repack as hex
2342
2343   $c->_from_hex($h);
2344   }
2345
2346 ##############################################################################
2347 # special modulus functions
2348
2349 sub _modinv
2350   {
2351   # modular inverse
2352   my ($c,$x,$y) = @_;
2353
2354   my $u = _zero($c); my $u1 = _one($c);
2355   my $a = _copy($c,$y); my $b = _copy($c,$x);
2356
2357   # Euclid's Algorithm for bgcd(), only that we calc bgcd() ($a) and the
2358   # result ($u) at the same time. See comments in BigInt for why this works.
2359   my $q;
2360   ($a, $q, $b) = ($b, _div($c,$a,$b));          # step 1
2361   my $sign = 1;
2362   while (!_is_zero($c,$b))
2363     {
2364     my $t = _add($c,                            # step 2:
2365        _mul($c,_copy($c,$u1), $q) ,             #  t =  u1 * q
2366        $u );                                    #     + u
2367     $u = $u1;                                   #  u = u1, u1 = t
2368     $u1 = $t;
2369     $sign = -$sign;
2370     ($a, $q, $b) = ($b, _div($c,$a,$b));        # step 1
2371     }
2372
2373   # if the gcd is not 1, then return NaN
2374   return (undef,undef) unless _is_one($c,$a);
2375
2376   ($u1, $sign == 1 ? '+' : '-');
2377   }
2378
2379 sub _modpow
2380   {
2381   # modulus of power ($x ** $y) % $z
2382   my ($c,$num,$exp,$mod) = @_;
2383
2384   # in the trivial case,
2385   if (_is_one($c,$mod))
2386     {
2387     splice @$num,0,1; $num->[0] = 0;
2388     return $num;
2389     }
2390   if ((scalar @$num == 1) && (($num->[0] == 0) || ($num->[0] == 1)))
2391     {
2392     $num->[0] = 1;
2393     return $num;
2394     }
2395
2396 #  $num = _mod($c,$num,$mod);   # this does not make it faster
2397
2398   my $acc = _copy($c,$num); my $t = _one();
2399
2400   my $expbin = _as_bin($c,$exp); $expbin =~ s/^0b//;
2401   my $len = length($expbin);
2402   while (--$len >= 0)
2403     {
2404     if ( substr($expbin,$len,1) eq '1')                 # is_odd
2405       {
2406       _mul($c,$t,$acc);
2407       $t = _mod($c,$t,$mod);
2408       }
2409     _mul($c,$acc,$acc);
2410     $acc = _mod($c,$acc,$mod);
2411     }
2412   @$num = @$t;
2413   $num;
2414   }
2415
2416 sub _gcd
2417   {
2418   # greatest common divisor
2419   my ($c,$x,$y) = @_;
2420
2421   while ( (scalar @$y != 1) || ($y->[0] != 0) )         # while ($y != 0)
2422     {
2423     my $t = _copy($c,$y);
2424     $y = _mod($c, $x, $y);
2425     $x = $t;
2426     }
2427   $x;
2428   }
2429
2430 ##############################################################################
2431 ##############################################################################
2432
2433 1;
2434 __END__
2435
2436 =head1 NAME
2437
2438 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
2439
2440 =head1 SYNOPSIS
2441
2442 Provides support for big integer calculations. Not intended to be used by other
2443 modules. Other modules which sport the same functions can also be used to support
2444 Math::BigInt, like Math::BigInt::GMP or Math::BigInt::Pari.
2445
2446 =head1 DESCRIPTION
2447
2448 In order to allow for multiple big integer libraries, Math::BigInt was
2449 rewritten to use library modules for core math routines. Any module which
2450 follows the same API as this can be used instead by using the following:
2451
2452         use Math::BigInt lib => 'libname';
2453
2454 'libname' is either the long name ('Math::BigInt::Pari'), or only the short
2455 version like 'Pari'.
2456
2457 =head1 STORAGE
2458
2459 =head1 METHODS
2460
2461 The following functions MUST be defined in order to support the use by
2462 Math::BigInt v1.70 or later:
2463
2464         api_version()   return API version, 1 for v1.70, 2 for v1.83
2465         _new(string)    return ref to new object from ref to decimal string
2466         _zero()         return a new object with value 0
2467         _one()          return a new object with value 1
2468         _two()          return a new object with value 2
2469         _ten()          return a new object with value 10
2470
2471         _str(obj)       return ref to a string representing the object
2472         _num(obj)       returns a Perl integer/floating point number
2473                         NOTE: because of Perl numeric notation defaults,
2474                         the _num'ified obj may lose accuracy due to
2475                         machine-dependent floating point size limitations
2476
2477         _add(obj,obj)   Simple addition of two objects
2478         _mul(obj,obj)   Multiplication of two objects
2479         _div(obj,obj)   Division of the 1st object by the 2nd
2480                         In list context, returns (result,remainder).
2481                         NOTE: this is integer math, so no
2482                         fractional part will be returned.
2483                         The second operand will be not be 0, so no need to
2484                         check for that.
2485         _sub(obj,obj)   Simple subtraction of 1 object from another
2486                         a third, optional parameter indicates that the params
2487                         are swapped. In this case, the first param needs to
2488                         be preserved, while you can destroy the second.
2489                         sub (x,y,1) => return x - y and keep x intact!
2490         _dec(obj)       decrement object by one (input is guaranteed to be > 0)
2491         _inc(obj)       increment object by one
2492
2493
2494         _acmp(obj,obj)  <=> operator for objects (return -1, 0 or 1)
2495
2496         _len(obj)       returns count of the decimal digits of the object
2497         _digit(obj,n)   returns the n'th decimal digit of object
2498
2499         _is_one(obj)    return true if argument is 1
2500         _is_two(obj)    return true if argument is 2
2501         _is_ten(obj)    return true if argument is 10
2502         _is_zero(obj)   return true if argument is 0
2503         _is_even(obj)   return true if argument is even (0,2,4,6..)
2504         _is_odd(obj)    return true if argument is odd (1,3,5,7..)
2505
2506         _copy           return a ref to a true copy of the object
2507
2508         _check(obj)     check whether internal representation is still intact
2509                         return 0 for ok, otherwise error message as string
2510
2511         _from_hex(str)  return new object from a hexadecimal string
2512         _from_bin(str)  return new object from a binary string
2513         _from_oct(str)  return new object from an octal string
2514
2515         _as_hex(str)    return string containing the value as
2516                         unsigned hex string, with the '0x' prepended.
2517                         Leading zeros must be stripped.
2518         _as_bin(str)    Like as_hex, only as binary string containing only
2519                         zeros and ones. Leading zeros must be stripped and a
2520                         '0b' must be prepended.
2521
2522         _rsft(obj,N,B)  shift object in base B by N 'digits' right
2523         _lsft(obj,N,B)  shift object in base B by N 'digits' left
2524
2525         _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
2526                         Note: XOR, AND and OR pad with zeros if size mismatches
2527         _and(obj1,obj2) AND (bit-wise) object 1 with object 2
2528         _or(obj1,obj2)  OR (bit-wise) object 1 with object 2
2529
2530         _mod(obj1,obj2) Return remainder of div of the 1st by the 2nd object
2531         _sqrt(obj)      return the square root of object (truncated to int)
2532         _root(obj)      return the n'th (n >= 3) root of obj (truncated to int)
2533         _fac(obj)       return factorial of object 1 (1*2*3*4..)
2534         _pow(obj1,obj2) return object 1 to the power of object 2
2535                         return undef for NaN
2536         _zeros(obj)     return number of trailing decimal zeros
2537         _modinv         return inverse modulus
2538         _modpow         return modulus of power ($x ** $y) % $z
2539         _log_int(X,N)   calculate integer log() of X in base N
2540                         X >= 0, N >= 0 (return undef for NaN)
2541                         returns (RESULT, EXACT) where EXACT is:
2542                          1     : result is exactly RESULT
2543                          0     : result was truncated to RESULT
2544                          undef : unknown whether result is exactly RESULT
2545         _gcd(obj,obj)   return Greatest Common Divisor of two objects
2546
2547 The following functions are REQUIRED for an api_version of 2 or greater:
2548
2549         _1ex($x)        create the number 1Ex where x >= 0
2550         _alen(obj)      returns approximate count of the decimal digits of the
2551                         object. This estimate MUST always be greater or equal
2552                         to what _len() returns.
2553         _nok(n,k)       calculate n over k (binomial coefficient)
2554
2555 The following functions are optional, and can be defined if the underlying lib
2556 has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
2557 slow) fallback routines to emulate these:
2558
2559         _signed_or
2560         _signed_and
2561         _signed_xor
2562
2563 Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
2564 or '0b1101').
2565
2566 So the library needs only to deal with unsigned big integers. Testing of input
2567 parameter validity is done by the caller, so you need not worry about
2568 underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by zero or similar
2569 cases.
2570
2571 The first parameter can be modified, that includes the possibility that you
2572 return a reference to a completely different object instead. Although keeping
2573 the reference and just changing its contents is preferred over creating and
2574 returning a different reference.
2575
2576 Return values are always references to objects, strings, or true/false for
2577 comparison routines.
2578
2579 =head1 WRAP YOUR OWN
2580
2581 If you want to port your own favourite c-lib for big numbers to the
2582 Math::BigInt interface, you can take any of the already existing modules as
2583 a rough guideline. You should really wrap up the latest BigInt and BigFloat
2584 testsuites with your module, and replace in them any of the following:
2585
2586         use Math::BigInt;
2587
2588 by this:
2589
2590         use Math::BigInt lib => 'yourlib';
2591
2592 This way you ensure that your library really works 100% within Math::BigInt.
2593
2594 =head1 LICENSE
2595
2596 This program is free software; you may redistribute it and/or modify it under
2597 the same terms as Perl itself.
2598
2599 =head1 AUTHORS
2600
2601 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
2602 in late 2000.
2603 Seperated from BigInt and shaped API with the help of John Peacock.
2604
2605 Fixed, speed-up, streamlined and enhanced by Tels 2001 - 2007.
2606
2607 =head1 SEE ALSO
2608
2609 L<Math::BigInt>, L<Math::BigFloat>,
2610 L<Math::BigInt::GMP>, L<Math::BigInt::FastCalc> and L<Math::BigInt::Pari>.
2611
2612 =cut