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