This is a live mirror of the Perl 5 development currently hosted at https://github.com/perl/perl5
squelch undef warnings
[perl5.git] / 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] || 0, $y[-1]);
262     while ($#x > $#y) {
263         ($u2,$u1,$u0) = ($x[-3] || 0, $x[-2] || 0, $x[-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