Introduce NVef, NVff, and NVgf, use the middle one.
[perl.git] / lib / bigint.pl
1 package bigint;
2 #
3 # This library is no longer being maintained, and is included for backward
4 # compatibility with Perl 4 programs which may require it.
5 #
6 # In particular, this should not be used as an example of modern Perl
7 # programming techniques.
8 #
9 # Suggested alternative:  Math::BigInt
10 #
11 # arbitrary size integer math package
12 #
13 # by Mark Biggar
14 #
15 # Canonical Big integer value are strings of the form
16 #       /^[+-]\d+$/ with leading zeros suppressed
17 # Input values to these routines may be strings of the form
18 #       /^\s*[+-]?[\d\s]+$/.
19 # Examples:
20 #   '+0'                            canonical zero value
21 #   '   -123 123 123'               canonical value '-123123123'
22 #   '1 23 456 7890'                 canonical value '+1234567890'
23 # Output values always always in canonical form
24 #
25 # Actual math is done in an internal format consisting of an array
26 #   whose first element is the sign (/^[+-]$/) and whose remaining 
27 #   elements are base 100000 digits with the least significant digit first.
28 # The string 'NaN' is used to represent the result when input arguments 
29 #   are not numbers, as well as the result of dividing by zero
30 #
31 # routines provided are:
32 #
33 #   bneg(BINT) return BINT              negation
34 #   babs(BINT) return BINT              absolute value
35 #   bcmp(BINT,BINT) return CODE         compare numbers (undef,<0,=0,>0)
36 #   badd(BINT,BINT) return BINT         addition
37 #   bsub(BINT,BINT) return BINT         subtraction
38 #   bmul(BINT,BINT) return BINT         multiplication
39 #   bdiv(BINT,BINT) return (BINT,BINT)  division (quo,rem) just quo if scalar
40 #   bmod(BINT,BINT) return BINT         modulus
41 #   bgcd(BINT,BINT) return BINT         greatest common divisor
42 #   bnorm(BINT) return BINT             normalization
43 #
44
45 $zero = 0;
46
47 \f
48 # normalize string form of number.   Strip leading zeros.  Strip any
49 #   white space and add a sign, if missing.
50 # Strings that are not numbers result the value 'NaN'.
51
52 sub main'bnorm { #(num_str) return num_str
53     local($_) = @_;
54     s/\s+//g;                           # strip white space
55     if (s/^([+-]?)0*(\d+)$/$1$2/) {     # test if number
56         substr($_,$[,0) = '+' unless $1; # Add missing sign
57         s/^-0/+0/;
58         $_;
59     } else {
60         'NaN';
61     }
62 }
63
64 # Convert a number from string format to internal base 100000 format.
65 #   Assumes normalized value as input.
66 sub internal { #(num_str) return int_num_array
67     local($d) = @_;
68     ($is,$il) = (substr($d,$[,1),length($d)-2);
69     substr($d,$[,1) = '';
70     ($is, reverse(unpack("a" . ($il%5+1) . ("a5" x ($il/5)), $d)));
71 }
72
73 # Convert a number from internal base 100000 format to string format.
74 #   This routine scribbles all over input array.
75 sub external { #(int_num_array) return num_str
76     $es = shift;
77     grep($_ > 9999 || ($_ = substr('0000'.$_,-5)), @_);   # zero pad
78     &'bnorm(join('', $es, reverse(@_)));    # reverse concat and normalize
79 }
80
81 # Negate input value.
82 sub main'bneg { #(num_str) return num_str
83     local($_) = &'bnorm(@_);
84     vec($_,0,8) ^= ord('+') ^ ord('-') unless $_ eq '+0';
85     s/^./N/ unless /^[-+]/; # works both in ASCII and EBCDIC
86     $_;
87 }
88
89 # Returns the absolute value of the input.
90 sub main'babs { #(num_str) return num_str
91     &abs(&'bnorm(@_));
92 }
93
94 sub abs { # post-normalized abs for internal use
95     local($_) = @_;
96     s/^-/+/;
97     $_;
98 }
99 \f
100 # Compares 2 values.  Returns one of undef, <0, =0, >0. (suitable for sort)
101 sub main'bcmp { #(num_str, num_str) return cond_code
102     local($x,$y) = (&'bnorm($_[$[]),&'bnorm($_[$[+1]));
103     if ($x eq 'NaN') {
104         undef;
105     } elsif ($y eq 'NaN') {
106         undef;
107     } else {
108         &cmp($x,$y);
109     }
110 }
111
112 sub cmp { # post-normalized compare for internal use
113     local($cx, $cy) = @_;
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
133 sub main'badd { #(num_str, num_str) return num_str
134     local(*x, *y); ($x, $y) = (&'bnorm($_[$[]),&'bnorm($_[$[+1]));
135     if ($x eq 'NaN') {
136         'NaN';
137     } elsif ($y eq 'NaN') {
138         'NaN';
139     } else {
140         @x = &internal($x);             # convert to internal form
141         @y = &internal($y);
142         local($sx, $sy) = (shift @x, shift @y); # get signs
143         if ($sx eq $sy) {
144             &external($sx, &add(*x, *y)); # if same sign add
145         } else {
146             ($x, $y) = (&abs($x),&abs($y)); # make abs
147             if (&cmp($y,$x) > 0) {
148                 &external($sy, &sub(*y, *x));
149             } else {
150                 &external($sx, &sub(*x, *y));
151             }
152         }
153     }
154 }
155
156 sub main'bsub { #(num_str, num_str) return num_str
157     &'badd($_[$[],&'bneg($_[$[+1]));    
158 }
159
160 # GCD -- Euclids algorithm Knuth Vol 2 pg 296
161 sub main'bgcd { #(num_str, num_str) return num_str
162     local($x,$y) = (&'bnorm($_[$[]),&'bnorm($_[$[+1]));
163     if ($x eq 'NaN' || $y eq 'NaN') {
164         'NaN';
165     } else {
166         ($x, $y) = ($y,&'bmod($x,$y)) while $y ne '+0';
167         $x;
168     }
169 }
170 \f
171 # routine to add two base 1e5 numbers
172 #   stolen from Knuth Vol 2 Algorithm A pg 231
173 #   there are separate routines to add and sub as per Kunth pg 233
174 sub add { #(int_num_array, int_num_array) return int_num_array
175     local(*x, *y) = @_;
176     $car = 0;
177     for $x (@x) {
178         last unless @y || $car;
179         $x -= 1e5 if $car = (($x += shift(@y) + $car) >= 1e5) ? 1 : 0;
180     }
181     for $y (@y) {
182         last unless $car;
183         $y -= 1e5 if $car = (($y += $car) >= 1e5) ? 1 : 0;
184     }
185     (@x, @y, $car);
186 }
187
188 # subtract base 1e5 numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
189 sub sub { #(int_num_array, int_num_array) return int_num_array
190     local(*sx, *sy) = @_;
191     $bar = 0;
192     for $sx (@sx) {
193         last unless @y || $bar;
194         $sx += 1e5 if $bar = (($sx -= shift(@sy) + $bar) < 0);
195     }
196     @sx;
197 }
198
199 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
200 sub main'bmul { #(num_str, num_str) return num_str
201     local(*x, *y); ($x, $y) = (&'bnorm($_[$[]), &'bnorm($_[$[+1]));
202     if ($x eq 'NaN') {
203         'NaN';
204     } elsif ($y eq 'NaN') {
205         'NaN';
206     } else {
207         @x = &internal($x);
208         @y = &internal($y);
209         local($signr) = (shift @x ne shift @y) ? '-' : '+';
210         @prod = ();
211         for $x (@x) {
212             ($car, $cty) = (0, $[);
213             for $y (@y) {
214                 $prod = $x * $y + $prod[$cty] + $car;
215                 $prod[$cty++] =
216                     $prod - ($car = int($prod * 1e-5)) * 1e5;
217             }
218             $prod[$cty] += $car if $car;
219             $x = shift @prod;
220         }
221         &external($signr, @x, @prod);
222     }
223 }
224
225 # modulus
226 sub main'bmod { #(num_str, num_str) return num_str
227     (&'bdiv(@_))[$[+1];
228 }
229 \f
230 sub main'bdiv { #(dividend: num_str, divisor: num_str) return num_str
231     local (*x, *y); ($x, $y) = (&'bnorm($_[$[]), &'bnorm($_[$[+1]));
232     return wantarray ? ('NaN','NaN') : 'NaN'
233         if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0');
234     return wantarray ? ('+0',$x) : '+0' if (&cmp(&abs($x),&abs($y)) < 0);
235     @x = &internal($x); @y = &internal($y);
236     $srem = $y[$[];
237     $sr = (shift @x ne shift @y) ? '-' : '+';
238     $car = $bar = $prd = 0;
239     if (($dd = int(1e5/($y[$#y]+1))) != 1) {
240         for $x (@x) {
241             $x = $x * $dd + $car;
242             $x -= ($car = int($x * 1e-5)) * 1e5;
243         }
244         push(@x, $car); $car = 0;
245         for $y (@y) {
246             $y = $y * $dd + $car;
247             $y -= ($car = int($y * 1e-5)) * 1e5;
248         }
249     }
250     else {
251         push(@x, 0);
252     }
253     @q = (); ($v2,$v1) = @y[-2,-1];
254     while ($#x > $#y) {
255         ($u2,$u1,$u0) = @x[-3..-1];
256         $q = (($u0 == $v1) ? 99999 : int(($u0*1e5+$u1)/$v1));
257         --$q while ($v2*$q > ($u0*1e5+$u1-$q*$v1)*1e5+$u2);
258         if ($q) {
259             ($car, $bar) = (0,0);
260             for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
261                 $prd = $q * $y[$y] + $car;
262                 $prd -= ($car = int($prd * 1e-5)) * 1e5;
263                 $x[$x] += 1e5 if ($bar = (($x[$x] -= $prd + $bar) < 0));
264             }
265             if ($x[$#x] < $car + $bar) {
266                 $car = 0; --$q;
267                 for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
268                     $x[$x] -= 1e5
269                         if ($car = (($x[$x] += $y[$y] + $car) > 1e5));
270                 }
271             }   
272         }
273         pop(@x); unshift(@q, $q);
274     }
275     if (wantarray) {
276         @d = ();
277         if ($dd != 1) {
278             $car = 0;
279             for $x (reverse @x) {
280                 $prd = $car * 1e5 + $x;
281                 $car = $prd - ($tmp = int($prd / $dd)) * $dd;
282                 unshift(@d, $tmp);
283             }
284         }
285         else {
286             @d = @x;
287         }
288         (&external($sr, @q), &external($srem, @d, $zero));
289     } else {
290         &external($sr, @q);
291     }
292 }
293 1;