home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / files / language / icon8_fx / exp.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-10-23  |  6.2 KB  |  221 lines
  1. /************************************************************************
  2.  *                                                                      *
  3.  * A replacement exp() routine for PML library.  An original algorithm  *
  4.  * is not changed but a table of fractionals power of 2 was recalcul-   *
  5.  * ated (believe it or not - this has an influence on a final result).  *
  6.  * Also divisions by power of 2 were replaced by applications of ldexp  *
  7.  * routine which is faster and better preserves precision               *
  8.  *                                                                      *
  9.  * Michal Jaegermann, December 1990                                     *    
  10.  *                                                                      *
  11.  ************************************************************************
  12.  */
  13.   /************************************************************************
  14.  *                                    *
  15.  *                N O T I C E                *
  16.  *                                    *
  17.  *            Copyright Abandoned, 1987, Fred Fish        *
  18.  *                                    *
  19.  *    This previously copyrighted work has been placed into the    *
  20.  *    public domain by the author (Fred Fish) and may be freely used    *
  21.  *    for any purpose, private or commercial.  I would appreciate    *
  22.  *    it, as a courtesy, if this notice is left in all copies and    *
  23.  *    derivative works.  Thank you, and enjoy...            *
  24.  *                                    *
  25.  *    The author makes no warranty of any kind with respect to this    *
  26.  *    product and explicitly disclaims any implied warranties of    *
  27.  *    merchantability or fitness for any particular purpose.        *
  28.  *                                    *
  29.  ************************************************************************
  30.  */
  31.  
  32.  
  33. /*
  34.  *  FUNCTION
  35.  *
  36.  *    exp   double precision exponential
  37.  *
  38.  *  KEY WORDS
  39.  *
  40.  *    exp
  41.  *    machine independent routines
  42.  *    math libraries
  43.  *
  44.  *  DESCRIPTION
  45.  *
  46.  *    Returns double precision exponential of double precision
  47.  *    floating point number.
  48.  *
  49.  *  USAGE
  50.  *
  51.  *    double exp (x)
  52.  *    double x;
  53.  *
  54.  *  REFERENCES
  55.  *
  56.  *    Fortran IV plus user's guide, Digital Equipment Corp. pp B-3
  57.  *
  58.  *    Computer Approximations, J.F. Hart et al, John Wiley & Sons,
  59.  *    1968, pp. 96-104.
  60.  *
  61.  *  RESTRICTIONS
  62.  *
  63.  *    Inputs greater than log(MAXDOUBLE) result in overflow.
  64.  *    Inputs less than log(MINDOUBLE) result in underflow.
  65.  *
  66.  *    The maximum relative error for the approximating polynomial
  67.  *    is 10**(-16.4).  However, this assumes exact arithmetic
  68.  *    in the polynomial evaluation.  Additional rounding and
  69.  *    truncation errors may occur as the argument is reduced
  70.  *    to the range over which the polynomial approximation
  71.  *    is valid, and as the polynomial is evaluated using
  72.  *    finite precision arithmetic.
  73.  *    
  74.  *  PROGRAMMER
  75.  *
  76.  *    Fred Fish
  77.  *
  78.  *  INTERNALS
  79.  *
  80.  *    Computes exponential from:
  81.  *
  82.  *        exp(x)    =    2**y  *  2**z  *  2**w
  83.  *
  84.  *    Where:
  85.  *
  86.  *        y    =    int ( x * log2(e) )
  87.  *
  88.  *        v    =    16 * frac ( x * log2(e))
  89.  *
  90.  *        z    =    (1/16) * int (v)
  91.  *
  92.  *        w    =    (1/16) * frac (v)
  93.  *
  94.  *    Note that:
  95.  *
  96.  *        0 =< v < 16
  97.  *
  98.  *        z = {0, 1/16, 2/16, ...15/16}
  99.  *
  100.  *        0 =< w < 1/16
  101.  *
  102.  *    Then:
  103.  *
  104.  *        2**z is looked up in a table of 2**0, 2**1/16, ...
  105.  *
  106.  *        2**w is computed from an approximation:
  107.  *
  108.  *            2**w    =  (A + B) / (A - B)
  109.  *
  110.  *            A    =  C + (D * w * w)
  111.  *
  112.  *            B    =  w * (E + (F * w * w))
  113.  *
  114.  *            C    =  20.8137711965230361973
  115.  *
  116.  *            D    =  1.0
  117.  *
  118.  *            E    =  7.2135034108448192083
  119.  *
  120.  *            F    =  0.057761135831801928
  121.  *
  122.  *        Coefficients are from HART, table #1121, pg 206.
  123.  *
  124.  *        Effective multiplication by 2**y is done by a
  125.  *        floating point scale with y as scale argument.
  126.  *
  127.  */
  128.  
  129. #include <stdio.h>
  130. #include <pmluser.h>
  131. #include "pml.h"
  132.  
  133. # define C  20.8137711965230361973    /* Polynomial approx coeff.    */
  134. /* # define D  1.0    */              /* Polynomial approx coeff.    */
  135. # define E  7.2135034108448192083    /* Polynomial approx coeff.    */
  136. # define F  0.057761135831801928    /* Polynomial approx coeff.    */
  137.  
  138.  
  139. /************************************************************************
  140.  *                                    *
  141.  *    This table is fixed in size and reasonably hardware        *
  142.  *    independent.  The given constants were computed on Atari ST     *
  143.  *      using integer arithmetic and decimal values for 2**(1/2),       *
  144.  *      2**(1/4) and 2**(1/8) taken from Appendix C of HART et al.      *
  145.  *                                    *
  146.  ************************************************************************
  147.  */
  148.  
  149. static double fpof2tbl[] = {
  150.     1.0/*0000000000000000000*/,        /*    2 ** 0/16    */
  151.     1.04427378242741384032,        /*    2 ** 1/16    */
  152.     1.09050773266525765920,        /*    2 ** 2/16    */
  153.     1.13878863475669165369,        /*    2 ** 3/16    */
  154.     1.18920711500272106671,        /*    2 ** 4/16    */
  155.     1.24185781207348404858,        /*    2 ** 5/16    */
  156.     1.29683955465100966592,        /*    2 ** 6/16    */
  157.     1.35425554693689272828,        /*    2 ** 7/16    */
  158.     1.41421356237309504880,        /*    2 ** 8/16    */
  159.     1.47682614593949931138,        /*    2 ** 9/16    */
  160.     1.54221082540794082360,        /*    2 ** 10/16    */
  161.     1.61049033194925430816,        /*    2 ** 11/16    */
  162.     1.68179283050742908605,        /*    2 ** 12/16    */
  163.     1.75625216037329948310,        /*    2 ** 13/16    */
  164.     1.83400808640934246346,        /*    2 ** 14/16    */
  165.     1.91520656139714729384        /*    2 ** 15/16    */
  166. };
  167.  
  168. static char funcname[] = "exp";
  169.  
  170.  
  171. double exp (x)
  172. double x;
  173. {
  174.     register int index;
  175.     auto int y;
  176.     auto double w;
  177.     auto double a;
  178.     auto double b;
  179.     auto double temp;
  180.     auto char * xcptstr;
  181.     extern double dabs ();
  182.     extern double ldexp ();
  183.     extern double modf ();
  184.     auto struct exception xcpt;
  185.  
  186.     DBUG_ENTER (funcname);
  187.     DBUG_3 ("expin", "arg %le", x);
  188.     if (x > LOGE_MAXDOUBLE) {
  189.     xcpt.type = OVERFLOW;
  190.     xcpt.retval = MAXDOUBLE;
  191.     xcptstr = "OVER";
  192.     goto eset;
  193.     }
  194.     if (x <= LOGE_MINDOUBLE) {
  195.     xcpt.type = UNDERFLOW;
  196.     xcpt.retval = 0.0;
  197.     xcptstr = "UNDER";
  198. eset:
  199.     xcpt.name = funcname;
  200.     xcpt.arg1 = x;
  201.     if (!matherr (&xcpt)) {
  202.             errno = ERANGE;
  203.         fprintf (stderr, "%s: %sFLOW error\n", funcname, xcptstr);
  204.     }
  205.     } else {
  206.     x *= LOG2E;
  207.     w = ldexp(modf(x,&a), 4);
  208.     y = a;
  209.     w = ldexp(modf(w, &temp), -4);
  210.     index = temp;
  211.     
  212.     b = w * w;
  213.     a = C + b;
  214.     b = w * (E + (F * b));
  215.     xcpt.retval = ldexp (((a + b) / (a - b)) *
  216.         (index < 0 ? ldexp(fpof2tbl[16 + index], -1) : fpof2tbl[index]), y);
  217.     }
  218.     DBUG_3 ("expout", "result %le", xcpt.retval);
  219.     DBUG_RETURN (xcpt.retval);
  220. }
  221.