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. (LICENSE)
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
new_sub_num(int length,int scale,char * value)53 static bc_num new_sub_num(int length, int scale, char *value)
54 {
55 bc_num temp;
56
57 temp = (bc_num) emalloc (sizeof(bc_struct));
58
59 temp->n_sign = PLUS;
60 temp->n_len = length;
61 temp->n_scale = scale;
62 temp->n_refs = 1;
63 temp->n_ptr = NULL;
64 temp->n_value = value;
65 return temp;
66 }
67
68 static void
_bc_simp_mul(bc_num n1,int n1len,bc_num n2,int n2len,bc_num * prod,int full_scale)69 _bc_simp_mul (bc_num n1, int n1len, bc_num n2, int n2len, bc_num *prod,
70 int full_scale)
71 {
72 char *n1ptr, *n2ptr, *pvptr;
73 char *n1end, *n2end; /* To the end of n1 and n2. */
74 int indx, sum, prodlen;
75
76 prodlen = n1len+n2len+1;
77
78 *prod = bc_new_num (prodlen, 0);
79
80 n1end = (char *) (n1->n_value + n1len - 1);
81 n2end = (char *) (n2->n_value + n2len - 1);
82 pvptr = (char *) ((*prod)->n_value + prodlen - 1);
83 sum = 0;
84
85 /* Here is the loop... */
86 for (indx = 0; indx < prodlen-1; indx++)
87 {
88 n1ptr = (char *) (n1end - MAX(0, indx-n2len+1));
89 n2ptr = (char *) (n2end - MIN(indx, n2len-1));
90 while ((n1ptr >= n1->n_value) && (n2ptr <= n2end))
91 sum += *n1ptr-- * *n2ptr++;
92 *pvptr-- = sum % BASE;
93 sum = sum / BASE;
94 }
95 *pvptr = sum;
96 }
97
98
99 /* A special adder/subtractor for the recursive divide and conquer
100 multiply algorithm. Note: if sub is called, accum must
101 be larger that what is being subtracted. Also, accum and val
102 must have n_scale = 0. (e.g. they must look like integers. *) */
103 static void
_bc_shift_addsub(bc_num accum,bc_num val,int shift,int sub)104 _bc_shift_addsub (bc_num accum, bc_num val, int shift, int sub)
105 {
106 signed char *accp, *valp;
107 int count, carry;
108
109 count = val->n_len;
110 if (val->n_value[0] == 0)
111 count--;
112 assert (accum->n_len+accum->n_scale >= shift+count);
113
114 /* Set up pointers and others */
115 accp = (signed char *)(accum->n_value +
116 accum->n_len + accum->n_scale - shift - 1);
117 valp = (signed char *)(val->n_value + val->n_len - 1);
118 carry = 0;
119
120 if (sub) {
121 /* Subtraction, carry is really borrow. */
122 while (count--) {
123 *accp -= *valp-- + carry;
124 if (*accp < 0) {
125 carry = 1;
126 *accp-- += BASE;
127 } else {
128 carry = 0;
129 accp--;
130 }
131 }
132 while (carry) {
133 *accp -= carry;
134 if (*accp < 0)
135 *accp-- += BASE;
136 else
137 carry = 0;
138 }
139 } else {
140 /* Addition */
141 while (count--) {
142 *accp += *valp-- + carry;
143 if (*accp > (BASE-1)) {
144 carry = 1;
145 *accp-- -= BASE;
146 } else {
147 carry = 0;
148 accp--;
149 }
150 }
151 while (carry) {
152 *accp += carry;
153 if (*accp > (BASE-1))
154 *accp-- -= BASE;
155 else
156 carry = 0;
157 }
158 }
159 }
160
161 /* Recursive divide and conquer multiply algorithm.
162 Based on
163 Let u = u0 + u1*(b^n)
164 Let v = v0 + v1*(b^n)
165 Then uv = (B^2n+B^n)*u1*v1 + B^n*(u1-u0)*(v0-v1) + (B^n+1)*u0*v0
166
167 B is the base of storage, number of digits in u1,u0 close to equal.
168 */
169 static void
_bc_rec_mul(bc_num u,int ulen,bc_num v,int vlen,bc_num * prod,int full_scale)170 _bc_rec_mul (bc_num u, int ulen, bc_num v, int vlen, bc_num *prod,
171 int full_scale)
172 {
173 bc_num u0, u1, v0, v1;
174 bc_num m1, m2, m3, d1, d2;
175 int n, prodlen, m1zero;
176 int d1len, d2len;
177
178 /* Base case? */
179 if ((ulen+vlen) < mul_base_digits
180 || ulen < MUL_SMALL_DIGITS
181 || vlen < MUL_SMALL_DIGITS ) {
182 _bc_simp_mul (u, ulen, v, vlen, prod, full_scale);
183 return;
184 }
185
186 /* Calculate n -- the u and v split point in digits. */
187 n = (MAX(ulen, vlen)+1) / 2;
188
189 /* Split u and v. */
190 if (ulen < n) {
191 u1 = bc_copy_num (BCG(_zero_));
192 u0 = new_sub_num (ulen,0, u->n_value);
193 } else {
194 u1 = new_sub_num (ulen-n, 0, u->n_value);
195 u0 = new_sub_num (n, 0, u->n_value+ulen-n);
196 }
197 if (vlen < n) {
198 v1 = bc_copy_num (BCG(_zero_));
199 v0 = new_sub_num (vlen,0, v->n_value);
200 } else {
201 v1 = new_sub_num (vlen-n, 0, v->n_value);
202 v0 = new_sub_num (n, 0, v->n_value+vlen-n);
203 }
204 _bc_rm_leading_zeros (u1);
205 _bc_rm_leading_zeros (u0);
206 _bc_rm_leading_zeros (v1);
207 _bc_rm_leading_zeros (v0);
208
209 m1zero = bc_is_zero(u1) || bc_is_zero(v1);
210
211 /* Calculate sub results ... */
212
213 bc_init_num(&d1);
214 bc_init_num(&d2);
215 bc_sub (u1, u0, &d1, 0);
216 d1len = d1->n_len;
217 bc_sub (v0, v1, &d2, 0);
218 d2len = d2->n_len;
219
220
221 /* Do recursive multiplies and shifted adds. */
222 if (m1zero)
223 m1 = bc_copy_num (BCG(_zero_));
224 else
225 _bc_rec_mul (u1, u1->n_len, v1, v1->n_len, &m1, 0);
226
227 if (bc_is_zero(d1) || bc_is_zero(d2))
228 m2 = bc_copy_num (BCG(_zero_));
229 else
230 _bc_rec_mul (d1, d1len, d2, d2len, &m2, 0);
231
232 if (bc_is_zero(u0) || bc_is_zero(v0))
233 m3 = bc_copy_num (BCG(_zero_));
234 else
235 _bc_rec_mul (u0, u0->n_len, v0, v0->n_len, &m3, 0);
236
237 /* Initialize product */
238 prodlen = ulen+vlen+1;
239 *prod = bc_new_num(prodlen, 0);
240
241 if (!m1zero) {
242 _bc_shift_addsub (*prod, m1, 2*n, 0);
243 _bc_shift_addsub (*prod, m1, n, 0);
244 }
245 _bc_shift_addsub (*prod, m3, n, 0);
246 _bc_shift_addsub (*prod, m3, 0, 0);
247 _bc_shift_addsub (*prod, m2, n, d1->n_sign != d2->n_sign);
248
249 /* Now clean up! */
250 bc_free_num (&u1);
251 bc_free_num (&u0);
252 bc_free_num (&v1);
253 bc_free_num (&m1);
254 bc_free_num (&v0);
255 bc_free_num (&m2);
256 bc_free_num (&m3);
257 bc_free_num (&d1);
258 bc_free_num (&d2);
259 }
260
261 /* The multiply routine. N2 times N1 is put int PROD with the scale of
262 the result being MIN(N2 scale+N1 scale, MAX (SCALE, N2 scale, N1 scale)).
263 */
264
265 void
bc_multiply(bc_num n1,bc_num n2,bc_num * prod,int scale)266 bc_multiply (bc_num n1, bc_num n2, bc_num *prod, int scale)
267 {
268 bc_num pval;
269 int len1, len2;
270 int full_scale, prod_scale;
271
272 /* Initialize things. */
273 len1 = n1->n_len + n1->n_scale;
274 len2 = n2->n_len + n2->n_scale;
275 full_scale = n1->n_scale + n2->n_scale;
276 prod_scale = MIN(full_scale,MAX(scale,MAX(n1->n_scale,n2->n_scale)));
277
278 /* Do the multiply */
279 _bc_rec_mul (n1, len1, n2, len2, &pval, full_scale);
280
281 /* Assign to prod and clean up the number. */
282 pval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS );
283 pval->n_value = pval->n_ptr;
284 pval->n_len = len2 + len1 + 1 - full_scale;
285 pval->n_scale = prod_scale;
286 _bc_rm_leading_zeros (pval);
287 if (bc_is_zero (pval))
288 pval->n_sign = PLUS;
289 bc_free_num (prod);
290 *prod = pval;
291 }
292
