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