1 #include <math.h>
2 #include <errno.h>
3
4
5 #define IBMPC 1
6 #define ANSIPROT 1
7 #define MINUSZERO 1
8 #define INFINITIES 1
9 #define NANS 1
10 #define DENORMAL 1
11 #define VOLATILE
12 #define mtherr(fname, code)
13 #define XPD 0,
14 #ifdef __x86_64__
15 #define XPD_SHORT 0, 0,
16 #define XPD_LONG 0,
17 #else
18 #define XPD_SHORT
19 #define XPD_LONG
20 #endif
21
22 #if UNK
23 typedef union uLD { long double ld; unsigned short sh[8]; long lo[4]; } uLD;
24 typedef union uD { double d; unsigned short sh[4]; } uD;
25 #elif IBMPC
26 typedef union uLD { unsigned short sh[8]; long double ld; long lo[4]; } uLD;
27 typedef union uD { unsigned short sh[4]; double d; } uD;
28 #elif MIEEE
29 typedef union uLD { long lo[4]; long double ld; unsigned short sh[8]; } uLD;
30 typedef union uD { unsigned short sh[4]; double d; } uD;
31 #else
32 #error Unknown uLD/uD type definition
33 #endif
34
35 #define _CEPHES_USE_ERRNO
36
37 #ifdef _CEPHES_USE_ERRNO
38 #define _SET_ERRNO(x) errno = (x)
39 #else
40 #define _SET_ERRNO(x)
41 #endif
42
43 /* constants used by cephes functions */
44
45 /* double */
46 #define MAXNUM 1.7976931348623158E308
47 #define MAXLOG 7.09782712893383996843E2
48 #define MINLOG -7.08396418532264106224E2
49 #define LOGE2 6.93147180559945309417E-1
50 #define LOG2E 1.44269504088896340736
51 #define PI 3.14159265358979323846
52 #define PIO2 1.57079632679489661923
53 #define PIO4 7.85398163397448309616E-1
54
55 #define NEGZERO (-0.0)
56 #undef NAN
57 #undef INFINITY
58 #if (__GNUC__ > 3 || (__GNUC__ == 3 && __GNUC_MINOR__ > 2))
59 #define INFINITY __builtin_huge_val()
60 #define NAN __builtin_nan("")
61 #else
62 extern double __INF;
63 #define INFINITY (__INF)
64 extern double __QNAN;
65 #define NAN (__QNAN)
66 #endif
67
68 /*long double*/
69 #if defined(__arm__) || defined(_ARM_) || defined(__aarch64__) || defined(_ARM64_)
70 #define MAXNUML 1.7976931348623158E308
71 #define MAXLOGL 7.09782712893383996843E2
72 #define MINLOGL -7.08396418532264106224E2
73 #define LOGE2L 6.93147180559945309417E-1
74 #define LOG2EL 1.44269504088896340736
75 #define PIL 3.14159265358979323846
76 #define PIO2L 1.57079632679489661923
77 #define PIO4L 7.85398163397448309616E-1
78 #else
79 #define MAXNUML 1.189731495357231765021263853E4932L
80 #define MAXLOGL 1.1356523406294143949492E4L
81 #define MINLOGL -1.13994985314888605586758E4L
82 #define LOGE2L 6.9314718055994530941723E-1L
83 #define LOG2EL 1.4426950408889634073599E0L
84 #define PIL 3.1415926535897932384626L
85 #define PIO2L 1.5707963267948966192313L
86 #define PIO4L 7.8539816339744830961566E-1L
87 #endif /* defined(__arm__) || defined(_ARM_) || defined(__aarch64__) || defined(_ARM64_) */
88
89 #define isfinitel isfinite
90 #define isinfl isinf
91 #define isnanl isnan
92 #define signbitl signbit
93
94 #define NEGZEROL (-0.0L)
95
96 #undef NANL
97 #undef INFINITYL
98 #if (__GNUC__ > 3 || (__GNUC__ == 3 && __GNUC_MINOR__ > 2))
99 #define INFINITYL __builtin_huge_vall()
100 #define NANL __builtin_nanl("")
101 #else
102 extern long double __INFL;
103 #define INFINITYL (__INFL)
104 extern long double __QNANL;
105 #define NANL (__QNANL)
106 #endif
107
108 /* float */
109
110 #define MAXNUMF 3.4028234663852885981170418348451692544e38F
111 #define MAXLOGF 88.72283905206835F
112 #define MINLOGF -103.278929903431851103F /* log(2^-149) */
113 #define LOG2EF 1.44269504088896341F
114 #define LOGE2F 0.693147180559945309F
115 #define PIF 3.141592653589793238F
116 #define PIO2F 1.5707963267948966192F
117 #define PIO4F 0.7853981633974483096F
118
119 #define isfinitef isfinite
120 #define isinff isinf
121 #define isnanf isnan
122 #define signbitf signbit
123
124 #define NEGZEROF (-0.0F)
125
126 #undef NANF
127 #undef INFINITYF
128 #if (__GNUC__ > 3 || (__GNUC__ == 3 && __GNUC_MINOR__ > 2))
129 #define INFINITYF __builtin_huge_valf()
130 #define NANF __builtin_nanf("")
131 #else
132 extern float __INFF;
133 #define INFINITYF (__INFF)
134 extern float __QNANF;
135 #define NANF (__QNANF)
136 #endif
137
138
139 /* double */
140
141 /*
142 Cephes Math Library Release 2.2: July, 1992
143 Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
144 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
145 */
146
147
148 /* polevl.c
149 * p1evl.c
150 *
151 * Evaluate polynomial
152 *
153 *
154 *
155 * SYNOPSIS:
156 *
157 * int N;
158 * double x, y, coef[N+1], polevl[];
159 *
160 * y = polevl( x, coef, N );
161 *
162 *
163 *
164 * DESCRIPTION:
165 *
166 * Evaluates polynomial of degree N:
167 *
168 * 2 N
169 * y = C + C x + C x +...+ C x
170 * 0 1 2 N
171 *
172 * Coefficients are stored in reverse order:
173 *
174 * coef[0] = C , ..., coef[N] = C .
175 * N 0
176 *
177 * The function p1evl() assumes that coef[N] = 1.0 and is
178 * omitted from the array. Its calling arguments are
179 * otherwise the same as polevl().
180 *
181 *
182 * SPEED:
183 *
184 * In the interest of speed, there are no checks for out
185 * of bounds arithmetic. This routine is used by most of
186 * the functions in the library. Depending on available
187 * equipment features, the user may wish to rewrite the
188 * program in microcode or assembly language.
189 *
190 */
191
192 /* Polynomial evaluator:
193 * P[0] x^n + P[1] x^(n-1) + ... + P[n]
194 */
polevl(double x,const uD * p,int n)195 static __inline__ double polevl(double x, const uD *p, int n)
196 {
197 register double y;
198
199 y = p->d;
200 p++;
201 do
202 {
203 y = y * x + p->d;
204 p++;
205 }
206 while (--n);
207 return (y);
208 }
209
210
211 /* Polynomial evaluator:
212 * x^n + P[0] x^(n-1) + P[1] x^(n-2) + ... + P[n]
213 */
p1evl(double x,const uD * p,int n)214 static __inline__ double p1evl(double x, const uD *p, int n)
215 {
216 register double y;
217
218 n -= 1;
219 y = x + p->d; p++;
220 do
221 {
222 y = y * x + p->d; p++;
223 }
224 while (--n);
225 return (y);
226 }
227
228
229 /* long double */
230 /*
231 Cephes Math Library Release 2.2: July, 1992
232 Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
233 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
234 */
235
236
237 /* polevll.c
238 * p1evll.c
239 *
240 * Evaluate polynomial
241 *
242 *
243 *
244 * SYNOPSIS:
245 *
246 * int N;
247 * long double x, y, coef[N+1], polevl[];
248 *
249 * y = polevll( x, coef, N );
250 *
251 *
252 *
253 * DESCRIPTION:
254 *
255 * Evaluates polynomial of degree N:
256 *
257 * 2 N
258 * y = C + C x + C x +...+ C x
259 * 0 1 2 N
260 *
261 * Coefficients are stored in reverse order:
262 *
263 * coef[0] = C , ..., coef[N] = C .
264 * N 0
265 *
266 * The function p1evll() assumes that coef[N] = 1.0 and is
267 * omitted from the array. Its calling arguments are
268 * otherwise the same as polevll().
269 *
270 *
271 * SPEED:
272 *
273 * In the interest of speed, there are no checks for out
274 * of bounds arithmetic. This routine is used by most of
275 * the functions in the library. Depending on available
276 * equipment features, the user may wish to rewrite the
277 * program in microcode or assembly language.
278 *
279 */
280
281 /* Polynomial evaluator:
282 * P[0] x^n + P[1] x^(n-1) + ... + P[n]
283 */
polevll(long double x,const uLD * p,int n)284 static __inline__ long double polevll(long double x, const uLD *p, int n)
285 {
286 register long double y;
287
288 y = p->ld;
289 p++;
290 do
291 {
292 y = y * x + p->ld;
293 p++;
294 }
295 while (--n);
296 return y;
297 }
298
299
300
301 /* Polynomial evaluator:
302 * x^n + P[0] x^(n-1) + P[1] x^(n-2) + ... + P[n]
303 */
p1evll(long double x,const uLD * p,int n)304 static __inline__ long double p1evll(long double x, const uLD *p, int n)
305 {
306 register long double y;
307
308 n -= 1;
309 y = x + p->ld;
310 p++;
311
312 do
313 {
314 y = y * x + p->ld;
315 p++;
316 }
317 while (--n);
318 return (y);
319 }
320
321 /* Float version */
322
323 /* polevlf.c
324 * p1evlf.c
325 *
326 * Evaluate polynomial
327 *
328 *
329 *
330 * SYNOPSIS:
331 *
332 * int N;
333 * float x, y, coef[N+1], polevlf[];
334 *
335 * y = polevlf( x, coef, N );
336 *
337 *
338 *
339 * DESCRIPTION:
340 *
341 * Evaluates polynomial of degree N:
342 *
343 * 2 N
344 * y = C + C x + C x +...+ C x
345 * 0 1 2 N
346 *
347 * Coefficients are stored in reverse order:
348 *
349 * coef[0] = C , ..., coef[N] = C .
350 * N 0
351 *
352 * The function p1evl() assumes that coef[N] = 1.0 and is
353 * omitted from the array. Its calling arguments are
354 * otherwise the same as polevl().
355 *
356 *
357 * SPEED:
358 *
359 * In the interest of speed, there are no checks for out
360 * of bounds arithmetic. This routine is used by most of
361 * the functions in the library. Depending on available
362 * equipment features, the user may wish to rewrite the
363 * program in microcode or assembly language.
364 *
365 */
366
367 /*
368 Cephes Math Library Release 2.1: December, 1988
369 Copyright 1984, 1987, 1988 by Stephen L. Moshier
370 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
371 */
372
polevlf(float x,const float * coef,int N)373 static __inline__ float polevlf(float x, const float* coef, int N)
374 {
375 float ans;
376 float *p;
377 int i;
378
379 p = (float*)coef;
380 ans = *p++;
381
382 /*
383 for (i = 0; i < N; i++)
384 ans = ans * x + *p++;
385 */
386
387 i = N;
388 do
389 ans = ans * x + *p++;
390 while (--i);
391
392 return (ans);
393 }
394
395 /* p1evl() */
396 /* N
397 * Evaluate polynomial when coefficient of x is 1.0.
398 * Otherwise same as polevl.
399 */
400
p1evlf(float x,const float * coef,int N)401 static __inline__ float p1evlf(float x, const float *coef, int N)
402 {
403 float ans;
404 float *p;
405 int i;
406
407 p = (float*)coef;
408 ans = x + *p++;
409 i = N - 1;
410
411 do
412 ans = ans * x + *p++;
413 while (--i);
414
415 return (ans);
416 }
417
418