xref: /PHP-5.4/ext/bcmath/libbcmath/src/recmul.c (revision 348b429b)
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
new_sub_num(length,scale,value)54 new_sub_num (length, scale, value)
55      int length, scale;
56      char *value;
57 {
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     temp = (bc_num) emalloc (sizeof(bc_struct));
67 #ifdef SANDER_0
68     if (temp == NULL) bc_out_of_memory ();
69   }
70 #endif
71   temp->n_sign = PLUS;
72   temp->n_len = length;
73   temp->n_scale = scale;
74   temp->n_refs = 1;
75   temp->n_ptr = NULL;
76   temp->n_value = value;
77   return temp;
78 }
79 
80 static void
_bc_simp_mul(bc_num n1,int n1len,bc_num n2,int n2len,bc_num * prod,int full_scale)81 _bc_simp_mul (bc_num n1, int n1len, bc_num n2, int n2len, bc_num *prod,
82 	      int full_scale)
83 {
84   char *n1ptr, *n2ptr, *pvptr;
85   char *n1end, *n2end;		/* To the end of n1 and n2. */
86   int indx, sum, prodlen;
87 
88   prodlen = n1len+n2len+1;
89 
90   *prod = bc_new_num (prodlen, 0);
91 
92   n1end = (char *) (n1->n_value + n1len - 1);
93   n2end = (char *) (n2->n_value + n2len - 1);
94   pvptr = (char *) ((*prod)->n_value + prodlen - 1);
95   sum = 0;
96 
97   /* Here is the loop... */
98   for (indx = 0; indx < prodlen-1; indx++)
99     {
100       n1ptr = (char *) (n1end - MAX(0, indx-n2len+1));
101       n2ptr = (char *) (n2end - MIN(indx, n2len-1));
102       while ((n1ptr >= n1->n_value) && (n2ptr <= n2end))
103 	sum += *n1ptr-- * *n2ptr++;
104       *pvptr-- = sum % BASE;
105       sum = sum / BASE;
106     }
107   *pvptr = sum;
108 }
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
_bc_shift_addsub(bc_num accum,bc_num val,int shift,int sub)116 _bc_shift_addsub (bc_num accum, bc_num val, int shift, int sub)
117 {
118   signed char *accp, *valp;
119   int  count, carry;
120 
121   count = val->n_len;
122   if (val->n_value[0] == 0)
123     count--;
124   assert (accum->n_len+accum->n_scale >= shift+count);
125 
126   /* Set up pointers and others */
127   accp = (signed char *)(accum->n_value +
128 			 accum->n_len + accum->n_scale - shift - 1);
129   valp = (signed char *)(val->n_value + val->n_len - 1);
130   carry = 0;
131 
132   if (sub) {
133     /* Subtraction, carry is really borrow. */
134     while (count--) {
135       *accp -= *valp-- + carry;
136       if (*accp < 0) {
137 	carry = 1;
138         *accp-- += BASE;
139       } else {
140 	carry = 0;
141 	accp--;
142       }
143     }
144     while (carry) {
145       *accp -= carry;
146       if (*accp < 0)
147 	*accp-- += BASE;
148       else
149 	carry = 0;
150     }
151   } else {
152     /* Addition */
153     while (count--) {
154       *accp += *valp-- + carry;
155       if (*accp > (BASE-1)) {
156 	carry = 1;
157         *accp-- -= BASE;
158       } else {
159 	carry = 0;
160 	accp--;
161       }
162     }
163     while (carry) {
164       *accp += carry;
165       if (*accp > (BASE-1))
166 	*accp-- -= BASE;
167       else
168 	carry = 0;
169     }
170   }
171 }
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
_bc_rec_mul(bc_num u,int ulen,bc_num v,int vlen,bc_num * prod,int full_scale TSRMLS_DC)182 _bc_rec_mul (bc_num u, int ulen, bc_num v, int vlen, bc_num *prod,
183 	     int full_scale TSRMLS_DC)
184 {
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   if ((ulen+vlen) < mul_base_digits
193       || ulen < MUL_SMALL_DIGITS
194       || vlen < MUL_SMALL_DIGITS ) {
195     _bc_simp_mul (u, ulen, v, vlen, prod, full_scale);
196     return;
197   }
198 
199   /* Calculate n -- the u and v split point in digits. */
200   n = (MAX(ulen, vlen)+1) / 2;
201 
202   /* Split u and v. */
203   if (ulen < n) {
204     u1 = bc_copy_num (BCG(_zero_));
205     u0 = new_sub_num (ulen,0, u->n_value);
206   } else {
207     u1 = new_sub_num (ulen-n, 0, u->n_value);
208     u0 = new_sub_num (n, 0, u->n_value+ulen-n);
209   }
210   if (vlen < n) {
211     v1 = bc_copy_num (BCG(_zero_));
212     v0 = new_sub_num (vlen,0, v->n_value);
213   } else {
214     v1 = new_sub_num (vlen-n, 0, v->n_value);
215     v0 = new_sub_num (n, 0, v->n_value+vlen-n);
216     }
217   _bc_rm_leading_zeros (u1);
218   _bc_rm_leading_zeros (u0);
219   u0len = u0->n_len;
220   _bc_rm_leading_zeros (v1);
221   _bc_rm_leading_zeros (v0);
222   v0len = v0->n_len;
223 
224   m1zero = bc_is_zero(u1 TSRMLS_CC) || bc_is_zero(v1 TSRMLS_CC);
225 
226   /* Calculate sub results ... */
227 
228   bc_init_num(&d1 TSRMLS_CC);
229   bc_init_num(&d2 TSRMLS_CC);
230   bc_sub (u1, u0, &d1, 0);
231   d1len = d1->n_len;
232   bc_sub (v0, v1, &d2, 0);
233   d2len = d2->n_len;
234 
235 
236   /* Do recursive multiplies and shifted adds. */
237   if (m1zero)
238     m1 = bc_copy_num (BCG(_zero_));
239   else
240     _bc_rec_mul (u1, u1->n_len, v1, v1->n_len, &m1, 0 TSRMLS_CC);
241 
242   if (bc_is_zero(d1 TSRMLS_CC) || bc_is_zero(d2 TSRMLS_CC))
243     m2 = bc_copy_num (BCG(_zero_));
244   else
245     _bc_rec_mul (d1, d1len, d2, d2len, &m2, 0 TSRMLS_CC);
246 
247   if (bc_is_zero(u0 TSRMLS_CC) || bc_is_zero(v0 TSRMLS_CC))
248     m3 = bc_copy_num (BCG(_zero_));
249   else
250     _bc_rec_mul (u0, u0->n_len, v0, v0->n_len, &m3, 0 TSRMLS_CC);
251 
252   /* Initialize product */
253   prodlen = ulen+vlen+1;
254   *prod = bc_new_num(prodlen, 0);
255 
256   if (!m1zero) {
257     _bc_shift_addsub (*prod, m1, 2*n, 0);
258     _bc_shift_addsub (*prod, m1, n, 0);
259   }
260   _bc_shift_addsub (*prod, m3, n, 0);
261   _bc_shift_addsub (*prod, m3, 0, 0);
262   _bc_shift_addsub (*prod, m2, n, d1->n_sign != d2->n_sign);
263 
264   /* Now clean up! */
265   bc_free_num (&u1);
266   bc_free_num (&u0);
267   bc_free_num (&v1);
268   bc_free_num (&m1);
269   bc_free_num (&v0);
270   bc_free_num (&m2);
271   bc_free_num (&m3);
272   bc_free_num (&d1);
273   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
bc_multiply(bc_num n1,bc_num n2,bc_num * prod,int scale TSRMLS_DC)281 bc_multiply (bc_num n1, bc_num n2, bc_num *prod, int scale TSRMLS_DC)
282 {
283   bc_num pval;
284   int len1, len2;
285   int full_scale, prod_scale;
286 
287   /* Initialize things. */
288   len1 = n1->n_len + n1->n_scale;
289   len2 = n2->n_len + n2->n_scale;
290   full_scale = n1->n_scale + n2->n_scale;
291   prod_scale = MIN(full_scale,MAX(scale,MAX(n1->n_scale,n2->n_scale)));
292 
293   /* Do the multiply */
294   _bc_rec_mul (n1, len1, n2, len2, &pval, full_scale TSRMLS_CC);
295 
296   /* Assign to prod and clean up the number. */
297   pval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS );
298   pval->n_value = pval->n_ptr;
299   pval->n_len = len2 + len1 + 1 - full_scale;
300   pval->n_scale = prod_scale;
301   _bc_rm_leading_zeros (pval);
302   if (bc_is_zero (pval TSRMLS_CC))
303     pval->n_sign = PLUS;
304   bc_free_num (prod);
305   *prod = pval;
306 }
307