1 /* md_exp2.c
2 *
3 * Base 2 exponential function
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * double x, y, md_exp2();
10 *
11 * y = md_exp2( x );
12 *
13 *
14 *
15 * DESCRIPTION:
16 *
17 * Returns 2 raised to the x power.
18 *
19 * Range reduction is accomplished by separating the argument
20 * into an integer k and fraction f such that
21 * x k f
22 * 2 = 2 2.
23 *
24 * A Pade' form
25 *
26 * 1 + 2x P(x**2) / (Q(x**2) - x P(x**2) )
27 *
28 * approximates 2**x in the basic range [-0.5, 0.5].
29 *
30 *
31 * ACCURACY:
32 *
33 * Relative error:
34 * arithmetic domain # trials peak rms
35 * IEEE -1022,+1024 30000 1.8e-16 5.4e-17
36 *
37 *
38 * See md_exp.c for comments on error amplification.
39 *
40 *
41 * ERROR MESSAGES:
42 *
43 * message condition value returned
44 * md_exp underflow x < -MAXL2 0.0
45 * md_exp overflow x > MAXL2 MAXNUM
46 *
47 * For DEC arithmetic, MAXL2 = 127.
48 * For IEEE arithmetic, MAXL2 = 1024.
49 */
50
51
52 /*
53 Cephes Math Library Release 2.8: June, 2000
54 Copyright 1984, 1995, 2000 by Stephen L. Moshier
55 */
56
57
58
59 #include "mconf.h"
60
61 #ifdef UNK
62 static double P[] = {
63 2.30933477057345225087E-2,
64 2.02020656693165307700E1,
65 1.51390680115615096133E3,
66 };
67 static double Q[] = {
68 /* 1.00000000000000000000E0,*/
69 2.33184211722314911771E2,
70 4.36821166879210612817E3,
71 };
72 #define MAXL2 1024.0
73 #define MINL2 -1024.0
74 #endif
75
76 #ifdef DEC
77 static unsigned short P[] = {
78 0036675,0027102,0122327,0053227,
79 0041241,0116724,0115412,0157355,
80 0042675,0036404,0101733,0132226,
81 };
82 static unsigned short Q[] = {
83 /*0040200,0000000,0000000,0000000,*/
84 0042151,0027450,0077732,0160744,
85 0043210,0100661,0077550,0056560,
86 };
87 #define MAXL2 127.0
88 #define MINL2 -127.0
89 #endif
90
91 #ifdef IBMPC
92 static unsigned short P[] = {
93 0xead3,0x549a,0xa5c8,0x3f97,
94 0x5bde,0x9361,0x33ba,0x4034,
95 0x7693,0x907b,0xa7a0,0x4097,
96 };
97 static unsigned short Q[] = {
98 /*0x0000,0x0000,0x0000,0x3ff0,*/
99 0x5c3c,0x0ffb,0x25e5,0x406d,
100 0x0bae,0x2fed,0x1036,0x40b1,
101 };
102 #define MAXL2 1024.0
103 #define MINL2 -1022.0
104 #endif
105
106 #ifdef MIEEE
107 static unsigned short P[] = {
108 0x3f97,0xa5c8,0x549a,0xead3,
109 0x4034,0x33ba,0x9361,0x5bde,
110 0x4097,0xa7a0,0x907b,0x7693,
111 };
112 static unsigned short Q[] = {
113 /*0x3ff0,0x0000,0x0000,0x0000,*/
114 0x406d,0x25e5,0x0ffb,0x5c3c,
115 0x40b1,0x1036,0x2fed,0x0bae,
116 };
117 #define MAXL2 1024.0
118 #define MINL2 -1022.0
119 #endif
120
121 #ifdef ANSIPROT
122 extern double polevl ( double, void *, int );
123 extern double p1evl ( double, void *, int );
124 extern double md_floor ( double );
125 extern double md_ldexp ( double, int );
126 extern int isnan ( double );
127 extern int isfinite ( double );
128 #else
129 double polevl(), p1evl(), md_floor(), md_ldexp();
130 int isnan(), isfinite();
131 #endif
132 #ifdef INFINITIES
133 extern double INFINITY;
134 #endif
135 extern double MAXNUM;
136
md_exp2(x)137 double md_exp2(x)
138 double x;
139 {
140 double px, xx;
141 short n;
142
143 #ifdef NANS
144 if( isnan(x) )
145 return(x);
146 #endif
147 if( x > MAXL2)
148 {
149 #ifdef INFINITIES
150 return( INFINITY );
151 #else
152 mtherr( "md_exp2", OVERFLOW );
153 return( MAXNUM );
154 #endif
155 }
156
157 if( x < MINL2 )
158 {
159 #ifndef INFINITIES
160 mtherr( "md_exp2", UNDERFLOW );
161 #endif
162 return(0.0);
163 }
164
165 xx = x; /* save x */
166 /* separate into integer and fractional parts */
167 px = md_floor(x+0.5);
168 n = px;
169 x = x - px;
170
171 /* rational approximation
172 * md_exp2(x) = 1 + 2xP(xx)/(Q(xx) - P(xx))
173 * where xx = x**2
174 */
175 xx = x * x;
176 px = x * polevl( xx, P, 2 );
177 x = px / ( p1evl( xx, Q, 2 ) - px );
178 x = 1.0 + md_ldexp( x, 1 );
179
180 /* scale by power of 2 */
181 x = md_ldexp( x, n );
182 return(x);
183 }
184