xref: /PHP-5.4/ext/bcmath/libbcmath/src/div.c (revision d0ee4977)
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
_one_mult(num,size,digit,result)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
bc_divide(bc_num n1,bc_num n2,bc_num * quot,int scale TSRMLS_DC)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 
275