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>