Return to div.c CVS log | Up to [ELWIX - Embedded LightWeight unIX -] / embedaddon / php / ext / bcmath / libbcmath / src |
1.1 misho 1: /* div.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: 42: /* Some utility routines for the divide: First a one digit multiply. 43: NUM (with SIZE digits) is multiplied by DIGIT and the result is 44: placed into RESULT. It is written so that NUM and RESULT can be 45: the same pointers. */ 46: 47: static void 48: _one_mult (num, size, digit, result) 49: unsigned char *num; 50: int size, digit; 51: unsigned char *result; 52: { 53: int carry, value; 54: unsigned char *nptr, *rptr; 55: 56: if (digit == 0) 57: memset (result, 0, size); 58: else 59: { 60: if (digit == 1) 61: memcpy (result, num, size); 62: else 63: { 64: /* Initialize */ 65: nptr = (unsigned char *) (num+size-1); 66: rptr = (unsigned char *) (result+size-1); 67: carry = 0; 68: 69: while (size-- > 0) 70: { 71: value = *nptr-- * digit + carry; 72: *rptr-- = value % BASE; 73: carry = value / BASE; 74: } 75: 76: if (carry != 0) *rptr = carry; 77: } 78: } 79: } 80: 81: 82: /* The full division routine. This computes N1 / N2. It returns 83: 0 if the division is ok and the result is in QUOT. The number of 84: digits after the decimal point is SCALE. It returns -1 if division 85: by zero is tried. The algorithm is found in Knuth Vol 2. p237. */ 86: 87: int 88: bc_divide (bc_num n1, bc_num n2, bc_num *quot, int scale TSRMLS_DC) 89: { 90: bc_num qval; 91: unsigned char *num1, *num2; 92: unsigned char *ptr1, *ptr2, *n2ptr, *qptr; 93: int scale1, val; 94: unsigned int len1, len2, scale2, qdigits, extra, count; 95: unsigned int qdig, qguess, borrow, carry; 96: unsigned char *mval; 97: char zero; 98: unsigned int norm; 99: 100: /* Test for divide by zero. */ 101: if (bc_is_zero (n2 TSRMLS_CC)) return -1; 102: 103: /* Test for divide by 1. If it is we must truncate. */ 104: if (n2->n_scale == 0) 105: { 106: if (n2->n_len == 1 && *n2->n_value == 1) 107: { 108: qval = bc_new_num (n1->n_len, scale); 109: qval->n_sign = (n1->n_sign == n2->n_sign ? PLUS : MINUS); 110: memset (&qval->n_value[n1->n_len],0,scale); 111: memcpy (qval->n_value, n1->n_value, 112: n1->n_len + MIN(n1->n_scale,scale)); 113: bc_free_num (quot); 114: *quot = qval; 115: } 116: } 117: 118: /* Set up the divide. Move the decimal point on n1 by n2's scale. 119: Remember, zeros on the end of num2 are wasted effort for dividing. */ 120: scale2 = n2->n_scale; 121: n2ptr = (unsigned char *) n2->n_value+n2->n_len+scale2-1; 122: while ((scale2 > 0) && (*n2ptr-- == 0)) scale2--; 123: 124: len1 = n1->n_len + scale2; 125: scale1 = n1->n_scale - scale2; 126: if (scale1 < scale) 127: extra = scale - scale1; 128: else 129: extra = 0; 130: num1 = (unsigned char *) safe_emalloc (1, n1->n_len+n1->n_scale, extra+2); 131: if (num1 == NULL) bc_out_of_memory(); 132: memset (num1, 0, n1->n_len+n1->n_scale+extra+2); 133: memcpy (num1+1, n1->n_value, n1->n_len+n1->n_scale); 134: 135: len2 = n2->n_len + scale2; 136: num2 = (unsigned char *) safe_emalloc (1, len2, 1); 137: if (num2 == NULL) bc_out_of_memory(); 138: memcpy (num2, n2->n_value, len2); 139: *(num2+len2) = 0; 140: n2ptr = num2; 141: while (*n2ptr == 0) 142: { 143: n2ptr++; 144: len2--; 145: } 146: 147: /* Calculate the number of quotient digits. */ 148: if (len2 > len1+scale) 149: { 150: qdigits = scale+1; 151: zero = TRUE; 152: } 153: else 154: { 155: zero = FALSE; 156: if (len2>len1) 157: qdigits = scale+1; /* One for the zero integer part. */ 158: else 159: qdigits = len1-len2+scale+1; 160: } 161: 162: /* Allocate and zero the storage for the quotient. */ 163: qval = bc_new_num (qdigits-scale,scale); 164: memset (qval->n_value, 0, qdigits); 165: 166: /* Allocate storage for the temporary storage mval. */ 167: mval = (unsigned char *) safe_emalloc (1, len2, 1); 168: if (mval == NULL) bc_out_of_memory (); 169: 170: /* Now for the full divide algorithm. */ 171: if (!zero) 172: { 173: /* Normalize */ 174: norm = 10 / ((int)*n2ptr + 1); 175: if (norm != 1) 176: { 177: _one_mult (num1, len1+scale1+extra+1, norm, num1); 178: _one_mult (n2ptr, len2, norm, n2ptr); 179: } 180: 181: /* Initialize divide loop. */ 182: qdig = 0; 183: if (len2 > len1) 184: qptr = (unsigned char *) qval->n_value+len2-len1; 185: else 186: qptr = (unsigned char *) qval->n_value; 187: 188: /* Loop */ 189: while (qdig <= len1+scale-len2) 190: { 191: /* Calculate the quotient digit guess. */ 192: if (*n2ptr == num1[qdig]) 193: qguess = 9; 194: else 195: qguess = (num1[qdig]*10 + num1[qdig+1]) / *n2ptr; 196: 197: /* Test qguess. */ 198: if (n2ptr[1]*qguess > 199: (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10 200: + num1[qdig+2]) 201: { 202: qguess--; 203: /* And again. */ 204: if (n2ptr[1]*qguess > 205: (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10 206: + num1[qdig+2]) 207: qguess--; 208: } 209: 210: /* Multiply and subtract. */ 211: borrow = 0; 212: if (qguess != 0) 213: { 214: *mval = 0; 215: _one_mult (n2ptr, len2, qguess, mval+1); 216: ptr1 = (unsigned char *) num1+qdig+len2; 217: ptr2 = (unsigned char *) mval+len2; 218: for (count = 0; count < len2+1; count++) 219: { 220: val = (int) *ptr1 - (int) *ptr2-- - borrow; 221: if (val < 0) 222: { 223: val += 10; 224: borrow = 1; 225: } 226: else 227: borrow = 0; 228: *ptr1-- = val; 229: } 230: } 231: 232: /* Test for negative result. */ 233: if (borrow == 1) 234: { 235: qguess--; 236: ptr1 = (unsigned char *) num1+qdig+len2; 237: ptr2 = (unsigned char *) n2ptr+len2-1; 238: carry = 0; 239: for (count = 0; count < len2; count++) 240: { 241: val = (int) *ptr1 + (int) *ptr2-- + carry; 242: if (val > 9) 243: { 244: val -= 10; 245: carry = 1; 246: } 247: else 248: carry = 0; 249: *ptr1-- = val; 250: } 251: if (carry == 1) *ptr1 = (*ptr1 + 1) % 10; 252: } 253: 254: /* We now know the quotient digit. */ 255: *qptr++ = qguess; 256: qdig++; 257: } 258: } 259: 260: /* Clean up and return the number. */ 261: qval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS ); 262: if (bc_is_zero (qval TSRMLS_CC)) qval->n_sign = PLUS; 263: _bc_rm_leading_zeros (qval); 264: bc_free_num (quot); 265: *quot = qval; 266: 267: /* Clean up temporary storage. */ 268: efree (mval); 269: efree (num1); 270: efree (num2); 271: 272: return 0; /* Everything is OK. */ 273: } 274: