math: add eval_as_float and eval_as_double
[musl] / src / internal / libm.h
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
12 union 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. */
22 union 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
31 union 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
45 union 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 /* Evaluate an expression as the specified type. With standard excess
63    precision handling a type cast or assignment is enough (with
64    -ffloat-store an assignment is required, in old compilers argument
65    passing and return statement may not drop excess precision).  */
66
67 static inline float eval_as_float(float x)
68 {
69         float y = x;
70         return y;
71 }
72
73 static inline double eval_as_double(double x)
74 {
75         double y = x;
76         return y;
77 }
78
79 /* fp_barrier returns its input, but limits code transformations
80    as if it had a side-effect (e.g. observable io) and returned
81    an arbitrary value.  */
82
83 #ifndef fp_barrierf
84 #define fp_barrierf fp_barrierf
85 static inline float fp_barrierf(float x)
86 {
87         volatile float y = x;
88         return y;
89 }
90 #endif
91
92 #ifndef fp_barrier
93 #define fp_barrier fp_barrier
94 static inline double fp_barrier(double x)
95 {
96         volatile double y = x;
97         return y;
98 }
99 #endif
100
101 #ifndef fp_barrierl
102 #define fp_barrierl fp_barrierl
103 static inline long double fp_barrierl(long double x)
104 {
105         volatile long double y = x;
106         return y;
107 }
108 #endif
109
110 /* fp_force_eval ensures that the input value is computed when that's
111    otherwise unused.  To prevent the constant folding of the input
112    expression, an additional fp_barrier may be needed or a compilation
113    mode that does so (e.g. -frounding-math in gcc). Then it can be
114    used to evaluate an expression for its fenv side-effects only.   */
115
116 #ifndef fp_force_evalf
117 #define fp_force_evalf fp_force_evalf
118 static inline void fp_force_evalf(float x)
119 {
120         volatile float y = x;
121 }
122 #endif
123
124 #ifndef fp_force_eval
125 #define fp_force_eval fp_force_eval
126 static inline void fp_force_eval(double x)
127 {
128         volatile double y = x;
129 }
130 #endif
131
132 #ifndef fp_force_evall
133 #define fp_force_evall fp_force_evall
134 static inline void fp_force_evall(long double x)
135 {
136         volatile long double y = x;
137 }
138 #endif
139
140 #define FORCE_EVAL(x) do {                        \
141         if (sizeof(x) == sizeof(float)) {         \
142                 fp_force_evalf(x);                \
143         } else if (sizeof(x) == sizeof(double)) { \
144                 fp_force_eval(x);                 \
145         } else {                                  \
146                 fp_force_evall(x);                \
147         }                                         \
148 } while(0)
149
150 #define asuint(f) ((union{float _f; uint32_t _i;}){f})._i
151 #define asfloat(i) ((union{uint32_t _i; float _f;}){i})._f
152 #define asuint64(f) ((union{double _f; uint64_t _i;}){f})._i
153 #define asdouble(i) ((union{uint64_t _i; double _f;}){i})._f
154
155 #define EXTRACT_WORDS(hi,lo,d)                    \
156 do {                                              \
157   uint64_t __u = asuint64(d);                     \
158   (hi) = __u >> 32;                               \
159   (lo) = (uint32_t)__u;                           \
160 } while (0)
161
162 #define GET_HIGH_WORD(hi,d)                       \
163 do {                                              \
164   (hi) = asuint64(d) >> 32;                       \
165 } while (0)
166
167 #define GET_LOW_WORD(lo,d)                        \
168 do {                                              \
169   (lo) = (uint32_t)asuint64(d);                   \
170 } while (0)
171
172 #define INSERT_WORDS(d,hi,lo)                     \
173 do {                                              \
174   (d) = asdouble(((uint64_t)(hi)<<32) | (uint32_t)(lo)); \
175 } while (0)
176
177 #define SET_HIGH_WORD(d,hi)                       \
178   INSERT_WORDS(d, hi, (uint32_t)asuint64(d))
179
180 #define SET_LOW_WORD(d,lo)                        \
181   INSERT_WORDS(d, asuint64(d)>>32, lo)
182
183 #define GET_FLOAT_WORD(w,d)                       \
184 do {                                              \
185   (w) = asuint(d);                                \
186 } while (0)
187
188 #define SET_FLOAT_WORD(d,w)                       \
189 do {                                              \
190   (d) = asfloat(w);                               \
191 } while (0)
192
193 hidden int    __rem_pio2_large(double*,double*,int,int,int);
194
195 hidden int    __rem_pio2(double,double*);
196 hidden double __sin(double,double,int);
197 hidden double __cos(double,double);
198 hidden double __tan(double,double,int);
199 hidden double __expo2(double);
200
201 hidden int    __rem_pio2f(float,double*);
202 hidden float  __sindf(double);
203 hidden float  __cosdf(double);
204 hidden float  __tandf(double,int);
205 hidden float  __expo2f(float);
206
207 hidden int __rem_pio2l(long double, long double *);
208 hidden long double __sinl(long double, long double, int);
209 hidden long double __cosl(long double, long double);
210 hidden long double __tanl(long double, long double, int);
211
212 hidden long double __polevll(long double, const long double *, int);
213 hidden long double __p1evll(long double, const long double *, int);
214
215 extern int __signgam;
216 hidden double __lgamma_r(double, int *);
217 hidden float __lgammaf_r(float, int *);
218
219 #endif