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