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