1 #ifndef _LIBM_H
2 #define _LIBM_H
3
4 #include <stdint.h>
5 #include <float.h>
6 #include <math.h>
7 #include <endian.h>
8 #include "../include/features.h"
9
10
11 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
12 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __LITTLE_ENDIAN
13 union ldshape {
14 long double f;
15 struct {
16 uint64_t m;
17 uint16_t se;
18 } i;
19 };
20 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __BIG_ENDIAN
21 /* This is the m68k variant of 80-bit long double, and this definition only works
22 * on archs where the alignment requirement of uint64_t is <= 4. */
23 union ldshape {
24 long double f;
25 struct {
26 uint16_t se;
27 uint16_t pad;
28 uint64_t m;
29 } i;
30 };
31 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __LITTLE_ENDIAN
32 union ldshape {
33 long double f;
34 struct {
35 uint64_t lo;
36 uint32_t mid;
37 uint16_t top;
38 uint16_t se;
39 } i;
40 struct {
41 uint64_t lo;
42 uint64_t hi;
43 } i2;
44 };
45 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __BIG_ENDIAN
46 union ldshape {
47 long double f;
48 struct {
49 uint16_t se;
50 uint16_t top;
51 uint32_t mid;
52 uint64_t lo;
53 } i;
54 struct {
55 uint64_t hi;
56 uint64_t lo;
57 } i2;
58 };
59 #else
60 #error Unsupported long double representation
61 #endif
62
63 /* Support non-nearest rounding mode. */
64 #define WANT_ROUNDING 1
65 /* Support signaling NaNs. */
66 #define WANT_SNAN 0
67
68 #if WANT_SNAN
69 #error SNaN is unsupported
70 #else
71 #define issignalingf_inline(x) 0
72 #define issignaling_inline(x) 0
73 #endif
74
75 #ifndef TOINT_INTRINSICS
76 #define TOINT_INTRINSICS 0
77 #endif
78
79 #if TOINT_INTRINSICS
80 /* Round x to nearest int in all rounding modes, ties have to be rounded
81 consistently with converttoint so the results match. If the result
82 would be outside of [-2^31, 2^31-1] then the semantics is unspecified. */
83 static double_t roundtoint(double_t);
84
85 /* Convert x to nearest int in all rounding modes, ties have to be rounded
86 consistently with roundtoint. If the result is not representible in an
87 int32_t then the semantics is unspecified. */
88 static int32_t converttoint(double_t);
89 #endif
90
91 /* Helps static branch prediction so hot path can be better optimized. */
92 #ifdef __GNUC__
93 #define predict_true(x) __builtin_expect(!!(x), 1)
94 #define predict_false(x) __builtin_expect(x, 0)
95 #else
96 #define predict_true(x) (x)
97 #define predict_false(x) (x)
98 #endif
99
100 /* Evaluate an expression as the specified type. With standard excess
101 precision handling a type cast or assignment is enough (with
102 -ffloat-store an assignment is required, in old compilers argument
103 passing and return statement may not drop excess precision). */
104
eval_as_float(float x)105 static inline float eval_as_float(float x)
106 {
107 float y = x;
108 return y;
109 }
110
eval_as_double(double x)111 static inline double eval_as_double(double x)
112 {
113 double y = x;
114 return y;
115 }
116
117 /* fp_barrier returns its input, but limits code transformations
118 as if it had a side-effect (e.g. observable io) and returned
119 an arbitrary value. */
120
121 #ifndef fp_barrierf
122 #define fp_barrierf fp_barrierf
fp_barrierf(float x)123 static inline float fp_barrierf(float x)
124 {
125 volatile float y = x;
126 return y;
127 }
128 #endif
129
130 #ifndef fp_barrier
131 #define fp_barrier fp_barrier
fp_barrier(double x)132 static inline double fp_barrier(double x)
133 {
134 volatile double y = x;
135 return y;
136 }
137 #endif
138
139 #ifndef fp_barrierl
140 #define fp_barrierl fp_barrierl
fp_barrierl(long double x)141 static inline long double fp_barrierl(long double x)
142 {
143 volatile long double y = x;
144 return y;
145 }
146 #endif
147
148 /* fp_force_eval ensures that the input value is computed when that's
149 otherwise unused. To prevent the constant folding of the input
150 expression, an additional fp_barrier may be needed or a compilation
151 mode that does so (e.g. -frounding-math in gcc). Then it can be
152 used to evaluate an expression for its fenv side-effects only. */
153
154 #ifndef fp_force_evalf
155 #define fp_force_evalf fp_force_evalf
fp_force_evalf(float x)156 static inline void fp_force_evalf(float x)
157 {
158 volatile float y;
159 y = x;
160 }
161 #endif
162
163 #ifndef fp_force_eval
164 #define fp_force_eval fp_force_eval
fp_force_eval(double x)165 static inline void fp_force_eval(double x)
166 {
167 volatile double y;
168 y = x;
169 }
170 #endif
171
172 #ifndef fp_force_evall
173 #define fp_force_evall fp_force_evall
fp_force_evall(long double x)174 static inline void fp_force_evall(long double x)
175 {
176 volatile long double y;
177 y = x;
178 }
179 #endif
180
181 #define FORCE_EVAL(x) do { \
182 if (sizeof(x) == sizeof(float)) { \
183 fp_force_evalf(x); \
184 } else if (sizeof(x) == sizeof(double)) { \
185 fp_force_eval(x); \
186 } else { \
187 fp_force_evall(x); \
188 } \
189 } while(0)
190
191 #define asuint(f) ((union{float _f; uint32_t _i;}){f})._i
192 #define asfloat(i) ((union{uint32_t _i; float _f;}){i})._f
193 #define asuint64(f) ((union{double _f; uint64_t _i;}){f})._i
194 #define asdouble(i) ((union{uint64_t _i; double _f;}){i})._f
195
196 #define EXTRACT_WORDS(hi,lo,d) \
197 do { \
198 uint64_t __u = asuint64(d); \
199 (hi) = __u >> 32; \
200 (lo) = (uint32_t)__u; \
201 } while (0)
202
203 #define GET_HIGH_WORD(hi,d) \
204 do { \
205 (hi) = asuint64(d) >> 32; \
206 } while (0)
207
208 #define GET_LOW_WORD(lo,d) \
209 do { \
210 (lo) = (uint32_t)asuint64(d); \
211 } while (0)
212
213 #define INSERT_WORDS(d,hi,lo) \
214 do { \
215 (d) = asdouble(((uint64_t)(hi)<<32) | (uint32_t)(lo)); \
216 } while (0)
217
218 #define SET_HIGH_WORD(d,hi) \
219 INSERT_WORDS(d, hi, (uint32_t)asuint64(d))
220
221 #define SET_LOW_WORD(d,lo) \
222 INSERT_WORDS(d, asuint64(d)>>32, lo)
223
224 #define GET_FLOAT_WORD(w,d) \
225 do { \
226 (w) = asuint(d); \
227 } while (0)
228
229 #define SET_FLOAT_WORD(d,w) \
230 do { \
231 (d) = asfloat(w); \
232 } while (0)
233
234 hidden int __rem_pio2_large(double*,double*,int,int,int);
235
236 hidden int __rem_pio2(double,double*);
237 hidden double __sin(double,double,int);
238 hidden double __cos(double,double);
239 hidden double __tan(double,double,int);
240 hidden double __expo2(double);
241
242 hidden int __rem_pio2f(float,double*);
243 hidden float __sindf(double);
244 hidden float __cosdf(double);
245 hidden float __tandf(double,int);
246 hidden float __expo2f(float);
247
248 hidden int __rem_pio2l(long double, long double *);
249 hidden long double __sinl(long double, long double, int);
250 hidden long double __cosl(long double, long double);
251 hidden long double __tanl(long double, long double, int);
252
253 hidden long double __polevll(long double, const long double *, int);
254 hidden long double __p1evll(long double, const long double *, int);
255
256 extern int __signgam;
257 hidden double __lgamma_r(double, int *);
258 hidden float __lgammaf_r(float, int *);
259
260 /* error handling functions */
261 hidden float __math_xflowf(uint32_t, float);
262 hidden float __math_uflowf(uint32_t);
263 hidden float __math_oflowf(uint32_t);
264 hidden float __math_divzerof(uint32_t);
265 hidden float __math_invalidf(float);
266 hidden double __math_xflow(uint32_t, double);
267 hidden double __math_uflow(uint32_t);
268 hidden double __math_oflow(uint32_t);
269 hidden double __math_divzero(uint32_t);
270 hidden double __math_invalid(double);
271
272 #endif
273