lib/Math/BigInt.pm

   1 package Math::BigInt;
   2 $VERSION='0.01';
   3
   4 use overload
   5 '+'     =>      sub {new Math::BigInt &badd},
   6 '-'     =>      sub {new Math::BigInt
   7                        $_[2]? bsub($_[1],${$_[0]}) : bsub(${$_[0]},$_[1])},
   8 '<=>'   =>      sub {$_[2]? bcmp($_[1],${$_[0]}) : bcmp(${$_[0]},$_[1])},
   9 'cmp'   =>      sub {$_[2]? ($_[1] cmp ${$_[0]}) : (${$_[0]} cmp $_[1])},
  10 '*'     =>      sub {new Math::BigInt &bmul},
  11 '/'     =>      sub {new Math::BigInt
  12                        $_[2]? scalar bdiv($_[1],${$_[0]}) :
  13                          scalar bdiv(${$_[0]},$_[1])},
  14 '%'     =>      sub {new Math::BigInt
  15                        $_[2]? bmod($_[1],${$_[0]}) : bmod(${$_[0]},$_[1])},
  16 '**'    =>      sub {new Math::BigInt
  17                        $_[2]? bpow($_[1],${$_[0]}) : bpow(${$_[0]},$_[1])},
  18 'neg'   =>      sub {new Math::BigInt &bneg},
  19 'abs'   =>      sub {new Math::BigInt &babs},
  20 '<<'    =>      sub {new Math::BigInt
  21                        $_[2]? blsft($_[1],${$_[0]}) : blsft(${$_[0]},$_[1])},
  22 '>>'    =>      sub {new Math::BigInt
  23                        $_[2]? brsft($_[1],${$_[0]}) : brsft(${$_[0]},$_[1])},
  24 '&'     =>      sub {new Math::BigInt &band},
  25 '|'     =>      sub {new Math::BigInt &bior},
  26 '^'     =>      sub {new Math::BigInt &bxor},
  27 '~'     =>      sub {new Math::BigInt &bnot},
  28
  29 qw(
  30 ""      stringify
  31 0+      numify)                 # Order of arguments unsignificant
  32 ;
  33
  34 $NaNOK=1;
  35
  36 sub new {
  37   my($class) = shift;
  38   my($foo) = bnorm(shift);
  39   die "Not a number initialized to Math::BigInt" if !$NaNOK && $foo eq "NaN";
  40   bless \$foo, $class;
  41 }
  42 sub stringify { "${$_[0]}" }
  43 sub numify { 0 + "${$_[0]}" }   # Not needed, additional overhead
  44                                 # comparing to direct compilation based on
  45                                 # stringify
  46 sub import {
  47   shift;
  48   return unless @_;
  49   die "unknown import: @_" unless @_ == 1 and $_[0] eq ':constant';
  50   overload::constant integer => sub {Math::BigInt->new(shift)};
  51 }
  52
  53 $zero = 0;
  54
  55
  56 # normalize string form of number.   Strip leading zeros.  Strip any
  57 #   white space and add a sign, if missing.
  58 # Strings that are not numbers result the value 'NaN'.
  59
  60 sub bnorm { #(num_str) return num_str
  61     local($_) = @_;
  62     s/\s+//g;                           # strip white space
  63     if (s/^([+-]?)0*(\d+)$/$1$2/) {     # test if number
  64         substr($_,$[,0) = '+' unless $1; # Add missing sign
  65         s/^-0/+0/;
  66         $_;
  67     } else {
  68         'NaN';
  69     }
  70 }
  71
  72 # Convert a number from string format to internal base 100000 format.
  73 #   Assumes normalized value as input.
  74 sub internal { #(num_str) return int_num_array
  75     local($d) = @_;
  76     ($is,$il) = (substr($d,$[,1),length($d)-2);
  77     substr($d,$[,1) = '';
  78     ($is, reverse(unpack("a" . ($il%5+1) . ("a5" x ($il/5)), $d)));
  79 }
  80
  81 # Convert a number from internal base 100000 format to string format.
  82 #   This routine scribbles all over input array.
  83 sub external { #(int_num_array) return num_str
  84     $es = shift;
  85     grep($_ > 9999 || ($_ = substr('0000'.$_,-5)), @_);   # zero pad
  86     &bnorm(join('', $es, reverse(@_)));    # reverse concat and normalize
  87 }
  88
  89 # Negate input value.
  90 sub bneg { #(num_str) return num_str
  91     local($_) = &bnorm(@_);
  92     return $_ if $_ eq '+0' or $_ eq 'NaN';
  93     vec($_,0,8) ^= ord('+') ^ ord('-');
  94     $_;
  95 }
  96
  97 # Returns the absolute value of the input.
  98 sub babs { #(num_str) return num_str
  99     &abs(&bnorm(@_));
 100 }
 101
 102 sub abs { # post-normalized abs for internal use
 103     local($_) = @_;
 104     s/^-/+/;
 105     $_;
 106 }
 107
 108 # Compares 2 values.  Returns one of undef, <0, =0, >0. (suitable for sort)
 109 sub bcmp { #(num_str, num_str) return cond_code
 110     local($x,$y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
 111     if ($x eq 'NaN') {
 112         undef;
 113     } elsif ($y eq 'NaN') {
 114         undef;
 115     } else {
 116         &cmp($x,$y) <=> 0;
 117     }
 118 }
 119
 120 sub cmp { # post-normalized compare for internal use
 121     local($cx, $cy) = @_;
 122
 123     return 0 if ($cx eq $cy);
 124
 125     local($sx, $sy) = (substr($cx, 0, 1), substr($cy, 0, 1));
 126     local($ld);
 127
 128     if ($sx eq '+') {
 129       return  1 if ($sy eq '-' || $cy eq '+0');
 130       $ld = length($cx) - length($cy);
 131       return $ld if ($ld);
 132       return $cx cmp $cy;
 133     } else { # $sx eq '-'
 134       return -1 if ($sy eq '+');
 135       $ld = length($cy) - length($cx);
 136       return $ld if ($ld);
 137       return $cy cmp $cx;
 138     }
 139 }
 140
 141 sub badd { #(num_str, num_str) return num_str
 142     local(*x, *y); ($x, $y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
 143     if ($x eq 'NaN') {
 144         'NaN';
 145     } elsif ($y eq 'NaN') {
 146         'NaN';
 147     } else {
 148         @x = &internal($x);             # convert to internal form
 149         @y = &internal($y);
 150         local($sx, $sy) = (shift @x, shift @y); # get signs
 151         if ($sx eq $sy) {
 152             &external($sx, &add(*x, *y)); # if same sign add
 153         } else {
 154             ($x, $y) = (&abs($x),&abs($y)); # make abs
 155             if (&cmp($y,$x) > 0) {
 156                 &external($sy, &sub(*y, *x));
 157             } else {
 158                 &external($sx, &sub(*x, *y));
 159             }
 160         }
 161     }
 162 }
 163
 164 sub bsub { #(num_str, num_str) return num_str
 165     &badd($_[$[],&bneg($_[$[+1]));
 166 }
 167
 168 # GCD -- Euclids algorithm Knuth Vol 2 pg 296
 169 sub bgcd { #(num_str, num_str) return num_str
 170     local($x,$y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
 171     if ($x eq 'NaN' || $y eq 'NaN') {
 172         'NaN';
 173     } else {
 174         ($x, $y) = ($y,&bmod($x,$y)) while $y ne '+0';
 175         $x;
 176     }
 177 }
 178
 179 # routine to add two base 1e5 numbers
 180 #   stolen from Knuth Vol 2 Algorithm A pg 231
 181 #   there are separate routines to add and sub as per Kunth pg 233
 182 sub add { #(int_num_array, int_num_array) return int_num_array
 183     local(*x, *y) = @_;
 184     $car = 0;
 185     for $x (@x) {
 186         last unless @y || $car;
 187         $x -= 1e5 if $car = (($x += (@y ? shift(@y) : 0) + $car) >= 1e5) ? 1 : 0;
 188     }
 189     for $y (@y) {
 190         last unless $car;
 191         $y -= 1e5 if $car = (($y += $car) >= 1e5) ? 1 : 0;
 192     }
 193     (@x, @y, $car);
 194 }
 195
 196 # subtract base 1e5 numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
 197 sub sub { #(int_num_array, int_num_array) return int_num_array
 198     local(*sx, *sy) = @_;
 199     $bar = 0;
 200     for $sx (@sx) {
 201         last unless @sy || $bar;
 202         $sx += 1e5 if $bar = (($sx -= (@sy ? shift(@sy) : 0) + $bar) < 0);
 203     }
 204     @sx;
 205 }
 206
 207 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
 208 sub bmul { #(num_str, num_str) return num_str
 209     local(*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
 210     if ($x eq 'NaN') {
 211         'NaN';
 212     } elsif ($y eq 'NaN') {
 213         'NaN';
 214     } else {
 215         @x = &internal($x);
 216         @y = &internal($y);
 217         &external(&mul(*x,*y));
 218     }
 219 }
 220
 221 # multiply two numbers in internal representation
 222 # destroys the arguments, supposes that two arguments are different
 223 sub mul { #(*int_num_array, *int_num_array) return int_num_array
 224     local(*x, *y) = (shift, shift);
 225     local($signr) = (shift @x ne shift @y) ? '-' : '+';
 226     @prod = ();
 227     for $x (@x) {
 228       ($car, $cty) = (0, $[);
 229       for $y (@y) {
 230         $prod = $x * $y + ($prod[$cty] || 0) + $car;
 231         $prod[$cty++] =
 232           $prod - ($car = int($prod * 1e-5)) * 1e5;
 233       }
 234       $prod[$cty] += $car if $car;
 235       $x = shift @prod;
 236     }
 237     ($signr, @x, @prod);
 238 }
 239
 240 # modulus
 241 sub bmod { #(num_str, num_str) return num_str
 242     (&bdiv(@_))[$[+1];
 243 }
 244
 245 sub bdiv { #(dividend: num_str, divisor: num_str) return num_str
 246     local (*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
 247     return wantarray ? ('NaN','NaN') : 'NaN'
 248         if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0');
 249     return wantarray ? ('+0',$x) : '+0' if (&cmp(&abs($x),&abs($y)) < 0);
 250     @x = &internal($x); @y = &internal($y);
 251     $srem = $y[$[];
 252     $sr = (shift @x ne shift @y) ? '-' : '+';
 253     $car = $bar = $prd = 0;
 254     if (($dd = int(1e5/($y[$#y]+1))) != 1) {
 255         for $x (@x) {
 256             $x = $x * $dd + $car;
 257             $x -= ($car = int($x * 1e-5)) * 1e5;
 258         }
 259         push(@x, $car); $car = 0;
 260         for $y (@y) {
 261             $y = $y * $dd + $car;
 262             $y -= ($car = int($y * 1e-5)) * 1e5;
 263         }
 264     }
 265     else {
 266         push(@x, 0);
 267     }
 268     @q = (); ($v2,$v1) = @y[-2,-1];
 269     $v2 = 0 unless $v2;
 270     while ($#x > $#y) {
 271         ($u2,$u1,$u0) = @x[-3..-1];
 272         $u2 = 0 unless $u2;
 273         $q = (($u0 == $v1) ? 99999 : int(($u0*1e5+$u1)/$v1));
 274         --$q while ($v2*$q > ($u0*1e5+$u1-$q*$v1)*1e5+$u2);
 275         if ($q) {
 276             ($car, $bar) = (0,0);
 277             for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
 278                 $prd = $q * $y[$y] + $car;
 279                 $prd -= ($car = int($prd * 1e-5)) * 1e5;
 280                 $x[$x] += 1e5 if ($bar = (($x[$x] -= $prd + $bar) < 0));
 281             }
 282             if ($x[$#x] < $car + $bar) {
 283                 $car = 0; --$q;
 284                 for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
 285                     $x[$x] -= 1e5
 286                         if ($car = (($x[$x] += $y[$y] + $car) > 1e5));
 287                 }
 288             }
 289         }
 290         pop(@x); unshift(@q, $q);
 291     }
 292     if (wantarray) {
 293         @d = ();
 294         if ($dd != 1) {
 295             $car = 0;
 296             for $x (reverse @x) {
 297                 $prd = $car * 1e5 + $x;
 298                 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
 299                 unshift(@d, $tmp);
 300             }
 301         }
 302         else {
 303             @d = @x;
 304         }
 305         (&external($sr, @q), &external($srem, @d, $zero));
 306     } else {
 307         &external($sr, @q);
 308     }
 309 }
 310
 311 # compute power of two numbers -- stolen from Knuth Vol 2 pg 233
 312 sub bpow { #(num_str, num_str) return num_str
 313     local(*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
 314     if ($x eq 'NaN') {
 315         'NaN';
 316     } elsif ($y eq 'NaN') {
 317         'NaN';
 318     } elsif ($x eq '+1') {
 319         '+1';
 320     } elsif ($x eq '-1') {
 321         &bmod($x,2) ? '-1': '+1';
 322     } elsif ($y =~ /^-/) {
 323         'NaN';
 324     } elsif ($x eq '+0' && $y eq '+0') {
 325         'NaN';
 326     } else {
 327         @x = &internal($x);
 328         local(@pow2)=@x;
 329         local(@pow)=&internal("+1");
 330         local($y1,$res,@tmp1,@tmp2)=(1); # need tmp to send to mul
 331         while ($y ne '+0') {
 332           ($y,$res)=&bdiv($y,2);
 333           if ($res ne '+0') {@tmp=@pow2; @pow=&mul(*pow,*tmp);}
 334           if ($y ne '+0') {@tmp=@pow2;@pow2=&mul(*pow2,*tmp);}
 335         }
 336         &external(@pow);
 337     }
 338 }
 339
 340 # compute x << y, y >= 0
 341 sub blsft { #(num_str, num_str) return num_str
 342     &bmul($_[$[], &bpow(2, $_[$[+1]));
 343 }
 344
 345 # compute x >> y, y >= 0
 346 sub brsft { #(num_str, num_str) return num_str
 347     &bdiv($_[$[], &bpow(2, $_[$[+1]));
 348 }
 349
 350 # compute x & y
 351 sub band { #(num_str, num_str) return num_str
 352     local($x,$y,$r,$m,$xr,$yr) = (&bnorm($_[$[]),&bnorm($_[$[+1]),0,1);
 353     if ($x eq 'NaN' || $y eq 'NaN') {
 354         'NaN';
 355     } else {
 356         while ($x ne '+0' && $y ne '+0') {
 357             ($x, $xr) = &bdiv($x, 0x10000);
 358             ($y, $yr) = &bdiv($y, 0x10000);
 359             $r = &badd(&bmul(int $xr & $yr, $m), $r);
 360             $m = &bmul($m, 0x10000);
 361         }
 362         $r;
 363     }
 364 }
 365
 366 # compute x | y
 367 sub bior { #(num_str, num_str) return num_str
 368     local($x,$y,$r,$m,$xr,$yr) = (&bnorm($_[$[]),&bnorm($_[$[+1]),0,1);
 369     if ($x eq 'NaN' || $y eq 'NaN') {
 370         'NaN';
 371     } else {
 372         while ($x ne '+0' || $y ne '+0') {
 373             ($x, $xr) = &bdiv($x, 0x10000);
 374             ($y, $yr) = &bdiv($y, 0x10000);
 375             $r = &badd(&bmul(int $xr | $yr, $m), $r);
 376             $m = &bmul($m, 0x10000);
 377         }
 378         $r;
 379     }
 380 }
 381
 382 # compute x ^ y
 383 sub bxor { #(num_str, num_str) return num_str
 384     local($x,$y,$r,$m,$xr,$yr) = (&bnorm($_[$[]),&bnorm($_[$[+1]),0,1);
 385     if ($x eq 'NaN' || $y eq 'NaN') {
 386         'NaN';
 387     } else {
 388         while ($x ne '+0' || $y ne '+0') {
 389             ($x, $xr) = &bdiv($x, 0x10000);
 390             ($y, $yr) = &bdiv($y, 0x10000);
 391             $r = &badd(&bmul(int $xr ^ $yr, $m), $r);
 392             $m = &bmul($m, 0x10000);
 393         }
 394         $r;
 395     }
 396 }
 397
 398 # represent ~x as twos-complement number
 399 sub bnot { #(num_str) return num_str
 400     &bsub(-1,$_[$[]);
 401 }
 402
 403 1;
 404 __END__
 405
 406 =head1 NAME
 407
 408 Math::BigInt - Arbitrary size integer math package
 409
 410 =head1 SYNOPSIS
 411
 412   use Math::BigInt;
 413   $i = Math::BigInt->new($string);
 414
 415   $i->bneg return BINT               negation
 416   $i->babs return BINT               absolute value
 417   $i->bcmp(BINT) return CODE         compare numbers (undef,<0,=0,>0)
 418   $i->badd(BINT) return BINT         addition
 419   $i->bsub(BINT) return BINT         subtraction
 420   $i->bmul(BINT) return BINT         multiplication
 421   $i->bdiv(BINT) return (BINT,BINT)  division (quo,rem) just quo if scalar
 422   $i->bmod(BINT) return BINT         modulus
 423   $i->bgcd(BINT) return BINT         greatest common divisor
 424   $i->bnorm return BINT              normalization
 425   $i->blsft(BINT) return BINT        left shift
 426   $i->brsft(BINT) return (BINT,BINT) right shift (quo,rem) just quo if scalar
 427   $i->band(BINT) return BINT         bit-wise and
 428   $i->bior(BINT) return BINT         bit-wise inclusive or
 429   $i->bxor(BINT) return BINT         bit-wise exclusive or
 430   $i->bnot return BINT               bit-wise not
 431
 432 =head1 DESCRIPTION
 433
 434 All basic math operations are overloaded if you declare your big
 435 integers as
 436
 437   $i = new Math::BigInt '123 456 789 123 456 789';
 438
 439
 440 =over 2
 441
 442 =item Canonical notation
 443
 444 Big integer value are strings of the form C</^[+-]\d+$/> with leading
 445 zeros suppressed.
 446
 447 =item Input
 448
 449 Input values to these routines may be strings of the form
 450 C</^\s*[+-]?[\d\s]+$/>.
 451
 452 =item Output
 453
 454 Output values always always in canonical form
 455
 456 =back
 457
 458 Actual math is done in an internal format consisting of an array
 459 whose first element is the sign (/^[+-]$/) and whose remaining
 460 elements are base 100000 digits with the least significant digit first.
 461 The string 'NaN' is used to represent the result when input arguments
 462 are not numbers, as well as the result of dividing by zero.
 463
 464 =head1 EXAMPLES
 465
 466    '+0'                            canonical zero value
 467    '   -123 123 123'               canonical value '-123123123'
 468    '1 23 456 7890'                 canonical value '+1234567890'
 469
 470
 471 =head1 Autocreating constants
 472
 473 After C<use Math::BigInt ':constant'> all the integer decimal constants
 474 in the given scope are converted to C<Math::BigInt>.  This conversion
 475 happens at compile time.
 476
 477 In particular
 478
 479   perl -MMath::BigInt=:constant -e 'print 2**100'
 480
 481 print the integer value of C<2**100>.  Note that without conversion of
 482 constants the expression 2**100 will be calculated as floating point number.
 483
 484 =head1 BUGS
 485
 486 The current version of this module is a preliminary version of the
 487 real thing that is currently (as of perl5.002) under development.
 488
 489 =head1 AUTHOR
 490
 491 Mark Biggar, overloaded interface by Ilya Zakharevich.
 492
 493 =cut