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
  20 qw(
  21 ""      stringify
  22 0+      numify)                 # Order of arguments unsignificant
  23 ;
  24
  25 $NaNOK=1;
  26
  27 sub new {
  28   my($class) = shift;
  29   my($foo) = bnorm(shift);
  30   die "Not a number initialized to Math::BigInt" if !$NaNOK && $foo eq "NaN";
  31   bless \$foo, $class;
  32 }
  33 sub stringify { "${$_[0]}" }
  34 sub numify { 0 + "${$_[0]}" }   # Not needed, additional overhead
  35                                 # comparing to direct compilation based on
  36                                 # stringify
  37 sub import {
  38   shift;
  39   return unless @_;
  40   die "unknown import: @_" unless @_ == 1 and $_[0] eq ':constant';
  41   overload::constant integer => sub {Math::BigInt->new(shift)};
  42 }
  43
  44 $zero = 0;
  45
  46
  47 # normalize string form of number.   Strip leading zeros.  Strip any
  48 #   white space and add a sign, if missing.
  49 # Strings that are not numbers result the value 'NaN'.
  50
  51 sub bnorm { #(num_str) return num_str
  52     local($_) = @_;
  53     s/\s+//g;                           # strip white space
  54     if (s/^([+-]?)0*(\d+)$/$1$2/) {     # test if number
  55         substr($_,$[,0) = '+' unless $1; # Add missing sign
  56         s/^-0/+0/;
  57         $_;
  58     } else {
  59         'NaN';
  60     }
  61 }
  62
  63 # Convert a number from string format to internal base 100000 format.
  64 #   Assumes normalized value as input.
  65 sub internal { #(num_str) return int_num_array
  66     local($d) = @_;
  67     ($is,$il) = (substr($d,$[,1),length($d)-2);
  68     substr($d,$[,1) = '';
  69     ($is, reverse(unpack("a" . ($il%5+1) . ("a5" x ($il/5)), $d)));
  70 }
  71
  72 # Convert a number from internal base 100000 format to string format.
  73 #   This routine scribbles all over input array.
  74 sub external { #(int_num_array) return num_str
  75     $es = shift;
  76     grep($_ > 9999 || ($_ = substr('0000'.$_,-5)), @_);   # zero pad
  77     &bnorm(join('', $es, reverse(@_)));    # reverse concat and normalize
  78 }
  79
  80 # Negate input value.
  81 sub bneg { #(num_str) return num_str
  82     local($_) = &bnorm(@_);
  83     return $_ if $_ eq '+0' or $_ eq 'NaN';
  84     vec($_,0,8) ^= ord('+') ^ ord('-');
  85     $_;
  86 }
  87
  88 # Returns the absolute value of the input.
  89 sub babs { #(num_str) return num_str
  90     &abs(&bnorm(@_));
  91 }
  92
  93 sub abs { # post-normalized abs for internal use
  94     local($_) = @_;
  95     s/^-/+/;
  96     $_;
  97 }
  98
  99 # Compares 2 values.  Returns one of undef, <0, =0, >0. (suitable for sort)
 100 sub bcmp { #(num_str, num_str) return cond_code
 101     local($x,$y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
 102     if ($x eq 'NaN') {
 103         undef;
 104     } elsif ($y eq 'NaN') {
 105         undef;
 106     } else {
 107         &cmp($x,$y) <=> 0;
 108     }
 109 }
 110
 111 sub cmp { # post-normalized compare for internal use
 112     local($cx, $cy) = @_;
 113
 114     return 0 if ($cx eq $cy);
 115
 116     local($sx, $sy) = (substr($cx, 0, 1), substr($cy, 0, 1));
 117     local($ld);
 118
 119     if ($sx eq '+') {
 120       return  1 if ($sy eq '-' || $cy eq '+0');
 121       $ld = length($cx) - length($cy);
 122       return $ld if ($ld);
 123       return $cx cmp $cy;
 124     } else { # $sx eq '-'
 125       return -1 if ($sy eq '+');
 126       $ld = length($cy) - length($cx);
 127       return $ld if ($ld);
 128       return $cy cmp $cx;
 129     }
 130 }
 131
 132 sub badd { #(num_str, num_str) return num_str
 133     local(*x, *y); ($x, $y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
 134     if ($x eq 'NaN') {
 135         'NaN';
 136     } elsif ($y eq 'NaN') {
 137         'NaN';
 138     } else {
 139         @x = &internal($x);             # convert to internal form
 140         @y = &internal($y);
 141         local($sx, $sy) = (shift @x, shift @y); # get signs
 142         if ($sx eq $sy) {
 143             &external($sx, &add(*x, *y)); # if same sign add
 144         } else {
 145             ($x, $y) = (&abs($x),&abs($y)); # make abs
 146             if (&cmp($y,$x) > 0) {
 147                 &external($sy, &sub(*y, *x));
 148             } else {
 149                 &external($sx, &sub(*x, *y));
 150             }
 151         }
 152     }
 153 }
 154
 155 sub bsub { #(num_str, num_str) return num_str
 156     &badd($_[$[],&bneg($_[$[+1]));
 157 }
 158
 159 # GCD -- Euclids algorithm Knuth Vol 2 pg 296
 160 sub bgcd { #(num_str, num_str) return num_str
 161     local($x,$y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
 162     if ($x eq 'NaN' || $y eq 'NaN') {
 163         'NaN';
 164     } else {
 165         ($x, $y) = ($y,&bmod($x,$y)) while $y ne '+0';
 166         $x;
 167     }
 168 }
 169
 170 # routine to add two base 1e5 numbers
 171 #   stolen from Knuth Vol 2 Algorithm A pg 231
 172 #   there are separate routines to add and sub as per Kunth pg 233
 173 sub add { #(int_num_array, int_num_array) return int_num_array
 174     local(*x, *y) = @_;
 175     $car = 0;
 176     for $x (@x) {
 177         last unless @y || $car;
 178         $x -= 1e5 if $car = (($x += (@y ? shift(@y) : 0) + $car) >= 1e5) ? 1 : 0;
 179     }
 180     for $y (@y) {
 181         last unless $car;
 182         $y -= 1e5 if $car = (($y += $car) >= 1e5) ? 1 : 0;
 183     }
 184     (@x, @y, $car);
 185 }
 186
 187 # subtract base 1e5 numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
 188 sub sub { #(int_num_array, int_num_array) return int_num_array
 189     local(*sx, *sy) = @_;
 190     $bar = 0;
 191     for $sx (@sx) {
 192         last unless @sy || $bar;
 193         $sx += 1e5 if $bar = (($sx -= (@sy ? shift(@sy) : 0) + $bar) < 0);
 194     }
 195     @sx;
 196 }
 197
 198 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
 199 sub bmul { #(num_str, num_str) return num_str
 200     local(*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
 201     if ($x eq 'NaN') {
 202         'NaN';
 203     } elsif ($y eq 'NaN') {
 204         'NaN';
 205     } else {
 206         @x = &internal($x);
 207         @y = &internal($y);
 208         &external(&mul(*x,*y));
 209     }
 210 }
 211
 212 # multiply two numbers in internal representation
 213 # destroys the arguments, supposes that two arguments are different
 214 sub mul { #(*int_num_array, *int_num_array) return int_num_array
 215     local(*x, *y) = (shift, shift);
 216     local($signr) = (shift @x ne shift @y) ? '-' : '+';
 217     @prod = ();
 218     for $x (@x) {
 219       ($car, $cty) = (0, $[);
 220       for $y (@y) {
 221         $prod = $x * $y + ($prod[$cty] || 0) + $car;
 222         $prod[$cty++] =
 223           $prod - ($car = int($prod * 1e-5)) * 1e5;
 224       }
 225       $prod[$cty] += $car if $car;
 226       $x = shift @prod;
 227     }
 228     ($signr, @x, @prod);
 229 }
 230
 231 # modulus
 232 sub bmod { #(num_str, num_str) return num_str
 233     (&bdiv(@_))[$[+1];
 234 }
 235
 236 sub bdiv { #(dividend: num_str, divisor: num_str) return num_str
 237     local (*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
 238     return wantarray ? ('NaN','NaN') : 'NaN'
 239         if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0');
 240     return wantarray ? ('+0',$x) : '+0' if (&cmp(&abs($x),&abs($y)) < 0);
 241     @x = &internal($x); @y = &internal($y);
 242     $srem = $y[$[];
 243     $sr = (shift @x ne shift @y) ? '-' : '+';
 244     $car = $bar = $prd = 0;
 245     if (($dd = int(1e5/($y[$#y]+1))) != 1) {
 246         for $x (@x) {
 247             $x = $x * $dd + $car;
 248             $x -= ($car = int($x * 1e-5)) * 1e5;
 249         }
 250         push(@x, $car); $car = 0;
 251         for $y (@y) {
 252             $y = $y * $dd + $car;
 253             $y -= ($car = int($y * 1e-5)) * 1e5;
 254         }
 255     }
 256     else {
 257         push(@x, 0);
 258     }
 259     @q = (); ($v2,$v1) = @y[-2,-1];
 260     $v2 = 0 unless $v2;
 261     while ($#x > $#y) {
 262         ($u2,$u1,$u0) = @x[-3..-1];
 263         $u2 = 0 unless $u2;
 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 conversion of
 404 constants the expression 2**100 will be calculated 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