| 1 | /*
|
|---|
| 2 | * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
|
|---|
| 3 | * See the copyright notice in the ACK home directory, in the file "Copyright".
|
|---|
| 4 | *
|
|---|
| 5 | * Author: Ceriel J.H. Jacobs
|
|---|
| 6 | */
|
|---|
| 7 |
|
|---|
| 8 | /* $Header: /cvsup/minix/src/lib/ack/libp/exp.c,v 1.1 2005/10/10 15:27:46 beng Exp $ */
|
|---|
| 9 | #define __NO_DEFS
|
|---|
| 10 | #include <math.h>
|
|---|
| 11 | #include <pc_err.h>
|
|---|
| 12 | extern _trp();
|
|---|
| 13 |
|
|---|
| 14 | #if __STDC__
|
|---|
| 15 | #include <float.h>
|
|---|
| 16 | #include <pc_math.h>
|
|---|
| 17 | #define M_MIN_D DBL_MIN
|
|---|
| 18 | #define M_MAX_D DBL_MAX
|
|---|
| 19 | #define M_DMINEXP DBL_MIN_EXP
|
|---|
| 20 | #endif
|
|---|
| 21 | #undef HUGE
|
|---|
| 22 | #define HUGE 1e1000
|
|---|
| 23 |
|
|---|
| 24 | static double
|
|---|
| 25 | Ldexp(fl,exp)
|
|---|
| 26 | double fl;
|
|---|
| 27 | int exp;
|
|---|
| 28 | {
|
|---|
| 29 | extern double _fef();
|
|---|
| 30 | int sign = 1;
|
|---|
| 31 | int currexp;
|
|---|
| 32 |
|
|---|
| 33 | if (fl<0) {
|
|---|
| 34 | fl = -fl;
|
|---|
| 35 | sign = -1;
|
|---|
| 36 | }
|
|---|
| 37 | fl = _fef(fl,&currexp);
|
|---|
| 38 | exp += currexp;
|
|---|
| 39 | if (exp > 0) {
|
|---|
| 40 | while (exp>30) {
|
|---|
| 41 | fl *= (double) (1L << 30);
|
|---|
| 42 | exp -= 30;
|
|---|
| 43 | }
|
|---|
| 44 | fl *= (double) (1L << exp);
|
|---|
| 45 | }
|
|---|
| 46 | else {
|
|---|
| 47 | while (exp<-30) {
|
|---|
| 48 | fl /= (double) (1L << 30);
|
|---|
| 49 | exp += 30;
|
|---|
| 50 | }
|
|---|
| 51 | fl /= (double) (1L << -exp);
|
|---|
| 52 | }
|
|---|
| 53 | return sign * fl;
|
|---|
| 54 | }
|
|---|
| 55 |
|
|---|
| 56 | double
|
|---|
| 57 | _exp(x)
|
|---|
| 58 | double x;
|
|---|
| 59 | {
|
|---|
| 60 | /* Algorithm and coefficients from:
|
|---|
| 61 | "Software manual for the elementary functions"
|
|---|
| 62 | by W.J. Cody and W. Waite, Prentice-Hall, 1980
|
|---|
| 63 | */
|
|---|
| 64 |
|
|---|
| 65 | static double p[] = {
|
|---|
| 66 | 0.25000000000000000000e+0,
|
|---|
| 67 | 0.75753180159422776666e-2,
|
|---|
| 68 | 0.31555192765684646356e-4
|
|---|
| 69 | };
|
|---|
| 70 |
|
|---|
| 71 | static double q[] = {
|
|---|
| 72 | 0.50000000000000000000e+0,
|
|---|
| 73 | 0.56817302698551221787e-1,
|
|---|
| 74 | 0.63121894374398503557e-3,
|
|---|
| 75 | 0.75104028399870046114e-6
|
|---|
| 76 | };
|
|---|
| 77 | double xn, g;
|
|---|
| 78 | int n;
|
|---|
| 79 | int negative = x < 0;
|
|---|
| 80 |
|
|---|
| 81 | if (x <= M_LN_MIN_D) {
|
|---|
| 82 | g = M_MIN_D/4.0;
|
|---|
| 83 |
|
|---|
| 84 | if (g != 0.0) {
|
|---|
| 85 | /* unnormalized numbers apparently exist */
|
|---|
| 86 | if (x < (M_LN2 * (M_DMINEXP - 53))) return 0.0;
|
|---|
| 87 | }
|
|---|
| 88 | else {
|
|---|
| 89 | if (x < M_LN_MIN_D) return 0.0;
|
|---|
| 90 | return M_MIN_D;
|
|---|
| 91 | }
|
|---|
| 92 | }
|
|---|
| 93 | if (x >= M_LN_MAX_D) {
|
|---|
| 94 | if (x > M_LN_MAX_D) {
|
|---|
| 95 | _trp(EEXP);
|
|---|
| 96 | return HUGE;
|
|---|
| 97 | }
|
|---|
| 98 | return M_MAX_D;
|
|---|
| 99 | }
|
|---|
| 100 | if (negative) x = -x;
|
|---|
| 101 |
|
|---|
| 102 | n = x * M_LOG2E + 0.5; /* 1/ln(2) = log2(e), 0.5 added for rounding */
|
|---|
| 103 | xn = n;
|
|---|
| 104 | {
|
|---|
| 105 | double x1 = (long) x;
|
|---|
| 106 | double x2 = x - x1;
|
|---|
| 107 |
|
|---|
| 108 | g = ((x1-xn*0.693359375)+x2) - xn*(-2.1219444005469058277e-4);
|
|---|
| 109 | }
|
|---|
| 110 | if (negative) {
|
|---|
| 111 | g = -g;
|
|---|
| 112 | n = -n;
|
|---|
| 113 | }
|
|---|
| 114 | xn = g * g;
|
|---|
| 115 | x = g * POLYNOM2(xn, p);
|
|---|
| 116 | n += 1;
|
|---|
| 117 | return (Ldexp(0.5 + x/(POLYNOM3(xn, q) - x), n));
|
|---|
| 118 | }
|
|---|
