lib/Math/BigInt/Calc.pm

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