lib/Math/BigInt.pm

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