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