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