1 /* 2 * Configuration for math routines. 3 * 4 * Copyright (c) 2017-2018, Arm Limited. 5 * SPDX-License-Identifier: MIT 6 */ 7 8 #ifndef _MATH_CONFIG_H 9 #define _MATH_CONFIG_H 10 11 #include <math.h> 12 #include <stdint.h> 13 14 #ifndef WANT_ROUNDING 15 /* Correct special case results in non-nearest rounding modes. */ 16 # define WANT_ROUNDING 1 17 #endif 18 #ifndef WANT_ERRNO 19 /* Set errno according to ISO C with (math_errhandling & MATH_ERRNO) != 0. */ 20 # define WANT_ERRNO 1 21 #endif 22 #ifndef WANT_ERRNO_UFLOW 23 /* Set errno to ERANGE if result underflows to 0 (in all rounding modes). */ 24 # define WANT_ERRNO_UFLOW (WANT_ROUNDING && WANT_ERRNO) 25 #endif 26 27 /* Compiler can inline round as a single instruction. */ 28 #ifndef HAVE_FAST_ROUND 29 # if __aarch64__ 30 # define HAVE_FAST_ROUND 1 31 # else 32 # define HAVE_FAST_ROUND 0 33 # endif 34 #endif 35 36 /* Compiler can inline lround, but not (long)round(x). */ 37 #ifndef HAVE_FAST_LROUND 38 # if __aarch64__ && (100*__GNUC__ + __GNUC_MINOR__) >= 408 && __NO_MATH_ERRNO__ 39 # define HAVE_FAST_LROUND 1 40 # else 41 # define HAVE_FAST_LROUND 0 42 # endif 43 #endif 44 45 /* Compiler can inline fma as a single instruction. */ 46 #ifndef HAVE_FAST_FMA 47 # if defined FP_FAST_FMA || __aarch64__ 48 # define HAVE_FAST_FMA 1 49 # else 50 # define HAVE_FAST_FMA 0 51 # endif 52 #endif 53 54 /* Provide *_finite symbols and some of the glibc hidden symbols 55 so libmathlib can be used with binaries compiled against glibc 56 to interpose math functions with both static and dynamic linking. */ 57 #ifndef USE_GLIBC_ABI 58 # if __GNUC__ 59 # define USE_GLIBC_ABI 1 60 # else 61 # define USE_GLIBC_ABI 0 62 # endif 63 #endif 64 65 /* Optionally used extensions. */ 66 #ifdef __GNUC__ 67 # define HIDDEN __attribute__ ((__visibility__ ("hidden"))) 68 # define NOINLINE __attribute__ ((noinline)) 69 # define UNUSED __attribute__ ((unused)) 70 # define likely(x) __builtin_expect (!!(x), 1) 71 # define unlikely(x) __builtin_expect (x, 0) 72 # if __GNUC__ >= 9 73 # define attribute_copy(f) __attribute__ ((copy (f))) 74 # else 75 # define attribute_copy(f) 76 # endif 77 # define strong_alias(f, a) \ 78 extern __typeof (f) a __attribute__ ((alias (#f))) attribute_copy (f); 79 # define hidden_alias(f, a) \ 80 extern __typeof (f) a __attribute__ ((alias (#f), visibility ("hidden"))) \ 81 attribute_copy (f); 82 #else 83 # define HIDDEN 84 # define NOINLINE 85 # define UNUSED 86 # define likely(x) (x) 87 # define unlikely(x) (x) 88 #endif 89 90 #if HAVE_FAST_ROUND 91 /* When set, the roundtoint and converttoint functions are provided with 92 the semantics documented below. */ 93 # define TOINT_INTRINSICS 1 94 95 /* Round x to nearest int in all rounding modes, ties have to be rounded 96 consistently with converttoint so the results match. If the result 97 would be outside of [-2^31, 2^31-1] then the semantics is unspecified. */ 98 static inline double_t 99 roundtoint (double_t x) 100 { 101 return round (x); 102 } 103 104 /* Convert x to nearest int in all rounding modes, ties have to be rounded 105 consistently with roundtoint. If the result is not representible in an 106 int32_t then the semantics is unspecified. */ 107 static inline int32_t 108 converttoint (double_t x) 109 { 110 # if HAVE_FAST_LROUND 111 return lround (x); 112 # else 113 return (long) round (x); 114 # endif 115 } 116 #endif 117 118 static inline uint32_t 119 asuint (float f) 120 { 121 union 122 { 123 float f; 124 uint32_t i; 125 } u = {f}; 126 return u.i; 127 } 128 129 static inline float 130 asfloat (uint32_t i) 131 { 132 union 133 { 134 uint32_t i; 135 float f; 136 } u = {i}; 137 return u.f; 138 } 139 140 static inline uint64_t 141 asuint64 (double f) 142 { 143 union 144 { 145 double f; 146 uint64_t i; 147 } u = {f}; 148 return u.i; 149 } 150 151 static inline double 152 asdouble (uint64_t i) 153 { 154 union 155 { 156 uint64_t i; 157 double f; 158 } u = {i}; 159 return u.f; 160 } 161 162 #ifndef IEEE_754_2008_SNAN 163 # define IEEE_754_2008_SNAN 1 164 #endif 165 static inline int 166 issignalingf_inline (float x) 167 { 168 uint32_t ix = asuint (x); 169 if (!IEEE_754_2008_SNAN) 170 return (ix & 0x7fc00000) == 0x7fc00000; 171 return 2 * (ix ^ 0x00400000) > 2u * 0x7fc00000; 172 } 173 174 static inline int 175 issignaling_inline (double x) 176 { 177 uint64_t ix = asuint64 (x); 178 if (!IEEE_754_2008_SNAN) 179 return (ix & 0x7ff8000000000000) == 0x7ff8000000000000; 180 return 2 * (ix ^ 0x0008000000000000) > 2 * 0x7ff8000000000000ULL; 181 } 182 183 #if __aarch64__ && __GNUC__ 184 /* Prevent the optimization of a floating-point expression. */ 185 static inline float 186 opt_barrier_float (float x) 187 { 188 __asm__ __volatile__ ("" : "+w" (x)); 189 return x; 190 } 191 static inline double 192 opt_barrier_double (double x) 193 { 194 __asm__ __volatile__ ("" : "+w" (x)); 195 return x; 196 } 197 /* Force the evaluation of a floating-point expression for its side-effect. */ 198 static inline void 199 force_eval_float (float x) 200 { 201 __asm__ __volatile__ ("" : "+w" (x)); 202 } 203 static inline void 204 force_eval_double (double x) 205 { 206 __asm__ __volatile__ ("" : "+w" (x)); 207 } 208 #else 209 static inline float 210 opt_barrier_float (float x) 211 { 212 volatile float y = x; 213 return y; 214 } 215 static inline double 216 opt_barrier_double (double x) 217 { 218 volatile double y = x; 219 return y; 220 } 221 static inline void 222 force_eval_float (float x) 223 { 224 volatile float y UNUSED = x; 225 } 226 static inline void 227 force_eval_double (double x) 228 { 229 volatile double y UNUSED = x; 230 } 231 #endif 232 233 /* Evaluate an expression as the specified type, normally a type 234 cast should be enough, but compilers implement non-standard 235 excess-precision handling, so when FLT_EVAL_METHOD != 0 then 236 these functions may need to be customized. */ 237 static inline float 238 eval_as_float (float x) 239 { 240 return x; 241 } 242 static inline double 243 eval_as_double (double x) 244 { 245 return x; 246 } 247 248 /* Error handling tail calls for special cases, with a sign argument. 249 The sign of the return value is set if the argument is non-zero. */ 250 251 /* The result overflows. */ 252 HIDDEN float __math_oflowf (uint32_t); 253 /* The result underflows to 0 in nearest rounding mode. */ 254 HIDDEN float __math_uflowf (uint32_t); 255 /* The result underflows to 0 in some directed rounding mode only. */ 256 HIDDEN float __math_may_uflowf (uint32_t); 257 /* Division by zero. */ 258 HIDDEN float __math_divzerof (uint32_t); 259 /* The result overflows. */ 260 HIDDEN double __math_oflow (uint32_t); 261 /* The result underflows to 0 in nearest rounding mode. */ 262 HIDDEN double __math_uflow (uint32_t); 263 /* The result underflows to 0 in some directed rounding mode only. */ 264 HIDDEN double __math_may_uflow (uint32_t); 265 /* Division by zero. */ 266 HIDDEN double __math_divzero (uint32_t); 267 268 /* Error handling using input checking. */ 269 270 /* Invalid input unless it is a quiet NaN. */ 271 HIDDEN float __math_invalidf (float); 272 /* Invalid input unless it is a quiet NaN. */ 273 HIDDEN double __math_invalid (double); 274 275 /* Error handling using output checking, only for errno setting. */ 276 277 /* Check if the result overflowed to infinity. */ 278 HIDDEN double __math_check_oflow (double); 279 /* Check if the result underflowed to 0. */ 280 HIDDEN double __math_check_uflow (double); 281 282 /* Check if the result overflowed to infinity. */ 283 static inline double 284 check_oflow (double x) 285 { 286 return WANT_ERRNO ? __math_check_oflow (x) : x; 287 } 288 289 /* Check if the result underflowed to 0. */ 290 static inline double 291 check_uflow (double x) 292 { 293 return WANT_ERRNO ? __math_check_uflow (x) : x; 294 } 295 296 297 /* Shared between expf, exp2f and powf. */ 298 #define EXP2F_TABLE_BITS 5 299 #define EXP2F_POLY_ORDER 3 300 extern const struct exp2f_data 301 { 302 uint64_t tab[1 << EXP2F_TABLE_BITS]; 303 double shift_scaled; 304 double poly[EXP2F_POLY_ORDER]; 305 double shift; 306 double invln2_scaled; 307 double poly_scaled[EXP2F_POLY_ORDER]; 308 } __exp2f_data HIDDEN; 309 310 #define LOGF_TABLE_BITS 4 311 #define LOGF_POLY_ORDER 4 312 extern const struct logf_data 313 { 314 struct 315 { 316 double invc, logc; 317 } tab[1 << LOGF_TABLE_BITS]; 318 double ln2; 319 double poly[LOGF_POLY_ORDER - 1]; /* First order coefficient is 1. */ 320 } __logf_data HIDDEN; 321 322 #define LOG2F_TABLE_BITS 4 323 #define LOG2F_POLY_ORDER 4 324 extern const struct log2f_data 325 { 326 struct 327 { 328 double invc, logc; 329 } tab[1 << LOG2F_TABLE_BITS]; 330 double poly[LOG2F_POLY_ORDER]; 331 } __log2f_data HIDDEN; 332 333 #define POWF_LOG2_TABLE_BITS 4 334 #define POWF_LOG2_POLY_ORDER 5 335 #if TOINT_INTRINSICS 336 # define POWF_SCALE_BITS EXP2F_TABLE_BITS 337 #else 338 # define POWF_SCALE_BITS 0 339 #endif 340 #define POWF_SCALE ((double) (1 << POWF_SCALE_BITS)) 341 extern const struct powf_log2_data 342 { 343 struct 344 { 345 double invc, logc; 346 } tab[1 << POWF_LOG2_TABLE_BITS]; 347 double poly[POWF_LOG2_POLY_ORDER]; 348 } __powf_log2_data HIDDEN; 349 350 351 #define EXP_TABLE_BITS 7 352 #define EXP_POLY_ORDER 5 353 /* Use polynomial that is optimized for a wider input range. This may be 354 needed for good precision in non-nearest rounding and !TOINT_INTRINSICS. */ 355 #define EXP_POLY_WIDE 0 356 /* Use close to nearest rounding toint when !TOINT_INTRINSICS. This may be 357 needed for good precision in non-nearest rouning and !EXP_POLY_WIDE. */ 358 #define EXP_USE_TOINT_NARROW 0 359 #define EXP2_POLY_ORDER 5 360 #define EXP2_POLY_WIDE 0 361 extern const struct exp_data 362 { 363 double invln2N; 364 double shift; 365 double negln2hiN; 366 double negln2loN; 367 double poly[4]; /* Last four coefficients. */ 368 double exp2_shift; 369 double exp2_poly[EXP2_POLY_ORDER]; 370 uint64_t tab[2*(1 << EXP_TABLE_BITS)]; 371 } __exp_data HIDDEN; 372 373 #define LOG_TABLE_BITS 7 374 #define LOG_POLY_ORDER 6 375 #define LOG_POLY1_ORDER 12 376 extern const struct log_data 377 { 378 double ln2hi; 379 double ln2lo; 380 double poly[LOG_POLY_ORDER - 1]; /* First coefficient is 1. */ 381 double poly1[LOG_POLY1_ORDER - 1]; 382 struct {double invc, logc;} tab[1 << LOG_TABLE_BITS]; 383 #if !HAVE_FAST_FMA 384 struct {double chi, clo;} tab2[1 << LOG_TABLE_BITS]; 385 #endif 386 } __log_data HIDDEN; 387 388 #define LOG2_TABLE_BITS 6 389 #define LOG2_POLY_ORDER 7 390 #define LOG2_POLY1_ORDER 11 391 extern const struct log2_data 392 { 393 double invln2hi; 394 double invln2lo; 395 double poly[LOG2_POLY_ORDER - 1]; 396 double poly1[LOG2_POLY1_ORDER - 1]; 397 struct {double invc, logc;} tab[1 << LOG2_TABLE_BITS]; 398 #if !HAVE_FAST_FMA 399 struct {double chi, clo;} tab2[1 << LOG2_TABLE_BITS]; 400 #endif 401 } __log2_data HIDDEN; 402 403 #define POW_LOG_TABLE_BITS 7 404 #define POW_LOG_POLY_ORDER 8 405 extern const struct pow_log_data 406 { 407 double ln2hi; 408 double ln2lo; 409 double poly[POW_LOG_POLY_ORDER - 1]; /* First coefficient is 1. */ 410 /* Note: the pad field is unused, but allows slightly faster indexing. */ 411 struct {double invc, pad, logc, logctail;} tab[1 << POW_LOG_TABLE_BITS]; 412 } __pow_log_data HIDDEN; 413 414 #endif 415