Return to exp.c CVS log

Up to [Research Unix] / researchv9 / libc / math

Annotation of researchv9/libc/math/exp.c, revision 1.1.1.1

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

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