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
|
/* mconf.h
*
* Common include file for math routines
*
*
*
* SYNOPSIS:
*
* #include "mconf.h"
*
*
*
* DESCRIPTION:
*
* This file contains definitions for error codes that are
* passed to the common error handling routine mtherr()
* (which see).
*
* The file also includes a conditional assembly definition
* for the type of computer arithmetic (IEEE, DEC, Motorola
* IEEE, or UNKnown).
*
* For Digital Equipment PDP-11 and VAX computers, certain
* IBM systems, and others that use numbers with a 56-bit
* significand, the symbol DEC should be defined. In this
* mode, most floating point constants are given as arrays
* of octal integers to eliminate decimal to binary conversion
* errors that might be introduced by the compiler.
*
* For little-endian computers, such as IBM PC, that follow the
* IEEE Standard for Binary Floating Point Arithmetic (ANSI/IEEE
* Std 754-1985), the symbol IBMPC should be defined. These
* numbers have 53-bit significands. In this mode, constants
* are provided as arrays of hexadecimal 16 bit integers.
*
* Big-endian IEEE format is denoted MIEEE. On some RISC
* systems such as Sun SPARC, double precision constants
* must be stored on 8-byte address boundaries. Since integer
* arrays may be aligned differently, the MIEEE configuration
* may fail on such machines.
*
* To accommodate other types of computer arithmetic, all
* constants are also provided in a normal decimal radix
* which one can hope are correctly converted to a suitable
* format by the available C language compiler. To invoke
* this mode, define the symbol UNK.
*
* An important difference among these modes is a predefined
* set of machine arithmetic constants for each. The numbers
* MACHEP (the machine roundoff error), MAXNUM (largest number
* represented), and several other parameters are preset by
* the configuration symbol. Check the file const.c to
* ensure that these values are correct for your computer.
*
* Configurations NANS, INFINITIES, MINUSZERO, and DENORMAL
* may fail on many systems. Verify that they are supposed
* to work on your computer.
*/
/*
Cephes Math Library Release 2.3: June, 1995
Copyright 1984, 1987, 1989, 1995 by Stephen L. Moshier
*/
#include "cephes_names.h"
/* Constant definitions for math error conditions
*/
#define DOMAIN 1 /* argument domain error */
#define SING 2 /* argument singularity */
#define OVERFLOW 3 /* overflow range error */
#define UNDERFLOW 4 /* underflow range error */
#define TLOSS 5 /* total loss of precision */
#define PLOSS 6 /* partial loss of precision */
#define TOOMANY 7 /* too many iterations */
#define MAXITER 500
#define EDOM 33
#define ERANGE 34
/* Complex numeral. */
typedef struct
{
double r;
double i;
} cmplx;
/* Long double complex numeral. */
/*
typedef struct
{
long double r;
long double i;
} cmplxl;
*/
/* Type of computer arithmetic */
/* This is kind of improper, as the byte-order of floats may not
* be the same as the byte-order of ints. However, it works.
*/
#include <pyconfig.h>
#ifdef WORDS_BIGENDIAN
# define MIEEE 1
# define BIGENDIAN 1
#else
# define IBMPC 1
# define BIGENDIAN 0
#endif
/* UNKnown arithmetic, invokes coefficients given in
* normal decimal format. Beware of range boundary
* problems (MACHEP, MAXLOG, etc. in const.c) and
* roundoff problems in pow.c:
* (Sun SPARCstation)
*/
/* #define UNK 1 */
/* Define this `volatile' if your compiler thinks
* that floating point arithmetic obeys the associative
* and distributive laws. It will defeat some optimizations
* (but probably not enough of them).
*
* #define VOLATILE volatile
*/
#define VOLATILE
/* For 12-byte long doubles on an i386, pad a 16-bit short 0
* to the end of real constants initialized by integer arrays.
*
* #define XPD 0,
*
* Otherwise, the type is 10 bytes long and XPD should be
* defined blank (e.g., Microsoft C).
*
* #define XPD
*/
#define XPD 0,
/* Define to support tiny denormal numbers, else undefine. */
#define DENORMAL 1
/* Define to ask for infinity support, else undefine. */
#define INFINITIES 1
#ifdef NOINFINITIES
#undef INFINITIES
#endif
/* Define to ask for support of numbers that are Not-a-Number,
else undefine. This may automatically define INFINITIES in some files. */
#define NANS 1
#ifdef NONANS
#undef NANS
#endif
/* Define to distinguish between -0.0 and +0.0. */
#define MINUSZERO 1
/* Define 1 for ANSI C atan2() function
See atan.c and clog.c. */
#define ANSIC 1
/* Get ANSI function prototypes, if you want them. */
#if defined(__STDC__) || defined(_MSC_EXTENSIONS)
#define ANSIPROT
#include "protos.h"
#else
int mtherr();
#endif
/* Variable for error reporting. See mtherr.c. */
extern int merror;
#define gamma Gamma
|
