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
67 static M_APM M_div_worka;
68 static M_APM M_div_workb;
69 static M_APM M_div_tmp7;
70 static M_APM M_div_tmp8;
71 static M_APM M_div_tmp9;
72
73 static int M_div_firsttime = TRUE;
74
75 /****************************************************************************/
M_free_all_div()76 void M_free_all_div()
77 {
78 if (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 /****************************************************************************/
m_apm_integer_div_rem(M_APM qq,M_APM rr,M_APM aa,M_APM bb)90 void m_apm_integer_div_rem(M_APM qq, M_APM rr, M_APM aa, M_APM bb)
91 {
92 m_apm_integer_divide(qq, aa, bb);
93 m_apm_multiply(M_div_tmp7, qq, bb);
94 m_apm_subtract(rr, aa, M_div_tmp7);
95 }
96 /****************************************************************************/
m_apm_integer_divide(M_APM rr,M_APM aa,M_APM bb)97 void 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
107 M_apm_sdivide(rr, 4, aa, bb);
108
109 if (rr->m_apm_exponent <= 0) /* result is 0 */
110 {
111 M_set_to_zero(rr);
112 }
113 else
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 /****************************************************************************/
M_apm_sdivide(M_APM r,int places,M_APM a,M_APM b)123 void M_apm_sdivide(M_APM r, int places, M_APM a, M_APM b)
124 {
125 int j, k, m, b0, sign, nexp, indexr, icompare, iterations;
126 long trial_numer;
127 void *vp;
128
129 if (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
140 sign = a->m_apm_sign * b->m_apm_sign;
141
142 if (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
158 if (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 }
163 else /* '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
177 b0 = 100 * (int)M_div_workb->m_apm_data[0];
178
179 if (M_div_workb->m_apm_datalength >= 3)
180 b0 += M_div_workb->m_apm_data[1];
181
182 nexp = M_div_worka->m_apm_exponent - M_div_workb->m_apm_exponent;
183
184 if (nexp > 0)
185 iterations = nexp + places + 1;
186 else
187 iterations = places + 1;
188
189 k = (iterations + 1) >> 1; /* required size of result, in bytes */
190
191 if (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
206 M_div_worka->m_apm_exponent = 0;
207 M_div_workb->m_apm_exponent = 0;
208
209 /* if numbers are equal, ratio == 1.00000... */
210
211 if ((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 }
217 else /* 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
331 r->m_apm_sign = sign;
332 r->m_apm_exponent = nexp;
333 r->m_apm_datalength = iterations;
334
335 M_apm_normalize(r);
336 }
337 /****************************************************************************/
338