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