Annotation of embedaddon/php/ext/bcmath/libbcmath/src/div.c, revision 1.1.1.1
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:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>