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