lib/Math/BigInt.pm

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