|
|
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.