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 <features.h> 9 #include "fp_arch.h" 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 105 static inline float eval_as_float(float x) 106 { 107 float y = x; 108 return y; 109 } 110 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 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 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 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 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 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 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