Return to exp.c CVS log

Up to [CSRG BSD Unix] / 43BSD / contrib / icon / rt

Annotation of 43BSD/contrib/icon/rt/exp.c, revision 1.1.1.1

1.1       root        1: /*
                      2:  * exp - calculate the exponential function of arg.
                      3:  *
                      4:  *   The coefficients are #1069 from Hart and Cheney. (22.35D)
                      5:  */
                      6: 
                      7: #include <errno.h>
                      8: #include <math.h>
                      9: 
                     10: int   errno;
                     11: static double   p0   = .2080384346694663001443843411e7;
                     12: static double   p1   = .3028697169744036299076048876e5;
                     13: static double   p2   = .6061485330061080841615584556e2;
                     14: static double   q0   = .6002720360238832528230907598e7;
                     15: static double   q1   = .3277251518082914423057964422e6;
                     16: static double   q2   = .1749287689093076403844945335e4;
                     17: static double   log2e   = 1.4426950408889634073599247;
                     18: static double   sqrt2   = 1.4142135623730950488016887;
                     19: static double   maxf   = 10000;
                     20: 
                     21: double
                     22: exp(arg)
                     23: double arg;
                     24:    {
                     25:    double fract;
                     26:    double temp1, temp2, xsq;
                     27:    int ent;
                     28: 
                     29:    if(arg == 0.)
                     30:       return(1.);
                     31:    if(arg < -maxf)
                     32:       return(0.);
                     33:    if(arg > maxf) {
                     34:       errno = ERANGE;
                     35:       return(HUGE);
                     36:    }
                     37:    arg *= log2e;
                     38:    ent = floor(arg);
                     39:    fract = (arg-ent) - 0.5;
                     40:    xsq = fract*fract;
                     41:    temp1 = ((p2*xsq+p1)*xsq+p0)*fract;
                     42:    temp2 = ((1.0*xsq+q2)*xsq+q1)*xsq + q0;
                     43:    return(ldexp(sqrt2*(temp2+temp1)/(temp2-temp1), ent));
                     44:    }

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.