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 /****************************************************************************/ 409