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
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
|
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* from: @(#)fdlibm.h 5.1 93/09/24
*/
#ifndef _MATH_PRIVATE_H_
#define _MATH_PRIVATE_H_
#include <endian.h>
#include <stdbool.h>
#include <stdint.h>
#include <sys/types.h>
/* Gather machine dependent _Floatn support. */
#include <bits/floatn.h>
/* The original fdlibm code used statements like:
n0 = ((*(int*)&one)>>29)^1; * index of high word *
ix0 = *(n0+(int*)&x); * high word of x *
ix1 = *((1-n0)+(int*)&x); * low word of x *
to dig two 32 bit words out of the 64 bit IEEE floating point
value. That is non-ANSI, and, moreover, the gcc instruction
scheduler gets it wrong. We instead use the following macros.
Unlike the original code, we determine the endianness at compile
time, not at run time; I don't see much benefit to selecting
endianness at run time. */
/* A union which permits us to convert between a double and two 32 bit
ints. */
#if __FLOAT_WORD_ORDER == __BIG_ENDIAN
typedef union
{
double value;
struct
{
uint32_t msw;
uint32_t lsw;
} parts;
uint64_t word;
} ieee_double_shape_type;
#endif
#if __FLOAT_WORD_ORDER == __LITTLE_ENDIAN
typedef union
{
double value;
struct
{
uint32_t lsw;
uint32_t msw;
} parts;
uint64_t word;
} ieee_double_shape_type;
#endif
/* Get two 32 bit ints from a double. */
#define EXTRACT_WORDS(ix0,ix1,d) \
do { \
ieee_double_shape_type ew_u; \
ew_u.value = (d); \
(ix0) = ew_u.parts.msw; \
(ix1) = ew_u.parts.lsw; \
} while (0)
/* Get the more significant 32 bit int from a double. */
#ifndef GET_HIGH_WORD
# define GET_HIGH_WORD(i,d) \
do { \
ieee_double_shape_type gh_u; \
gh_u.value = (d); \
(i) = gh_u.parts.msw; \
} while (0)
#endif
/* Get the less significant 32 bit int from a double. */
#ifndef GET_LOW_WORD
# define GET_LOW_WORD(i,d) \
do { \
ieee_double_shape_type gl_u; \
gl_u.value = (d); \
(i) = gl_u.parts.lsw; \
} while (0)
#endif
/* Get all in one, efficient on 64-bit machines. */
#ifndef EXTRACT_WORDS64
# define EXTRACT_WORDS64(i,d) \
do { \
ieee_double_shape_type gh_u; \
gh_u.value = (d); \
(i) = gh_u.word; \
} while (0)
#endif
/* Set a double from two 32 bit ints. */
#ifndef INSERT_WORDS
# define INSERT_WORDS(d,ix0,ix1) \
do { \
ieee_double_shape_type iw_u; \
iw_u.parts.msw = (ix0); \
iw_u.parts.lsw = (ix1); \
(d) = iw_u.value; \
} while (0)
#endif
/* Get all in one, efficient on 64-bit machines. */
#ifndef INSERT_WORDS64
# define INSERT_WORDS64(d,i) \
do { \
ieee_double_shape_type iw_u; \
iw_u.word = (i); \
(d) = iw_u.value; \
} while (0)
#endif
/* Set the more significant 32 bits of a double from an int. */
#ifndef SET_HIGH_WORD
#define SET_HIGH_WORD(d,v) \
do { \
ieee_double_shape_type sh_u; \
sh_u.value = (d); \
sh_u.parts.msw = (v); \
(d) = sh_u.value; \
} while (0)
#endif
/* Set the less significant 32 bits of a double from an int. */
#ifndef SET_LOW_WORD
# define SET_LOW_WORD(d,v) \
do { \
ieee_double_shape_type sl_u; \
sl_u.value = (d); \
sl_u.parts.lsw = (v); \
(d) = sl_u.value; \
} while (0)
#endif
/* A union which permits us to convert between a float and a 32 bit
int. */
typedef union
{
float value;
uint32_t word;
} ieee_float_shape_type;
/* Get a 32 bit int from a float. */
#ifndef GET_FLOAT_WORD
# define GET_FLOAT_WORD(i,d) \
do { \
ieee_float_shape_type gf_u; \
gf_u.value = (d); \
(i) = gf_u.word; \
} while (0)
#endif
/* Set a float from a 32 bit int. */
#ifndef SET_FLOAT_WORD
# define SET_FLOAT_WORD(d,i) \
do { \
ieee_float_shape_type sf_u; \
sf_u.word = (i); \
(d) = sf_u.value; \
} while (0)
#endif
/* We need to guarantee an expansion of name when building
ldbl-128 files as another type (e.g _Float128). */
#define mathx_hidden_def(name) hidden_def(name)
/* Get long double macros from a separate header. */
#include <math_ldbl.h>
/* Include function declarations for each floating-point. */
#define _Mdouble_ double
#define _MSUF_
#include <math_private_calls.h>
#undef _MSUF_
#undef _Mdouble_
#define _Mdouble_ float
#define _MSUF_ f
#define __MATH_DECLARING_FLOAT
#include <math_private_calls.h>
#undef __MATH_DECLARING_FLOAT
#undef _MSUF_
#undef _Mdouble_
#define _Mdouble_ long double
#define _MSUF_ l
#define __MATH_DECLARING_LONG_DOUBLE
#include <math_private_calls.h>
#undef __MATH_DECLARING_LONG_DOUBLE
#undef _MSUF_
#undef _Mdouble_
#if __HAVE_DISTINCT_FLOAT128
# define _Mdouble_ _Float128
# define _MSUF_ f128
# define __MATH_DECLARING_FLOATN
# include <math_private_calls.h>
# undef __MATH_DECLARING_FLOATN
# undef _MSUF_
# undef _Mdouble_
#endif
#if __HAVE_DISTINCT_FLOAT128
/* __builtin_isinf_sign is broken in GCC < 7 for float128. */
# if ! __GNUC_PREREQ (7, 0)
# include <ieee754_float128.h>
extern inline int
__isinff128 (_Float128 x)
{
int64_t hx, lx;
GET_FLOAT128_WORDS64 (hx, lx, x);
lx |= (hx & 0x7fffffffffffffffLL) ^ 0x7fff000000000000LL;
lx |= -lx;
return ~(lx >> 63) & (hx >> 62);
}
# endif
extern inline _Float128
fabsf128 (_Float128 x)
{
return __builtin_fabsf128 (x);
}
#endif
/* Prototypes for functions of the IBM Accurate Mathematical Library. */
extern double __exp1 (double __x, double __xx);
extern double __sin (double __x);
extern double __cos (double __x);
extern int __branred (double __x, double *__a, double *__aa);
extern void __doasin (double __x, double __dx, double __v[]);
extern void __dubsin (double __x, double __dx, double __v[]);
extern void __dubcos (double __x, double __dx, double __v[]);
extern double __sin32 (double __x, double __res, double __res1);
extern double __cos32 (double __x, double __res, double __res1);
extern double __mpsin (double __x, double __dx, bool __range_reduce);
extern double __mpcos (double __x, double __dx, bool __range_reduce);
extern void __docos (double __x, double __dx, double __v[]);
#endif /* _MATH_PRIVATE_H_ */
|
