Back to index

php5  5.3.10
div.c
Go to the documentation of this file.
00001 /* div.c: bcmath library file. */
00002 /*
00003     Copyright (C) 1991, 1992, 1993, 1994, 1997 Free Software Foundation, Inc.
00004     Copyright (C) 2000 Philip A. Nelson
00005 
00006     This library is free software; you can redistribute it and/or
00007     modify it under the terms of the GNU Lesser General Public
00008     License as published by the Free Software Foundation; either
00009     version 2 of the License, or (at your option) any later version.
00010 
00011     This library is distributed in the hope that it will be useful,
00012     but WITHOUT ANY WARRANTY; without even the implied warranty of
00013     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014     Lesser General Public License for more details.  (COPYING.LIB)
00015 
00016     You should have received a copy of the GNU Lesser General Public
00017     License along with this library; if not, write to:
00018 
00019       The Free Software Foundation, Inc.
00020       59 Temple Place, Suite 330
00021       Boston, MA 02111-1307 USA.
00022 
00023     You may contact the author by:
00024        e-mail:  philnelson@acm.org
00025       us-mail:  Philip A. Nelson
00026                 Computer Science Department, 9062
00027                 Western Washington University
00028                 Bellingham, WA 98226-9062
00029        
00030 *************************************************************************/
00031 
00032 #include <config.h>
00033 #include <stdio.h>
00034 #include <assert.h>
00035 #include <stdlib.h>
00036 #include <ctype.h>
00037 #include <stdarg.h>
00038 #include "bcmath.h"
00039 #include "private.h"
00040 
00041 
00042 /* Some utility routines for the divide:  First a one digit multiply.
00043    NUM (with SIZE digits) is multiplied by DIGIT and the result is
00044    placed into RESULT.  It is written so that NUM and RESULT can be
00045    the same pointers.  */
00046 
00047 static void
00048 _one_mult (num, size, digit, result)
00049      unsigned char *num;
00050      int size, digit;
00051      unsigned char *result;
00052 {
00053   int carry, value;
00054   unsigned char *nptr, *rptr;
00055 
00056   if (digit == 0)
00057     memset (result, 0, size);
00058   else
00059     {
00060       if (digit == 1)
00061        memcpy (result, num, size);
00062       else
00063        {
00064          /* Initialize */
00065          nptr = (unsigned char *) (num+size-1);
00066          rptr = (unsigned char *) (result+size-1);
00067          carry = 0;
00068 
00069          while (size-- > 0)
00070            {
00071              value = *nptr-- * digit + carry;
00072              *rptr-- = value % BASE;
00073              carry = value / BASE;
00074            }
00075 
00076          if (carry != 0) *rptr = carry;
00077        }
00078     }
00079 }
00080 
00081 
00082 /* The full division routine. This computes N1 / N2.  It returns
00083    0 if the division is ok and the result is in QUOT.  The number of
00084    digits after the decimal point is SCALE. It returns -1 if division
00085    by zero is tried.  The algorithm is found in Knuth Vol 2. p237. */
00086 
00087 int
00088 bc_divide (bc_num n1, bc_num n2, bc_num *quot, int scale TSRMLS_DC)
00089 {
00090   bc_num qval;
00091   unsigned char *num1, *num2;
00092   unsigned char *ptr1, *ptr2, *n2ptr, *qptr;
00093   int  scale1, val;
00094   unsigned int  len1, len2, scale2, qdigits, extra, count;
00095   unsigned int  qdig, qguess, borrow, carry;
00096   unsigned char *mval;
00097   char zero;
00098   unsigned int  norm;
00099 
00100   /* Test for divide by zero. */
00101   if (bc_is_zero (n2 TSRMLS_CC)) return -1;
00102 
00103   /* Test for divide by 1.  If it is we must truncate. */
00104   if (n2->n_scale == 0)
00105     {
00106       if (n2->n_len == 1 && *n2->n_value == 1)
00107        {
00108          qval = bc_new_num (n1->n_len, scale);
00109          qval->n_sign = (n1->n_sign == n2->n_sign ? PLUS : MINUS);
00110          memset (&qval->n_value[n1->n_len],0,scale);
00111          memcpy (qval->n_value, n1->n_value,
00112                 n1->n_len + MIN(n1->n_scale,scale));
00113          bc_free_num (quot);
00114          *quot = qval;
00115        }
00116     }
00117 
00118   /* Set up the divide.  Move the decimal point on n1 by n2's scale.
00119      Remember, zeros on the end of num2 are wasted effort for dividing. */
00120   scale2 = n2->n_scale;
00121   n2ptr = (unsigned char *) n2->n_value+n2->n_len+scale2-1;
00122   while ((scale2 > 0) && (*n2ptr-- == 0)) scale2--;
00123 
00124   len1 = n1->n_len + scale2;
00125   scale1 = n1->n_scale - scale2;
00126   if (scale1 < scale)
00127     extra = scale - scale1;
00128   else
00129     extra = 0;
00130   num1 = (unsigned char *) safe_emalloc (1, n1->n_len+n1->n_scale, extra+2);
00131   if (num1 == NULL) bc_out_of_memory();
00132   memset (num1, 0, n1->n_len+n1->n_scale+extra+2);
00133   memcpy (num1+1, n1->n_value, n1->n_len+n1->n_scale);
00134 
00135   len2 = n2->n_len + scale2;
00136   num2 = (unsigned char *) safe_emalloc (1, len2, 1);
00137   if (num2 == NULL) bc_out_of_memory();
00138   memcpy (num2, n2->n_value, len2);
00139   *(num2+len2) = 0;
00140   n2ptr = num2;
00141   while (*n2ptr == 0)
00142     {
00143       n2ptr++;
00144       len2--;
00145     }
00146 
00147   /* Calculate the number of quotient digits. */
00148   if (len2 > len1+scale)
00149     {
00150       qdigits = scale+1;
00151       zero = TRUE;
00152     }
00153   else
00154     {
00155       zero = FALSE;
00156       if (len2>len1)
00157        qdigits = scale+1;   /* One for the zero integer part. */
00158       else
00159        qdigits = len1-len2+scale+1;
00160     }
00161 
00162   /* Allocate and zero the storage for the quotient. */
00163   qval = bc_new_num (qdigits-scale,scale);
00164   memset (qval->n_value, 0, qdigits);
00165 
00166   /* Allocate storage for the temporary storage mval. */
00167   mval = (unsigned char *) safe_emalloc (1, len2, 1);
00168   if (mval == NULL) bc_out_of_memory ();
00169 
00170   /* Now for the full divide algorithm. */
00171   if (!zero)
00172     {
00173       /* Normalize */
00174       norm =  10 / ((int)*n2ptr + 1);
00175       if (norm != 1)
00176        {
00177          _one_mult (num1, len1+scale1+extra+1, norm, num1);
00178          _one_mult (n2ptr, len2, norm, n2ptr);
00179        }
00180 
00181       /* Initialize divide loop. */
00182       qdig = 0;
00183       if (len2 > len1)
00184        qptr = (unsigned char *) qval->n_value+len2-len1;
00185       else
00186        qptr = (unsigned char *) qval->n_value;
00187 
00188       /* Loop */
00189       while (qdig <= len1+scale-len2)
00190        {
00191          /* Calculate the quotient digit guess. */
00192          if (*n2ptr == num1[qdig])
00193            qguess = 9;
00194          else
00195            qguess = (num1[qdig]*10 + num1[qdig+1]) / *n2ptr;
00196 
00197          /* Test qguess. */
00198          if (n2ptr[1]*qguess >
00199              (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10
00200               + num1[qdig+2])
00201            {
00202              qguess--;
00203              /* And again. */
00204              if (n2ptr[1]*qguess >
00205                 (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10
00206                 + num1[qdig+2])
00207               qguess--;
00208            }
00209 
00210          /* Multiply and subtract. */
00211          borrow = 0;
00212          if (qguess != 0)
00213            {
00214              *mval = 0;
00215              _one_mult (n2ptr, len2, qguess, mval+1);
00216              ptr1 = (unsigned char *) num1+qdig+len2;
00217              ptr2 = (unsigned char *) mval+len2;
00218              for (count = 0; count < len2+1; count++)
00219               {
00220                 val = (int) *ptr1 - (int) *ptr2-- - borrow;
00221                 if (val < 0)
00222                   {
00223                     val += 10;
00224                     borrow = 1;
00225                   }
00226                 else
00227                   borrow = 0;
00228                 *ptr1-- = val;
00229               }
00230            }
00231 
00232          /* Test for negative result. */
00233          if (borrow == 1)
00234            {
00235              qguess--;
00236              ptr1 = (unsigned char *) num1+qdig+len2;
00237              ptr2 = (unsigned char *) n2ptr+len2-1;
00238              carry = 0;
00239              for (count = 0; count < len2; count++)
00240               {
00241                 val = (int) *ptr1 + (int) *ptr2-- + carry;
00242                 if (val > 9)
00243                   {
00244                     val -= 10;
00245                     carry = 1;
00246                   }
00247                 else
00248                   carry = 0;
00249                 *ptr1-- = val;
00250               }
00251              if (carry == 1) *ptr1 = (*ptr1 + 1) % 10;
00252            }
00253 
00254          /* We now know the quotient digit. */
00255          *qptr++ =  qguess;
00256          qdig++;
00257        }
00258     }
00259 
00260   /* Clean up and return the number. */
00261   qval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS );
00262   if (bc_is_zero (qval TSRMLS_CC)) qval->n_sign = PLUS;
00263   _bc_rm_leading_zeros (qval);
00264   bc_free_num (quot);
00265   *quot = qval;
00266 
00267   /* Clean up temporary storage. */
00268   efree (mval);
00269   efree (num1);
00270   efree (num2);
00271 
00272   return 0;   /* Everything is OK. */
00273 }
00274