blob: 98bf5c68ce4d187e79180b5c8c360d3e8abb7793 [file] [log] [blame]
Rich Felkerb69f6952012-03-13 01:17:53 -04001#ifndef _LIBM_H
2#define _LIBM_H
3
4#include <stdint.h>
5#include <float.h>
6#include <math.h>
Szabolcs Nagyaf5f6d92013-09-02 00:38:51 +00007#include <endian.h>
Szabolcs Nagyb50d3152018-11-26 23:30:00 +00008#include "fp_arch.h"
Rich Felkerb69f6952012-03-13 01:17:53 -04009
Szabolcs Nagyaf5f6d92013-09-02 00:38:51 +000010#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};
Rich Felker18f02c42018-06-14 11:00:58 -040019#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};
Szabolcs Nagyaf5f6d92013-09-02 00:38:51 +000030#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};
Szabolcs Nagyf4e46322015-03-10 20:01:20 +000044#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};
Szabolcs Nagyaf5f6d92013-09-02 00:38:51 +000058#else
59#error Unsupported long double representation
60#endif
61
Szabolcs Nagy169fc002018-12-02 18:53:37 +000062/* Support non-nearest rounding mode. */
63#define WANT_ROUNDING 1
64/* Support signaling NaNs. */
65#define WANT_SNAN 0
66
Szabolcs Nagy3f94c642017-10-22 18:06:00 +000067#ifndef TOINT_INTRINSICS
68#define TOINT_INTRINSICS 0
69#endif
70
71#if TOINT_INTRINSICS
72/* Round x to nearest int in all rounding modes, ties have to be rounded
73 consistently with converttoint so the results match. If the result
74 would be outside of [-2^31, 2^31-1] then the semantics is unspecified. */
75static double_t roundtoint(double_t);
76
77/* Convert x to nearest int in all rounding modes, ties have to be rounded
78 consistently with roundtoint. If the result is not representible in an
79 int32_t then the semantics is unspecified. */
80static int32_t converttoint(double_t);
81#endif
82
Szabolcs Nagy3c3763f2018-11-29 23:33:40 +000083/* Helps static branch prediction so hot path can be better optimized. */
84#ifdef __GNUC__
85#define predict_true(x) __builtin_expect(!!(x), 1)
86#define predict_false(x) __builtin_expect(x, 0)
87#else
88#define predict_true(x) (x)
89#define predict_false(x) (x)
90#endif
91
Szabolcs Nagyfe545442018-12-01 23:52:34 +000092/* Evaluate an expression as the specified type. With standard excess
93 precision handling a type cast or assignment is enough (with
94 -ffloat-store an assignment is required, in old compilers argument
95 passing and return statement may not drop excess precision). */
96
97static inline float eval_as_float(float x)
98{
99 float y = x;
100 return y;
101}
102
103static inline double eval_as_double(double x)
104{
105 double y = x;
106 return y;
107}
108
Szabolcs Nagyb50d3152018-11-26 23:30:00 +0000109/* fp_barrier returns its input, but limits code transformations
110 as if it had a side-effect (e.g. observable io) and returned
111 an arbitrary value. */
112
113#ifndef fp_barrierf
114#define fp_barrierf fp_barrierf
115static inline float fp_barrierf(float x)
116{
117 volatile float y = x;
118 return y;
119}
120#endif
121
122#ifndef fp_barrier
123#define fp_barrier fp_barrier
124static inline double fp_barrier(double x)
125{
126 volatile double y = x;
127 return y;
128}
129#endif
130
131#ifndef fp_barrierl
132#define fp_barrierl fp_barrierl
133static inline long double fp_barrierl(long double x)
134{
135 volatile long double y = x;
136 return y;
137}
138#endif
139
140/* fp_force_eval ensures that the input value is computed when that's
141 otherwise unused. To prevent the constant folding of the input
142 expression, an additional fp_barrier may be needed or a compilation
143 mode that does so (e.g. -frounding-math in gcc). Then it can be
144 used to evaluate an expression for its fenv side-effects only. */
145
146#ifndef fp_force_evalf
147#define fp_force_evalf fp_force_evalf
148static inline void fp_force_evalf(float x)
149{
150 volatile float y = x;
151}
152#endif
153
154#ifndef fp_force_eval
155#define fp_force_eval fp_force_eval
156static inline void fp_force_eval(double x)
157{
158 volatile double y = x;
159}
160#endif
161
162#ifndef fp_force_evall
163#define fp_force_evall fp_force_evall
164static inline void fp_force_evall(long double x)
165{
166 volatile long double y = x;
167}
168#endif
169
Szabolcs Nagy8dba5482013-09-04 17:36:00 +0000170#define FORCE_EVAL(x) do { \
171 if (sizeof(x) == sizeof(float)) { \
Szabolcs Nagyb50d3152018-11-26 23:30:00 +0000172 fp_force_evalf(x); \
Szabolcs Nagy8dba5482013-09-04 17:36:00 +0000173 } else if (sizeof(x) == sizeof(double)) { \
Szabolcs Nagyb50d3152018-11-26 23:30:00 +0000174 fp_force_eval(x); \
Szabolcs Nagy8dba5482013-09-04 17:36:00 +0000175 } else { \
Szabolcs Nagyb50d3152018-11-26 23:30:00 +0000176 fp_force_evall(x); \
Szabolcs Nagy8dba5482013-09-04 17:36:00 +0000177 } \
nsz6ab81362012-05-06 21:24:28 +0200178} while(0)
179
Szabolcs Nagyd59e5042017-10-21 21:09:02 +0000180#define asuint(f) ((union{float _f; uint32_t _i;}){f})._i
181#define asfloat(i) ((union{uint32_t _i; float _f;}){i})._f
182#define asuint64(f) ((union{double _f; uint64_t _i;}){f})._i
183#define asdouble(i) ((union{uint64_t _i; double _f;}){i})._f
184
Szabolcs Nagy8dba5482013-09-04 17:36:00 +0000185#define EXTRACT_WORDS(hi,lo,d) \
186do { \
Szabolcs Nagyd59e5042017-10-21 21:09:02 +0000187 uint64_t __u = asuint64(d); \
188 (hi) = __u >> 32; \
189 (lo) = (uint32_t)__u; \
Rich Felkerb69f6952012-03-13 01:17:53 -0400190} while (0)
191
Szabolcs Nagy8dba5482013-09-04 17:36:00 +0000192#define GET_HIGH_WORD(hi,d) \
193do { \
Szabolcs Nagyd59e5042017-10-21 21:09:02 +0000194 (hi) = asuint64(d) >> 32; \
Rich Felkerb69f6952012-03-13 01:17:53 -0400195} while (0)
196
Szabolcs Nagy8dba5482013-09-04 17:36:00 +0000197#define GET_LOW_WORD(lo,d) \
198do { \
Szabolcs Nagyd59e5042017-10-21 21:09:02 +0000199 (lo) = (uint32_t)asuint64(d); \
Rich Felkerb69f6952012-03-13 01:17:53 -0400200} while (0)
201
Szabolcs Nagy8dba5482013-09-04 17:36:00 +0000202#define INSERT_WORDS(d,hi,lo) \
203do { \
Szabolcs Nagyd59e5042017-10-21 21:09:02 +0000204 (d) = asdouble(((uint64_t)(hi)<<32) | (uint32_t)(lo)); \
Rich Felkerb69f6952012-03-13 01:17:53 -0400205} while (0)
206
Szabolcs Nagy8dba5482013-09-04 17:36:00 +0000207#define SET_HIGH_WORD(d,hi) \
Szabolcs Nagyd59e5042017-10-21 21:09:02 +0000208 INSERT_WORDS(d, hi, (uint32_t)asuint64(d))
Rich Felkerb69f6952012-03-13 01:17:53 -0400209
Szabolcs Nagy8dba5482013-09-04 17:36:00 +0000210#define SET_LOW_WORD(d,lo) \
Szabolcs Nagyd59e5042017-10-21 21:09:02 +0000211 INSERT_WORDS(d, asuint64(d)>>32, lo)
Rich Felkerb69f6952012-03-13 01:17:53 -0400212
Szabolcs Nagy8dba5482013-09-04 17:36:00 +0000213#define GET_FLOAT_WORD(w,d) \
214do { \
Szabolcs Nagyd59e5042017-10-21 21:09:02 +0000215 (w) = asuint(d); \
Rich Felkerb69f6952012-03-13 01:17:53 -0400216} while (0)
217
Szabolcs Nagy8dba5482013-09-04 17:36:00 +0000218#define SET_FLOAT_WORD(d,w) \
219do { \
Szabolcs Nagyd59e5042017-10-21 21:09:02 +0000220 (d) = asfloat(w); \
Rich Felkerb69f6952012-03-13 01:17:53 -0400221} while (0)
222
Rich Felker46e38952018-09-11 11:04:35 -0400223hidden int __rem_pio2_large(double*,double*,int,int,int);
Rich Felkerb69f6952012-03-13 01:17:53 -0400224
Rich Felker46e38952018-09-11 11:04:35 -0400225hidden int __rem_pio2(double,double*);
226hidden double __sin(double,double,int);
227hidden double __cos(double,double);
228hidden double __tan(double,double,int);
229hidden double __expo2(double);
Rich Felkerb69f6952012-03-13 01:17:53 -0400230
Rich Felker46e38952018-09-11 11:04:35 -0400231hidden int __rem_pio2f(float,double*);
232hidden float __sindf(double);
233hidden float __cosdf(double);
234hidden float __tandf(double,int);
235hidden float __expo2f(float);
Rich Felkerb69f6952012-03-13 01:17:53 -0400236
Rich Felker46e38952018-09-11 11:04:35 -0400237hidden int __rem_pio2l(long double, long double *);
238hidden long double __sinl(long double, long double, int);
239hidden long double __cosl(long double, long double);
240hidden long double __tanl(long double, long double, int);
Rich Felkerb69f6952012-03-13 01:17:53 -0400241
Rich Felker46e38952018-09-11 11:04:35 -0400242hidden long double __polevll(long double, const long double *, int);
243hidden long double __p1evll(long double, const long double *, int);
Rich Felkerb69f6952012-03-13 01:17:53 -0400244
Rich Felker59d88942018-09-07 12:32:12 -0400245extern int __signgam;
Rich Felker46e38952018-09-11 11:04:35 -0400246hidden double __lgamma_r(double, int *);
247hidden float __lgammaf_r(float, int *);
Rich Felker59d88942018-09-07 12:32:12 -0400248
Szabolcs Nagy9ef6ca42017-10-22 13:51:35 +0000249/* error handling functions */
250hidden float __math_xflowf(uint32_t, float);
251hidden float __math_uflowf(uint32_t);
252hidden float __math_oflowf(uint32_t);
253hidden float __math_divzerof(uint32_t);
254hidden float __math_invalidf(float);
Szabolcs Nagy4f8acf92018-11-30 21:15:23 +0000255hidden double __math_xflow(uint32_t, double);
256hidden double __math_uflow(uint32_t);
257hidden double __math_oflow(uint32_t);
258hidden double __math_divzero(uint32_t);
259hidden double __math_invalid(double);
Szabolcs Nagy9ef6ca42017-10-22 13:51:35 +0000260
Rich Felkerb69f6952012-03-13 01:17:53 -0400261#endif