root/ext/bcmath/libbcmath/src/recmul.c

DEFINITIONS

1. new_sub_num
2. _bc_simp_mul
3. _bc_shift_addsub
4. _bc_rec_mul
5. bc_multiply

   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 {
  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
  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
 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
 182 _bc_rec_mul (bc_num u, int ulen, bc_num v, int vlen, bc_num *prod,
 183              int full_scale)
 184 {
 185   bc_num u0, u1, v0, v1;
 186   bc_num m1, m2, m3, d1, d2;
 187   int n, prodlen, m1zero;
 188   int d1len, d2len;
 189 
 190   /* Base case? */
 191   if ((ulen+vlen) < mul_base_digits
 192       || ulen < MUL_SMALL_DIGITS
 193       || vlen < MUL_SMALL_DIGITS ) {
 194     _bc_simp_mul (u, ulen, v, vlen, prod, full_scale);
 195     return;
 196   }
 197 
 198   /* Calculate n -- the u and v split point in digits. */
 199   n = (MAX(ulen, vlen)+1) / 2;
 200 
 201   /* Split u and v. */
 202   if (ulen < n) {
 203     u1 = bc_copy_num (BCG(_zero_));
 204     u0 = new_sub_num (ulen,0, u->n_value);
 205   } else {
 206     u1 = new_sub_num (ulen-n, 0, u->n_value);
 207     u0 = new_sub_num (n, 0, u->n_value+ulen-n);
 208   }
 209   if (vlen < n) {
 210     v1 = bc_copy_num (BCG(_zero_));
 211     v0 = new_sub_num (vlen,0, v->n_value);
 212   } else {
 213     v1 = new_sub_num (vlen-n, 0, v->n_value);
 214     v0 = new_sub_num (n, 0, v->n_value+vlen-n);
 215     }
 216   _bc_rm_leading_zeros (u1);
 217   _bc_rm_leading_zeros (u0);
 218   _bc_rm_leading_zeros (v1);
 219   _bc_rm_leading_zeros (v0);
 220 
 221   m1zero = bc_is_zero(u1) || bc_is_zero(v1);
 222 
 223   /* Calculate sub results ... */
 224 
 225   bc_init_num(&d1);
 226   bc_init_num(&d2);
 227   bc_sub (u1, u0, &d1, 0);
 228   d1len = d1->n_len;
 229   bc_sub (v0, v1, &d2, 0);
 230   d2len = d2->n_len;
 231 
 232 
 233   /* Do recursive multiplies and shifted adds. */
 234   if (m1zero)
 235     m1 = bc_copy_num (BCG(_zero_));
 236   else
 237     _bc_rec_mul (u1, u1->n_len, v1, v1->n_len, &m1, 0);
 238 
 239   if (bc_is_zero(d1) || bc_is_zero(d2))
 240     m2 = bc_copy_num (BCG(_zero_));
 241   else
 242     _bc_rec_mul (d1, d1len, d2, d2len, &m2, 0);
 243 
 244   if (bc_is_zero(u0) || bc_is_zero(v0))
 245     m3 = bc_copy_num (BCG(_zero_));
 246   else
 247     _bc_rec_mul (u0, u0->n_len, v0, v0->n_len, &m3, 0);
 248 
 249   /* Initialize product */
 250   prodlen = ulen+vlen+1;
 251   *prod = bc_new_num(prodlen, 0);
 252 
 253   if (!m1zero) {
 254     _bc_shift_addsub (*prod, m1, 2*n, 0);
 255     _bc_shift_addsub (*prod, m1, n, 0);
 256   }
 257   _bc_shift_addsub (*prod, m3, n, 0);
 258   _bc_shift_addsub (*prod, m3, 0, 0);
 259   _bc_shift_addsub (*prod, m2, n, d1->n_sign != d2->n_sign);
 260 
 261   /* Now clean up! */
 262   bc_free_num (&u1);
 263   bc_free_num (&u0);
 264   bc_free_num (&v1);
 265   bc_free_num (&m1);
 266   bc_free_num (&v0);
 267   bc_free_num (&m2);
 268   bc_free_num (&m3);
 269   bc_free_num (&d1);
 270   bc_free_num (&d2);
 271 }
 272 
 273 /* The multiply routine.  N2 times N1 is put int PROD with the scale of
 274    the result being MIN(N2 scale+N1 scale, MAX (SCALE, N2 scale, N1 scale)).
 275    */
 276 
 277 void
 278 bc_multiply (bc_num n1, bc_num n2, bc_num *prod, int scale)
 279 {
 280   bc_num pval;
 281   int len1, len2;
 282   int full_scale, prod_scale;
 283 
 284   /* Initialize things. */
 285   len1 = n1->n_len + n1->n_scale;
 286   len2 = n2->n_len + n2->n_scale;
 287   full_scale = n1->n_scale + n2->n_scale;
 288   prod_scale = MIN(full_scale,MAX(scale,MAX(n1->n_scale,n2->n_scale)));
 289 
 290   /* Do the multiply */
 291   _bc_rec_mul (n1, len1, n2, len2, &pval, full_scale);
 292 
 293   /* Assign to prod and clean up the number. */
 294   pval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS );
 295   pval->n_value = pval->n_ptr;
 296   pval->n_len = len2 + len1 + 1 - full_scale;
 297   pval->n_scale = prod_scale;
 298   _bc_rm_leading_zeros (pval);
 299   if (bc_is_zero (pval))
 300     pval->n_sign = PLUS;
 301   bc_free_num (prod);
 302   *prod = pval;
 303 }