Annotation of embedaddon/php/ext/bcmath/libbcmath/src/recmul.c, revision 1.1
1.1 ! misho 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 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
! 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: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to recmul.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [ELWIX - Embedded LightWeight unIX -] / embedaddon / php / ext / bcmath / libbcmath / src