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