1 : /* recmul.c: bcmath library file. */
2 : /*
3 : Copyright (C) 1991, 1992, 1993, 1994, 1997 Free Software Foundation, Inc.
4 : Copyright (C) 2000 Philip A. Nelson
5 :
6 : This library is free software; you can redistribute it and/or
7 : modify it under the terms of the GNU Lesser General Public
8 : License as published by the Free Software Foundation; either
9 : version 2 of the License, or (at your option) any later version.
10 :
11 : This library is distributed in the hope that it will be useful,
12 : but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 : Lesser General Public License for more details. (COPYING.LIB)
15 :
16 : You should have received a copy of the GNU Lesser General Public
17 : License along with this library; if not, write to:
18 :
19 : The Free Software Foundation, Inc.
20 : 59 Temple Place, Suite 330
21 : Boston, MA 02111-1307 USA.
22 :
23 : You may contact the author by:
24 : e-mail: philnelson@acm.org
25 : us-mail: Philip A. Nelson
26 : Computer Science Department, 9062
27 : Western Washington University
28 : Bellingham, WA 98226-9062
29 :
30 : *************************************************************************/
31 :
32 : #include <config.h>
33 : #include <stdio.h>
34 : #include <assert.h>
35 : #include <stdlib.h>
36 : #include <ctype.h>
37 : #include <stdarg.h>
38 : #include "bcmath.h"
39 : #include "private.h"
40 :
41 : /* Recursive vs non-recursive multiply crossover ranges. */
42 : #if defined(MULDIGITS)
43 : #include "muldigits.h"
44 : #else
45 : #define MUL_BASE_DIGITS 80
46 : #endif
47 :
48 : int mul_base_digits = MUL_BASE_DIGITS;
49 : #define MUL_SMALL_DIGITS mul_base_digits/4
50 :
51 : /* Multiply utility routines */
52 :
53 : static bc_num
54 : new_sub_num (length, scale, value)
55 : int length, scale;
56 : char *value;
57 0 : {
58 : bc_num temp;
59 :
60 : #ifdef SANDER_0
61 : if (_bc_Free_list != NULL) {
62 : temp = _bc_Free_list;
63 : _bc_Free_list = temp->n_next;
64 : } else {
65 : #endif
66 0 : temp = (bc_num) emalloc (sizeof(bc_struct));
67 : #ifdef SANDER_0
68 : if (temp == NULL) bc_out_of_memory ();
69 : }
70 : #endif
71 0 : temp->n_sign = PLUS;
72 0 : temp->n_len = length;
73 0 : temp->n_scale = scale;
74 0 : temp->n_refs = 1;
75 0 : temp->n_ptr = NULL;
76 0 : temp->n_value = value;
77 0 : return temp;
78 : }
79 :
80 : static void
81 : _bc_simp_mul (bc_num n1, int n1len, bc_num n2, int n2len, bc_num *prod,
82 : int full_scale)
83 38 : {
84 : char *n1ptr, *n2ptr, *pvptr;
85 : char *n1end, *n2end; /* To the end of n1 and n2. */
86 : int indx, sum, prodlen;
87 :
88 38 : prodlen = n1len+n2len+1;
89 :
90 38 : *prod = bc_new_num (prodlen, 0);
91 :
92 38 : n1end = (char *) (n1->n_value + n1len - 1);
93 38 : n2end = (char *) (n2->n_value + n2len - 1);
94 38 : pvptr = (char *) ((*prod)->n_value + prodlen - 1);
95 38 : sum = 0;
96 :
97 : /* Here is the loop... */
98 382 : for (indx = 0; indx < prodlen-1; indx++)
99 : {
100 344 : n1ptr = (char *) (n1end - MAX(0, indx-n2len+1));
101 344 : n2ptr = (char *) (n2end - MIN(indx, n2len-1));
102 1555 : while ((n1ptr >= n1->n_value) && (n2ptr <= n2end))
103 867 : sum += *n1ptr-- * *n2ptr++;
104 344 : *pvptr-- = sum % BASE;
105 344 : sum = sum / BASE;
106 : }
107 38 : *pvptr = sum;
108 38 : }
109 :
110 :
111 : /* A special adder/subtractor for the recursive divide and conquer
112 : multiply algorithm. Note: if sub is called, accum must
113 : be larger that what is being subtracted. Also, accum and val
114 : must have n_scale = 0. (e.g. they must look like integers. *) */
115 : static void
116 : _bc_shift_addsub (bc_num accum, bc_num val, int shift, int sub)
117 0 : {
118 : signed char *accp, *valp;
119 : int count, carry;
120 :
121 0 : count = val->n_len;
122 0 : if (val->n_value[0] == 0)
123 0 : count--;
124 : assert (accum->n_len+accum->n_scale >= shift+count);
125 :
126 : /* Set up pointers and others */
127 0 : accp = (signed char *)(accum->n_value +
128 : accum->n_len + accum->n_scale - shift - 1);
129 0 : valp = (signed char *)(val->n_value + val->n_len - 1);
130 0 : carry = 0;
131 :
132 0 : if (sub) {
133 : /* Subtraction, carry is really borrow. */
134 0 : while (count--) {
135 0 : *accp -= *valp-- + carry;
136 0 : if (*accp < 0) {
137 0 : carry = 1;
138 0 : *accp-- += BASE;
139 : } else {
140 0 : carry = 0;
141 0 : accp--;
142 : }
143 : }
144 0 : while (carry) {
145 0 : *accp -= carry;
146 0 : if (*accp < 0)
147 0 : *accp-- += BASE;
148 : else
149 0 : carry = 0;
150 : }
151 : } else {
152 : /* Addition */
153 0 : while (count--) {
154 0 : *accp += *valp-- + carry;
155 0 : if (*accp > (BASE-1)) {
156 0 : carry = 1;
157 0 : *accp-- -= BASE;
158 : } else {
159 0 : carry = 0;
160 0 : accp--;
161 : }
162 : }
163 0 : while (carry) {
164 0 : *accp += carry;
165 0 : if (*accp > (BASE-1))
166 0 : *accp-- -= BASE;
167 : else
168 0 : carry = 0;
169 : }
170 : }
171 0 : }
172 :
173 : /* Recursive divide and conquer multiply algorithm.
174 : Based on
175 : Let u = u0 + u1*(b^n)
176 : Let v = v0 + v1*(b^n)
177 : Then uv = (B^2n+B^n)*u1*v1 + B^n*(u1-u0)*(v0-v1) + (B^n+1)*u0*v0
178 :
179 : B is the base of storage, number of digits in u1,u0 close to equal.
180 : */
181 : static void
182 : _bc_rec_mul (bc_num u, int ulen, bc_num v, int vlen, bc_num *prod,
183 : int full_scale TSRMLS_DC)
184 38 : {
185 : bc_num u0, u1, v0, v1;
186 : int u0len, v0len;
187 : bc_num m1, m2, m3, d1, d2;
188 : int n, prodlen, m1zero;
189 : int d1len, d2len;
190 :
191 : /* Base case? */
192 38 : if ((ulen+vlen) < mul_base_digits
193 : || ulen < MUL_SMALL_DIGITS
194 : || vlen < MUL_SMALL_DIGITS ) {
195 38 : _bc_simp_mul (u, ulen, v, vlen, prod, full_scale);
196 38 : return;
197 : }
198 :
199 : /* Calculate n -- the u and v split point in digits. */
200 0 : n = (MAX(ulen, vlen)+1) / 2;
201 :
202 : /* Split u and v. */
203 0 : if (ulen < n) {
204 0 : u1 = bc_copy_num (BCG(_zero_));
205 0 : u0 = new_sub_num (ulen,0, u->n_value);
206 : } else {
207 0 : u1 = new_sub_num (ulen-n, 0, u->n_value);
208 0 : u0 = new_sub_num (n, 0, u->n_value+ulen-n);
209 : }
210 0 : if (vlen < n) {
211 0 : v1 = bc_copy_num (BCG(_zero_));
212 0 : v0 = new_sub_num (vlen,0, v->n_value);
213 : } else {
214 0 : v1 = new_sub_num (vlen-n, 0, v->n_value);
215 0 : v0 = new_sub_num (n, 0, v->n_value+vlen-n);
216 : }
217 0 : _bc_rm_leading_zeros (u1);
218 0 : _bc_rm_leading_zeros (u0);
219 0 : u0len = u0->n_len;
220 0 : _bc_rm_leading_zeros (v1);
221 0 : _bc_rm_leading_zeros (v0);
222 0 : v0len = v0->n_len;
223 :
224 0 : m1zero = bc_is_zero(u1 TSRMLS_CC) || bc_is_zero(v1 TSRMLS_CC);
225 :
226 : /* Calculate sub results ... */
227 :
228 0 : bc_init_num(&d1 TSRMLS_CC);
229 0 : bc_init_num(&d2 TSRMLS_CC);
230 0 : bc_sub (u1, u0, &d1, 0);
231 0 : d1len = d1->n_len;
232 0 : bc_sub (v0, v1, &d2, 0);
233 0 : d2len = d2->n_len;
234 :
235 :
236 : /* Do recursive multiplies and shifted adds. */
237 0 : if (m1zero)
238 0 : m1 = bc_copy_num (BCG(_zero_));
239 : else
240 0 : _bc_rec_mul (u1, u1->n_len, v1, v1->n_len, &m1, 0 TSRMLS_CC);
241 :
242 0 : if (bc_is_zero(d1 TSRMLS_CC) || bc_is_zero(d2 TSRMLS_CC))
243 0 : m2 = bc_copy_num (BCG(_zero_));
244 : else
245 0 : _bc_rec_mul (d1, d1len, d2, d2len, &m2, 0 TSRMLS_CC);
246 :
247 0 : if (bc_is_zero(u0 TSRMLS_CC) || bc_is_zero(v0 TSRMLS_CC))
248 0 : m3 = bc_copy_num (BCG(_zero_));
249 : else
250 0 : _bc_rec_mul (u0, u0->n_len, v0, v0->n_len, &m3, 0 TSRMLS_CC);
251 :
252 : /* Initialize product */
253 0 : prodlen = ulen+vlen+1;
254 0 : *prod = bc_new_num(prodlen, 0);
255 :
256 0 : if (!m1zero) {
257 0 : _bc_shift_addsub (*prod, m1, 2*n, 0);
258 0 : _bc_shift_addsub (*prod, m1, n, 0);
259 : }
260 0 : _bc_shift_addsub (*prod, m3, n, 0);
261 0 : _bc_shift_addsub (*prod, m3, 0, 0);
262 0 : _bc_shift_addsub (*prod, m2, n, d1->n_sign != d2->n_sign);
263 :
264 : /* Now clean up! */
265 0 : bc_free_num (&u1);
266 0 : bc_free_num (&u0);
267 0 : bc_free_num (&v1);
268 0 : bc_free_num (&m1);
269 0 : bc_free_num (&v0);
270 0 : bc_free_num (&m2);
271 0 : bc_free_num (&m3);
272 0 : bc_free_num (&d1);
273 0 : bc_free_num (&d2);
274 : }
275 :
276 : /* The multiply routine. N2 times N1 is put int PROD with the scale of
277 : the result being MIN(N2 scale+N1 scale, MAX (SCALE, N2 scale, N1 scale)).
278 : */
279 :
280 : void
281 : bc_multiply (bc_num n1, bc_num n2, bc_num *prod, int scale TSRMLS_DC)
282 38 : {
283 : bc_num pval;
284 : int len1, len2;
285 : int full_scale, prod_scale;
286 :
287 : /* Initialize things. */
288 38 : len1 = n1->n_len + n1->n_scale;
289 38 : len2 = n2->n_len + n2->n_scale;
290 38 : full_scale = n1->n_scale + n2->n_scale;
291 38 : prod_scale = MIN(full_scale,MAX(scale,MAX(n1->n_scale,n2->n_scale)));
292 :
293 : /* Do the multiply */
294 38 : _bc_rec_mul (n1, len1, n2, len2, &pval, full_scale TSRMLS_CC);
295 :
296 : /* Assign to prod and clean up the number. */
297 38 : pval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS );
298 38 : pval->n_value = pval->n_ptr;
299 38 : pval->n_len = len2 + len1 + 1 - full_scale;
300 38 : pval->n_scale = prod_scale;
301 38 : _bc_rm_leading_zeros (pval);
302 38 : if (bc_is_zero (pval TSRMLS_CC))
303 1 : pval->n_sign = PLUS;
304 38 : bc_free_num (prod);
305 38 : *prod = pval;
306 38 : }
|
