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