Back to index

php5  5.3.10
recmul.c
Go to the documentation of this file.
00001 /* recmul.c: bcmath library file. */
00002 /*
00003     Copyright (C) 1991, 1992, 1993, 1994, 1997 Free Software Foundation, Inc.
00004     Copyright (C) 2000 Philip A. Nelson
00005 
00006     This library is free software; you can redistribute it and/or
00007     modify it under the terms of the GNU Lesser General Public
00008     License as published by the Free Software Foundation; either
00009     version 2 of the License, or (at your option) any later version.
00010 
00011     This library is distributed in the hope that it will be useful,
00012     but WITHOUT ANY WARRANTY; without even the implied warranty of
00013     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014     Lesser General Public License for more details.  (COPYING.LIB)
00015 
00016     You should have received a copy of the GNU Lesser General Public
00017     License along with this library; if not, write to:
00018 
00019       The Free Software Foundation, Inc.
00020       59 Temple Place, Suite 330
00021       Boston, MA 02111-1307 USA.
00022 
00023     You may contact the author by:
00024        e-mail:  philnelson@acm.org
00025       us-mail:  Philip A. Nelson
00026                 Computer Science Department, 9062
00027                 Western Washington University
00028                 Bellingham, WA 98226-9062
00029        
00030 *************************************************************************/
00031 
00032 #include <config.h>
00033 #include <stdio.h>
00034 #include <assert.h>
00035 #include <stdlib.h>
00036 #include <ctype.h>
00037 #include <stdarg.h>
00038 #include "bcmath.h"
00039 #include "private.h"
00040 
00041 /* Recursive vs non-recursive multiply crossover ranges. */
00042 #if defined(MULDIGITS)
00043 #include "muldigits.h"
00044 #else
00045 #define MUL_BASE_DIGITS 80
00046 #endif
00047 
00048 int mul_base_digits = MUL_BASE_DIGITS;
00049 #define MUL_SMALL_DIGITS mul_base_digits/4
00050 
00051 /* Multiply utility routines */
00052 
00053 static bc_num
00054 new_sub_num (length, scale, value)
00055      int length, scale;
00056      char *value;
00057 {
00058   bc_num temp;
00059 
00060 #ifdef SANDER_0
00061   if (_bc_Free_list != NULL) {
00062     temp = _bc_Free_list;
00063     _bc_Free_list = temp->n_next;
00064   } else {
00065 #endif
00066     temp = (bc_num) emalloc (sizeof(bc_struct));
00067 #ifdef SANDER_0
00068     if (temp == NULL) bc_out_of_memory ();
00069   }
00070 #endif
00071   temp->n_sign = PLUS;
00072   temp->n_len = length;
00073   temp->n_scale = scale;
00074   temp->n_refs = 1;
00075   temp->n_ptr = NULL;
00076   temp->n_value = value;
00077   return temp;
00078 }
00079 
00080 static void
00081 _bc_simp_mul (bc_num n1, int n1len, bc_num n2, int n2len, bc_num *prod,
00082              int full_scale)
00083 {
00084   char *n1ptr, *n2ptr, *pvptr;
00085   char *n1end, *n2end;             /* To the end of n1 and n2. */
00086   int indx, sum, prodlen;
00087 
00088   prodlen = n1len+n2len+1;
00089 
00090   *prod = bc_new_num (prodlen, 0);
00091 
00092   n1end = (char *) (n1->n_value + n1len - 1);
00093   n2end = (char *) (n2->n_value + n2len - 1);
00094   pvptr = (char *) ((*prod)->n_value + prodlen - 1);
00095   sum = 0;
00096 
00097   /* Here is the loop... */
00098   for (indx = 0; indx < prodlen-1; indx++)
00099     {
00100       n1ptr = (char *) (n1end - MAX(0, indx-n2len+1));
00101       n2ptr = (char *) (n2end - MIN(indx, n2len-1));
00102       while ((n1ptr >= n1->n_value) && (n2ptr <= n2end))
00103        sum += *n1ptr-- * *n2ptr++;
00104       *pvptr-- = sum % BASE;
00105       sum = sum / BASE;
00106     }
00107   *pvptr = sum;
00108 }
00109 
00110 
00111 /* A special adder/subtractor for the recursive divide and conquer
00112    multiply algorithm.  Note: if sub is called, accum must
00113    be larger that what is being subtracted.  Also, accum and val
00114    must have n_scale = 0.  (e.g. they must look like integers. *) */
00115 static void
00116 _bc_shift_addsub (bc_num accum, bc_num val, int shift, int sub)
00117 {
00118   signed char *accp, *valp;
00119   int  count, carry;
00120 
00121   count = val->n_len;
00122   if (val->n_value[0] == 0)
00123     count--;
00124   assert (accum->n_len+accum->n_scale >= shift+count);
00125   
00126   /* Set up pointers and others */
00127   accp = (signed char *)(accum->n_value +
00128                       accum->n_len + accum->n_scale - shift - 1);
00129   valp = (signed char *)(val->n_value + val->n_len - 1);
00130   carry = 0;
00131 
00132   if (sub) {
00133     /* Subtraction, carry is really borrow. */
00134     while (count--) {
00135       *accp -= *valp-- + carry;
00136       if (*accp < 0) {
00137        carry = 1;
00138         *accp-- += BASE;
00139       } else {
00140        carry = 0;
00141        accp--;
00142       }
00143     }
00144     while (carry) {
00145       *accp -= carry;
00146       if (*accp < 0)
00147        *accp-- += BASE;
00148       else
00149        carry = 0;
00150     }
00151   } else {
00152     /* Addition */
00153     while (count--) {
00154       *accp += *valp-- + carry;
00155       if (*accp > (BASE-1)) {
00156        carry = 1;
00157         *accp-- -= BASE;
00158       } else {
00159        carry = 0;
00160        accp--;
00161       }
00162     }
00163     while (carry) {
00164       *accp += carry;
00165       if (*accp > (BASE-1))
00166        *accp-- -= BASE;
00167       else
00168        carry = 0;
00169     }
00170   }
00171 }
00172 
00173 /* Recursive divide and conquer multiply algorithm.  
00174    Based on 
00175    Let u = u0 + u1*(b^n)
00176    Let v = v0 + v1*(b^n)
00177    Then uv = (B^2n+B^n)*u1*v1 + B^n*(u1-u0)*(v0-v1) + (B^n+1)*u0*v0
00178 
00179    B is the base of storage, number of digits in u1,u0 close to equal.
00180 */
00181 static void
00182 _bc_rec_mul (bc_num u, int ulen, bc_num v, int vlen, bc_num *prod,
00183             int full_scale TSRMLS_DC)
00184 { 
00185   bc_num u0, u1, v0, v1;
00186   int u0len, v0len;
00187   bc_num m1, m2, m3, d1, d2;
00188   int n, prodlen, m1zero;
00189   int d1len, d2len;
00190 
00191   /* Base case? */
00192   if ((ulen+vlen) < mul_base_digits
00193       || ulen < MUL_SMALL_DIGITS
00194       || vlen < MUL_SMALL_DIGITS ) {
00195     _bc_simp_mul (u, ulen, v, vlen, prod, full_scale);
00196     return;
00197   }
00198 
00199   /* Calculate n -- the u and v split point in digits. */
00200   n = (MAX(ulen, vlen)+1) / 2;
00201 
00202   /* Split u and v. */
00203   if (ulen < n) {
00204     u1 = bc_copy_num (BCG(_zero_));
00205     u0 = new_sub_num (ulen,0, u->n_value);
00206   } else {
00207     u1 = new_sub_num (ulen-n, 0, u->n_value);
00208     u0 = new_sub_num (n, 0, u->n_value+ulen-n);
00209   }
00210   if (vlen < n) {
00211     v1 = bc_copy_num (BCG(_zero_));
00212     v0 = new_sub_num (vlen,0, v->n_value);
00213   } else {
00214     v1 = new_sub_num (vlen-n, 0, v->n_value);
00215     v0 = new_sub_num (n, 0, v->n_value+vlen-n);
00216     }
00217   _bc_rm_leading_zeros (u1);
00218   _bc_rm_leading_zeros (u0);
00219   u0len = u0->n_len;
00220   _bc_rm_leading_zeros (v1);
00221   _bc_rm_leading_zeros (v0);
00222   v0len = v0->n_len;
00223 
00224   m1zero = bc_is_zero(u1 TSRMLS_CC) || bc_is_zero(v1 TSRMLS_CC);
00225 
00226   /* Calculate sub results ... */
00227 
00228   bc_init_num(&d1 TSRMLS_CC);
00229   bc_init_num(&d2 TSRMLS_CC);
00230   bc_sub (u1, u0, &d1, 0);
00231   d1len = d1->n_len;
00232   bc_sub (v0, v1, &d2, 0);
00233   d2len = d2->n_len;
00234 
00235 
00236   /* Do recursive multiplies and shifted adds. */
00237   if (m1zero)
00238     m1 = bc_copy_num (BCG(_zero_));
00239   else
00240     _bc_rec_mul (u1, u1->n_len, v1, v1->n_len, &m1, 0 TSRMLS_CC);
00241 
00242   if (bc_is_zero(d1 TSRMLS_CC) || bc_is_zero(d2 TSRMLS_CC))
00243     m2 = bc_copy_num (BCG(_zero_));
00244   else
00245     _bc_rec_mul (d1, d1len, d2, d2len, &m2, 0 TSRMLS_CC);
00246 
00247   if (bc_is_zero(u0 TSRMLS_CC) || bc_is_zero(v0 TSRMLS_CC))
00248     m3 = bc_copy_num (BCG(_zero_));
00249   else
00250     _bc_rec_mul (u0, u0->n_len, v0, v0->n_len, &m3, 0 TSRMLS_CC);
00251 
00252   /* Initialize product */
00253   prodlen = ulen+vlen+1;
00254   *prod = bc_new_num(prodlen, 0);
00255 
00256   if (!m1zero) {
00257     _bc_shift_addsub (*prod, m1, 2*n, 0);
00258     _bc_shift_addsub (*prod, m1, n, 0);
00259   }
00260   _bc_shift_addsub (*prod, m3, n, 0);
00261   _bc_shift_addsub (*prod, m3, 0, 0);
00262   _bc_shift_addsub (*prod, m2, n, d1->n_sign != d2->n_sign);
00263 
00264   /* Now clean up! */
00265   bc_free_num (&u1);
00266   bc_free_num (&u0);
00267   bc_free_num (&v1);
00268   bc_free_num (&m1);
00269   bc_free_num (&v0);
00270   bc_free_num (&m2);
00271   bc_free_num (&m3);
00272   bc_free_num (&d1);
00273   bc_free_num (&d2);
00274 }
00275 
00276 /* The multiply routine.  N2 times N1 is put int PROD with the scale of
00277    the result being MIN(N2 scale+N1 scale, MAX (SCALE, N2 scale, N1 scale)).
00278    */
00279 
00280 void
00281 bc_multiply (bc_num n1, bc_num n2, bc_num *prod, int scale TSRMLS_DC)
00282 {
00283   bc_num pval; 
00284   int len1, len2;
00285   int full_scale, prod_scale;
00286 
00287   /* Initialize things. */
00288   len1 = n1->n_len + n1->n_scale;
00289   len2 = n2->n_len + n2->n_scale;
00290   full_scale = n1->n_scale + n2->n_scale;
00291   prod_scale = MIN(full_scale,MAX(scale,MAX(n1->n_scale,n2->n_scale)));
00292 
00293   /* Do the multiply */
00294   _bc_rec_mul (n1, len1, n2, len2, &pval, full_scale TSRMLS_CC);
00295 
00296   /* Assign to prod and clean up the number. */
00297   pval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS );
00298   pval->n_value = pval->n_ptr;
00299   pval->n_len = len2 + len1 + 1 - full_scale;
00300   pval->n_scale = prod_scale;
00301   _bc_rm_leading_zeros (pval);
00302   if (bc_is_zero (pval TSRMLS_CC))
00303     pval->n_sign = PLUS;
00304   bc_free_num (prod);
00305   *prod = pval;
00306 }