"Fossies" - the Fresh Open Source Software Archive 
Member "mapm_4.9.5a/mapm_exp.c" (21 Feb 2010, 10480 Bytes) of package /linux/misc/old/mapm-4.9.5a.tar.gz:
As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style:
standard) with prefixed line numbers and
code folding option.
Alternatively you can here
view or
download the uninterpreted source code file.
For more information about "mapm_exp.c" see the
Fossies "Dox" file reference documentation.
1
2 /*
3 * M_APM - mapm_exp.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_exp.c,v 1.37 2007/12/03 01:36:00 mike Exp $
23 *
24 * This file contains the EXP function.
25 *
26 * $Log: mapm_exp.c,v $
27 * Revision 1.37 2007/12/03 01:36:00 mike
28 * Update license
29 *
30 * Revision 1.36 2004/06/02 00:29:03 mike
31 * simplify logic in compute_nn
32 *
33 * Revision 1.35 2004/05/29 18:29:44 mike
34 * move exp nn calculation into its own function
35 *
36 * Revision 1.34 2004/05/21 20:41:01 mike
37 * fix potential buffer overflow
38 *
39 * Revision 1.33 2004/02/18 02:46:45 mike
40 * fix comment
41 *
42 * Revision 1.32 2004/02/18 02:41:35 mike
43 * check to make sure 'nn' does not overflow
44 *
45 * Revision 1.31 2004/01/01 00:06:38 mike
46 * make a comment more clear
47 *
48 * Revision 1.30 2003/12/31 21:44:57 mike
49 * simplify logic in _exp
50 *
51 * Revision 1.29 2003/12/28 00:03:27 mike
52 * dont allow 'tmp7' to get too small prior to divide by 256
53 *
54 * Revision 1.28 2003/12/27 22:53:04 mike
55 * change 1024 to 512, if input is already small, call
56 * raw_exp directly and return
57 *
58 * Revision 1.27 2003/03/30 21:16:37 mike
59 * use global version of log(2) instead of local copy
60 *
61 * Revision 1.26 2002/11/03 22:30:18 mike
62 * Updated function parameters to use the modern style
63 *
64 * Revision 1.25 2001/08/06 22:07:00 mike
65 * round the 'nn' calculation to the nearest int
66 *
67 * Revision 1.24 2001/08/05 21:58:59 mike
68 * make 1 / log(2) constant shorter
69 *
70 * Revision 1.23 2001/07/16 19:10:23 mike
71 * add function M_free_all_exp
72 *
73 * Revision 1.22 2001/07/10 21:55:36 mike
74 * optimize raw_exp by using fewer digits as the
75 * subsequent terms get smaller
76 *
77 * Revision 1.21 2001/02/06 21:20:31 mike
78 * optimize 'nn' calculation (for the future)
79 *
80 * Revision 1.20 2001/02/05 21:55:12 mike
81 * minor optimization, use a multiply instead
82 * of a divide
83 *
84 * Revision 1.19 2000/08/22 21:34:41 mike
85 * increase local precision
86 *
87 * Revision 1.18 2000/05/18 22:05:22 mike
88 * move _pow to a separate file
89 *
90 * Revision 1.17 2000/05/04 23:21:01 mike
91 * use global var 256R
92 *
93 * Revision 1.16 2000/03/30 21:33:30 mike
94 * change termination of raw_exp to use ints, not MAPM numbers
95 *
96 * Revision 1.15 2000/02/05 22:59:46 mike
97 * adjust decimal places on calculation
98 *
99 * Revision 1.14 2000/02/04 20:45:21 mike
100 * re-compute log(2) on the fly if we need a more precise value
101 *
102 * Revision 1.13 2000/02/04 19:35:14 mike
103 * use just an approx log(2) for the integer divide
104 *
105 * Revision 1.12 2000/02/04 16:47:32 mike
106 * use the real algorithm for EXP
107 *
108 * Revision 1.11 1999/09/18 01:27:40 mike
109 * if X is 0 on the pow function, return 0 right away
110 *
111 * Revision 1.10 1999/06/19 20:54:07 mike
112 * changed local static MAPM to stack variables
113 *
114 * Revision 1.9 1999/06/01 22:37:44 mike
115 * adjust decimal places passed to raw function
116 *
117 * Revision 1.8 1999/06/01 01:44:03 mike
118 * change threshold from 1000 to 100 for 65536 divisor
119 *
120 * Revision 1.7 1999/06/01 01:03:31 mike
121 * vary 'q' instead of checking input against +/- 10 and +/- 40
122 *
123 * Revision 1.6 1999/05/15 01:54:27 mike
124 * add check for number of decimal places
125 *
126 * Revision 1.5 1999/05/15 01:09:49 mike
127 * minor tweak to POW decimal places
128 *
129 * Revision 1.4 1999/05/13 00:14:00 mike
130 * added more comments
131 *
132 * Revision 1.3 1999/05/12 23:39:05 mike
133 * change #places passed to sub functions
134 *
135 * Revision 1.2 1999/05/10 21:35:13 mike
136 * added some comments
137 *
138 * Revision 1.1 1999/05/10 20:56:31 mike
139 * Initial revision
140 */
141
142 #include "m_apm_lc.h"
143
144 static M_APM MM_exp_log2R;
145 static M_APM MM_exp_512R;
146 static int MM_firsttime1 = TRUE;
147
148 /****************************************************************************/
149 void M_free_all_exp()
150 {
151 if (MM_firsttime1 == FALSE)
152 {
153 m_apm_free(MM_exp_log2R);
154 m_apm_free(MM_exp_512R);
155
156 MM_firsttime1 = TRUE;
157 }
158 }
159 /****************************************************************************/
160 void m_apm_exp(M_APM r, int places, M_APM x)
161 {
162 M_APM tmp7, tmp8, tmp9;
163 int dplaces, nn, ii;
164
165 if (MM_firsttime1)
166 {
167 MM_firsttime1 = FALSE;
168
169 MM_exp_log2R = m_apm_init();
170 MM_exp_512R = m_apm_init();
171
172 m_apm_set_string(MM_exp_log2R, "1.44269504089"); /* ~ 1 / log(2) */
173 m_apm_set_string(MM_exp_512R, "1.953125E-3"); /* 1 / 512 */
174 }
175
176 tmp7 = M_get_stack_var();
177 tmp8 = M_get_stack_var();
178 tmp9 = M_get_stack_var();
179
180 if (x->m_apm_sign == 0) /* if input == 0, return '1' */
181 {
182 m_apm_copy(r, MM_One);
183 M_restore_stack(3);
184 return;
185 }
186
187 if (x->m_apm_exponent <= -3) /* already small enough so call _raw directly */
188 {
189 M_raw_exp(tmp9, (places + 6), x);
190 m_apm_round(r, places, tmp9);
191 M_restore_stack(3);
192 return;
193 }
194
195 /*
196 From David H. Bailey's MPFUN Fortran package :
197
198 exp (t) = (1 + r + r^2 / 2! + r^3 / 3! + r^4 / 4! ...) ^ q * 2 ^ n
199
200 where q = 256, r = t' / q, t' = t - n Log(2) and where n is chosen so
201 that -0.5 Log(2) < t' <= 0.5 Log(2). Reducing t mod Log(2) and
202 dividing by 256 insures that -0.001 < r <= 0.001, which accelerates
203 convergence in the above series.
204
205 I use q = 512 and also limit how small 'r' can become. The 'r' used
206 here is limited in magnitude from 1.95E-4 < |r| < 1.35E-3. Forcing
207 'r' into a narrow range keeps the algorithm 'well behaved'.
208
209 ( the range is [0.1 / 512] to [log(2) / 512] )
210 */
211
212 if (M_exp_compute_nn(&nn, tmp7, x) != 0)
213 {
214 M_apm_log_error_msg(M_APM_RETURN,
215 "\'m_apm_exp\', Input too large, Overflow");
216
217 M_set_to_zero(r);
218 M_restore_stack(3);
219 return;
220 }
221
222 dplaces = places + 8;
223
224 /* check to make sure our log(2) is accurate enough */
225
226 M_check_log_places(dplaces);
227
228 m_apm_multiply(tmp8, tmp7, MM_lc_log2);
229 m_apm_subtract(tmp7, x, tmp8);
230
231 /*
232 * guarantee that |tmp7| is between 0.1 and 0.9999999....
233 * (in practice, the upper limit only reaches log(2), 0.693... )
234 */
235
236 while (TRUE)
237 {
238 if (tmp7->m_apm_sign != 0)
239 {
240 if (tmp7->m_apm_exponent == 0)
241 break;
242 }
243
244 if (tmp7->m_apm_sign >= 0)
245 {
246 nn++;
247 m_apm_subtract(tmp8, tmp7, MM_lc_log2);
248 m_apm_copy(tmp7, tmp8);
249 }
250 else
251 {
252 nn--;
253 m_apm_add(tmp8, tmp7, MM_lc_log2);
254 m_apm_copy(tmp7, tmp8);
255 }
256 }
257
258 m_apm_multiply(tmp9, tmp7, MM_exp_512R);
259
260 /* perform the series expansion ... */
261
262 M_raw_exp(tmp8, dplaces, tmp9);
263
264 /*
265 * raise result to the 512 power
266 *
267 * note : x ^ 512 = (((x ^ 2) ^ 2) ^ 2) ... 9 times
268 */
269
270 ii = 9;
271
272 while (TRUE)
273 {
274 m_apm_multiply(tmp9, tmp8, tmp8);
275 m_apm_round(tmp8, dplaces, tmp9);
276
277 if (--ii == 0)
278 break;
279 }
280
281 /* now compute 2 ^ N */
282
283 m_apm_integer_pow(tmp7, dplaces, MM_Two, nn);
284
285 m_apm_multiply(tmp9, tmp7, tmp8);
286 m_apm_round(r, places, tmp9);
287
288 M_restore_stack(3); /* restore the 3 locals we used here */
289 }
290 /****************************************************************************/
291 /*
292 compute int *n = round_to_nearest_int(a / log(2))
293 M_APM b = MAPM version of *n
294
295 returns 0: OK
296 -1, 1: failure
297 */
298 int M_exp_compute_nn(int *n, M_APM b, M_APM a)
299 {
300 M_APM tmp0, tmp1;
301 void *vp;
302 char *cp, sbuf[48];
303 int kk;
304
305 *n = 0;
306 vp = NULL;
307 cp = sbuf;
308 tmp0 = M_get_stack_var();
309 tmp1 = M_get_stack_var();
310
311 /* find 'n' and convert it to a normal C int */
312 /* we just need an approx 1/log(2) for this calculation */
313
314 m_apm_multiply(tmp1, a, MM_exp_log2R);
315
316 /* round to the nearest int */
317
318 if (tmp1->m_apm_sign >= 0)
319 {
320 m_apm_add(tmp0, tmp1, MM_0_5);
321 m_apm_floor(tmp1, tmp0);
322 }
323 else
324 {
325 m_apm_subtract(tmp0, tmp1, MM_0_5);
326 m_apm_ceil(tmp1, tmp0);
327 }
328
329 kk = tmp1->m_apm_exponent;
330 if (kk >= 42)
331 {
332 if ((vp = (void *)MAPM_MALLOC((kk + 16) * sizeof(char))) == NULL)
333 {
334 /* fatal, this does not return */
335
336 M_apm_log_error_msg(M_APM_FATAL, "\'M_exp_compute_nn\', Out of memory");
337 }
338
339 cp = (char *)vp;
340 }
341
342 m_apm_to_integer_string(cp, tmp1);
343 *n = atoi(cp);
344
345 m_apm_set_long(b, (long)(*n));
346
347 kk = m_apm_compare(b, tmp1);
348
349 if (vp != NULL)
350 MAPM_FREE(vp);
351
352 M_restore_stack(2);
353 return(kk);
354 }
355 /****************************************************************************/
356 /*
357 calculate the exponential function using the following series :
358
359 x^2 x^3 x^4 x^5
360 exp(x) == 1 + x + --- + --- + --- + --- ...
361 2! 3! 4! 5!
362
363 */
364 void M_raw_exp(M_APM rr, int places, M_APM xx)
365 {
366 M_APM tmp0, digit, term;
367 int tolerance, local_precision, prev_exp;
368 long m1;
369
370 tmp0 = M_get_stack_var();
371 term = M_get_stack_var();
372 digit = M_get_stack_var();
373
374 local_precision = places + 8;
375 tolerance = -(places + 4);
376 prev_exp = 0;
377
378 m_apm_add(rr, MM_One, xx);
379 m_apm_copy(term, xx);
380
381 m1 = 2L;
382
383 while (TRUE)
384 {
385 m_apm_set_long(digit, m1);
386 m_apm_multiply(tmp0, term, xx);
387 m_apm_divide(term, local_precision, tmp0, digit);
388 m_apm_add(tmp0, rr, term);
389 m_apm_copy(rr, tmp0);
390
391 if ((term->m_apm_exponent < tolerance) || (term->m_apm_sign == 0))
392 break;
393
394 if (m1 != 2L)
395 {
396 local_precision = local_precision + term->m_apm_exponent - prev_exp;
397
398 if (local_precision < 20)
399 local_precision = 20;
400 }
401
402 prev_exp = term->m_apm_exponent;
403 m1++;
404 }
405
406 M_restore_stack(3); /* restore the 3 locals we used here */
407 }
408 /****************************************************************************/