Annotation of embedaddon/php/ext/bcmath/libbcmath/src/div.c, revision 1.1.1.1

1.1       misho       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 TSRMLS_DC)
                     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 TSRMLS_CC)) 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 TSRMLS_CC)) 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: 

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>