lib/Math/BigInt/Calc.pm

   1 package Math::BigInt::Calc;
   2
   3 use 5.005;
   4 use strict;
   5 # use warnings; # dont use warnings for older Perls
   6
   7 require Exporter;
   8 use vars qw/@ISA $VERSION/;
   9 @ISA = qw(Exporter);
  10
  11 $VERSION = '0.32';
  12
  13 # Package to store unsigned big integers in decimal and do math with them
  14
  15 # Internally the numbers are stored in an array with at least 1 element, no
  16 # leading zero parts (except the first) and in base 1eX where X is determined
  17 # automatically at loading time to be the maximum possible value
  18
  19 # todo:
  20 # - fully remove funky $# stuff (maybe)
  21
  22 # USE_MUL: due to problems on certain os (os390, posix-bc) "* 1e-5" is used
  23 # instead of "/ 1e5" at some places, (marked with USE_MUL). Other platforms
  24 # BS2000, some Crays need USE_DIV instead.
  25 # The BEGIN block is used to determine which of the two variants gives the
  26 # correct result.
  27
  28 ##############################################################################
  29 # global constants, flags and accessory
  30
  31 # constants for easier life
  32 my $nan = 'NaN';
  33 my ($MBASE,$BASE,$RBASE,$BASE_LEN,$MAX_VAL,$BASE_LEN2,$BASE_LEN_SMALL);
  34 my ($AND_BITS,$XOR_BITS,$OR_BITS);
  35 my ($AND_MASK,$XOR_MASK,$OR_MASK);
  36 my ($LEN_CONVERT);
  37
  38 sub _base_len
  39   {
  40   # set/get the BASE_LEN and assorted other, connected values
  41   # used only be the testsuite, set is used only by the BEGIN block below
  42   shift;
  43
  44   my $b = shift;
  45   if (defined $b)
  46     {
  47     # find whether we can use mul or div or none in mul()/div()
  48     # (in last case reduce BASE_LEN_SMALL)
  49     $BASE_LEN_SMALL = $b+1;
  50     my $caught = 0;
  51     while (--$BASE_LEN_SMALL > 5)
  52       {
  53       $MBASE = int("1e".$BASE_LEN_SMALL);
  54       $RBASE = abs('1e-'.$BASE_LEN_SMALL);              # see USE_MUL
  55       $caught = 0;
  56       $caught += 1 if (int($MBASE * $RBASE) != 1);      # should be 1
  57       $caught += 2 if (int($MBASE / $MBASE) != 1);      # should be 1
  58       last if $caught != 3;
  59       }
  60     # BASE_LEN is used for anything else than mul()/div()
  61     $BASE_LEN = $BASE_LEN_SMALL;
  62     $BASE_LEN = shift if (defined $_[0]);               # one more arg?
  63     $BASE = int("1e".$BASE_LEN);
  64
  65     $BASE_LEN2 = int($BASE_LEN_SMALL / 2);              # for mul shortcut
  66     $MBASE = int("1e".$BASE_LEN_SMALL);
  67     $RBASE = abs('1e-'.$BASE_LEN_SMALL);                # see USE_MUL
  68     $MAX_VAL = $MBASE-1;
  69     $LEN_CONVERT = 0;
  70     $LEN_CONVERT = 1 if $BASE_LEN_SMALL != $BASE_LEN;
  71
  72     #print "BASE_LEN: $BASE_LEN MAX_VAL: $MAX_VAL BASE: $BASE RBASE: $RBASE ";
  73     #print "BASE_LEN_SMALL: $BASE_LEN_SMALL MBASE: $MBASE\n";
  74
  75     undef &_mul;
  76     undef &_div;
  77
  78     if ($caught & 1 != 0)
  79       {
  80       # must USE_MUL
  81       *{_mul} = \&_mul_use_mul;
  82       *{_div} = \&_div_use_mul;
  83       }
  84     else                # $caught must be 2, since it can't be 1 nor 3
  85       {
  86       # can USE_DIV instead
  87       *{_mul} = \&_mul_use_div;
  88       *{_div} = \&_div_use_div;
  89       }
  90     }
  91   return $BASE_LEN unless wantarray;
  92   return ($BASE_LEN, $AND_BITS, $XOR_BITS, $OR_BITS, $BASE_LEN_SMALL, $MAX_VAL);
  93   }
  94
  95 BEGIN
  96   {
  97   # from Daniel Pfeiffer: determine largest group of digits that is precisely
  98   # multipliable with itself plus carry
  99   # Test now changed to expect the proper pattern, not a result off by 1 or 2
 100   my ($e, $num) = 3;    # lowest value we will use is 3+1-1 = 3
 101   do
 102     {
 103     $num = ('9' x ++$e) + 0;
 104     $num *= $num + 1.0;
 105     } while ("$num" =~ /9{$e}0{$e}/);   # must be a certain pattern
 106   $e--;                                 # last test failed, so retract one step
 107   # the limits below brush the problems with the test above under the rug:
 108   # the test should be able to find the proper $e automatically
 109   $e = 5 if $^O =~ /^uts/;      # UTS get's some special treatment
 110   $e = 5 if $^O =~ /^unicos/;   # unicos is also problematic (6 seems to work
 111                                 # there, but we play safe)
 112   $e = 5 if $] < 5.006;         # cap, for older Perls
 113   $e = 7 if $e > 7;             # cap, for VMS, OS/390 and other 64 bit systems
 114                                 # 8 fails inside random testsuite, so take 7
 115
 116   # determine how many digits fit into an integer and can be safely added
 117   # together plus carry w/o causing an overflow
 118
 119   # this below detects 15 on a 64 bit system, because after that it becomes
 120   # 1e16  and not 1000000 :/ I can make it detect 18, but then I get a lot of
 121   # test failures. Ugh! (Tomake detect 18: uncomment lines marked with *)
 122   use integer;
 123   my $bi = 5;                   # approx. 16 bit
 124   $num = int('9' x $bi);
 125   # $num = 99999; # *
 126   # while ( ($num+$num+1) eq '1' . '9' x $bi)   # *
 127   while ( int($num+$num+1) eq '1' . '9' x $bi)
 128     {
 129     $bi++; $num = int('9' x $bi);
 130     # $bi++; $num *= 10; $num += 9;     # *
 131     }
 132   $bi--;                                # back off one step
 133   # by setting them equal, we ignore the findings and use the default
 134   # one-size-fits-all approach from former versions
 135   $bi = $e;                             # XXX, this should work always
 136
 137   __PACKAGE__->_base_len($e,$bi);       # set and store
 138
 139   # find out how many bits _and, _or and _xor can take (old default = 16)
 140   # I don't think anybody has yet 128 bit scalars, so let's play safe.
 141   local $^W = 0;        # don't warn about 'nonportable number'
 142   $AND_BITS = 15; $XOR_BITS = 15; $OR_BITS  = 15;
 143
 144   # find max bits, we will not go higher than numberofbits that fit into $BASE
 145   # to make _and etc simpler (and faster for smaller, slower for large numbers)
 146   my $max = 16;
 147   while (2 ** $max < $BASE) { $max++; }
 148   {
 149     no integer;
 150     $max = 16 if $] < 5.006;    # older Perls might not take >16 too well
 151   }
 152   my ($x,$y,$z);
 153   do {
 154     $AND_BITS++;
 155     $x = oct('0b' . '1' x $AND_BITS); $y = $x & $x;
 156     $z = (2 ** $AND_BITS) - 1;
 157     } while ($AND_BITS < $max && $x == $z && $y == $x);
 158   $AND_BITS --;                                         # retreat one step
 159   do {
 160     $XOR_BITS++;
 161     $x = oct('0b' . '1' x $XOR_BITS); $y = $x ^ 0;
 162     $z = (2 ** $XOR_BITS) - 1;
 163     } while ($XOR_BITS < $max && $x == $z && $y == $x);
 164   $XOR_BITS --;                                         # retreat one step
 165   do {
 166     $OR_BITS++;
 167     $x = oct('0b' . '1' x $OR_BITS); $y = $x | $x;
 168     $z = (2 ** $OR_BITS) - 1;
 169     } while ($OR_BITS < $max && $x == $z && $y == $x);
 170   $OR_BITS --;                                          # retreat one step
 171
 172   }
 173
 174 ##############################################################################
 175 # convert between the "small" and the "large" representation
 176
 177 sub _to_large
 178   {
 179   # take an array in base $BASE_LEN_SMALL and convert it in-place to $BASE_LEN
 180   my ($c,$x) = @_;
 181
 182 #  print "_to_large $BASE_LEN_SMALL => $BASE_LEN\n";
 183
 184   return $x if $LEN_CONVERT == 0 ||             # nothing to converconvertor
 185          @$x == 1;                              # only one element => early out
 186
 187   #     12345    67890    12345    67890   contents
 188   # to      3        2        1        0   index
 189   #             123456  7890123  4567890   contents
 190
 191 #  # faster variant
 192 #  my @d; my $str = '';
 193 #  my $z = '0' x $BASE_LEN_SMALL;
 194 #  foreach (@$x)
 195 #    {
 196 #    # ... . 04321 . 000321
 197 #    $str = substr($z.$_,-$BASE_LEN_SMALL,$BASE_LEN_SMALL) . $str;
 198 #    if (length($str) > $BASE_LEN)
 199 #      {
 200 #      push @d, substr($str,-$BASE_LEN,$BASE_LEN);      # extract one piece
 201 #      substr($str,-$BASE_LEN,$BASE_LEN) = '';          # remove it
 202 #      }
 203 #    }
 204 #  push @d, $str if $str !~ /^0*$/;                     # extract last piece
 205 #  @$x = @d;
 206 #  $x->[-1] = int($x->[-1]);                    # strip leading zero
 207 #  $x;
 208
 209   my $ret = "";
 210   my $l = scalar @$x;           # number of parts
 211   $l --; $ret .= int($x->[$l]); $l--;
 212   my $z = '0' x ($BASE_LEN_SMALL-1);
 213   while ($l >= 0)
 214     {
 215     $ret .= substr($z.$x->[$l],-$BASE_LEN_SMALL);
 216     $l--;
 217     }
 218   my $str = _new($c,\$ret);                     # make array
 219   @$x = @$str;                                  # clobber contents of $x
 220   $x->[-1] = int($x->[-1]);                     # strip leading zero
 221   }
 222
 223 sub _to_small
 224   {
 225   # take an array in base $BASE_LEN and convert it in-place to $BASE_LEN_SMALL
 226   my ($c,$x) = @_;
 227
 228   return $x if $LEN_CONVERT == 0;               # nothing to do
 229   return $x if @$x == 1 && length(int($x->[0])) <= $BASE_LEN_SMALL;
 230
 231   my $d = _str($c,$x);
 232   my $il = length($$d)-1;
 233   ## this leaves '00000' instead of int 0 and will be corrected after any op
 234   # clobber contents of $x
 235   @$x = reverse(unpack("a" . ($il % $BASE_LEN_SMALL+1)
 236       . ("a$BASE_LEN_SMALL" x ($il / $BASE_LEN_SMALL)), $$d));
 237
 238   $x->[-1] = int($x->[-1]);                     # strip leading zero
 239   }
 240
 241 ###############################################################################
 242
 243 sub _new
 244   {
 245   # (ref to string) return ref to num_array
 246   # Convert a number from string format (without sign) to internal base
 247   # 1ex format. Assumes normalized value as input.
 248   my $d = $_[1];
 249   my $il = length($$d)-1;
 250   # this leaves '00000' instead of int 0 and will be corrected after any op
 251   [ reverse(unpack("a" . ($il % $BASE_LEN+1)
 252     . ("a$BASE_LEN" x ($il / $BASE_LEN)), $$d)) ];
 253   }
 254
 255 BEGIN
 256   {
 257   $AND_MASK = __PACKAGE__->_new( \( 2 ** $AND_BITS ));
 258   $XOR_MASK = __PACKAGE__->_new( \( 2 ** $XOR_BITS ));
 259   $OR_MASK = __PACKAGE__->_new( \( 2 ** $OR_BITS ));
 260   }
 261
 262 sub _zero
 263   {
 264   # create a zero
 265   [ 0 ];
 266   }
 267
 268 sub _one
 269   {
 270   # create a one
 271   [ 1 ];
 272   }
 273
 274 sub _two
 275   {
 276   # create a two (used internally for shifting)
 277   [ 2 ];
 278   }
 279
 280 sub _copy
 281   {
 282   [ @{$_[1]} ];
 283   }
 284
 285 # catch and throw away
 286 sub import { }
 287
 288 ##############################################################################
 289 # convert back to string and number
 290
 291 sub _str
 292   {
 293   # (ref to BINT) return num_str
 294   # Convert number from internal base 100000 format to string format.
 295   # internal format is always normalized (no leading zeros, "-0" => "+0")
 296   my $ar = $_[1];
 297   my $ret = "";
 298
 299   my $l = scalar @$ar;          # number of parts
 300   return $nan if $l < 1;        # should not happen
 301
 302   # handle first one different to strip leading zeros from it (there are no
 303   # leading zero parts in internal representation)
 304   $l --; $ret .= int($ar->[$l]); $l--;
 305   # Interestingly, the pre-padd method uses more time
 306   # the old grep variant takes longer (14 to 10 sec)
 307   my $z = '0' x ($BASE_LEN-1);
 308   while ($l >= 0)
 309     {
 310     $ret .= substr($z.$ar->[$l],-$BASE_LEN); # fastest way I could think of
 311     $l--;
 312     }
 313   \$ret;
 314   }
 315
 316 sub _num
 317   {
 318   # Make a number (scalar int/float) from a BigInt object
 319   my $x = $_[1];
 320   return $x->[0] if scalar @$x == 1;  # below $BASE
 321   my $fac = 1;
 322   my $num = 0;
 323   foreach (@$x)
 324     {
 325     $num += $fac*$_; $fac *= $BASE;
 326     }
 327   $num;
 328   }
 329
 330 ##############################################################################
 331 # actual math code
 332
 333 sub _add
 334   {
 335   # (ref to int_num_array, ref to int_num_array)
 336   # routine to add two base 1eX numbers
 337   # stolen from Knuth Vol 2 Algorithm A pg 231
 338   # there are separate routines to add and sub as per Knuth pg 233
 339   # This routine clobbers up array x, but not y.
 340
 341   my ($c,$x,$y) = @_;
 342
 343   return $x if (@$y == 1) && $y->[0] == 0;              # $x + 0 => $x
 344   if ((@$x == 1) && $x->[0] == 0)                       # 0 + $y => $y->copy
 345     {
 346     # twice as slow as $x = [ @$y ], but necc. to retain $x as ref :(
 347     @$x = @$y; return $x;
 348     }
 349
 350   # for each in Y, add Y to X and carry. If after that, something is left in
 351   # X, foreach in X add carry to X and then return X, carry
 352   # Trades one "$j++" for having to shift arrays, $j could be made integer
 353   # but this would impose a limit to number-length of 2**32.
 354   my $i; my $car = 0; my $j = 0;
 355   for $i (@$y)
 356     {
 357     $x->[$j] -= $BASE if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
 358     $j++;
 359     }
 360   while ($car != 0)
 361     {
 362     $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
 363     }
 364   $x;
 365   }
 366
 367 sub _inc
 368   {
 369   # (ref to int_num_array, ref to int_num_array)
 370   # routine to add 1 to a base 1eX numbers
 371   # This routine modifies array x
 372   my ($c,$x) = @_;
 373
 374   for my $i (@$x)
 375     {
 376     return $x if (($i += 1) < $BASE);           # early out
 377     $i = 0;                                     # overflow, next
 378     }
 379   push @$x,1 if ($x->[-1] == 0);                # last overflowed, so extend
 380   $x;
 381   }
 382
 383 sub _dec
 384   {
 385   # (ref to int_num_array, ref to int_num_array)
 386   # routine to add 1 to a base 1eX numbers
 387   # This routine modifies array x
 388   my ($c,$x) = @_;
 389
 390   my $MAX = $BASE-1;                            # since MAX_VAL based on MBASE
 391   for my $i (@$x)
 392     {
 393     last if (($i -= 1) >= 0);                   # early out
 394     $i = $MAX;                                  # overflow, next
 395     }
 396   pop @$x if $x->[-1] == 0 && @$x > 1;          # last overflowed (but leave 0)
 397   $x;
 398   }
 399
 400 sub _sub
 401   {
 402   # (ref to int_num_array, ref to int_num_array, swap)
 403   # subtract base 1eX numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
 404   # subtract Y from X by modifying x in place
 405   my ($c,$sx,$sy,$s) = @_;
 406
 407   my $car = 0; my $i; my $j = 0;
 408   if (!$s)
 409     {
 410     #print "case 2\n";
 411     for $i (@$sx)
 412       {
 413       last unless defined $sy->[$j] || $car;
 414       $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
 415       }
 416     # might leave leading zeros, so fix that
 417     return __strip_zeros($sx);
 418     }
 419   #print "case 1 (swap)\n";
 420   for $i (@$sx)
 421     {
 422     # we can't do an early out if $x is < than $y, since we
 423     # need to copy the high chunks from $y. Found by Bob Mathews.
 424     #last unless defined $sy->[$j] || $car;
 425     $sy->[$j] += $BASE
 426      if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
 427     $j++;
 428     }
 429   # might leave leading zeros, so fix that
 430   __strip_zeros($sy);
 431   }
 432
 433 sub _mul_use_mul
 434   {
 435   # (ref to int_num_array, ref to int_num_array)
 436   # multiply two numbers in internal representation
 437   # modifies first arg, second need not be different from first
 438   my ($c,$xv,$yv) = @_;
 439
 440   # shortcut for two very short numbers (improved by Nathan Zook)
 441   # works also if xv and yv are the same reference
 442   if ((@$xv == 1) && (@$yv == 1))
 443     {
 444     if (($xv->[0] *= $yv->[0]) >= $MBASE)
 445        {
 446        $xv->[0] = $xv->[0] - ($xv->[1] = int($xv->[0] * $RBASE)) * $MBASE;
 447        };
 448     return $xv;
 449     }
 450   # shortcut for result == 0
 451   if ( ((@$xv == 1) && ($xv->[0] == 0)) ||
 452        ((@$yv == 1) && ($yv->[0] == 0)) )
 453     {
 454     @$xv = (0);
 455     return $xv;
 456     }
 457
 458   # since multiplying $x with $x fails, make copy in this case
 459   $yv = [@$xv] if $xv == $yv;   # same references?
 460
 461   if ($LEN_CONVERT != 0)
 462     {
 463     $c->_to_small($xv); $c->_to_small($yv);
 464     }
 465
 466   my @prod = (); my ($prod,$car,$cty,$xi,$yi);
 467
 468   for $xi (@$xv)
 469     {
 470     $car = 0; $cty = 0;
 471
 472     # slow variant
 473 #    for $yi (@$yv)
 474 #      {
 475 #      $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
 476 #      $prod[$cty++] =
 477 #       $prod - ($car = int($prod * RBASE)) * $MBASE;  # see USE_MUL
 478 #      }
 479 #    $prod[$cty] += $car if $car; # need really to check for 0?
 480 #    $xi = shift @prod;
 481
 482     # faster variant
 483     # looping through this if $xi == 0 is silly - so optimize it away!
 484     $xi = (shift @prod || 0), next if $xi == 0;
 485     for $yi (@$yv)
 486       {
 487       $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
 488 ##     this is actually a tad slower
 489 ##        $prod = $prod[$cty]; $prod += ($car + $xi * $yi);     # no ||0 here
 490       $prod[$cty++] =
 491        $prod - ($car = int($prod * $RBASE)) * $MBASE;  # see USE_MUL
 492       }
 493     $prod[$cty] += $car if $car; # need really to check for 0?
 494     $xi = shift @prod || 0;     # || 0 makes v5.005_3 happy
 495     }
 496   push @$xv, @prod;
 497   if ($LEN_CONVERT != 0)
 498     {
 499     $c->_to_large($yv);
 500     $c->_to_large($xv);
 501     }
 502   else
 503     {
 504     __strip_zeros($xv);
 505     }
 506   $xv;
 507   }
 508
 509 sub _mul_use_div
 510   {
 511   # (ref to int_num_array, ref to int_num_array)
 512   # multiply two numbers in internal representation
 513   # modifies first arg, second need not be different from first
 514   my ($c,$xv,$yv) = @_;
 515
 516   # shortcut for two very short numbers (improved by Nathan Zook)
 517   # works also if xv and yv are the same reference
 518   if ((@$xv == 1) && (@$yv == 1))
 519     {
 520     if (($xv->[0] *= $yv->[0]) >= $MBASE)
 521         {
 522         $xv->[0] =
 523             $xv->[0] - ($xv->[1] = int($xv->[0] / $MBASE)) * $MBASE;
 524         };
 525     return $xv;
 526     }
 527   # shortcut for result == 0
 528   if ( ((@$xv == 1) && ($xv->[0] == 0)) ||
 529        ((@$yv == 1) && ($yv->[0] == 0)) )
 530     {
 531     @$xv = (0);
 532     return $xv;
 533     }
 534
 535
 536   # since multiplying $x with $x fails, make copy in this case
 537   $yv = [@$xv] if $xv == $yv;   # same references?
 538
 539   if ($LEN_CONVERT != 0)
 540     {
 541     $c->_to_small($xv); $c->_to_small($yv);
 542     }
 543
 544   my @prod = (); my ($prod,$car,$cty,$xi,$yi);
 545   for $xi (@$xv)
 546     {
 547     $car = 0; $cty = 0;
 548     # looping through this if $xi == 0 is silly - so optimize it away!
 549     $xi = (shift @prod || 0), next if $xi == 0;
 550     for $yi (@$yv)
 551       {
 552       $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
 553       $prod[$cty++] =
 554        $prod - ($car = int($prod / $MBASE)) * $MBASE;
 555       }
 556     $prod[$cty] += $car if $car; # need really to check for 0?
 557     $xi = shift @prod || 0;     # || 0 makes v5.005_3 happy
 558     }
 559   push @$xv, @prod;
 560   if ($LEN_CONVERT != 0)
 561     {
 562     $c->_to_large($yv);
 563     $c->_to_large($xv);
 564     }
 565   else
 566     {
 567     __strip_zeros($xv);
 568     }
 569   $xv;
 570   }
 571
 572 sub _div_use_mul
 573   {
 574   # ref to array, ref to array, modify first array and return remainder if
 575   # in list context
 576   my ($c,$x,$yorg) = @_;
 577
 578   if (@$x == 1 && @$yorg == 1)
 579     {
 580     # shortcut, $yorg and $x are two small numbers
 581     if (wantarray)
 582       {
 583       my $r = [ $x->[0] % $yorg->[0] ];
 584       $x->[0] = int($x->[0] / $yorg->[0]);
 585       return ($x,$r);
 586       }
 587     else
 588       {
 589       $x->[0] = int($x->[0] / $yorg->[0]);
 590       return $x;
 591       }
 592     }
 593   if (@$yorg == 1)
 594     {
 595     my $rem;
 596     $rem = _mod($c,[ @$x ],$yorg) if wantarray;
 597
 598     # shortcut, $y is < $BASE
 599     my $j = scalar @$x; my $r = 0;
 600     my $y = $yorg->[0]; my $b;
 601     while ($j-- > 0)
 602       {
 603       $b = $r * $MBASE + $x->[$j];
 604       $x->[$j] = int($b/$y);
 605       $r = $b % $y;
 606       }
 607     pop @$x if @$x > 1 && $x->[-1] == 0;        # splice up a leading zero
 608     return ($x,$rem) if wantarray;
 609     return $x;
 610     }
 611
 612   my $y = [ @$yorg ];                           # always make copy to preserve
 613   if ($LEN_CONVERT != 0)
 614     {
 615     $c->_to_small($x); $c->_to_small($y);
 616     }
 617
 618   my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
 619
 620   $car = $bar = $prd = 0;
 621   if (($dd = int($MBASE/($y->[-1]+1))) != 1)
 622     {
 623     for $xi (@$x)
 624       {
 625       $xi = $xi * $dd + $car;
 626       $xi -= ($car = int($xi * $RBASE)) * $MBASE;       # see USE_MUL
 627       }
 628     push(@$x, $car); $car = 0;
 629     for $yi (@$y)
 630       {
 631       $yi = $yi * $dd + $car;
 632       $yi -= ($car = int($yi * $RBASE)) * $MBASE;       # see USE_MUL
 633       }
 634     }
 635   else
 636     {
 637     push(@$x, 0);
 638     }
 639   @q = (); ($v2,$v1) = @$y[-2,-1];
 640   $v2 = 0 unless $v2;
 641   while ($#$x > $#$y)
 642     {
 643     ($u2,$u1,$u0) = @$x[-3..-1];
 644     $u2 = 0 unless $u2;
 645     #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
 646     # if $v1 == 0;
 647      $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
 648     --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
 649     if ($q)
 650       {
 651       ($car, $bar) = (0,0);
 652       for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
 653         {
 654         $prd = $q * $y->[$yi] + $car;
 655         $prd -= ($car = int($prd * $RBASE)) * $MBASE;   # see USE_MUL
 656         $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
 657         }
 658       if ($x->[-1] < $car + $bar)
 659         {
 660         $car = 0; --$q;
 661         for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
 662           {
 663           $x->[$xi] -= $MBASE
 664            if ($car = (($x->[$xi] += $y->[$yi] + $car) > $MBASE));
 665           }
 666         }
 667       }
 668       pop(@$x); unshift(@q, $q);
 669     }
 670   if (wantarray)
 671     {
 672     @d = ();
 673     if ($dd != 1)
 674       {
 675       $car = 0;
 676       for $xi (reverse @$x)
 677         {
 678         $prd = $car * $MBASE + $xi;
 679         $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
 680         unshift(@d, $tmp);
 681         }
 682       }
 683     else
 684       {
 685       @d = @$x;
 686       }
 687     @$x = @q;
 688     my $d = \@d;
 689     if ($LEN_CONVERT != 0)
 690       {
 691       $c->_to_large($x); $c->_to_large($d);
 692       }
 693     else
 694       {
 695       __strip_zeros($x);
 696       __strip_zeros($d);
 697       }
 698     return ($x,$d);
 699     }
 700   @$x = @q;
 701   if ($LEN_CONVERT != 0)
 702     {
 703     $c->_to_large($x);
 704     }
 705   else
 706     {
 707     __strip_zeros($x);
 708     }
 709   $x;
 710   }
 711
 712 sub _div_use_div
 713   {
 714   # ref to array, ref to array, modify first array and return remainder if
 715   # in list context
 716   my ($c,$x,$yorg) = @_;
 717
 718   if (@$x == 1 && @$yorg == 1)
 719     {
 720     # shortcut, $yorg and $x are two small numbers
 721     if (wantarray)
 722       {
 723       my $r = [ $x->[0] % $yorg->[0] ];
 724       $x->[0] = int($x->[0] / $yorg->[0]);
 725       return ($x,$r);
 726       }
 727     else
 728       {
 729       $x->[0] = int($x->[0] / $yorg->[0]);
 730       return $x;
 731       }
 732     }
 733   if (@$yorg == 1)
 734     {
 735     my $rem;
 736     $rem = _mod($c,[ @$x ],$yorg) if wantarray;
 737
 738     # shortcut, $y is < $BASE
 739     my $j = scalar @$x; my $r = 0;
 740     my $y = $yorg->[0]; my $b;
 741     while ($j-- > 0)
 742       {
 743       $b = $r * $MBASE + $x->[$j];
 744       $x->[$j] = int($b/$y);
 745       $r = $b % $y;
 746       }
 747     pop @$x if @$x > 1 && $x->[-1] == 0;        # splice up a leading zero
 748     return ($x,$rem) if wantarray;
 749     return $x;
 750     }
 751
 752   my $y = [ @$yorg ];                           # always make copy to preserve
 753   if ($LEN_CONVERT != 0)
 754     {
 755     $c->_to_small($x); $c->_to_small($y);
 756     }
 757
 758   my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
 759
 760   $car = $bar = $prd = 0;
 761   if (($dd = int($MBASE/($y->[-1]+1))) != 1)
 762     {
 763     for $xi (@$x)
 764       {
 765       $xi = $xi * $dd + $car;
 766       $xi -= ($car = int($xi / $MBASE)) * $MBASE;
 767       }
 768     push(@$x, $car); $car = 0;
 769     for $yi (@$y)
 770       {
 771       $yi = $yi * $dd + $car;
 772       $yi -= ($car = int($yi / $MBASE)) * $MBASE;
 773       }
 774     }
 775   else
 776     {
 777     push(@$x, 0);
 778     }
 779   @q = (); ($v2,$v1) = @$y[-2,-1];
 780   $v2 = 0 unless $v2;
 781   while ($#$x > $#$y)
 782     {
 783     ($u2,$u1,$u0) = @$x[-3..-1];
 784     $u2 = 0 unless $u2;
 785     #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
 786     # if $v1 == 0;
 787      $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
 788     --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
 789     if ($q)
 790       {
 791       ($car, $bar) = (0,0);
 792       for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
 793         {
 794         $prd = $q * $y->[$yi] + $car;
 795         $prd -= ($car = int($prd / $MBASE)) * $MBASE;
 796         $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
 797         }
 798       if ($x->[-1] < $car + $bar)
 799         {
 800         $car = 0; --$q;
 801         for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
 802           {
 803           $x->[$xi] -= $MBASE
 804            if ($car = (($x->[$xi] += $y->[$yi] + $car) > $MBASE));
 805           }
 806         }
 807       }
 808     pop(@$x); unshift(@q, $q);
 809     }
 810   if (wantarray)
 811     {
 812     @d = ();
 813     if ($dd != 1)
 814       {
 815       $car = 0;
 816       for $xi (reverse @$x)
 817         {
 818         $prd = $car * $MBASE + $xi;
 819         $car = $prd - ($tmp = int($prd / $dd)) * $dd;
 820         unshift(@d, $tmp);
 821         }
 822       }
 823     else
 824       {
 825       @d = @$x;
 826       }
 827     @$x = @q;
 828     my $d = \@d;
 829     if ($LEN_CONVERT != 0)
 830       {
 831       $c->_to_large($x); $c->_to_large($d);
 832       }
 833     else
 834       {
 835       __strip_zeros($x);
 836       __strip_zeros($d);
 837       }
 838     return ($x,$d);
 839     }
 840   @$x = @q;
 841   if ($LEN_CONVERT != 0)
 842     {
 843     $c->_to_large($x);
 844     }
 845   else
 846     {
 847     __strip_zeros($x);
 848     }
 849   $x;
 850   }
 851
 852 ##############################################################################
 853 # testing
 854
 855 sub _acmp
 856   {
 857   # internal absolute post-normalized compare (ignore signs)
 858   # ref to array, ref to array, return <0, 0, >0
 859   # arrays must have at least one entry; this is not checked for
 860
 861   my ($c,$cx,$cy) = @_;
 862
 863   # fast comp based on number of array elements (aka pseudo-length)
 864   my $lxy = scalar @$cx - scalar @$cy;
 865   return -1 if $lxy < 0;                                # already differs, ret
 866   return 1 if $lxy > 0;                                 # ditto
 867
 868   # now calculate length based on digits, not parts
 869   # we need only the length of the last element, since both array have the
 870   # same number of parts
 871   $lxy = length(int($cx->[-1])) - length(int($cy->[-1]));
 872   return -1 if $lxy < 0;
 873   return 1 if $lxy > 0;
 874
 875   # hm, same lengths,  but same contents? So we need to check all parts:
 876   my $a; my $j = scalar @$cx - 1;
 877   # manual way (abort if unequal, good for early ne)
 878   while ($j >= 0)
 879     {
 880     last if ($a = $cx->[$j] - $cy->[$j]); $j--;
 881     }
 882   return 1 if $a > 0;
 883   return -1 if $a < 0;
 884   0;                                            # numbers are equal
 885   }
 886
 887 sub _len
 888   {
 889   # compute number of digits
 890
 891   # int() because add/sub sometimes leaves strings (like '00005') instead of
 892   # '5' in this place, thus causing length() to report wrong length
 893   my $cx = $_[1];
 894
 895   (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
 896   }
 897
 898 sub _digit
 899   {
 900   # return the nth digit, negative values count backward
 901   # zero is rightmost, so _digit(123,0) will give 3
 902   my ($c,$x,$n) = @_;
 903
 904   my $len = _len('',$x);
 905
 906   $n = $len+$n if $n < 0;               # -1 last, -2 second-to-last
 907   $n = abs($n);                         # if negative was too big
 908   $len--; $n = $len if $n > $len;       # n to big?
 909
 910   my $elem = int($n / $BASE_LEN);       # which array element
 911   my $digit = $n % $BASE_LEN;           # which digit in this element
 912   $elem = '0000'.@$x[$elem];            # get element padded with 0's
 913   substr($elem,-$digit-1,1);
 914   }
 915
 916 sub _zeros
 917   {
 918   # return amount of trailing zeros in decimal
 919   # check each array elem in _m for having 0 at end as long as elem == 0
 920   # Upon finding a elem != 0, stop
 921   my $x = $_[1];
 922   my $zeros = 0; my $elem;
 923   foreach my $e (@$x)
 924     {
 925     if ($e != 0)
 926       {
 927       $elem = "$e";                             # preserve x
 928       $elem =~ s/.*?(0*$)/$1/;                  # strip anything not zero
 929       $zeros *= $BASE_LEN;                      # elems * 5
 930       $zeros += length($elem);                  # count trailing zeros
 931       last;                                     # early out
 932       }
 933     $zeros ++;                                  # real else branch: 50% slower!
 934     }
 935   $zeros;
 936   }
 937
 938 ##############################################################################
 939 # _is_* routines
 940
 941 sub _is_zero
 942   {
 943   # return true if arg (BINT or num_str) is zero (array '+', '0')
 944   my $x = $_[1];
 945
 946   (((scalar @$x == 1) && ($x->[0] == 0))) <=> 0;
 947   }
 948
 949 sub _is_even
 950   {
 951   # return true if arg (BINT or num_str) is even
 952   my $x = $_[1];
 953   (!($x->[0] & 1)) <=> 0;
 954   }
 955
 956 sub _is_odd
 957   {
 958   # return true if arg (BINT or num_str) is even
 959   my $x = $_[1];
 960
 961   (($x->[0] & 1)) <=> 0;
 962   }
 963
 964 sub _is_one
 965   {
 966   # return true if arg (BINT or num_str) is one (array '+', '1')
 967   my $x = $_[1];
 968
 969   (scalar @$x == 1) && ($x->[0] == 1) <=> 0;
 970   }
 971
 972 sub __strip_zeros
 973   {
 974   # internal normalization function that strips leading zeros from the array
 975   # args: ref to array
 976   my $s = shift;
 977
 978   my $cnt = scalar @$s; # get count of parts
 979   my $i = $cnt-1;
 980   push @$s,0 if $i < 0;         # div might return empty results, so fix it
 981
 982   return $s if @$s == 1;                # early out
 983
 984   #print "strip: cnt $cnt i $i\n";
 985   # '0', '3', '4', '0', '0',
 986   #  0    1    2    3    4
 987   # cnt = 5, i = 4
 988   # i = 4
 989   # i = 3
 990   # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
 991   # >= 1: skip first part (this can be zero)
 992   while ($i > 0) { last if $s->[$i] != 0; $i--; }
 993   $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
 994   $s;
 995   }
 996
 997 ###############################################################################
 998 # check routine to test internal state of corruptions
 999
1000 sub _check
1001   {
1002   # used by the test suite
1003   my $x = $_[1];
1004
1005   return "$x is not a reference" if !ref($x);
1006
1007   # are all parts are valid?
1008   my $i = 0; my $j = scalar @$x; my ($e,$try);
1009   while ($i < $j)
1010     {
1011     $e = $x->[$i]; $e = 'undef' unless defined $e;
1012     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
1013     last if $e !~ /^[+]?[0-9]+$/;
1014     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
1015     last if "$e" !~ /^[+]?[0-9]+$/;
1016     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
1017     last if '' . "$e" !~ /^[+]?[0-9]+$/;
1018     $try = ' < 0 || >= $BASE; '."($x, $e)";
1019     last if $e <0 || $e >= $BASE;
1020     # this test is disabled, since new/bnorm and certain ops (like early out
1021     # in add/sub) are allowed/expected to leave '00000' in some elements
1022     #$try = '=~ /^00+/; '."($x, $e)";
1023     #last if $e =~ /^00+/;
1024     $i++;
1025     }
1026   return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
1027   return 0;
1028   }
1029
1030
1031 ###############################################################################
1032 ###############################################################################
1033 # some optional routines to make BigInt faster
1034
1035 sub _mod
1036   {
1037   # if possible, use mod shortcut
1038   my ($c,$x,$yo) = @_;
1039
1040   # slow way since $y to big
1041   if (scalar @$yo > 1)
1042     {
1043     my ($xo,$rem) = _div($c,$x,$yo);
1044     return $rem;
1045     }
1046   my $y = $yo->[0];
1047   # both are single element arrays
1048   if (scalar @$x == 1)
1049     {
1050     $x->[0] %= $y;
1051     return $x;
1052     }
1053
1054   # @y is single element, but @x has more than one
1055   my $b = $BASE % $y;
1056   if ($b == 0)
1057     {
1058     # when BASE % Y == 0 then (B * BASE) % Y == 0
1059     # (B * BASE) % $y + A % Y => A % Y
1060     # so need to consider only last element: O(1)
1061     $x->[0] %= $y;
1062     }
1063   elsif ($b == 1)
1064     {
1065     # else need to go trough all elements: O(N), but loop is a bit simplified
1066     my $r = 0;
1067     foreach (@$x)
1068       {
1069       $r = ($r + $_) % $y;              # not much faster, but heh...
1070       #$r += $_ % $y; $r %= $y;
1071       }
1072     $r = 0 if $r == $y;
1073     $x->[0] = $r;
1074     }
1075   else
1076     {
1077     # else need to go trough all elements: O(N)
1078     my $r = 0; my $bm = 1;
1079     foreach (@$x)
1080       {
1081       $r = ($_ * $bm + $r) % $y;
1082       $bm = ($bm * $b) % $y;
1083
1084       #$r += ($_ % $y) * $bm;
1085       #$bm *= $b;
1086       #$bm %= $y;
1087       #$r %= $y;
1088       }
1089     $r = 0 if $r == $y;
1090     $x->[0] = $r;
1091     }
1092   splice (@$x,1);
1093   $x;
1094   }
1095
1096 ##############################################################################
1097 # shifts
1098
1099 sub _rsft
1100   {
1101   my ($c,$x,$y,$n) = @_;
1102
1103   if ($n != 10)
1104     {
1105     $n = _new($c,\$n); return _div($c,$x, _pow($c,$n,$y));
1106     }
1107
1108   # shortcut (faster) for shifting by 10)
1109   # multiples of $BASE_LEN
1110   my $dst = 0;                          # destination
1111   my $src = _num($c,$y);                # as normal int
1112   my $xlen = (@$x-1)*$BASE_LEN+length(int($x->[-1]));  # len of x in digits
1113   if ($src > $xlen)
1114     {
1115     # 12345 67890 shifted right by more than 10 digits => 0
1116     splice (@$x,1);                    # leave only one element
1117     $x->[0] = 0;                       # set to zero
1118     return $x;
1119     }
1120   my $rem = $src % $BASE_LEN;           # remainder to shift
1121   $src = int($src / $BASE_LEN);         # source
1122   if ($rem == 0)
1123     {
1124     splice (@$x,0,$src);                # even faster, 38.4 => 39.3
1125     }
1126   else
1127     {
1128     my $len = scalar @$x - $src;        # elems to go
1129     my $vd; my $z = '0'x $BASE_LEN;
1130     $x->[scalar @$x] = 0;               # avoid || 0 test inside loop
1131     while ($dst < $len)
1132       {
1133       $vd = $z.$x->[$src];
1134       $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
1135       $src++;
1136       $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
1137       $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1138       $x->[$dst] = int($vd);
1139       $dst++;
1140       }
1141     splice (@$x,$dst) if $dst > 0;              # kill left-over array elems
1142     pop @$x if $x->[-1] == 0 && @$x > 1;        # kill last element if 0
1143     } # else rem == 0
1144   $x;
1145   }
1146
1147 sub _lsft
1148   {
1149   my ($c,$x,$y,$n) = @_;
1150
1151   if ($n != 10)
1152     {
1153     $n = _new($c,\$n); return _mul($c,$x, _pow($c,$n,$y));
1154     }
1155
1156   # shortcut (faster) for shifting by 10) since we are in base 10eX
1157   # multiples of $BASE_LEN:
1158   my $src = scalar @$x;                 # source
1159   my $len = _num($c,$y);                # shift-len as normal int
1160   my $rem = $len % $BASE_LEN;           # remainder to shift
1161   my $dst = $src + int($len/$BASE_LEN); # destination
1162   my $vd;                               # further speedup
1163   $x->[$src] = 0;                       # avoid first ||0 for speed
1164   my $z = '0' x $BASE_LEN;
1165   while ($src >= 0)
1166     {
1167     $vd = $x->[$src]; $vd = $z.$vd;
1168     $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
1169     $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
1170     $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1171     $x->[$dst] = int($vd);
1172     $dst--; $src--;
1173     }
1174   # set lowest parts to 0
1175   while ($dst >= 0) { $x->[$dst--] = 0; }
1176   # fix spurios last zero element
1177   splice @$x,-1 if $x->[-1] == 0;
1178   $x;
1179   }
1180
1181 sub _pow
1182   {
1183   # power of $x to $y
1184   # ref to array, ref to array, return ref to array
1185   my ($c,$cx,$cy) = @_;
1186
1187   my $pow2 = _one();
1188
1189   my $y_bin = ${_as_bin($c,$cy)}; $y_bin =~ s/^0b//;
1190   my $len = length($y_bin);
1191   while (--$len > 0)
1192     {
1193     _mul($c,$pow2,$cx) if substr($y_bin,$len,1) eq '1';         # is odd?
1194     _mul($c,$cx,$cx);
1195     }
1196
1197   _mul($c,$cx,$pow2);
1198   $cx;
1199   }
1200
1201 sub _fac
1202   {
1203   # factorial of $x
1204   # ref to array, return ref to array
1205   my ($c,$cx) = @_;
1206
1207   if ((@$cx == 1) && ($cx->[0] <= 2))
1208     {
1209     $cx->[0] = 1 * ($cx->[0]||1); # 0,1 => 1, 2 => 2
1210     return $cx;
1211     }
1212
1213   # go forward until $base is exceeded
1214   # limit is either $x or $base (x == 100 means as result too high)
1215   my $steps = 100; $steps = $cx->[0] if @$cx == 1;
1216   my $r = 2; my $cf = 3; my $step = 1; my $last = $r;
1217   while ($r < $BASE && $step < $steps)
1218     {
1219     $last = $r; $r *= $cf++; $step++;
1220     }
1221   if ((@$cx == 1) && ($step == $cx->[0]))
1222     {
1223     # completely done
1224     $cx = [$last];
1225     return $cx;
1226     }
1227   my $n = _copy($c,$cx);
1228   $cx = [$last];
1229
1230   while (!(@$n == 1 && $n->[0] == $step))
1231     {
1232     _mul($c,$cx,$n); _dec($c,$n);
1233     }
1234   $cx;
1235   }
1236
1237 # for debugging:
1238   use constant DEBUG => 0;
1239   my $steps = 0;
1240   sub steps { $steps };
1241
1242 sub _sqrt
1243   {
1244   # square-root of $x in place
1245   # Compute a guess of the result (rule of thumb), then improve it via
1246   # Newton's method.
1247   my ($c,$x) = @_;
1248
1249   if (scalar @$x == 1)
1250     {
1251     # fit's into one Perl scalar, so result can be computed directly
1252     $x->[0] = int(sqrt($x->[0]));
1253     return $x;
1254     }
1255   my $y = _copy($c,$x);
1256   # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
1257   # since our guess will "grow"
1258   my $l = int((_len($c,$x)-1) / 2);
1259
1260   my $lastelem = $x->[-1];                                      # for guess
1261   my $elems = scalar @$x - 1;
1262   # not enough digits, but could have more?
1263   if ((length($lastelem) <= 3) && ($elems > 1))
1264     {
1265     # right-align with zero pad
1266     my $len = length($lastelem) & 1;
1267     print "$lastelem => " if DEBUG;
1268     $lastelem .= substr($x->[-2] . '0' x $BASE_LEN,0,$BASE_LEN);
1269     # former odd => make odd again, or former even to even again
1270     $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
1271     print "$lastelem\n" if DEBUG;
1272     }
1273
1274   # construct $x (instead of _lsft($c,$x,$l,10)
1275   my $r = $l % $BASE_LEN;       # 10000 00000 00000 00000 ($BASE_LEN=5)
1276   $l = int($l / $BASE_LEN);
1277   print "l =  $l " if DEBUG;
1278
1279   splice @$x,$l;                # keep ref($x), but modify it
1280
1281   # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
1282   # that gives us:
1283   # 14400 00000 => sqrt(14400) => guess first digits to be 120
1284   # 144000 000000 => sqrt(144000) => guess 379
1285
1286   print "$lastelem (elems $elems) => " if DEBUG;
1287   $lastelem = $lastelem / 10 if ($elems & 1 == 1);              # odd or even?
1288   my $g = sqrt($lastelem); $g =~ s/\.//;                        # 2.345 => 2345
1289   $r -= 1 if $elems & 1 == 0;                                   # 70 => 7
1290
1291   # padd with zeros if result is too short
1292   $x->[$l--] = int(substr($g . '0' x $r,0,$r+1));
1293   print "now ",$x->[-1] if DEBUG;
1294   print " would have been ", int('1' . '0' x $r),"\n" if DEBUG;
1295
1296   # If @$x > 1, we could compute the second elem of the guess, too, to create
1297   # an even better guess. Not implemented yet. Does it improve performance?
1298   $x->[$l--] = 0 while ($l >= 0);       # all other digits of guess are zero
1299
1300   print "start x= ",${_str($c,$x)},"\n" if DEBUG;
1301   my $two = _two();
1302   my $last = _zero();
1303   my $lastlast = _zero();
1304   $steps = 0 if DEBUG;
1305   while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
1306     {
1307     $steps++ if DEBUG;
1308     $lastlast = _copy($c,$last);
1309     $last = _copy($c,$x);
1310     _add($c,$x, _div($c,_copy($c,$y),$x));
1311     _div($c,$x, $two );
1312     print " x= ",${_str($c,$x)},"\n" if DEBUG;
1313     }
1314   print "\nsteps in sqrt: $steps, " if DEBUG;
1315   _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0;     # overshot?
1316   print " final ",$x->[-1],"\n" if DEBUG;
1317   $x;
1318   }
1319
1320 ##############################################################################
1321 # binary stuff
1322
1323 sub _and
1324   {
1325   my ($c,$x,$y) = @_;
1326
1327   # the shortcut makes equal, large numbers _really_ fast, and makes only a
1328   # very small performance drop for small numbers (e.g. something with less
1329   # than 32 bit) Since we optimize for large numbers, this is enabled.
1330   return $x if _acmp($c,$x,$y) == 0;            # shortcut
1331
1332   my $m = _one(); my ($xr,$yr);
1333   my $mask = $AND_MASK;
1334
1335   my $x1 = $x;
1336   my $y1 = _copy($c,$y);                        # make copy
1337   $x = _zero();
1338   my ($b,$xrr,$yrr);
1339   use integer;
1340   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1341     {
1342     ($x1, $xr) = _div($c,$x1,$mask);
1343     ($y1, $yr) = _div($c,$y1,$mask);
1344
1345     # make ints() from $xr, $yr
1346     # this is when the AND_BITS are greater tahn $BASE and is slower for
1347     # small (<256 bits) numbers, but faster for large numbers. Disabled
1348     # due to KISS principle
1349
1350 #    $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1351 #    $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1352 #    _add($c,$x, _mul($c, _new( $c, \($xrr & $yrr) ), $m) );
1353
1354     # 0+ due to '&' doesn't work in strings
1355     _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
1356     _mul($c,$m,$mask);
1357     }
1358   $x;
1359   }
1360
1361 sub _xor
1362   {
1363   my ($c,$x,$y) = @_;
1364
1365   return _zero() if _acmp($c,$x,$y) == 0;       # shortcut (see -and)
1366
1367   my $m = _one(); my ($xr,$yr);
1368   my $mask = $XOR_MASK;
1369
1370   my $x1 = $x;
1371   my $y1 = _copy($c,$y);                        # make copy
1372   $x = _zero();
1373   my ($b,$xrr,$yrr);
1374   use integer;
1375   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1376     {
1377     ($x1, $xr) = _div($c,$x1,$mask);
1378     ($y1, $yr) = _div($c,$y1,$mask);
1379     # make ints() from $xr, $yr (see _and())
1380     #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1381     #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1382     #_add($c,$x, _mul($c, _new( $c, \($xrr ^ $yrr) ), $m) );
1383
1384     # 0+ due to '^' doesn't work in strings
1385     _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
1386     _mul($c,$m,$mask);
1387     }
1388   # the loop stops when the shorter of the two numbers is exhausted
1389   # the remainder of the longer one will survive bit-by-bit, so we simple
1390   # multiply-add it in
1391   _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1392   _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1393
1394   $x;
1395   }
1396
1397 sub _or
1398   {
1399   my ($c,$x,$y) = @_;
1400
1401   return $x if _acmp($c,$x,$y) == 0;            # shortcut (see _and)
1402
1403   my $m = _one(); my ($xr,$yr);
1404   my $mask = $OR_MASK;
1405
1406   my $x1 = $x;
1407   my $y1 = _copy($c,$y);                        # make copy
1408   $x = _zero();
1409   my ($b,$xrr,$yrr);
1410   use integer;
1411   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1412     {
1413     ($x1, $xr) = _div($c,$x1,$mask);
1414     ($y1, $yr) = _div($c,$y1,$mask);
1415     # make ints() from $xr, $yr (see _and())
1416 #    $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1417 #    $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1418 #    _add($c,$x, _mul($c, _new( $c, \($xrr | $yrr) ), $m) );
1419
1420     # 0+ due to '|' doesn't work in strings
1421     _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
1422     _mul($c,$m,$mask);
1423     }
1424   # the loop stops when the shorter of the two numbers is exhausted
1425   # the remainder of the longer one will survive bit-by-bit, so we simple
1426   # multiply-add it in
1427   _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1428   _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1429
1430   $x;
1431   }
1432
1433 sub _as_hex
1434   {
1435   # convert a decimal number to hex (ref to array, return ref to string)
1436   my ($c,$x) = @_;
1437
1438   my $x1 = _copy($c,$x);
1439
1440   my $es = '';
1441   my ($xr, $h, $x10000);
1442   if ($] >= 5.006)
1443     {
1444     $x10000 = [ 0x10000 ]; $h = 'h4';
1445     }
1446   else
1447     {
1448     $x10000 = [ 0x1000 ]; $h = 'h3';
1449     }
1450   while (! _is_zero($c,$x1))
1451     {
1452     ($x1, $xr) = _div($c,$x1,$x10000);
1453     $es .= unpack($h,pack('v',$xr->[0]));
1454     }
1455   $es = reverse $es;
1456   $es =~ s/^[0]+//;   # strip leading zeros
1457   $es = '0x' . $es;
1458   \$es;
1459   }
1460
1461 sub _as_bin
1462   {
1463   # convert a decimal number to bin (ref to array, return ref to string)
1464   my ($c,$x) = @_;
1465
1466   my $x1 = _copy($c,$x);
1467
1468   my $es = '';
1469   my ($xr, $b, $x10000);
1470   if ($] >= 5.006)
1471     {
1472     $x10000 = [ 0x10000 ]; $b = 'b16';
1473     }
1474   else
1475     {
1476     $x10000 = [ 0x1000 ]; $b = 'b12';
1477     }
1478   while (! _is_zero($c,$x1))
1479     {
1480     ($x1, $xr) = _div($c,$x1,$x10000);
1481     $es .= unpack($b,pack('v',$xr->[0]));
1482     }
1483   $es = reverse $es;
1484   $es =~ s/^[0]+//;   # strip leading zeros
1485   $es = '0b' . $es;
1486   \$es;
1487   }
1488
1489 sub _from_hex
1490   {
1491   # convert a hex number to decimal (ref to string, return ref to array)
1492   my ($c,$hs) = @_;
1493
1494   my $mul = _one();
1495   my $m = [ 0x10000 ];                          # 16 bit at a time
1496   my $x = _zero();
1497
1498   my $len = length($$hs)-2;
1499   $len = int($len/4);                           # 4-digit parts, w/o '0x'
1500   my $val; my $i = -4;
1501   while ($len >= 0)
1502     {
1503     $val = substr($$hs,$i,4);
1504     $val =~ s/^[+-]?0x// if $len == 0;          # for last part only because
1505     $val = hex($val);                           # hex does not like wrong chars
1506     $i -= 4; $len --;
1507     _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1508     _mul ($c, $mul, $m ) if $len >= 0;          # skip last mul
1509     }
1510   $x;
1511   }
1512
1513 sub _from_bin
1514   {
1515   # convert a hex number to decimal (ref to string, return ref to array)
1516   my ($c,$bs) = @_;
1517
1518   # instead of converting 8 bit at a time, it is faster to convert the
1519   # number to hex, and then call _from_hex.
1520
1521   my $hs = $$bs;
1522   $hs =~ s/^[+-]?0b//;                                  # remove sign and 0b
1523   my $l = length($hs);                                  # bits
1524   $hs = '0' x (8-($l % 8)) . $hs if ($l % 8) != 0;      # padd left side w/ 0
1525   my $h = unpack('H*', pack ('B*', $hs));               # repack as hex
1526   return $c->_from_hex(\('0x'.$h));
1527
1528   my $mul = _one();
1529   my $m = [ 0x100 ];                            # 8 bit at a time
1530   my $x = _zero();
1531
1532   my $len = length($$bs)-2;
1533   $len = int($len/8);                           # 4-digit parts, w/o '0x'
1534   my $val; my $i = -8;
1535   while ($len >= 0)
1536     {
1537     $val = substr($$bs,$i,8);
1538     $val =~ s/^[+-]?0b// if $len == 0;          # for last part only
1539
1540     $val = ord(pack('B8',substr('00000000'.$val,-8,8)));
1541
1542     $i -= 8; $len --;
1543     _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1544     _mul ($c, $mul, $m ) if $len >= 0;          # skip last mul
1545     }
1546   $x;
1547   }
1548
1549 ##############################################################################
1550 # special modulus functions
1551
1552 sub _modinv
1553   {
1554   # modular inverse
1555   my ($c,$x,$y) = @_;
1556
1557   my $u = _zero($c); my $u1 = _one($c);
1558   my $a = _copy($c,$y); my $b = _copy($c,$x);
1559
1560   # Euclid's Algorithm for bgcd(), only that we calc bgcd() ($a) and the
1561   # result ($u) at the same time. See comments in BigInt for why this works.
1562   my $q;
1563   ($a, $q, $b) = ($b, _div($c,$a,$b));          # step 1
1564   my $sign = 1;
1565   while (!_is_zero($c,$b))
1566     {
1567     my $t = _add($c,                            # step 2:
1568        _mul($c,_copy($c,$u1), $q) ,             #  t =  u1 * q
1569        $u );                                    #     + u
1570     $u = $u1;                                   #  u = u1, u1 = t
1571     $u1 = $t;
1572     $sign = -$sign;
1573     ($a, $q, $b) = ($b, _div($c,$a,$b));        # step 1
1574     }
1575
1576   # if the gcd is not 1, then return NaN
1577   return (undef,undef) unless _is_one($c,$a);
1578
1579   $sign = $sign == 1 ? '+' : '-';
1580   ($u1,$sign);
1581   }
1582
1583 sub _modpow
1584   {
1585   # modulus of power ($x ** $y) % $z
1586   my ($c,$num,$exp,$mod) = @_;
1587
1588   # in the trivial case,
1589   if (_is_one($c,$mod))
1590     {
1591     splice @$num,0,1; $num->[0] = 0;
1592     return $num;
1593     }
1594   if ((scalar @$num == 1) && (($num->[0] == 0) || ($num->[0] == 1)))
1595     {
1596     $num->[0] = 1;
1597     return $num;
1598     }
1599
1600 #  $num = _mod($c,$num,$mod);   # this does not make it faster
1601
1602   my $acc = _copy($c,$num); my $t = _one();
1603
1604   my $expbin = ${_as_bin($c,$exp)}; $expbin =~ s/^0b//;
1605   my $len = length($expbin);
1606   while (--$len >= 0)
1607     {
1608     if ( substr($expbin,$len,1) eq '1')                 # is_odd
1609       {
1610       _mul($c,$t,$acc);
1611       $t = _mod($c,$t,$mod);
1612       }
1613     _mul($c,$acc,$acc);
1614     $acc = _mod($c,$acc,$mod);
1615     }
1616   @$num = @$t;
1617   $num;
1618   }
1619
1620 ##############################################################################
1621 ##############################################################################
1622
1623 1;
1624 __END__
1625
1626 =head1 NAME
1627
1628 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
1629
1630 =head1 SYNOPSIS
1631
1632 Provides support for big integer calculations. Not intended to be used by other
1633 modules (except Math::BigInt::Cached). Other modules which sport the same
1634 functions can also be used to support Math::Bigint, like Math::BigInt::Pari.
1635
1636 =head1 DESCRIPTION
1637
1638 In order to allow for multiple big integer libraries, Math::BigInt was
1639 rewritten to use library modules for core math routines. Any module which
1640 follows the same API as this can be used instead by using the following:
1641
1642         use Math::BigInt lib => 'libname';
1643
1644 'libname' is either the long name ('Math::BigInt::Pari'), or only the short
1645 version like 'Pari'.
1646
1647 =head1 EXPORT
1648
1649 The following functions MUST be defined in order to support the use by
1650 Math::BigInt:
1651
1652         _new(string)    return ref to new object from ref to decimal string
1653         _zero()         return a new object with value 0
1654         _one()          return a new object with value 1
1655
1656         _str(obj)       return ref to a string representing the object
1657         _num(obj)       returns a Perl integer/floating point number
1658                         NOTE: because of Perl numeric notation defaults,
1659                         the _num'ified obj may lose accuracy due to
1660                         machine-dependend floating point size limitations
1661
1662         _add(obj,obj)   Simple addition of two objects
1663         _mul(obj,obj)   Multiplication of two objects
1664         _div(obj,obj)   Division of the 1st object by the 2nd
1665                         In list context, returns (result,remainder).
1666                         NOTE: this is integer math, so no
1667                         fractional part will be returned.
1668         _sub(obj,obj)   Simple subtraction of 1 object from another
1669                         a third, optional parameter indicates that the params
1670                         are swapped. In this case, the first param needs to
1671                         be preserved, while you can destroy the second.
1672                         sub (x,y,1) => return x - y and keep x intact!
1673         _dec(obj)       decrement object by one (input is garant. to be > 0)
1674         _inc(obj)       increment object by one
1675
1676
1677         _acmp(obj,obj)  <=> operator for objects (return -1, 0 or 1)
1678
1679         _len(obj)       returns count of the decimal digits of the object
1680         _digit(obj,n)   returns the n'th decimal digit of object
1681
1682         _is_one(obj)    return true if argument is +1
1683         _is_zero(obj)   return true if argument is 0
1684         _is_even(obj)   return true if argument is even (0,2,4,6..)
1685         _is_odd(obj)    return true if argument is odd (1,3,5,7..)
1686
1687         _copy           return a ref to a true copy of the object
1688
1689         _check(obj)     check whether internal representation is still intact
1690                         return 0 for ok, otherwise error message as string
1691
1692 The following functions are optional, and can be defined if the underlying lib
1693 has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
1694 slow) fallback routines to emulate these:
1695
1696         _from_hex(str)  return ref to new object from ref to hexadecimal string
1697         _from_bin(str)  return ref to new object from ref to binary string
1698
1699         _as_hex(str)    return ref to scalar string containing the value as
1700                         unsigned hex string, with the '0x' prepended.
1701                         Leading zeros must be stripped.
1702         _as_bin(str)    Like as_hex, only as binary string containing only
1703                         zeros and ones. Leading zeros must be stripped and a
1704                         '0b' must be prepended.
1705
1706         _rsft(obj,N,B)  shift object in base B by N 'digits' right
1707                         For unsupported bases B, return undef to signal failure
1708         _lsft(obj,N,B)  shift object in base B by N 'digits' left
1709                         For unsupported bases B, return undef to signal failure
1710
1711         _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
1712                         Note: XOR, AND and OR pad with zeros if size mismatches
1713         _and(obj1,obj2) AND (bit-wise) object 1 with object 2
1714         _or(obj1,obj2)  OR (bit-wise) object 1 with object 2
1715
1716         _mod(obj,obj)   Return remainder of div of the 1st by the 2nd object
1717         _sqrt(obj)      return the square root of object (truncate to int)
1718         _fac(obj)       return factorial of object 1 (1*2*3*4..)
1719         _pow(obj,obj)   return object 1 to the power of object 2
1720         _gcd(obj,obj)   return Greatest Common Divisor of two objects
1721
1722         _zeros(obj)     return number of trailing decimal zeros
1723         _modinv         return inverse modulus
1724         _modpow         return modulus of power ($x ** $y) % $z
1725
1726 Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
1727 or '0b1101').
1728
1729 Testing of input parameter validity is done by the caller, so you need not
1730 worry about underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by
1731 zero or similar cases.
1732
1733 The first parameter can be modified, that includes the possibility that you
1734 return a reference to a completely different object instead. Although keeping
1735 the reference and just changing it's contents is prefered over creating and
1736 returning a different reference.
1737
1738 Return values are always references to objects or strings. Exceptions are
1739 C<_lsft()> and C<_rsft()>, which return undef if they can not shift the
1740 argument. This is used to delegate shifting of bases different than the one
1741 you can support back to Math::BigInt, which will use some generic code to
1742 calculate the result.
1743
1744 =head1 WRAP YOUR OWN
1745
1746 If you want to port your own favourite c-lib for big numbers to the
1747 Math::BigInt interface, you can take any of the already existing modules as
1748 a rough guideline. You should really wrap up the latest BigInt and BigFloat
1749 testsuites with your module, and replace in them any of the following:
1750
1751         use Math::BigInt;
1752
1753 by this:
1754
1755         use Math::BigInt lib => 'yourlib';
1756
1757 This way you ensure that your library really works 100% within Math::BigInt.
1758
1759 =head1 LICENSE
1760
1761 This program is free software; you may redistribute it and/or modify it under
1762 the same terms as Perl itself.
1763
1764 =head1 AUTHORS
1765
1766 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
1767 in late 2000, 2001.
1768 Seperated from BigInt and shaped API with the help of John Peacock.
1769
1770 =head1 SEE ALSO
1771
1772 L<Math::BigInt>, L<Math::BigFloat>, L<Math::BigInt::BitVect>,
1773 L<Math::BigInt::GMP>, L<Math::BigInt::Cached> and L<Math::BigInt::Pari>.
1774
1775 =cut