Return to expd.c CVS log

Up to [CSRG BSD Unix] / 43BSDTahoe / cci / libM

Annotation of 43BSDTahoe/cci/libM/expd.c, revision 1.1.1.1

1.1       root        1: /*     @(#)exp.c       4.1/4.2         10/31/84        CCI-CPG */
                      2: 
                      3: 
                      4: /*
                      5:  * The double-precision 'exp' returns the exponential
                      6:  * function of its floating-point argument.
                      7:  *
                      8:  * New version by Les Powers (3/23/85).
                      9:  */
                     10: 
                     11: 
                     12: #include <errno.h>
                     13: 
                     14: 
                     15: 
                     16: /*
                     17:  * The following number 'forces' the hex value of 0x7fffffff 0xffffffff
                     18:  * to be used for HUGE.
                     19:  */
                     20: #define        HUGE            1.701411834604692350e+40
                     21: 
                     22: #define EXPONENT28     0x0e000000
                     23: #define EXP_SIZE0      12
                     24: #define EXP_SIZE1      128
                     25: #define EXP_SIZE2      256
                     26: #define EXP_SIZE3      256
                     27: #define EXP_SIZE4      256
                     28: 
                     29: int    errno;
                     30: 
                     31: /*
                     32:  * This is the natural log of the smallest number that can
                     33:  * be represented with double precision format (2^-128).
                     34:  */
                     35: double minf = -88.7228391116729996054057115466466007136640;
                     36: 
                     37: /*
                     38:  * This is the natural log of the biggest number that can
                     39:  * be represented with double precision format (2^127 (1-2^-56)).
                     40:  */
                     41: double maxf = 88.0296919311130542821106916173739672939966;
                     42: 
                     43: double _ep0[EXP_SIZE0];
                     44: double _ep1[EXP_SIZE1];
                     45: double _ep2[EXP_SIZE2];
                     46: double _ep3[EXP_SIZE3];
                     47: double _ep4[EXP_SIZE4];
                     48: double _en0[EXP_SIZE0];
                     49: double _en1[EXP_SIZE1];
                     50: double _en2[EXP_SIZE2];
                     51: double _en3[EXP_SIZE3];
                     52: double _en4[EXP_SIZE4];
                     53: 
                     54: 
                     55: double
                     56: exp(arg)
                     57: double arg;
                     58: {
                     59:        int a0;
                     60:        union {
                     61:                int i;
                     62:                struct {
                     63:                        unsigned char b0;
                     64:                        unsigned char b1;
                     65:                        unsigned char b2;
                     66:                        unsigned char b3;
                     67:                } b;
                     68:        } u;
                     69:        register union {
                     70:                double d;
                     71:                int i;
                     72:        } abs_arg;
                     73: 
                     74:        if (arg == 0.0)
                     75:                return(1.0);
                     76: 
                     77:        if (arg >= 0.0) {
                     78:                abs_arg.d = arg;
                     79:                if (abs_arg.i >= 0x42000000 ) { /* if (abs_arg.i >= 8.0) */
                     80:                        if (abs_arg.d > maxf) {
                     81:                                errno = ERANGE;
                     82:                                return(HUGE);
                     83:                        }
                     84:                        a0 = abs_arg.d;
                     85:                        return(exp(arg - (a0&0x78)) * _ep0[a0>>3]);
                     86:                }
                     87: 
                     88:                abs_arg.i += EXPONENT28;    /* multiply by 2 to 28th power */
                     89:                u.i = abs_arg.d;
                     90:                abs_arg.d = u.i;
                     91:                abs_arg.i -= EXPONENT28;    /* divide by 2 to 28th power */
                     92:                return(((arg-abs_arg.d)+1.)*
                     93:                        _ep1[u.b.b0]*_ep2[u.b.b1]*_ep3[u.b.b2]*_ep4[u.b.b3]);
                     94:        } else {
                     95:                abs_arg.d = -arg;
                     96:                if (abs_arg.i >= 0x42000000) {
                     97:                        if (arg < minf)
                     98:                                return(0.);
                     99:                        a0 = -arg;
                    100:                        return(exp(arg + (a0&0x78)) * _en0[a0>>3]);
                    101:                }
                    102: 
                    103:                abs_arg.i += EXPONENT28;    /* multiply by 2 to 28th power */
                    104:                u.i = abs_arg.d;
                    105:                abs_arg.d = u.i;
                    106:                abs_arg.i -= EXPONENT28;    /* divide by 2 to 28th power */
                    107:                return(((arg+abs_arg.d)+1.)*
                    108:                      _en1[u.b.b0]*_en2[u.b.b1]*_en3[u.b.b2]*_en4[u.b.b3]);
                    109:        }
                    110: }

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