lib/Math/BigInt/Calc.pm

   1 package Math::BigInt::Calc;
   2
   3 use 5.005;
   4 use strict;
   5 use warnings;
   6
   7 require Exporter;
   8
   9 use vars qw/ @ISA @EXPORT $VERSION/;
  10 @ISA = qw(Exporter);
  11
  12 @EXPORT = qw(
  13         _add _mul _div _mod _sub
  14         _new
  15         _str _num _acmp _len
  16         _digit
  17         _is_zero _is_one
  18         _is_even _is_odd
  19         _check _zero _one _copy _zeros
  20 );
  21 $VERSION = '0.06';
  22
  23 # Package to store unsigned big integers in decimal and do math with them
  24
  25 # Internally the numbers are stored in an array with at least 1 element, no
  26 # leading zero parts (except the first) and in base 100000
  27
  28 # todo:
  29 # - fully remove funky $# stuff (maybe)
  30 # - use integer; vs 1e7 as base
  31
  32 # USE_MUL: due to problems on certain os (os390, posix-bc) "* 1e-5" is used
  33 # instead of "/ 1e5" at some places, (marked with USE_MUL). But instead of
  34 # using the reverse only on problematic machines, I used it everytime to avoid
  35 # the costly comparisons. This _should_ work everywhere. Thanx Peter Prymmer
  36
  37 ##############################################################################
  38 # global constants, flags and accessory
  39
  40 # constants for easier life
  41 my $nan = 'NaN';
  42 my $BASE_LEN = 5;
  43 my $BASE = int("1e".$BASE_LEN);         # var for trying to change it to 1e7
  44 my $RBASE = 1e-5;                       # see USE_MUL
  45 my $class = 'Math::BigInt::Calc';
  46
  47 ##############################################################################
  48 # create objects from various representations
  49
  50 sub _new
  51   {
  52   # (string) return ref to num_array
  53   # Convert a number from string format to internal base 100000 format.
  54   # Assumes normalized value as input.
  55   shift @_ if $_[0] eq $class;
  56   my $d = shift;
  57   # print "_new $d $$d\n";
  58   my $il = CORE::length($$d)-1;
  59   # these leaves '00000' instead of int 0 and will be corrected after any op
  60   return [ reverse(unpack("a" . ($il%5+1) . ("a5" x ($il/5)), $$d)) ];
  61   }
  62
  63 sub _zero
  64   {
  65   # create a zero
  66   return [ 0 ];
  67   }
  68
  69 sub _one
  70   {
  71   # create a one
  72   return [ 1 ];
  73   }
  74
  75 sub _copy
  76   {
  77   shift @_ if $_[0] eq $class;
  78   my $x = shift;
  79   return [ @$x ];
  80   }
  81
  82 ##############################################################################
  83 # convert back to string and number
  84
  85 sub _str
  86   {
  87   # (ref to BINT) return num_str
  88   # Convert number from internal base 100000 format to string format.
  89   # internal format is always normalized (no leading zeros, "-0" => "+0")
  90   shift @_ if $_[0] eq $class;
  91   my $ar = shift;
  92   my $ret = "";
  93   my $l = scalar @$ar;         # number of parts
  94   return $nan if $l < 1;       # should not happen
  95   # handle first one different to strip leading zeros from it (there are no
  96   # leading zero parts in internal representation)
  97   $l --; $ret .= $ar->[$l]; $l--;
  98   # Interestingly, the pre-padd method uses more time
  99   # the old grep variant takes longer (14 to 10 sec)
 100   while ($l >= 0)
 101     {
 102     $ret .= substr('0000'.$ar->[$l],-5);   # fastest way I could think of
 103     $l--;
 104     }
 105   return \$ret;
 106   }
 107
 108 sub _num
 109   {
 110   # Make a number (scalar int/float) from a BigInt object
 111   shift @_ if $_[0] eq $class;
 112   my $x = shift;
 113   return $x->[0] if scalar @$x == 1;  # below $BASE
 114   my $fac = 1;
 115   my $num = 0;
 116   foreach (@$x)
 117     {
 118     $num += $fac*$_; $fac *= $BASE;
 119     }
 120   return $num;
 121   }
 122
 123 ##############################################################################
 124 # actual math code
 125
 126 sub _add
 127   {
 128   # (ref to int_num_array, ref to int_num_array)
 129   # routine to add two base 1e5 numbers
 130   # stolen from Knuth Vol 2 Algorithm A pg 231
 131   # there are separate routines to add and sub as per Knuth pg 233
 132   # This routine clobbers up array x, but not y.
 133
 134   shift @_ if $_[0] eq $class;
 135   my ($x,$y) = @_;
 136
 137   # for each in Y, add Y to X and carry. If after that, something is left in
 138   # X, foreach in X add carry to X and then return X, carry
 139   # Trades one "$j++" for having to shift arrays, $j could be made integer
 140   # but this would impose a limit to number-length of 2**32.
 141   my $i; my $car = 0; my $j = 0;
 142   for $i (@$y)
 143     {
 144     $x->[$j] -= $BASE
 145       if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
 146     $j++;
 147     }
 148   while ($car != 0)
 149     {
 150     $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
 151     }
 152     return $x;
 153   }
 154
 155 sub _sub
 156   {
 157   # (ref to int_num_array, ref to int_num_array)
 158   # subtract base 1e5 numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
 159   # subtract Y from X (X is always greater/equal!) by modifying x in place
 160   shift @_ if $_[0] eq $class;
 161   my ($sx,$sy,$s) = @_;
 162
 163   my $car = 0; my $i; my $j = 0;
 164   if (!$s)
 165     {
 166     #print "case 2\n";
 167     for $i (@$sx)
 168       {
 169       last unless defined $sy->[$j] || $car;
 170       #print "x: $i y: $sy->[$j] c: $car\n";
 171       $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
 172       #print "x: $i y: $sy->[$j-1] c: $car\n";
 173       }
 174     # might leave leading zeros, so fix that
 175     __strip_zeros($sx);
 176     return $sx;
 177     }
 178   else
 179     {
 180     #print "case 1 (swap)\n";
 181     for $i (@$sx)
 182       {
 183       last unless defined $sy->[$j] || $car;
 184       #print "$sy->[$j] $i $car => $sx->[$j]\n";
 185       $sy->[$j] += $BASE
 186        if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
 187       #print "$sy->[$j] $i $car => $sy->[$j]\n";
 188       $j++;
 189       }
 190     # might leave leading zeros, so fix that
 191     __strip_zeros($sy);
 192     return $sy;
 193     }
 194   }
 195
 196 sub _mul
 197   {
 198   # (BINT, BINT) return nothing
 199   # multiply two numbers in internal representation
 200   # modifies first arg, second need not be different from first
 201   shift @_ if $_[0] eq $class;
 202   my ($xv,$yv) = @_;
 203
 204   my @prod = (); my ($prod,$car,$cty,$xi,$yi);
 205   # since multiplying $x with $x fails, make copy in this case
 206   $yv = [@$xv] if "$xv" eq "$yv";
 207     # looping through @$y if $xi == 0 is silly! optimize it!
 208   for $xi (@$xv)
 209     {
 210     $car = 0; $cty = 0;
 211     for $yi (@$yv)
 212       {
 213       $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
 214       $prod[$cty++] =
 215        $prod - ($car = int($prod * 1e-5)) * $BASE;  # see USE_MUL
 216       }
 217     $prod[$cty] += $car if $car; # need really to check for 0?
 218     $xi = shift @prod;
 219     }
 220 #  for $xi (@$xv)
 221 #    {
 222 #    $car = 0; $cty = 0;
 223 #    # looping through this if $xi == 0 is silly! optimize it!
 224 #    if (($xi||0) != 0)
 225 #      {
 226 #      for $yi (@$yv)
 227 #        {
 228 #        $prod = $prod[$cty]; $prod += ($car + $xi * $yi);      # no ||0 here
 229 #        $prod[$cty++] =
 230 #         $prod - ($car = int($prod * 1e-5)) * $BASE;  # see USE_MUL
 231 #        }
 232 #      }
 233 #    $prod[$cty] += $car if $car; # need really to check for 0?
 234 #    $xi = shift @prod;
 235 #    }
 236   push @$xv, @prod;
 237   __strip_zeros($xv);
 238   # normalize (handled last to save check for $y->is_zero()
 239   return $xv;
 240   }
 241
 242 sub _div
 243   {
 244   # ref to array, ref to array, modify first array and return remainder if
 245   # in list context
 246   # no longer handles sign
 247   shift @_ if $_[0] eq $class;
 248   my ($x,$yorg) = @_;
 249   my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1);
 250
 251   my (@d,$tmp,$q,$u2,$u1,$u0);
 252
 253   $car = $bar = $prd = 0;
 254
 255   my $y = [ @$yorg ];
 256   if (($dd = int($BASE/($y->[-1]+1))) != 1)
 257     {
 258     for $xi (@$x)
 259       {
 260       $xi = $xi * $dd + $car;
 261       $xi -= ($car = int($xi * $RBASE)) * $BASE;        # see USE_MUL
 262       }
 263     push(@$x, $car); $car = 0;
 264     for $yi (@$y)
 265       {
 266       $yi = $yi * $dd + $car;
 267       $yi -= ($car = int($yi * $RBASE)) * $BASE;        # see USE_MUL
 268       }
 269     }
 270   else
 271     {
 272     push(@$x, 0);
 273     }
 274   @q = (); ($v2,$v1) = @$y[-2,-1];
 275   $v2 = 0 unless $v2;
 276   while ($#$x > $#$y)
 277     {
 278     ($u2,$u1,$u0) = @$x[-3..-1];
 279     $u2 = 0 unless $u2;
 280     #warn "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
 281     # if $v1 == 0;
 282     $q = (($u0 == $v1) ? 99999 : int(($u0*$BASE+$u1)/$v1));
 283     --$q while ($v2*$q > ($u0*1e5+$u1-$q*$v1)*$BASE+$u2);
 284     if ($q)
 285       {
 286       ($car, $bar) = (0,0);
 287       for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
 288         {
 289         $prd = $q * $y->[$yi] + $car;
 290         $prd -= ($car = int($prd * $RBASE)) * $BASE;    # see USE_MUL
 291         $x->[$xi] += 1e5 if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
 292         }
 293       if ($x->[-1] < $car + $bar)
 294         {
 295         $car = 0; --$q;
 296         for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
 297           {
 298           $x->[$xi] -= 1e5
 299            if ($car = (($x->[$xi] += $y->[$yi] + $car) > $BASE));
 300           }
 301         }
 302       }
 303       pop(@$x); unshift(@q, $q);
 304     }
 305   if (wantarray)
 306     {
 307     @d = ();
 308     if ($dd != 1)
 309       {
 310       $car = 0;
 311       for $xi (reverse @$x)
 312         {
 313         $prd = $car * $BASE + $xi;
 314         $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
 315         unshift(@d, $tmp);
 316         }
 317       }
 318     else
 319       {
 320       @d = @$x;
 321       }
 322     @$x = @q;
 323     __strip_zeros($x);
 324     __strip_zeros(\@d);
 325     return ($x,\@d);
 326     }
 327   @$x = @q;
 328   __strip_zeros($x);
 329   return $x;
 330   }
 331
 332 ##############################################################################
 333 # testing
 334
 335 sub _acmp
 336   {
 337   # internal absolute post-normalized compare (ignore signs)
 338   # ref to array, ref to array, return <0, 0, >0
 339   # arrays must have at least one entry; this is not checked for
 340
 341   shift @_ if $_[0] eq $class;
 342   my ($cx, $cy) = @_;
 343
 344   #print "$cx $cy\n";
 345   my ($i,$a,$x,$y,$k);
 346   # calculate length based on digits, not parts
 347   $x = _len($cx); $y = _len($cy);
 348   # print "length: ",($x-$y),"\n";
 349   return $x-$y if ($x - $y);              # if different in length
 350   #print "full compare\n";
 351   $i = 0; $a = 0;
 352   # first way takes 5.49 sec instead of 4.87, but has the early out advantage
 353   # so grep is slightly faster, but more inflexible. hm. $_ instead of $k
 354   # yields 5.6 instead of 5.5 sec huh?
 355   # manual way (abort if unequal, good for early ne)
 356   my $j = scalar @$cx - 1;
 357   while ($j >= 0)
 358    {
 359    # print "$cx->[$j] $cy->[$j] $a",$cx->[$j]-$cy->[$j],"\n";
 360    last if ($a = $cx->[$j] - $cy->[$j]); $j--;
 361    }
 362   return $a;
 363   # while it early aborts, it is even slower than the manual variant
 364   #grep { return $a if ($a = $_ - $cy->[$i++]); } @$cx;
 365   # grep way, go trough all (bad for early ne)
 366   #grep { $a = $_ - $cy->[$i++]; } @$cx;
 367   #return $a;
 368   }
 369
 370 sub _len
 371   {
 372   # computer number of digits in bigint, minus the sign
 373   # int() because add/sub sometimes leaves strings (like '00005') instead of
 374   # int ('5') in this place, causing length to fail
 375   shift @_ if $_[0] eq $class;
 376   my $cx = shift;
 377
 378   return (@$cx-1)*5+length(int($cx->[-1]));
 379   }
 380
 381 sub _digit
 382   {
 383   # return the nth digit, negative values count backward
 384   # zero is rightmost, so _digit(123,0) will give 3
 385   shift @_ if $_[0] eq $class;
 386   my $x = shift;
 387   my $n = shift || 0;
 388
 389   my $len = _len($x);
 390
 391   $n = $len+$n if $n < 0;               # -1 last, -2 second-to-last
 392   $n = abs($n);                         # if negative was too big
 393   $len--; $n = $len if $n > $len;       # n to big?
 394
 395   my $elem = int($n / 5);               # which array element
 396   my $digit = $n % 5;                   # which digit in this element
 397   $elem = '0000'.@$x[$elem];            # get element padded with 0's
 398   return substr($elem,-$digit-1,1);
 399   }
 400
 401 sub _zeros
 402   {
 403   # return amount of trailing zeros in decimal
 404   # check each array elem in _m for having 0 at end as long as elem == 0
 405   # Upon finding a elem != 0, stop
 406   shift @_ if $_[0] eq $class;
 407   my $x = shift;
 408   my $zeros = 0; my $elem;
 409   foreach my $e (@$x)
 410     {
 411     if ($e != 0)
 412       {
 413       $elem = "$e";                            # preserve x
 414       $elem =~ s/.*?(0*$)/$1/;                 # strip anything not zero
 415       $zeros *= 5;                             # elems * 5
 416       $zeros += CORE::length($elem);           # count trailing zeros
 417       last;                                    # early out
 418       }
 419     $zeros ++;                                 # real else branch: 50% slower!
 420     }
 421   return $zeros;
 422   }
 423
 424 ##############################################################################
 425 # _is_* routines
 426
 427 sub _is_zero
 428   {
 429   # return true if arg (BINT or num_str) is zero (array '+', '0')
 430   shift @_ if $_[0] eq $class;
 431   my ($x) = shift;
 432   return (((scalar @$x == 1) && ($x->[0] == 0))) <=> 0;
 433   }
 434
 435 sub _is_even
 436   {
 437   # return true if arg (BINT or num_str) is even
 438   shift @_ if $_[0] eq $class;
 439   my ($x) = shift;
 440   return (!($x->[0] & 1)) <=> 0;
 441   }
 442
 443 sub _is_odd
 444   {
 445   # return true if arg (BINT or num_str) is even
 446   shift @_ if $_[0] eq $class;
 447   my ($x) = shift;
 448   return (($x->[0] & 1)) <=> 0;
 449   }
 450
 451 sub _is_one
 452   {
 453   # return true if arg (BINT or num_str) is one (array '+', '1')
 454   shift @_ if $_[0] eq $class;
 455   my ($x) = shift;
 456   return (scalar @$x == 1) && ($x->[0] == 1) <=> 0;
 457   }
 458
 459 sub __strip_zeros
 460   {
 461   # internal normalization function that strips leading zeros from the array
 462   # args: ref to array
 463   #trace(@_);
 464   shift @_ if $_[0] eq $class;
 465   my $s = shift;
 466
 467   my $cnt = scalar @$s; # get count of parts
 468   my $i = $cnt-1;
 469   #print "strip: cnt $cnt i $i\n";
 470   # '0', '3', '4', '0', '0',
 471   #  0    1    2    3    4
 472   # cnt = 5, i = 4
 473   # i = 4
 474   # i = 3
 475   # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
 476   # >= 1: skip first part (this can be zero)
 477   while ($i > 0) { last if $s->[$i] != 0; $i--; }
 478   $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
 479   return $s;
 480   }
 481
 482 ###############################################################################
 483 # check routine to test internal state of corruptions
 484
 485 sub _check
 486   {
 487   # no checks yet, pull it out from the test suite
 488   shift @_ if $_[0] eq $class;
 489
 490   my ($x) = shift;
 491   return "$x is not a reference" if !ref($x);
 492
 493   # are all parts are valid?
 494   my $i = 0; my $j = scalar @$x; my ($e,$try);
 495   while ($i < $j)
 496     {
 497     $e = $x->[$i]; $e = 'undef' unless defined $e;
 498     $try = '=~ /^[\+]?[0-9]+\$/; '."($x, $e)";
 499     last if $e !~ /^[+]?[0-9]+$/;
 500     $try = ' < 0 || >= $BASE; '."($x, $e)";
 501     last if $e <0 || $e >= $BASE;
 502     # this test is disabled, since new/bnorm and certain ops (like early out
 503     # in add/sub) are allowed/expected to leave '00000' in some elements
 504     #$try = '=~ /^00+/; '."($x, $e)";
 505     #last if $e =~ /^00+/;
 506     $i++;
 507     }
 508   return "Illegal part '$e' at pos $i (tested: $try)" if $i < $j;
 509   return 0;
 510   }
 511
 512 1;
 513 __END__
 514
 515 =head1 NAME
 516
 517 Math::BigInt::Calc - Pure Perl module to support Math::BigInt
 518
 519 =head1 SYNOPSIS
 520
 521 Provides support for big integer calculations. Not intended
 522 to be used by other modules. Other modules which export the
 523 same functions can also be used to support Math::Bigint
 524
 525 =head1 DESCRIPTION
 526
 527 In order to allow for multiple big integer libraries, Math::BigInt
 528 was rewritten to use library modules for core math routines. Any
 529 module which follows the same API as this can be used instead by
 530 using the following call:
 531
 532         use Math::BigInt Calc => BigNum;
 533
 534 =head1 EXPORT
 535
 536 The following functions MUST be exported in order to support
 537 the use by Math::BigInt:
 538
 539         _new(string)    return ref to new object from ref to decimal string
 540         _zero()         return a new object with value 0
 541         _one()          return a new object with value 1
 542
 543         _str(obj)       return ref to a string representing the object
 544         _num(obj)       returns a Perl integer/floating point number
 545                         NOTE: because of Perl numeric notation defaults,
 546                         the _num'ified obj may lose accuracy due to
 547                         machine-dependend floating point size limitations
 548
 549         _add(obj,obj)   Simple addition of two objects
 550         _mul(obj,obj)   Multiplication of two objects
 551         _div(obj,obj)   Division of the 1st object by the 2nd
 552                         In list context, returns (result,remainder).
 553                         NOTE: this is integer math, so no
 554                         fractional part will be returned.
 555         _sub(obj,obj)   Simple subtraction of 1 object from another
 556                         a third, optional parameter indicates that the params
 557                         are swapped. In this case, the first param needs to
 558                         be preserved, while you can destroy the second.
 559                         sub (x,y,1) => return x - y and keep x intact!
 560
 561         _acmp(obj,obj)  <=> operator for objects (return -1, 0 or 1)
 562
 563         _len(obj)       returns count of the decimal digits of the object
 564         _digit(obj,n)   returns the n'th decimal digit of object
 565
 566         _is_one(obj)    return true if argument is +1
 567         _is_zero(obj)   return true if argument is 0
 568         _is_even(obj)   return true if argument is even (0,2,4,6..)
 569         _is_odd(obj)    return true if argument is odd (1,3,5,7..)
 570
 571         _copy           return a ref to a true copy of the object
 572
 573         _check(obj)     check whether internal representation is still intact
 574                         return 0 for ok, otherwise error message as string
 575
 576 The following functions are optional, and can be exported if the underlying lib
 577 has a fast way to do them. If not defined, Math::BigInt will use a pure, but
 578 slow, Perl function as fallback to emulate these:
 579
 580         _from_hex(str)  return ref to new object from ref to hexadecimal string
 581         _from_bin(str)  return ref to new object from ref to binary string
 582
 583         _rsft(obj,N,B)  shift object in base B by N 'digits' right
 584         _lsft(obj,N,B)  shift object in base B by N 'digits' left
 585
 586         _xor(obj1,obj2) XOR (bit-wise) object 1 with object 2
 587                         Mote: XOR, AND and OR pad with zeros if size mismatches
 588         _and(obj1,obj2) AND (bit-wise) object 1 with object 2
 589         _or(obj1,obj2)  OR (bit-wise) object 1 with object 2
 590
 591         _sqrt(obj)      return the square root of object
 592         _pow(obj,obj)   return object 1 to the power of object 2
 593         _gcd(obj,obj)   return Greatest Common Divisor of two objects
 594
 595         _zeros(obj)     return number of trailing decimal zeros
 596
 597         _dec(obj)       decrement object by one (input is >= 1)
 598         _inc(obj)       increment object by one
 599
 600 Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
 601 or '0b1101').
 602
 603 Testing of input parameter validity is done by the caller, so you need not
 604 worry about underflow (C<_sub()>, C<_dec()>) nor about division by zero or
 605 similar cases.
 606
 607 =head1 LICENSE
 608
 609 This program is free software; you may redistribute it and/or modify it under
 610 the same terms as Perl itself.
 611
 612 =head1 AUTHORS
 613
 614 Original math code by Mark Biggar, rewritten by Tels L<http://bloodgate.com/>
 615 in late 2000, 2001.
 616 Seperated from BigInt and shaped API with the help of John Peacock.
 617
 618 =head1 SEE ALSO
 619
 620 L<Math::BigInt>, L<Math::BigFloat>.
 621
 622 =cut