1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
|
/* exp10l.c
*
* Base 10 exponential function, long double precision
* (Common antilogarithm)
*
*
*
* SYNOPSIS:
*
* long double x, y, exp10l()
*
* y = exp10l( x );
*
*
*
* DESCRIPTION:
*
* Returns 10 raised to the x power.
*
* Range reduction is accomplished by expressing the argument
* as 10**x = 2**n 10**f, with |f| < 0.5 log10(2).
* The Pade' form
*
* 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
*
* is used to approximate 10**f.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE +-4900 30000 1.0e-19 2.7e-20
*
* ERROR MESSAGES:
*
* message condition value returned
* exp10l underflow x < -MAXL10 0.0
* exp10l overflow x > MAXL10 MAXNUM
*
* IEEE arithmetic: MAXL10 = 4932.0754489586679023819
*
*/
/*
Cephes Math Library Release 2.2: January, 1991
Copyright 1984, 1991 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
#include <math.h>
#ifdef UNK
static long double P[] = {
3.1341179396892496811523E1L,
4.5618283154904699073999E3L,
1.3433113468542797218610E5L,
7.6025447914440301593592E5L,
};
static long double Q[] = {
/* 1.0000000000000000000000E0,*/
4.7705440288425157637739E2L,
2.9732606548049614870598E4L,
4.0843697951001026189583E5L,
6.6034865026929015925608E5L,
};
/*static long double LOG102 = 3.0102999566398119521373889e-1L;*/
static long double LOG210 = 3.3219280948873623478703L;
static long double LG102A = 3.01025390625e-1L;
static long double LG102B = 4.6050389811952137388947e-6L;
#endif
#ifdef IBMPC
static short P[] = {
0x399a,0x7dc7,0xbc43,0xfaba,0x4003, XPD
0xb526,0xdf32,0xa063,0x8e8e,0x400b, XPD
0x18da,0xafa1,0xc89e,0x832e,0x4010, XPD
0x503d,0x9352,0xe7aa,0xb99b,0x4012, XPD
};
static short Q[] = {
/*0x0000,0x0000,0x0000,0x8000,0x3fff,*/
0x947d,0x7855,0xf6ac,0xee86,0x4007, XPD
0x18cf,0x7749,0x368d,0xe849,0x400d, XPD
0x85be,0x2560,0x9f58,0xc76e,0x4011, XPD
0x6d3c,0x80c5,0xca67,0xa137,0x4012, XPD
};
/*
static short L102[] = {0xf799,0xfbcf,0x9a84,0x9a20,0x3ffd, XPD};
#define LOG102 *(long double *)L102
*/
static short L210[] = {0x8afe,0xcd1b,0x784b,0xd49a,0x4000, XPD};
#define LOG210 *(long double *)L210
static short L102A[] = {0x0000,0x0000,0x0000,0x9a20,0x3ffd, XPD};
#define LG102A *(long double *)L102A
static short L102B[] = {0x8f89,0xf798,0xfbcf,0x9a84,0x3fed, XPD};
#define LG102B *(long double *)L102B
#endif
#ifdef MIEEE
static long P[] = {
0x40030000,0xfababc43,0x7dc7399a,
0x400b0000,0x8e8ea063,0xdf32b526,
0x40100000,0x832ec89e,0xafa118da,
0x40120000,0xb99be7aa,0x9352503d,
};
static long Q[] = {
/* 0x3fff0000,0x80000000,0x00000000, */
0x40070000,0xee86f6ac,0x7855947d,
0x400d0000,0xe849368d,0x774918cf,
0x40110000,0xc76e9f58,0x256085be,
0x40120000,0xa137ca67,0x80c56d3c,
};
/*
static long L102[] = {0x3ffd0000,0x9a209a84,0xfbcff799};
#define LOG102 *(long double *)L102
*/
static long L210[] = {0x40000000,0xd49a784b,0xcd1b8afe};
#define LOG210 *(long double *)L210
static long L102A[] = {0x3ffd0000,0x9a200000,0x00000000};
#define LG102A *(long double *)L102A
static long L102B[] = {0x3fed0000,0x9a84fbcf,0xf7988f89};
#define LG102B *(long double *)L102B
#endif
static long double MAXL10 = 4.9320754489586679023819e3L;
extern long double MAXNUML;
#ifdef ANSIPROT
extern long double floorl ( long double );
extern long double ldexpl ( long double, int );
extern long double polevll ( long double, void *, int );
extern long double p1evll ( long double, void *, int );
extern int isnanl ( long double );
#else
long double floorl(), ldexpl(), polevll(), p1evll(), isnanl();
#endif
#ifdef INFINITIES
extern long double INFINITYL;
#endif
long double exp10l(x)
long double x;
{
long double px, xx;
short n;
#ifdef NANS
if( isnanl(x) )
return(x);
#endif
if( x > MAXL10 )
{
#ifdef INFINITIES
return( INFINITYL );
#else
mtherr( "exp10l", OVERFLOW );
return( MAXNUML );
#endif
}
if( x < -MAXL10 ) /* Would like to use MINLOG but can't */
{
#ifndef INFINITIES
mtherr( "exp10l", UNDERFLOW );
#endif
return(0.0L);
}
/* Express 10**x = 10**g 2**n
* = 10**g 10**( n log10(2) )
* = 10**( g + n log10(2) )
*/
px = floorl( LOG210 * x + 0.5L );
n = px;
x -= px * LG102A;
x -= px * LG102B;
/* rational approximation for exponential
* of the fractional part:
* 10**x = 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
*/
xx = x * x;
px = x * polevll( xx, P, 3 );
x = px/( p1evll( xx, Q, 4 ) - px );
x = 1.0L + ldexpl( x, 1 );
/* multiply by power of 2 */
x = ldexpl( x, n );
return(x);
}
|
