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