1
2/*
3 *  M_APM  -  mapm_div.c
4 *
5 *  Copyright (C) 1999 - 2007   Michael C. Ring
6 *
7 *  Permission to use, copy, and distribute this software and its
8 *  documentation for any purpose with or without fee is hereby granted,
9 *  provided that the above copyright notice appear in all copies and
10 *  that both that copyright notice and this permission notice appear
11 *  in supporting documentation.
12 *
13 *  Permission to modify the software is granted. Permission to distribute
14 *  the modified code is granted. Modifications are to be distributed by
15 *  using the file 'license.txt' as a template to modify the file header.
16 *  'license.txt' is available in the official MAPM distribution.
17 *
18 *  This software is provided "as is" without express or implied warranty.
19 */
20
21/*
22 *      $Id: mapm_div.c,v 1.12 2007/12/03 01:35:18 mike Exp $
23 *
24 *      This file contains the basic division functions
25 *
26 *      $Log: mapm_div.c,v $
27 *      Revision 1.12  2007/12/03 01:35:18  mike
28 *      Update license
29 *
30 *      Revision 1.11  2003/07/21 20:07:13  mike
31 *      Modify error messages to be in a consistent format.
32 *
33 *      Revision 1.10  2003/03/31 22:09:18  mike
34 *      call generic error handling function
35 *
36 *      Revision 1.9  2002/11/03 22:06:50  mike
37 *      Updated function parameters to use the modern style
38 *
39 *      Revision 1.8  2001/07/16 19:03:22  mike
40 *      add function M_free_all_div
41 *
42 *      Revision 1.7  2001/02/11 22:30:42  mike
43 *      modify parameters to REALLOC
44 *
45 *      Revision 1.6  2000/09/23 19:07:17  mike
46 *      change _divide to M_apm_sdivide function name
47 *
48 *      Revision 1.5  2000/04/11 18:38:55  mike
49 *      use new algorithm to determine q-hat. uses more digits of
50 *      the numerator and denominator.
51 *
52 *      Revision 1.4  2000/02/03 22:45:08  mike
53 *      use MAPM_* generic memory function
54 *
55 *      Revision 1.3  1999/06/23 01:10:49  mike
56 *      use predefined constant for '15'
57 *
58 *      Revision 1.2  1999/06/23 00:55:09  mike
59 *      change mult factor to 15
60 *
61 *      Revision 1.1  1999/05/10 20:56:31  mike
62 *      Initial revision
63 */
64
65#include "m_apm_lc.h"
66
67static	M_APM	M_div_worka;
68static	M_APM	M_div_workb;
69static	M_APM	M_div_tmp7;
70static	M_APM	M_div_tmp8;
71static	M_APM	M_div_tmp9;
72
73static	int	M_div_firsttime = TRUE;
74
75/****************************************************************************/
76void	M_free_all_div()
77{
78if (M_div_firsttime == FALSE)
79  {
80   m_apm_free(M_div_worka);
81   m_apm_free(M_div_workb);
82   m_apm_free(M_div_tmp7);
83   m_apm_free(M_div_tmp8);
84   m_apm_free(M_div_tmp9);
85
86   M_div_firsttime = TRUE;
87  }
88}
89/****************************************************************************/
90void	m_apm_integer_div_rem(M_APM qq, M_APM rr, M_APM aa, M_APM bb)
91{
92m_apm_integer_divide(qq, aa, bb);
93m_apm_multiply(M_div_tmp7, qq, bb);
94m_apm_subtract(rr, aa, M_div_tmp7);
95}
96/****************************************************************************/
97void	m_apm_integer_divide(M_APM rr, M_APM aa, M_APM bb)
98{
99/*
100 *    we must use this divide function since the
101 *    faster divide function using the reciprocal
102 *    will round the result (possibly changing
103 *    nnm.999999...  -->  nn(m+1).0000 which would
104 *    invalidate the 'integer_divide' goal).
105 */
106
107M_apm_sdivide(rr, 4, aa, bb);
108
109if (rr->m_apm_exponent <= 0)        /* result is 0 */
110  {
111   M_set_to_zero(rr);
112  }
113else
114  {
115   if (rr->m_apm_datalength > rr->m_apm_exponent)
116     {
117      rr->m_apm_datalength = rr->m_apm_exponent;
118      M_apm_normalize(rr);
119     }
120  }
121}
122/****************************************************************************/
123void	M_apm_sdivide(M_APM r, int places, M_APM a, M_APM b)
124{
125int	j, k, m, b0, sign, nexp, indexr, icompare, iterations;
126long    trial_numer;
127void	*vp;
128
129if (M_div_firsttime)
130  {
131   M_div_firsttime = FALSE;
132
133   M_div_worka = m_apm_init();
134   M_div_workb = m_apm_init();
135   M_div_tmp7  = m_apm_init();
136   M_div_tmp8  = m_apm_init();
137   M_div_tmp9  = m_apm_init();
138  }
139
140sign = a->m_apm_sign * b->m_apm_sign;
141
142if (sign == 0)      /* one number is zero, result is zero */
143  {
144   if (b->m_apm_sign == 0)
145     {
146      M_apm_log_error_msg(M_APM_RETURN, "\'M_apm_sdivide\', Divide by 0");
147     }
148
149   M_set_to_zero(r);
150   return;
151  }
152
153/*
154 *  Knuth step D1. Since base = 100, base / 2 = 50.
155 *  (also make the working copies positive)
156 */
157
158if (b->m_apm_data[0] >= 50)
159  {
160   m_apm_absolute_value(M_div_worka, a);
161   m_apm_absolute_value(M_div_workb, b);
162  }
163else       /* 'normal' step D1 */
164  {
165   k = 100 / (b->m_apm_data[0] + 1);
166   m_apm_set_long(M_div_tmp9, (long)k);
167
168   m_apm_multiply(M_div_worka, M_div_tmp9, a);
169   m_apm_multiply(M_div_workb, M_div_tmp9, b);
170
171   M_div_worka->m_apm_sign = 1;
172   M_div_workb->m_apm_sign = 1;
173  }
174
175/* setup trial denominator for step D3 */
176
177b0 = 100 * (int)M_div_workb->m_apm_data[0];
178
179if (M_div_workb->m_apm_datalength >= 3)
180  b0 += M_div_workb->m_apm_data[1];
181
182nexp = M_div_worka->m_apm_exponent - M_div_workb->m_apm_exponent;
183
184if (nexp > 0)
185  iterations = nexp + places + 1;
186else
187  iterations = places + 1;
188
189k = (iterations + 1) >> 1;     /* required size of result, in bytes */
190
191if (k > r->m_apm_malloclength)
192  {
193   if ((vp = MAPM_REALLOC(r->m_apm_data, (k + 32))) == NULL)
194     {
195      /* fatal, this does not return */
196
197      M_apm_log_error_msg(M_APM_FATAL, "\'M_apm_sdivide\', Out of memory");
198     }
199
200   r->m_apm_malloclength = k + 28;
201   r->m_apm_data = (UCHAR *)vp;
202  }
203
204/* clear the exponent in the working copies */
205
206M_div_worka->m_apm_exponent = 0;
207M_div_workb->m_apm_exponent = 0;
208
209/* if numbers are equal, ratio == 1.00000... */
210
211if ((icompare = m_apm_compare(M_div_worka, M_div_workb)) == 0)
212  {
213   iterations = 1;
214   r->m_apm_data[0] = 10;
215   nexp++;
216  }
217else			           /* ratio not 1, do the real division */
218  {
219   if (icompare == 1)                        /* numerator > denominator */
220     {
221      nexp++;                           /* to adjust the final exponent */
222      M_div_worka->m_apm_exponent += 1;     /* multiply numerator by 10 */
223     }
224   else                                      /* numerator < denominator */
225     {
226      M_div_worka->m_apm_exponent += 2;    /* multiply numerator by 100 */
227     }
228
229   indexr = 0;
230   m      = 0;
231
232   while (TRUE)
233     {
234      /*
235       *  Knuth step D3. Only use the 3rd -> 6th digits if the number
236       *  actually has that many digits.
237       */
238
239      trial_numer = 10000L * (long)M_div_worka->m_apm_data[0];
240
241      if (M_div_worka->m_apm_datalength >= 5)
242        {
243         trial_numer += 100 * M_div_worka->m_apm_data[1]
244                            + M_div_worka->m_apm_data[2];
245	}
246      else
247        {
248         if (M_div_worka->m_apm_datalength >= 3)
249           trial_numer += 100 * M_div_worka->m_apm_data[1];
250        }
251
252      j = (int)(trial_numer / b0);
253
254      /*
255       *    Since the library 'normalizes' all the results, we need
256       *    to look at the exponent of the number to decide if we
257       *    have a lead in 0n or 00.
258       */
259
260      if ((k = 2 - M_div_worka->m_apm_exponent) > 0)
261        {
262	 while (TRUE)
263	   {
264	    j /= 10;
265	    if (--k == 0)
266	      break;
267	   }
268	}
269
270      if (j == 100)     /* qhat == base ??      */
271        j = 99;         /* if so, decrease by 1 */
272
273      m_apm_set_long(M_div_tmp8, (long)j);
274      m_apm_multiply(M_div_tmp7, M_div_tmp8, M_div_workb);
275
276      /*
277       *    Compare our q-hat (j) against the desired number.
278       *    j is either correct, 1 too large, or 2 too large
279       *    per Theorem B on pg 272 of Art of Compter Programming,
280       *    Volume 2, 3rd Edition.
281       *
282       *    The above statement is only true if using the 2 leading
283       *    digits of the numerator and the leading digit of the
284       *    denominator. Since we are using the (3) leading digits
285       *    of the numerator and the (2) leading digits of the
286       *    denominator, we eliminate the case where our q-hat is
287       *    2 too large, (and q-hat being 1 too large is quite remote).
288       */
289
290      if (m_apm_compare(M_div_tmp7, M_div_worka) == 1)
291        {
292	 j--;
293         m_apm_subtract(M_div_tmp8, M_div_tmp7, M_div_workb);
294         m_apm_copy(M_div_tmp7, M_div_tmp8);
295	}
296
297      /*
298       *  Since we know q-hat is correct, step D6 is unnecessary.
299       *
300       *  Store q-hat, step D5. Since D6 is unnecessary, we can
301       *  do D5 before D4 and decide if we are done.
302       */
303
304      r->m_apm_data[indexr++] = (UCHAR)j;    /* j == 'qhat' */
305      m += 2;
306
307      if (m >= iterations)
308        break;
309
310      /* step D4 */
311
312      m_apm_subtract(M_div_tmp9, M_div_worka, M_div_tmp7);
313
314      /*
315       *  if the subtraction yields zero, the division is exact
316       *  and we are done early.
317       */
318
319      if (M_div_tmp9->m_apm_sign == 0)
320        {
321	 iterations = m;
322	 break;
323	}
324
325      /* multiply by 100 and re-save */
326      M_div_tmp9->m_apm_exponent += 2;
327      m_apm_copy(M_div_worka, M_div_tmp9);
328     }
329  }
330
331r->m_apm_sign       = sign;
332r->m_apm_exponent   = nexp;
333r->m_apm_datalength = iterations;
334
335M_apm_normalize(r);
336}
337/****************************************************************************/
338