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