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