ext/bcmath/libbcmath/src/div.c

/* [<][>][^][v][top][bottom][index][help] */
This source file includes following definitions.
_one_mult
bc_divide
   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
  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
  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   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)) 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
/* [<][>][^][v][top][bottom][index][help] */
root/ext/bcmath/libbcmath/src/div.c

DEFINITIONS
