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
23typedef union uLD { long double ld; unsigned short sh[8]; long lo[4]; } uLD;
24typedef union uD { double d; unsigned short sh[4]; } uD;
25#elif IBMPC
26typedef union uLD { unsigned short sh[8]; long double ld; long lo[4]; } uLD;
27typedef union uD { unsigned short sh[4]; double d; } uD;
28#elif MIEEE
29typedef union uLD { long lo[4]; long double ld; unsigned short sh[8]; } uLD;
30typedef 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
62extern double __INF;
63#define INFINITY (__INF)
64extern double __QNAN;
65#define NAN (__QNAN)
66#endif
67
68/*long double*/
69#if __SIZEOF_LONG_DOUBLE__ == __SIZEOF_DOUBLE__
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 /* __SIZEOF_LONG_DOUBLE__ == __SIZEOF_DOUBLE__ */
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
102extern long double __INFL;
103#define INFINITYL (__INFL)
104extern 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
132extern float __INFF;
133#define INFINITYF (__INFF)
134extern float __QNANF;
135#define NANF (__QNANF)
136#endif
137
138
139/* double */
140
141/*
142Cephes Math Library Release 2.2: July, 1992
143Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
144Direct 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 */
195static __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 */
214static __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/*
231Cephes Math Library Release 2.2: July, 1992
232Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
233Direct 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 */
284static __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 */
304static __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/*
368Cephes Math Library Release 2.1: December, 1988
369Copyright 1984, 1987, 1988 by Stephen L. Moshier
370Direct inquiries to 30 Frost Street, Cambridge, MA 02140
371*/
372
373static __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
401static __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