xref: /PHP-7.2/ext/bcmath/libbcmath/src/div.c (revision 60a69dae)
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)88 bc_divide (bc_num n1, bc_num n2, bc_num *quot, int scale)
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)) 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   memset (num1, 0, n1->n_len+n1->n_scale+extra+2);
132   memcpy (num1+1, n1->n_value, n1->n_len+n1->n_scale);
133 
134   len2 = n2->n_len + scale2;
135   num2 = (unsigned char *) safe_emalloc (1, len2, 1);
136   memcpy (num2, n2->n_value, len2);
137   *(num2+len2) = 0;
138   n2ptr = num2;
139   while (*n2ptr == 0)
140     {
141       n2ptr++;
142       len2--;
143     }
144 
145   /* Calculate the number of quotient digits. */
146   if (len2 > len1+scale)
147     {
148       qdigits = scale+1;
149       zero = TRUE;
150     }
151   else
152     {
153       zero = FALSE;
154       if (len2>len1)
155 	qdigits = scale+1;  	/* One for the zero integer part. */
156       else
157 	qdigits = len1-len2+scale+1;
158     }
159 
160   /* Allocate and zero the storage for the quotient. */
161   qval = bc_new_num (qdigits-scale,scale);
162   memset (qval->n_value, 0, qdigits);
163 
164   /* Allocate storage for the temporary storage mval. */
165   mval = (unsigned char *) safe_emalloc (1, len2, 1);
166 
167   /* Now for the full divide algorithm. */
168   if (!zero)
169     {
170       /* Normalize */
171       norm =  10 / ((int)*n2ptr + 1);
172       if (norm != 1)
173 	{
174 	  _one_mult (num1, len1+scale1+extra+1, norm, num1);
175 	  _one_mult (n2ptr, len2, norm, n2ptr);
176 	}
177 
178       /* Initialize divide loop. */
179       qdig = 0;
180       if (len2 > len1)
181 	qptr = (unsigned char *) qval->n_value+len2-len1;
182       else
183 	qptr = (unsigned char *) qval->n_value;
184 
185       /* Loop */
186       while (qdig <= len1+scale-len2)
187 	{
188 	  /* Calculate the quotient digit guess. */
189 	  if (*n2ptr == num1[qdig])
190 	    qguess = 9;
191 	  else
192 	    qguess = (num1[qdig]*10 + num1[qdig+1]) / *n2ptr;
193 
194 	  /* Test qguess. */
195 	  if (n2ptr[1]*qguess >
196 	      (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10
197 	       + num1[qdig+2])
198 	    {
199 	      qguess--;
200 	      /* And again. */
201 	      if (n2ptr[1]*qguess >
202 		  (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10
203 		  + num1[qdig+2])
204 		qguess--;
205 	    }
206 
207 	  /* Multiply and subtract. */
208 	  borrow = 0;
209 	  if (qguess != 0)
210 	    {
211 	      *mval = 0;
212 	      _one_mult (n2ptr, len2, qguess, mval+1);
213 	      ptr1 = (unsigned char *) num1+qdig+len2;
214 	      ptr2 = (unsigned char *) mval+len2;
215 	      for (count = 0; count < len2+1; count++)
216 		{
217 		  val = (int) *ptr1 - (int) *ptr2-- - borrow;
218 		  if (val < 0)
219 		    {
220 		      val += 10;
221 		      borrow = 1;
222 		    }
223 		  else
224 		    borrow = 0;
225 		  *ptr1-- = val;
226 		}
227 	    }
228 
229 	  /* Test for negative result. */
230 	  if (borrow == 1)
231 	    {
232 	      qguess--;
233 	      ptr1 = (unsigned char *) num1+qdig+len2;
234 	      ptr2 = (unsigned char *) n2ptr+len2-1;
235 	      carry = 0;
236 	      for (count = 0; count < len2; count++)
237 		{
238 		  val = (int) *ptr1 + (int) *ptr2-- + carry;
239 		  if (val > 9)
240 		    {
241 		      val -= 10;
242 		      carry = 1;
243 		    }
244 		  else
245 		    carry = 0;
246 		  *ptr1-- = val;
247 		}
248 	      if (carry == 1) *ptr1 = (*ptr1 + 1) % 10;
249 	    }
250 
251 	  /* We now know the quotient digit. */
252 	  *qptr++ =  qguess;
253 	  qdig++;
254 	}
255     }
256 
257   /* Clean up and return the number. */
258   qval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS );
259   if (bc_is_zero (qval)) qval->n_sign = PLUS;
260   _bc_rm_leading_zeros (qval);
261   bc_free_num (quot);
262   *quot = qval;
263 
264   /* Clean up temporary storage. */
265   efree (mval);
266   efree (num1);
267   efree (num2);
268 
269   return 0;	/* Everything is OK. */
270 }
271