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