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