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.26';
  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 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     # we can't do an early out if $x is than $y, since we
 414     # need to copy the high chunks from $y. Found by Bob Mathews.
 415     #last unless defined $sy->[$j] || $car;
 416     $sy->[$j] += $BASE
 417      if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
 418     $j++;
 419     }
 420   # might leave leading zeros, so fix that
 421   __strip_zeros($sy);
 422   }
 423
 424 sub _square_use_mul
 425   {
 426   # compute $x ** 2 or $x * $x in-place and return $x
 427   my ($c,$x) = @_;
 428
 429   # From: Handbook of Applied Cryptography by A. Menezes, P. van Oorschot and
 430   #       S. Vanstone., Chapter 14
 431
 432   #14.16 Algorithm Multiple-precision squaring
 433   #INPUT: positive integer x = (xt 1 xt 2 ... x1 x0)b.
 434   #OUTPUT: x * x = x ** 2 in radix b representation.
 435   #1. For i from 0 to (2t - 1) do: wi <- 0.
 436   #2.  For i from 0 to (t - 1) do the following:
 437   # 2.1 (uv)b w2i + xi * xi, w2i v, c u.
 438   # 2.2 For j from (i + 1)to (t - 1) do the following:
 439   #      (uv)b <- wi+j + 2*xj * xi + c, wi+j <- v, c <- u.
 440   # 2.3 wi+t <- u.
 441   #3. Return((w2t-1 w2t-2 ... w1 w0)b).
 442
 443 #  # Note: That description is crap. Half of the symbols are not explained or
 444 #  # used with out beeing set.
 445 #  my $t = scalar @$x;          # count
 446 #  my ($c,$i,$j);
 447 #  for ($i = 0; $i < $t; $i++)
 448 #    {
 449 #    $x->[$i] = $x->[$i*2] + $x[$i]*$x[$i];
 450 #    $x->[$i*2] = $x[$i]; $c = $x[$i];
 451 #    for ($j = $i+1; $j < $t; $j++)
 452 #      {
 453 #      $x->[$i] = $x->[$i+$j] + 2 * $x->[$i] * $x->[$j];
 454 #      $x->[$i+$j] = $x[$j]; $c = $x[$i];
 455 #      }
 456 #    $x->[$i+$t] = $x[$i];
 457 #    }
 458   $x;
 459   }
 460
 461 sub _mul_use_mul
 462   {
 463   # (ref to int_num_array, ref to int_num_array)
 464   # multiply two numbers in internal representation
 465   # modifies first arg, second need not be different from first
 466   my ($c,$xv,$yv) = @_;
 467
 468   # shortcut for two very short numbers (improved by Nathan Zook)
 469   # works also if xv and yv are the same reference
 470   if ((@$xv == 1) && (@$yv == 1))
 471     {
 472     if (($xv->[0] *= $yv->[0]) >= $MBASE)
 473        {
 474        $xv->[0] = $xv->[0] - ($xv->[1] = int($xv->[0] * $RBASE)) * $MBASE;
 475        };
 476     return $xv;
 477     }
 478   # shortcut for result == 0
 479   if ( ((@$xv == 1) && ($xv->[0] == 0)) ||
 480        ((@$yv == 1) && ($yv->[0] == 0)) )
 481     {
 482     @$xv = (0);
 483     return $xv;
 484     }
 485
 486   # since multiplying $x with $x fails, make copy in this case
 487   $yv = [@$xv] if "$xv" eq "$yv";       # same references?
 488   # since multiplying $x with $x would fail here, use the faster squaring
 489 #  return _square($c,$xv) if "$xv" eq "$yv";    # same reference?
 490
 491   if ($LEN_CONVERT != 0)
 492     {
 493     $c->_to_small($xv); $c->_to_small($yv);
 494     }
 495
 496   my @prod = (); my ($prod,$car,$cty,$xi,$yi);
 497
 498   for $xi (@$xv)
 499     {
 500     $car = 0; $cty = 0;
 501
 502     # slow variant
 503 #    for $yi (@$yv)
 504 #      {
 505 #      $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
 506 #      $prod[$cty++] =
 507 #       $prod - ($car = int($prod * RBASE)) * $MBASE;  # see USE_MUL
 508 #      }
 509 #    $prod[$cty] += $car if $car; # need really to check for 0?
 510 #    $xi = shift @prod;
 511
 512     # faster variant
 513     # looping through this if $xi == 0 is silly - so optimize it away!
 514     $xi = (shift @prod || 0), next if $xi == 0;
 515     for $yi (@$yv)
 516       {
 517       $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
 518 ##     this is actually a tad slower
 519 ##        $prod = $prod[$cty]; $prod += ($car + $xi * $yi);     # no ||0 here
 520       $prod[$cty++] =
 521        $prod - ($car = int($prod * $RBASE)) * $MBASE;  # see USE_MUL
 522       }
 523     $prod[$cty] += $car if $car; # need really to check for 0?
 524     $xi = shift @prod || 0;     # || 0 makes v5.005_3 happy
 525     }
 526   push @$xv, @prod;
 527   if ($LEN_CONVERT != 0)
 528     {
 529     $c->_to_large($yv);
 530     $c->_to_large($xv);
 531     }
 532   else
 533     {
 534     __strip_zeros($xv);
 535     }
 536   $xv;
 537   }
 538
 539 sub _mul_use_div
 540   {
 541   # (ref to int_num_array, ref to int_num_array)
 542   # multiply two numbers in internal representation
 543   # modifies first arg, second need not be different from first
 544   my ($c,$xv,$yv) = @_;
 545
 546   # shortcut for two very short numbers (improved by Nathan Zook)
 547   # works also if xv and yv are the same reference
 548   if ((@$xv == 1) && (@$yv == 1))
 549     {
 550     if (($xv->[0] *= $yv->[0]) >= $MBASE)
 551         {
 552         $xv->[0] =
 553             $xv->[0] - ($xv->[1] = int($xv->[0] / $MBASE)) * $MBASE;
 554         };
 555     return $xv;
 556     }
 557   # shortcut for result == 0
 558   if ( ((@$xv == 1) && ($xv->[0] == 0)) ||
 559        ((@$yv == 1) && ($yv->[0] == 0)) )
 560     {
 561     @$xv = (0);
 562     return $xv;
 563     }
 564
 565
 566   # since multiplying $x with $x fails, make copy in this case
 567   $yv = [@$xv] if "$xv" eq "$yv";       # same references?
 568   # since multiplying $x with $x would fail here, use the faster squaring
 569 #  return _square($c,$xv) if "$xv" eq "$yv";    # same reference?
 570
 571   if ($LEN_CONVERT != 0)
 572     {
 573     $c->_to_small($xv); $c->_to_small($yv);
 574     }
 575
 576   my @prod = (); my ($prod,$car,$cty,$xi,$yi);
 577   for $xi (@$xv)
 578     {
 579     $car = 0; $cty = 0;
 580     # looping through this if $xi == 0 is silly - so optimize it away!
 581     $xi = (shift @prod || 0), next if $xi == 0;
 582     for $yi (@$yv)
 583       {
 584       $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
 585       $prod[$cty++] =
 586        $prod - ($car = int($prod / $MBASE)) * $MBASE;
 587       }
 588     $prod[$cty] += $car if $car; # need really to check for 0?
 589     $xi = shift @prod || 0;     # || 0 makes v5.005_3 happy
 590     }
 591   push @$xv, @prod;
 592   if ($LEN_CONVERT != 0)
 593     {
 594     $c->_to_large($yv);
 595     $c->_to_large($xv);
 596     }
 597   else
 598     {
 599     __strip_zeros($xv);
 600     }
 601   $xv;
 602   }
 603
 604 sub _div_use_mul
 605   {
 606   # ref to array, ref to array, modify first array and return remainder if
 607   # in list context
 608   my ($c,$x,$yorg) = @_;
 609
 610   if (@$x == 1 && @$yorg == 1)
 611     {
 612     # shortcut, $yorg and $x are two small numbers
 613     if (wantarray)
 614       {
 615       my $r = [ $x->[0] % $yorg->[0] ];
 616       $x->[0] = int($x->[0] / $yorg->[0]);
 617       return ($x,$r);
 618       }
 619     else
 620       {
 621       $x->[0] = int($x->[0] / $yorg->[0]);
 622       return $x;
 623       }
 624     }
 625   if (@$yorg == 1)
 626     {
 627     my $rem;
 628     $rem = _mod($c,[ @$x ],$yorg) if wantarray;
 629
 630     # shortcut, $y is < $BASE
 631     my $j = scalar @$x; my $r = 0;
 632     my $y = $yorg->[0]; my $b;
 633     while ($j-- > 0)
 634       {
 635       $b = $r * $MBASE + $x->[$j];
 636       $x->[$j] = int($b/$y);
 637       $r = $b % $y;
 638       }
 639     pop @$x if @$x > 1 && $x->[-1] == 0;        # splice up a leading zero
 640     return ($x,$rem) if wantarray;
 641     return $x;
 642     }
 643
 644   my $y = [ @$yorg ];
 645   if ($LEN_CONVERT != 0)
 646     {
 647     $c->_to_small($x); $c->_to_small($y);
 648     }
 649
 650   my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
 651
 652   $car = $bar = $prd = 0;
 653   if (($dd = int($MBASE/($y->[-1]+1))) != 1)
 654     {
 655     for $xi (@$x)
 656       {
 657       $xi = $xi * $dd + $car;
 658       $xi -= ($car = int($xi * $RBASE)) * $MBASE;       # see USE_MUL
 659       }
 660     push(@$x, $car); $car = 0;
 661     for $yi (@$y)
 662       {
 663       $yi = $yi * $dd + $car;
 664       $yi -= ($car = int($yi * $RBASE)) * $MBASE;       # see USE_MUL
 665       }
 666     }
 667   else
 668     {
 669     push(@$x, 0);
 670     }
 671   @q = (); ($v2,$v1) = @$y[-2,-1];
 672   $v2 = 0 unless $v2;
 673   while ($#$x > $#$y)
 674     {
 675     ($u2,$u1,$u0) = @$x[-3..-1];
 676     $u2 = 0 unless $u2;
 677     #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
 678     # if $v1 == 0;
 679      $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
 680     --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
 681     if ($q)
 682       {
 683       ($car, $bar) = (0,0);
 684       for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
 685         {
 686         $prd = $q * $y->[$yi] + $car;
 687         $prd -= ($car = int($prd * $RBASE)) * $MBASE;   # see USE_MUL
 688         $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
 689         }
 690       if ($x->[-1] < $car + $bar)
 691         {
 692         $car = 0; --$q;
 693         for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
 694           {
 695           $x->[$xi] -= $MBASE
 696            if ($car = (($x->[$xi] += $y->[$yi] + $car) > $MBASE));
 697           }
 698         }
 699       }
 700       pop(@$x); unshift(@q, $q);
 701     }
 702   if (wantarray)
 703     {
 704     @d = ();
 705     if ($dd != 1)
 706       {
 707       $car = 0;
 708       for $xi (reverse @$x)
 709         {
 710         $prd = $car * $MBASE + $xi;
 711         $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
 712         unshift(@d, $tmp);
 713         }
 714       }
 715     else
 716       {
 717       @d = @$x;
 718       }
 719     @$x = @q;
 720     my $d = \@d;
 721     if ($LEN_CONVERT != 0)
 722       {
 723       $c->_to_large($x); $c->_to_large($d);
 724       }
 725     else
 726       {
 727       __strip_zeros($x);
 728       __strip_zeros($d);
 729       }
 730     return ($x,$d);
 731     }
 732   @$x = @q;
 733   if ($LEN_CONVERT != 0)
 734     {
 735     $c->_to_large($x);
 736     }
 737   else
 738     {
 739     __strip_zeros($x);
 740     }
 741   $x;
 742   }
 743
 744 sub _div_use_div
 745   {
 746   # ref to array, ref to array, modify first array and return remainder if
 747   # in list context
 748   my ($c,$x,$yorg) = @_;
 749
 750   if (@$x == 1 && @$yorg == 1)
 751     {
 752     # shortcut, $yorg and $x are two small numbers
 753     if (wantarray)
 754       {
 755       my $r = [ $x->[0] % $yorg->[0] ];
 756       $x->[0] = int($x->[0] / $yorg->[0]);
 757       return ($x,$r);
 758       }
 759     else
 760       {
 761       $x->[0] = int($x->[0] / $yorg->[0]);
 762       return $x;
 763       }
 764     }
 765   if (@$yorg == 1)
 766     {
 767     my $rem;
 768     $rem = _mod($c,[ @$x ],$yorg) if wantarray;
 769
 770     # shortcut, $y is < $BASE
 771     my $j = scalar @$x; my $r = 0;
 772     my $y = $yorg->[0]; my $b;
 773     while ($j-- > 0)
 774       {
 775       $b = $r * $MBASE + $x->[$j];
 776       $x->[$j] = int($b/$y);
 777       $r = $b % $y;
 778       }
 779     pop @$x if @$x > 1 && $x->[-1] == 0;        # splice up a leading zero
 780     return ($x,$rem) if wantarray;
 781     return $x;
 782     }
 783
 784   my $y = [ @$yorg ];
 785   if ($LEN_CONVERT != 0)
 786     {
 787     $c->_to_small($x); $c->_to_small($y);
 788     }
 789
 790   my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1,@d,$tmp,$q,$u2,$u1,$u0);
 791
 792   $car = $bar = $prd = 0;
 793   if (($dd = int($MBASE/($y->[-1]+1))) != 1)
 794     {
 795     for $xi (@$x)
 796       {
 797       $xi = $xi * $dd + $car;
 798       $xi -= ($car = int($xi / $MBASE)) * $MBASE;
 799       }
 800     push(@$x, $car); $car = 0;
 801     for $yi (@$y)
 802       {
 803       $yi = $yi * $dd + $car;
 804       $yi -= ($car = int($yi / $MBASE)) * $MBASE;
 805       }
 806     }
 807   else
 808     {
 809     push(@$x, 0);
 810     }
 811   @q = (); ($v2,$v1) = @$y[-2,-1];
 812   $v2 = 0 unless $v2;
 813   while ($#$x > $#$y)
 814     {
 815     ($u2,$u1,$u0) = @$x[-3..-1];
 816     $u2 = 0 unless $u2;
 817     #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
 818     # if $v1 == 0;
 819      $q = (($u0 == $v1) ? $MAX_VAL : int(($u0*$MBASE+$u1)/$v1));
 820     --$q while ($v2*$q > ($u0*$MBASE+$u1-$q*$v1)*$MBASE+$u2);
 821     if ($q)
 822       {
 823       ($car, $bar) = (0,0);
 824       for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
 825         {
 826         $prd = $q * $y->[$yi] + $car;
 827         $prd -= ($car = int($prd / $MBASE)) * $MBASE;
 828         $x->[$xi] += $MBASE if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
 829         }
 830       if ($x->[-1] < $car + $bar)
 831         {
 832         $car = 0; --$q;
 833         for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
 834           {
 835           $x->[$xi] -= $MBASE
 836            if ($car = (($x->[$xi] += $y->[$yi] + $car) > $MBASE));
 837           }
 838         }
 839       }
 840     pop(@$x); unshift(@q, $q);
 841     }
 842   if (wantarray)
 843     {
 844     @d = ();
 845     if ($dd != 1)
 846       {
 847       $car = 0;
 848       for $xi (reverse @$x)
 849         {
 850         $prd = $car * $MBASE + $xi;
 851         $car = $prd - ($tmp = int($prd / $dd)) * $dd;
 852         unshift(@d, $tmp);
 853         }
 854       }
 855     else
 856       {
 857       @d = @$x;
 858       }
 859     @$x = @q;
 860     my $d = \@d;
 861     if ($LEN_CONVERT != 0)
 862       {
 863       $c->_to_large($x); $c->_to_large($d);
 864       }
 865     else
 866       {
 867       __strip_zeros($x);
 868       __strip_zeros($d);
 869       }
 870     return ($x,$d);
 871     }
 872   @$x = @q;
 873   if ($LEN_CONVERT != 0)
 874     {
 875     $c->_to_large($x);
 876     }
 877   else
 878     {
 879     __strip_zeros($x);
 880     }
 881   $x;
 882   }
 883
 884 ##############################################################################
 885 # testing
 886
 887 sub _acmp
 888   {
 889   # internal absolute post-normalized compare (ignore signs)
 890   # ref to array, ref to array, return <0, 0, >0
 891   # arrays must have at least one entry; this is not checked for
 892
 893   my ($c,$cx,$cy) = @_;
 894
 895   # fast comp based on array elements
 896   my $lxy = scalar @$cx - scalar @$cy;
 897   return -1 if $lxy < 0;                                # already differs, ret
 898   return 1 if $lxy > 0;                                 # ditto
 899
 900   # now calculate length based on digits, not parts
 901   $lxy = _len($c,$cx) - _len($c,$cy);                   # difference
 902   return -1 if $lxy < 0;
 903   return 1 if $lxy > 0;
 904
 905   # hm, same lengths,  but same contents?
 906   my $i = 0; my $a;
 907   # first way takes 5.49 sec instead of 4.87, but has the early out advantage
 908   # so grep is slightly faster, but more inflexible. hm. $_ instead of $k
 909   # yields 5.6 instead of 5.5 sec huh?
 910   # manual way (abort if unequal, good for early ne)
 911   my $j = scalar @$cx - 1;
 912   while ($j >= 0)
 913     {
 914     last if ($a = $cx->[$j] - $cy->[$j]); $j--;
 915     }
 916 #  my $j = scalar @$cx;
 917 #  while (--$j >= 0)
 918 #    {
 919 #    last if ($a = $cx->[$j] - $cy->[$j]);
 920 #    }
 921   return 1 if $a > 0;
 922   return -1 if $a < 0;
 923   0;                                    # equal
 924
 925   # while it early aborts, it is even slower than the manual variant
 926   #grep { return $a if ($a = $_ - $cy->[$i++]); } @$cx;
 927   # grep way, go trough all (bad for early ne)
 928   #grep { $a = $_ - $cy->[$i++]; } @$cx;
 929   #return $a;
 930   }
 931
 932 sub _len
 933   {
 934   # compute number of digits in bigint, minus the sign
 935
 936   # int() because add/sub sometimes leaves strings (like '00005') instead of
 937   # '5' in this place, thus causing length() to report wrong length
 938   my $cx = $_[1];
 939
 940   return (@$cx-1)*$BASE_LEN+length(int($cx->[-1]));
 941   }
 942
 943 sub _digit
 944   {
 945   # return the nth digit, negative values count backward
 946   # zero is rightmost, so _digit(123,0) will give 3
 947   my ($c,$x,$n) = @_;
 948
 949   my $len = _len('',$x);
 950
 951   $n = $len+$n if $n < 0;               # -1 last, -2 second-to-last
 952   $n = abs($n);                         # if negative was too big
 953   $len--; $n = $len if $n > $len;       # n to big?
 954
 955   my $elem = int($n / $BASE_LEN);       # which array element
 956   my $digit = $n % $BASE_LEN;           # which digit in this element
 957   $elem = '0000'.@$x[$elem];            # get element padded with 0's
 958   return substr($elem,-$digit-1,1);
 959   }
 960
 961 sub _zeros
 962   {
 963   # return amount of trailing zeros in decimal
 964   # check each array elem in _m for having 0 at end as long as elem == 0
 965   # Upon finding a elem != 0, stop
 966   my $x = $_[1];
 967   my $zeros = 0; my $elem;
 968   foreach my $e (@$x)
 969     {
 970     if ($e != 0)
 971       {
 972       $elem = "$e";                             # preserve x
 973       $elem =~ s/.*?(0*$)/$1/;                  # strip anything not zero
 974       $zeros *= $BASE_LEN;                      # elems * 5
 975       $zeros += length($elem);                  # count trailing zeros
 976       last;                                     # early out
 977       }
 978     $zeros ++;                                  # real else branch: 50% slower!
 979     }
 980   $zeros;
 981   }
 982
 983 ##############################################################################
 984 # _is_* routines
 985
 986 sub _is_zero
 987   {
 988   # return true if arg (BINT or num_str) is zero (array '+', '0')
 989   my $x = $_[1];
 990
 991   (((scalar @$x == 1) && ($x->[0] == 0))) <=> 0;
 992   }
 993
 994 sub _is_even
 995   {
 996   # return true if arg (BINT or num_str) is even
 997   my $x = $_[1];
 998   (!($x->[0] & 1)) <=> 0;
 999   }
1000
1001 sub _is_odd
1002   {
1003   # return true if arg (BINT or num_str) is even
1004   my $x = $_[1];
1005
1006   (($x->[0] & 1)) <=> 0;
1007   }
1008
1009 sub _is_one
1010   {
1011   # return true if arg (BINT or num_str) is one (array '+', '1')
1012   my $x = $_[1];
1013
1014   (scalar @$x == 1) && ($x->[0] == 1) <=> 0;
1015   }
1016
1017 sub __strip_zeros
1018   {
1019   # internal normalization function that strips leading zeros from the array
1020   # args: ref to array
1021   my $s = shift;
1022
1023   my $cnt = scalar @$s; # get count of parts
1024   my $i = $cnt-1;
1025   push @$s,0 if $i < 0;         # div might return empty results, so fix it
1026
1027   return $s if @$s == 1;                # early out
1028
1029   #print "strip: cnt $cnt i $i\n";
1030   # '0', '3', '4', '0', '0',
1031   #  0    1    2    3    4
1032   # cnt = 5, i = 4
1033   # i = 4
1034   # i = 3
1035   # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
1036   # >= 1: skip first part (this can be zero)
1037   while ($i > 0) { last if $s->[$i] != 0; $i--; }
1038   $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
1039   $s;
1040   }
1041
1042 ###############################################################################
1043 # check routine to test internal state of corruptions
1044
1045 sub _check
1046   {
1047   # used by the test suite
1048   my $x = $_[1];
1049
1050   return "$x is not a reference" if !ref($x);
1051
1052   # are all parts are valid?
1053   my $i = 0; my $j = scalar @$x; my ($e,$try);
1054   while ($i < $j)
1055     {
1056     $e = $x->[$i]; $e = 'undef' unless defined $e;
1057     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
1058     last if $e !~ /^[+]?[0-9]+$/;
1059     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (stringify)";
1060     last if "$e" !~ /^[+]?[0-9]+$/;
1061     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e) (cat-stringify)";
1062     last if '' . "$e" !~ /^[+]?[0-9]+$/;
1063     $try = ' < 0 || >= $BASE; '."($x, $e)";
1064     last if $e <0 || $e >= $BASE;
1065     # this test is disabled, since new/bnorm and certain ops (like early out
1066     # in add/sub) are allowed/expected to leave '00000' in some elements
1067     #$try = '=~ /^00+/; '."($x, $e)";
1068     #last if $e =~ /^00+/;
1069     $i++;
1070     }
1071   return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
1072   return 0;
1073   }
1074
1075
1076 ###############################################################################
1077 ###############################################################################
1078 # some optional routines to make BigInt faster
1079
1080 sub _mod
1081   {
1082   # if possible, use mod shortcut
1083   my ($c,$x,$yo) = @_;
1084
1085   # slow way since $y to big
1086   if (scalar @$yo > 1)
1087     {
1088     my ($xo,$rem) = _div($c,$x,$yo);
1089     return $rem;
1090     }
1091   my $y = $yo->[0];
1092   # both are single element arrays
1093   if (scalar @$x == 1)
1094     {
1095     $x->[0] %= $y;
1096     return $x;
1097     }
1098
1099   # @y is single element, but @x has more than one
1100   my $b = $BASE % $y;
1101   if ($b == 0)
1102     {
1103     # when BASE % Y == 0 then (B * BASE) % Y == 0
1104     # (B * BASE) % $y + A % Y => A % Y
1105     # so need to consider only last element: O(1)
1106     $x->[0] %= $y;
1107     }
1108   elsif ($b == 1)
1109     {
1110     # else need to go trough all elements: O(N), but loop is a bit simplified
1111     my $r = 0;
1112     foreach (@$x)
1113       {
1114       $r = ($r + $_) % $y;              # not much faster, but heh...
1115       #$r += $_ % $y; $r %= $y;
1116       }
1117     $r = 0 if $r == $y;
1118     $x->[0] = $r;
1119     }
1120   else
1121     {
1122     # else need to go trough all elements: O(N)
1123     my $r = 0; my $bm = 1;
1124     foreach (@$x)
1125       {
1126       $r = ($_ * $bm + $r) % $y;
1127       $bm = ($bm * $b) % $y;
1128
1129       #$r += ($_ % $y) * $bm;
1130       #$bm *= $b;
1131       #$bm %= $y;
1132       #$r %= $y;
1133       }
1134     $r = 0 if $r == $y;
1135     $x->[0] = $r;
1136     }
1137   splice (@$x,1);
1138   $x;
1139   }
1140
1141 ##############################################################################
1142 # shifts
1143
1144 sub _rsft
1145   {
1146   my ($c,$x,$y,$n) = @_;
1147
1148   if ($n != 10)
1149     {
1150     $n = _new($c,\$n); return _div($c,$x, _pow($c,$n,$y));
1151     }
1152
1153   # shortcut (faster) for shifting by 10)
1154   # multiples of $BASE_LEN
1155   my $dst = 0;                          # destination
1156   my $src = _num($c,$y);                # as normal int
1157   my $rem = $src % $BASE_LEN;           # remainder to shift
1158   $src = int($src / $BASE_LEN);         # source
1159   if ($rem == 0)
1160     {
1161     splice (@$x,0,$src);                # even faster, 38.4 => 39.3
1162     }
1163   else
1164     {
1165     my $len = scalar @$x - $src;        # elems to go
1166     my $vd; my $z = '0'x $BASE_LEN;
1167     $x->[scalar @$x] = 0;               # avoid || 0 test inside loop
1168     while ($dst < $len)
1169       {
1170       $vd = $z.$x->[$src];
1171       $vd = substr($vd,-$BASE_LEN,$BASE_LEN-$rem);
1172       $src++;
1173       $vd = substr($z.$x->[$src],-$rem,$rem) . $vd;
1174       $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1175       $x->[$dst] = int($vd);
1176       $dst++;
1177       }
1178     splice (@$x,$dst) if $dst > 0;              # kill left-over array elems
1179     pop @$x if $x->[-1] == 0 && @$x > 1;        # kill last element if 0
1180     } # else rem == 0
1181   $x;
1182   }
1183
1184 sub _lsft
1185   {
1186   my ($c,$x,$y,$n) = @_;
1187
1188   if ($n != 10)
1189     {
1190     $n = _new($c,\$n); return _mul($c,$x, _pow($c,$n,$y));
1191     }
1192
1193   # shortcut (faster) for shifting by 10) since we are in base 10eX
1194   # multiples of $BASE_LEN:
1195   my $src = scalar @$x;                 # source
1196   my $len = _num($c,$y);                # shift-len as normal int
1197   my $rem = $len % $BASE_LEN;           # remainder to shift
1198   my $dst = $src + int($len/$BASE_LEN); # destination
1199   my $vd;                               # further speedup
1200   $x->[$src] = 0;                       # avoid first ||0 for speed
1201   my $z = '0' x $BASE_LEN;
1202   while ($src >= 0)
1203     {
1204     $vd = $x->[$src]; $vd = $z.$vd;
1205     $vd = substr($vd,-$BASE_LEN+$rem,$BASE_LEN-$rem);
1206     $vd .= $src > 0 ? substr($z.$x->[$src-1],-$BASE_LEN,$rem) : '0' x $rem;
1207     $vd = substr($vd,-$BASE_LEN,$BASE_LEN) if length($vd) > $BASE_LEN;
1208     $x->[$dst] = int($vd);
1209     $dst--; $src--;
1210     }
1211   # set lowest parts to 0
1212   while ($dst >= 0) { $x->[$dst--] = 0; }
1213   # fix spurios last zero element
1214   splice @$x,-1 if $x->[-1] == 0;
1215   $x;
1216   }
1217
1218 sub _pow
1219   {
1220   # power of $x to $y
1221   # ref to array, ref to array, return ref to array
1222   my ($c,$cx,$cy) = @_;
1223
1224   my $pow2 = _one();
1225   my $two = _two();
1226   my $y1 = _copy($c,$cy);
1227   while (!_is_one($c,$y1))
1228     {
1229     _mul($c,$pow2,$cx) if _is_odd($c,$y1);
1230     _div($c,$y1,$two);
1231     _mul($c,$cx,$cx);
1232     }
1233   _mul($c,$cx,$pow2) unless _is_one($c,$pow2);
1234   $cx;
1235   }
1236
1237 sub _fac
1238   {
1239   # factorial of $x
1240   # ref to array, return ref to array
1241   my ($c,$cx) = @_;
1242
1243   if ((@$cx == 1) && ($cx->[0] <= 2))
1244     {
1245     $cx->[0] = 1 * ($cx->[0]||1); # 0,1 => 1, 2 => 2
1246     return $cx;
1247     }
1248
1249   # go forward until $base is exceeded
1250   # limit is either $x or $base (x == 100 means as result too high)
1251   my $steps = 100; $steps = $cx->[0] if @$cx == 1;
1252   my $r = 2; my $cf = 3; my $step = 1; my $last = $r;
1253   while ($r < $BASE && $step < $steps)
1254     {
1255     $last = $r; $r *= $cf++; $step++;
1256     }
1257   if ((@$cx == 1) && ($step == $cx->[0]))
1258     {
1259     # completely done
1260     $cx = [$last];
1261     return $cx;
1262     }
1263   my $n = _copy($c,$cx);
1264   $cx = [$last];
1265
1266   #$cx = _one();
1267   while (!(@$n == 1 && $n->[0] == $step))
1268     {
1269     _mul($c,$cx,$n); _dec($c,$n);
1270     }
1271   $cx;
1272   }
1273
1274 use constant DEBUG => 0;
1275
1276 my $steps = 0;
1277
1278 sub steps { $steps };
1279
1280 sub _sqrt
1281   {
1282   # square-root of $x
1283   # ref to array, return ref to array
1284   my ($c,$x) = @_;
1285
1286   if (scalar @$x == 1)
1287     {
1288     # fit's into one Perl scalar
1289     $x->[0] = int(sqrt($x->[0]));
1290     return $x;
1291     }
1292   my $y = _copy($c,$x);
1293   # hopefully _len/2 is < $BASE, the -1 is to always undershot the guess
1294   # since our guess will "grow"
1295   my $l = int((_len($c,$x)-1) / 2);
1296
1297   my $lastelem = $x->[-1];      # for guess
1298   my $elems = scalar @$x - 1;
1299   # not enough digits, but could have more?
1300   if ((length($lastelem) <= 3) && ($elems > 1))
1301     {
1302     # right-align with zero pad
1303     my $len = length($lastelem) & 1;
1304     print "$lastelem => " if DEBUG;
1305     $lastelem .= substr($x->[-2] . '0' x $BASE_LEN,0,$BASE_LEN);
1306     # former odd => make odd again, or former even to even again
1307     $lastelem = $lastelem / 10 if (length($lastelem) & 1) != $len;
1308     print "$lastelem\n" if DEBUG;
1309     }
1310
1311   # construct $x (instead of _lsft($c,$x,$l,10)
1312   my $r = $l % $BASE_LEN;       # 10000 00000 00000 00000 ($BASE_LEN=5)
1313   $l = int($l / $BASE_LEN);
1314   print "l =  $l " if DEBUG;
1315
1316   splice @$x,$l;                # keep ref($x), but modify it
1317
1318   # we make the first part of the guess not '1000...0' but int(sqrt($lastelem))
1319   # that gives us:
1320   # 14400 00000 => sqrt(14400) => 120
1321   # 144000 000000 => sqrt(144000) => 379
1322
1323   # $x->[$l--] = int('1' . '0' x $r);                   # old way of guessing
1324   print "$lastelem (elems $elems) => " if DEBUG;
1325   $lastelem = $lastelem / 10 if ($elems & 1 == 1);              # odd or even?
1326   my $g = sqrt($lastelem); $g =~ s/\.//;                        # 2.345 => 2345
1327   $r -= 1 if $elems & 1 == 0;                                   # 70 => 7
1328
1329   # padd with zeros if result is too short
1330   $x->[$l--] = int(substr($g . '0' x $r,0,$r+1));
1331   print "now ",$x->[-1] if DEBUG;
1332   print " would have been ", int('1' . '0' x $r),"\n" if DEBUG;
1333
1334   # If @$x > 1, we could compute the second elem of the guess, too, to create
1335   # an even better guess. Not implemented yet.
1336   $x->[$l--] = 0 while ($l >= 0);       # all other digits of guess are zero
1337
1338   print "start x= ",${_str($c,$x)},"\n" if DEBUG;
1339   my $two = _two();
1340   my $last = _zero();
1341   my $lastlast = _zero();
1342   $steps = 0 if DEBUG;
1343   while (_acmp($c,$last,$x) != 0 && _acmp($c,$lastlast,$x) != 0)
1344     {
1345     $steps++ if DEBUG;
1346     $lastlast = _copy($c,$last);
1347     $last = _copy($c,$x);
1348     _add($c,$x, _div($c,_copy($c,$y),$x));
1349     _div($c,$x, $two );
1350     print "      x= ",${_str($c,$x)},"\n" if DEBUG;
1351     }
1352   print "\nsteps in sqrt: $steps, " if DEBUG;
1353   _dec($c,$x) if _acmp($c,$y,_mul($c,_copy($c,$x),$x)) < 0;     # overshot?
1354   print " final ",$x->[-1],"\n" if DEBUG;
1355   $x;
1356   }
1357
1358 ##############################################################################
1359 # binary stuff
1360
1361 sub _and
1362   {
1363   my ($c,$x,$y) = @_;
1364
1365   # the shortcut makes equal, large numbers _really_ fast, and makes only a
1366   # very small performance drop for small numbers (e.g. something with less
1367   # than 32 bit) Since we optimize for large numbers, this is enabled.
1368   return $x if _acmp($c,$x,$y) == 0;            # shortcut
1369
1370   my $m = _one(); my ($xr,$yr);
1371   my $mask = $AND_MASK;
1372
1373   my $x1 = $x;
1374   my $y1 = _copy($c,$y);                        # make copy
1375   $x = _zero();
1376   my ($b,$xrr,$yrr);
1377   use integer;
1378   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1379     {
1380     ($x1, $xr) = _div($c,$x1,$mask);
1381     ($y1, $yr) = _div($c,$y1,$mask);
1382
1383     # make ints() from $xr, $yr
1384     # this is when the AND_BITS are greater tahn $BASE and is slower for
1385     # small (<256 bits) numbers, but faster for large numbers. Disabled
1386     # due to KISS principle
1387
1388 #    $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1389 #    $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1390 #    _add($c,$x, _mul($c, _new( $c, \($xrr & $yrr) ), $m) );
1391
1392     # 0+ due to '&' doesn't work in strings
1393     _add($c,$x, _mul($c, [ 0+$xr->[0] & 0+$yr->[0] ], $m) );
1394     _mul($c,$m,$mask);
1395     }
1396   $x;
1397   }
1398
1399 sub _xor
1400   {
1401   my ($c,$x,$y) = @_;
1402
1403   return _zero() if _acmp($c,$x,$y) == 0;       # shortcut (see -and)
1404
1405   my $m = _one(); my ($xr,$yr);
1406   my $mask = $XOR_MASK;
1407
1408   my $x1 = $x;
1409   my $y1 = _copy($c,$y);                        # make copy
1410   $x = _zero();
1411   my ($b,$xrr,$yrr);
1412   use integer;
1413   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1414     {
1415     ($x1, $xr) = _div($c,$x1,$mask);
1416     ($y1, $yr) = _div($c,$y1,$mask);
1417     # make ints() from $xr, $yr (see _and())
1418     #$b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1419     #$b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1420     #_add($c,$x, _mul($c, _new( $c, \($xrr ^ $yrr) ), $m) );
1421
1422     # 0+ due to '^' doesn't work in strings
1423     _add($c,$x, _mul($c, [ 0+$xr->[0] ^ 0+$yr->[0] ], $m) );
1424     _mul($c,$m,$mask);
1425     }
1426   # the loop stops when the shorter of the two numbers is exhausted
1427   # the remainder of the longer one will survive bit-by-bit, so we simple
1428   # multiply-add it in
1429   _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1430   _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1431
1432   $x;
1433   }
1434
1435 sub _or
1436   {
1437   my ($c,$x,$y) = @_;
1438
1439   return $x if _acmp($c,$x,$y) == 0;            # shortcut (see _and)
1440
1441   my $m = _one(); my ($xr,$yr);
1442   my $mask = $OR_MASK;
1443
1444   my $x1 = $x;
1445   my $y1 = _copy($c,$y);                        # make copy
1446   $x = _zero();
1447   my ($b,$xrr,$yrr);
1448   use integer;
1449   while (!_is_zero($c,$x1) && !_is_zero($c,$y1))
1450     {
1451     ($x1, $xr) = _div($c,$x1,$mask);
1452     ($y1, $yr) = _div($c,$y1,$mask);
1453     # make ints() from $xr, $yr (see _and())
1454 #    $b = 1; $xrr = 0; foreach (@$xr) { $xrr += $_ * $b; $b *= $BASE; }
1455 #    $b = 1; $yrr = 0; foreach (@$yr) { $yrr += $_ * $b; $b *= $BASE; }
1456 #    _add($c,$x, _mul($c, _new( $c, \($xrr | $yrr) ), $m) );
1457
1458     # 0+ due to '|' doesn't work in strings
1459     _add($c,$x, _mul($c, [ 0+$xr->[0] | 0+$yr->[0] ], $m) );
1460     _mul($c,$m,$mask);
1461     }
1462   # the loop stops when the shorter of the two numbers is exhausted
1463   # the remainder of the longer one will survive bit-by-bit, so we simple
1464   # multiply-add it in
1465   _add($c,$x, _mul($c, $x1, $m) ) if !_is_zero($c,$x1);
1466   _add($c,$x, _mul($c, $y1, $m) ) if !_is_zero($c,$y1);
1467
1468   $x;
1469   }
1470
1471 sub _as_hex
1472   {
1473   # convert a decimal number to hex (ref to array, return ref to string)
1474   my ($c,$x) = @_;
1475
1476   my $x1 = _copy($c,$x);
1477
1478   my $es = '';
1479   my $xr;
1480   my $x10000 = [ 0x10000 ];
1481   while (! _is_zero($c,$x1))
1482     {
1483     ($x1, $xr) = _div($c,$x1,$x10000);
1484     $es .= unpack('h4',pack('v',$xr->[0]));
1485     }
1486   $es = reverse $es;
1487   $es =~ s/^[0]+//;   # strip leading zeros
1488   $es = '0x' . $es;
1489   \$es;
1490   }
1491
1492 sub _as_bin
1493   {
1494   # convert a decimal number to bin (ref to array, return ref to string)
1495   my ($c,$x) = @_;
1496
1497   my $x1 = _copy($c,$x);
1498
1499   my $es = '';
1500   my $xr;
1501   my $x10000 = [ 0x10000 ];
1502   while (! _is_zero($c,$x1))
1503     {
1504     ($x1, $xr) = _div($c,$x1,$x10000);
1505     $es .= unpack('b16',pack('v',$xr->[0]));
1506     }
1507   $es = reverse $es;
1508   $es =~ s/^[0]+//;   # strip leading zeros
1509   $es = '0b' . $es;
1510   \$es;
1511   }
1512
1513 sub _from_hex
1514   {
1515   # convert a hex number to decimal (ref to string, return ref to array)
1516   my ($c,$hs) = @_;
1517
1518   my $mul = _one();
1519   my $m = [ 0x10000 ];                          # 16 bit at a time
1520   my $x = _zero();
1521
1522   my $len = length($$hs)-2;
1523   $len = int($len/4);                           # 4-digit parts, w/o '0x'
1524   my $val; my $i = -4;
1525   while ($len >= 0)
1526     {
1527     $val = substr($$hs,$i,4);
1528     $val =~ s/^[+-]?0x// if $len == 0;          # for last part only because
1529     $val = hex($val);                           # hex does not like wrong chars
1530     $i -= 4; $len --;
1531     _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1532     _mul ($c, $mul, $m ) if $len >= 0;          # skip last mul
1533     }
1534   $x;
1535   }
1536
1537 sub _from_bin
1538   {
1539   # convert a hex number to decimal (ref to string, return ref to array)
1540   my ($c,$bs) = @_;
1541
1542   # instead of converting 8 bit at a time, it is faster to convert the
1543   # number to hex, and then call _from_hex.
1544
1545   my $hs = $$bs;
1546   $hs =~ s/^[+-]?0b//;                                  # remove sign and 0b
1547   my $l = length($hs);                                  # bits
1548   $hs = '0' x (8-($l % 8)) . $hs if ($l % 8) != 0;      # padd left side w/ 0
1549   my $h = unpack('H*', pack ('B*', $hs));               # repack as hex
1550   return $c->_from_hex(\('0x'.$h));
1551
1552   my $mul = _one();
1553   my $m = [ 0x100 ];                            # 8 bit at a time
1554   my $x = _zero();
1555
1556   my $len = length($$bs)-2;
1557   $len = int($len/8);                           # 4-digit parts, w/o '0x'
1558   my $val; my $i = -8;
1559   while ($len >= 0)
1560     {
1561     $val = substr($$bs,$i,8);
1562     $val =~ s/^[+-]?0b// if $len == 0;          # for last part only
1563
1564     $val = ord(pack('B8',substr('00000000'.$val,-8,8)));
1565
1566     $i -= 8; $len --;
1567     _add ($c, $x, _mul ($c, [ $val ], $mul ) ) if $val != 0;
1568     _mul ($c, $mul, $m ) if $len >= 0;          # skip last mul
1569     }
1570   $x;
1571   }
1572
1573 ##############################################################################
1574 ##############################################################################
1575
1576 1;
1577 __END__
1578
1579 =head1 NAME
1580
1581 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
1582
1583 =head1 SYNOPSIS
1584
1585 Provides support for big integer calculations. Not intended to be used by other
1586 modules (except Math::BigInt::Cached). Other modules which sport the same
1587 functions can also be used to support Math::Bigint, like Math::BigInt::Pari.
1588
1589 =head1 DESCRIPTION
1590
1591 In order to allow for multiple big integer libraries, Math::BigInt was
1592 rewritten to use library modules for core math routines. Any module which
1593 follows the same API as this can be used instead by using the following:
1594
1595         use Math::BigInt lib => 'libname';
1596
1597 'libname' is either the long name ('Math::BigInt::Pari'), or only the short
1598 version like 'Pari'.
1599
1600 =head1 EXPORT
1601
1602 The following functions MUST be defined in order to support the use by
1603 Math::BigInt:
1604
1605         _new(string)    return ref to new object from ref to decimal string
1606         _zero()         return a new object with value 0
1607         _one()          return a new object with value 1
1608
1609         _str(obj)       return ref to a string representing the object
1610         _num(obj)       returns a Perl integer/floating point number
1611                         NOTE: because of Perl numeric notation defaults,
1612                         the _num'ified obj may lose accuracy due to
1613                         machine-dependend floating point size limitations
1614
1615         _add(obj,obj)   Simple addition of two objects
1616         _mul(obj,obj)   Multiplication of two objects
1617         _div(obj,obj)   Division of the 1st object by the 2nd
1618                         In list context, returns (result,remainder).
1619                         NOTE: this is integer math, so no
1620                         fractional part will be returned.
1621         _sub(obj,obj)   Simple subtraction of 1 object from another
1622                         a third, optional parameter indicates that the params
1623                         are swapped. In this case, the first param needs to
1624                         be preserved, while you can destroy the second.
1625                         sub (x,y,1) => return x - y and keep x intact!
1626         _dec(obj)       decrement object by one (input is garant. to be > 0)
1627         _inc(obj)       increment object by one
1628
1629
1630         _acmp(obj,obj)  <=> operator for objects (return -1, 0 or 1)
1631
1632         _len(obj)       returns count of the decimal digits of the object
1633         _digit(obj,n)   returns the n'th decimal digit of object
1634
1635         _is_one(obj)    return true if argument is +1
1636         _is_zero(obj)   return true if argument is 0
1637         _is_even(obj)   return true if argument is even (0,2,4,6..)
1638         _is_odd(obj)    return true if argument is odd (1,3,5,7..)
1639
1640         _copy           return a ref to a true copy of the object
1641
1642         _check(obj)     check whether internal representation is still intact
1643                         return 0 for ok, otherwise error message as string
1644
1645 The following functions are optional, and can be defined if the underlying lib
1646 has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
1647 slow) fallback routines to emulate these:
1648
1649         _from_hex(str)  return ref to new object from ref to hexadecimal string
1650         _from_bin(str)  return ref to new object from ref to binary string
1651
1652         _as_hex(str)    return ref to scalar string containing the value as
1653                         unsigned hex string, with the '0x' prepended.
1654                         Leading zeros must be stripped.
1655         _as_bin(str)    Like as_hex, only as binary string containing only
1656                         zeros and ones. Leading zeros must be stripped and a
1657                         '0b' must be prepended.
1658
1659         _rsft(obj,N,B)  shift object in base B by N 'digits' right
1660                         For unsupported bases B, return undef to signal failure
1661         _lsft(obj,N,B)  shift object in base B by N 'digits' left
1662                         For unsupported bases B, return undef to signal failure
1663
1664         _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
1665                         Note: XOR, AND and OR pad with zeros if size mismatches
1666         _and(obj1,obj2) AND (bit-wise) object 1 with object 2
1667         _or(obj1,obj2)  OR (bit-wise) object 1 with object 2
1668
1669         _mod(obj,obj)   Return remainder of div of the 1st by the 2nd object
1670         _sqrt(obj)      return the square root of object (truncate to int)
1671         _fac(obj)       return factorial of object 1 (1*2*3*4..)
1672         _pow(obj,obj)   return object 1 to the power of object 2
1673         _gcd(obj,obj)   return Greatest Common Divisor of two objects
1674
1675         _zeros(obj)     return number of trailing decimal zeros
1676
1677 Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
1678 or '0b1101').
1679
1680 Testing of input parameter validity is done by the caller, so you need not
1681 worry about underflow (f.i. in C<_sub()>, C<_dec()>) nor about division by
1682 zero or similar cases.
1683
1684 The first parameter can be modified, that includes the possibility that you
1685 return a reference to a completely different object instead. Although keeping
1686 the reference and just changing it's contents is prefered over creating and
1687 returning a different reference.
1688
1689 Return values are always references to objects or strings. Exceptions are
1690 C<_lsft()> and C<_rsft()>, which return undef if they can not shift the
1691 argument. This is used to delegate shifting of bases different than the one
1692 you can support back to Math::BigInt, which will use some generic code to
1693 calculate the result.
1694
1695 =head1 WRAP YOUR OWN
1696
1697 If you want to port your own favourite c-lib for big numbers to the
1698 Math::BigInt interface, you can take any of the already existing modules as
1699 a rough guideline. You should really wrap up the latest BigInt and BigFloat
1700 testsuites with your module, and replace in them any of the following:
1701
1702         use Math::BigInt;
1703
1704 by this:
1705
1706         use Math::BigInt lib => 'yourlib';
1707
1708 This way you ensure that your library really works 100% within Math::BigInt.
1709
1710 =head1 LICENSE
1711
1712 This program is free software; you may redistribute it and/or modify it under
1713 the same terms as Perl itself.
1714
1715 =head1 AUTHORS
1716
1717 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
1718 in late 2000, 2001.
1719 Seperated from BigInt and shaped API with the help of John Peacock.
1720
1721 =head1 SEE ALSO
1722
1723 L<Math::BigInt>, L<Math::BigFloat>, L<Math::BigInt::BitVect>,
1724 L<Math::BigInt::GMP>, L<Math::BigInt::Cached> and L<Math::BigInt::Pari>.
1725
1726 =cut