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     while ($#x > $#y) {
 261         ($u2,$u1,$u0) = @x[-3..-1];
 262         $q = (($u0 == $v1) ? 99999 : int(($u0*1e5+$u1)/$v1));
 263         --$q while ($v2*$q > ($u0*1e5+$u1-$q*$v1)*1e5+$u2);
 264         if ($q) {
 265             ($car, $bar) = (0,0);
 266             for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
 267                 $prd = $q * $y[$y] + $car;
 268                 $prd -= ($car = int($prd * 1e-5)) * 1e5;
 269                 $x[$x] += 1e5 if ($bar = (($x[$x] -= $prd + $bar) < 0));
 270             }
 271             if ($x[$#x] < $car + $bar) {
 272                 $car = 0; --$q;
 273                 for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
 274                     $x[$x] -= 1e5
 275                         if ($car = (($x[$x] += $y[$y] + $car) > 1e5));
 276                 }
 277             }
 278         }
 279         pop(@x); unshift(@q, $q);
 280     }
 281     if (wantarray) {
 282         @d = ();
 283         if ($dd != 1) {
 284             $car = 0;
 285             for $x (reverse @x) {
 286                 $prd = $car * 1e5 + $x;
 287                 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
 288                 unshift(@d, $tmp);
 289             }
 290         }
 291         else {
 292             @d = @x;
 293         }
 294         (&external($sr, @q), &external($srem, @d, $zero));
 295     } else {
 296         &external($sr, @q);
 297     }
 298 }
 299
 300 # compute power of two numbers -- stolen from Knuth Vol 2 pg 233
 301 sub bpow { #(num_str, num_str) return num_str
 302     local(*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
 303     if ($x eq 'NaN') {
 304         'NaN';
 305     } elsif ($y eq 'NaN') {
 306         'NaN';
 307     } elsif ($x eq '+1') {
 308         '+1';
 309     } elsif ($x eq '-1') {
 310         &bmod($x,2) ? '-1': '+1';
 311     } elsif ($y =~ /^-/) {
 312         'NaN';
 313     } elsif ($x eq '+0' && $y eq '+0') {
 314         'NaN';
 315     } else {
 316         @x = &internal($x);
 317         local(@pow2)=@x;
 318         local(@pow)=&internal("+1");
 319         local($y1,$res,@tmp1,@tmp2)=(1); # need tmp to send to mul
 320         while ($y ne '+0') {
 321           ($y,$res)=&bdiv($y,2);
 322           if ($res ne '+0') {@tmp=@pow2; @pow=&mul(*pow,*tmp);}
 323           if ($y ne '+0') {@tmp=@pow2;@pow2=&mul(*pow2,*tmp);}
 324         }
 325         &external(@pow);
 326     }
 327 }
 328
 329 1;
 330 __END__
 331
 332 =head1 NAME
 333
 334 Math::BigInt - Arbitrary size integer math package
 335
 336 =head1 SYNOPSIS
 337
 338   use Math::BigInt;
 339   $i = Math::BigInt->new($string);
 340
 341   $i->bneg return BINT               negation
 342   $i->babs return BINT               absolute value
 343   $i->bcmp(BINT) return CODE         compare numbers (undef,<0,=0,>0)
 344   $i->badd(BINT) return BINT         addition
 345   $i->bsub(BINT) return BINT         subtraction
 346   $i->bmul(BINT) return BINT         multiplication
 347   $i->bdiv(BINT) return (BINT,BINT)  division (quo,rem) just quo if scalar
 348   $i->bmod(BINT) return BINT         modulus
 349   $i->bgcd(BINT) return BINT         greatest common divisor
 350   $i->bnorm return BINT              normalization
 351
 352 =head1 DESCRIPTION
 353
 354 All basic math operations are overloaded if you declare your big
 355 integers as
 356
 357   $i = new Math::BigInt '123 456 789 123 456 789';
 358
 359
 360 =over 2
 361
 362 =item Canonical notation
 363
 364 Big integer value are strings of the form C</^[+-]\d+$/> with leading
 365 zeros suppressed.
 366
 367 =item Input
 368
 369 Input values to these routines may be strings of the form
 370 C</^\s*[+-]?[\d\s]+$/>.
 371
 372 =item Output
 373
 374 Output values always always in canonical form
 375
 376 =back
 377
 378 Actual math is done in an internal format consisting of an array
 379 whose first element is the sign (/^[+-]$/) and whose remaining
 380 elements are base 100000 digits with the least significant digit first.
 381 The string 'NaN' is used to represent the result when input arguments
 382 are not numbers, as well as the result of dividing by zero.
 383
 384 =head1 EXAMPLES
 385
 386    '+0'                            canonical zero value
 387    '   -123 123 123'               canonical value '-123123123'
 388    '1 23 456 7890'                 canonical value '+1234567890'
 389
 390
 391 =head1 Autocreating constants
 392
 393 After C<use Math::BigInt ':constant'> all the integer decimal constants
 394 in the given scope are converted to C<Math::BigInt>.  This conversion
 395 happens at compile time.
 396
 397 In particular
 398
 399   perl -MMath::BigInt=:constant -e 'print 2**100'
 400
 401 print the integer value of C<2**100>.  Note that without conversion of
 402 constants the expression 2**100 will be calculated as floating point number.
 403
 404 =head1 BUGS
 405
 406 The current version of this module is a preliminary version of the
 407 real thing that is currently (as of perl5.002) under development.
 408
 409 =head1 AUTHOR
 410
 411 Mark Biggar, overloaded interface by Ilya Zakharevich.
 412
 413 =cut