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
144static  M_APM  MM_exp_log2R;
145static  M_APM  MM_exp_512R;
146static	int    MM_firsttime1 = TRUE;
147
148/****************************************************************************/
149void	M_free_all_exp()
150{
151if (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/****************************************************************************/
160void	m_apm_exp(M_APM r, int places, M_APM x)
161{
162M_APM   tmp7, tmp8, tmp9;
163int	dplaces, nn, ii;
164
165if (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
176tmp7 = M_get_stack_var();
177tmp8 = M_get_stack_var();
178tmp9 = M_get_stack_var();
179
180if (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
187if (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
212if (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
222dplaces = places + 8;
223
224/* check to make sure our log(2) is accurate enough */
225
226M_check_log_places(dplaces);
227
228m_apm_multiply(tmp8, tmp7, MM_lc_log2);
229m_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
236while (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
258m_apm_multiply(tmp9, tmp7, MM_exp_512R);
259
260/* perform the series expansion ... */
261
262M_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
270ii = 9;
271
272while (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
283m_apm_integer_pow(tmp7, dplaces, MM_Two, nn);
284
285m_apm_multiply(tmp9, tmp7, tmp8);
286m_apm_round(r, places, tmp9);
287
288M_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*/
298int	M_exp_compute_nn(int *n, M_APM b, M_APM a)
299{
300M_APM	tmp0, tmp1;
301void	*vp;
302char    *cp, sbuf[48];
303int	kk;
304
305*n   = 0;
306vp   = NULL;
307cp   = sbuf;
308tmp0 = M_get_stack_var();
309tmp1 = 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
314m_apm_multiply(tmp1, a, MM_exp_log2R);
315
316/* round to the nearest int */
317
318if (tmp1->m_apm_sign >= 0)
319  {
320   m_apm_add(tmp0, tmp1, MM_0_5);
321   m_apm_floor(tmp1, tmp0);
322  }
323else
324  {
325   m_apm_subtract(tmp0, tmp1, MM_0_5);
326   m_apm_ceil(tmp1, tmp0);
327  }
328
329kk = tmp1->m_apm_exponent;
330if (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
342m_apm_to_integer_string(cp, tmp1);
343*n = atoi(cp);
344
345m_apm_set_long(b, (long)(*n));
346
347kk = m_apm_compare(b, tmp1);
348
349if (vp != NULL)
350  MAPM_FREE(vp);
351
352M_restore_stack(2);
353return(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*/
364void	M_raw_exp(M_APM rr, int places, M_APM xx)
365{
366M_APM   tmp0, digit, term;
367int	tolerance,  local_precision, prev_exp;
368long    m1;
369
370tmp0  = M_get_stack_var();
371term  = M_get_stack_var();
372digit = M_get_stack_var();
373
374local_precision = places + 8;
375tolerance       = -(places + 4);
376prev_exp        = 0;
377
378m_apm_add(rr, MM_One, xx);
379m_apm_copy(term, xx);
380
381m1 = 2L;
382
383while (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
406M_restore_stack(3);                    /* restore the 3 locals we used here */
407}
408/****************************************************************************/
409