src/libs/mapm/mapm_exp.c

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