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 13union 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. */ 23union 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 32union 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 46union 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. */ 83static 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. */ 88static 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 105static inline float eval_as_float(float x) 106{ 107 float y = x; 108 return y; 109} 110 111static 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 123static 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 132static 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 141static 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 156static 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 165static 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 174static 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) \ 197do { \ 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) \ 204do { \ 205 (hi) = asuint64(d) >> 32; \ 206} while (0) 207 208#define GET_LOW_WORD(lo,d) \ 209do { \ 210 (lo) = (uint32_t)asuint64(d); \ 211} while (0) 212 213#define INSERT_WORDS(d,hi,lo) \ 214do { \ 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) \ 225do { \ 226 (w) = asuint(d); \ 227} while (0) 228 229#define SET_FLOAT_WORD(d,w) \ 230do { \ 231 (d) = asfloat(w); \ 232} while (0) 233 234hidden int __rem_pio2_large(double*,double*,int,int,int); 235 236hidden int __rem_pio2(double,double*); 237hidden double __sin(double,double,int); 238hidden double __cos(double,double); 239hidden double __tan(double,double,int); 240hidden double __expo2(double); 241 242hidden int __rem_pio2f(float,double*); 243hidden float __sindf(double); 244hidden float __cosdf(double); 245hidden float __tandf(double,int); 246hidden float __expo2f(float); 247 248hidden int __rem_pio2l(long double, long double *); 249hidden long double __sinl(long double, long double, int); 250hidden long double __cosl(long double, long double); 251hidden long double __tanl(long double, long double, int); 252 253hidden long double __polevll(long double, const long double *, int); 254hidden long double __p1evll(long double, const long double *, int); 255 256extern int __signgam; 257hidden double __lgamma_r(double, int *); 258hidden float __lgammaf_r(float, int *); 259 260/* error handling functions */ 261hidden float __math_xflowf(uint32_t, float); 262hidden float __math_uflowf(uint32_t); 263hidden float __math_oflowf(uint32_t); 264hidden float __math_divzerof(uint32_t); 265hidden float __math_invalidf(float); 266hidden double __math_xflow(uint32_t, double); 267hidden double __math_uflow(uint32_t); 268hidden double __math_oflow(uint32_t); 269hidden double __math_divzero(uint32_t); 270hidden double __math_invalid(double); 271 272#endif 273