diff options
Diffstat (limited to 'src/internal/libm.h')
| -rw-r--r-- | src/internal/libm.h | 233 |
1 files changed, 154 insertions, 79 deletions
diff --git a/src/internal/libm.h b/src/internal/libm.h index fd916277..72ad17d8 100644 --- a/src/internal/libm.h +++ b/src/internal/libm.h @@ -1,23 +1,11 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/math_private.h */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - #ifndef _LIBM_H #define _LIBM_H #include <stdint.h> #include <float.h> #include <math.h> -#include <complex.h> #include <endian.h> +#include "fp_arch.h" #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __LITTLE_ENDIAN @@ -71,124 +59,196 @@ union ldshape { #error Unsupported long double representation #endif +/* Support non-nearest rounding mode. */ +#define WANT_ROUNDING 1 +/* Support signaling NaNs. */ +#define WANT_SNAN 0 + +#if WANT_SNAN +#error SNaN is unsupported +#else +#define issignalingf_inline(x) 0 +#define issignaling_inline(x) 0 +#endif + +#ifndef TOINT_INTRINSICS +#define TOINT_INTRINSICS 0 +#endif + +#if TOINT_INTRINSICS +/* Round x to nearest int in all rounding modes, ties have to be rounded + consistently with converttoint so the results match. If the result + would be outside of [-2^31, 2^31-1] then the semantics is unspecified. */ +static double_t roundtoint(double_t); + +/* Convert x to nearest int in all rounding modes, ties have to be rounded + consistently with roundtoint. If the result is not representible in an + int32_t then the semantics is unspecified. */ +static int32_t converttoint(double_t); +#endif + +/* Helps static branch prediction so hot path can be better optimized. */ +#ifdef __GNUC__ +#define predict_true(x) __builtin_expect(!!(x), 1) +#define predict_false(x) __builtin_expect(x, 0) +#else +#define predict_true(x) (x) +#define predict_false(x) (x) +#endif + +/* Evaluate an expression as the specified type. With standard excess + precision handling a type cast or assignment is enough (with + -ffloat-store an assignment is required, in old compilers argument + passing and return statement may not drop excess precision). */ + +static inline float eval_as_float(float x) +{ + float y = x; + return y; +} + +static inline double eval_as_double(double x) +{ + double y = x; + return y; +} + +/* fp_barrier returns its input, but limits code transformations + as if it had a side-effect (e.g. observable io) and returned + an arbitrary value. */ + +#ifndef fp_barrierf +#define fp_barrierf fp_barrierf +static inline float fp_barrierf(float x) +{ + volatile float y = x; + return y; +} +#endif + +#ifndef fp_barrier +#define fp_barrier fp_barrier +static inline double fp_barrier(double x) +{ + volatile double y = x; + return y; +} +#endif + +#ifndef fp_barrierl +#define fp_barrierl fp_barrierl +static inline long double fp_barrierl(long double x) +{ + volatile long double y = x; + return y; +} +#endif + +/* fp_force_eval ensures that the input value is computed when that's + otherwise unused. To prevent the constant folding of the input + expression, an additional fp_barrier may be needed or a compilation + mode that does so (e.g. -frounding-math in gcc). Then it can be + used to evaluate an expression for its fenv side-effects only. */ + +#ifndef fp_force_evalf +#define fp_force_evalf fp_force_evalf +static inline void fp_force_evalf(float x) +{ + volatile float y; + y = x; +} +#endif + +#ifndef fp_force_eval +#define fp_force_eval fp_force_eval +static inline void fp_force_eval(double x) +{ + volatile double y; + y = x; +} +#endif + +#ifndef fp_force_evall +#define fp_force_evall fp_force_evall +static inline void fp_force_evall(long double x) +{ + volatile long double y; + y = x; +} +#endif + #define FORCE_EVAL(x) do { \ if (sizeof(x) == sizeof(float)) { \ - volatile float __x; \ - __x = (x); \ + fp_force_evalf(x); \ } else if (sizeof(x) == sizeof(double)) { \ - volatile double __x; \ - __x = (x); \ + fp_force_eval(x); \ } else { \ - volatile long double __x; \ - __x = (x); \ + fp_force_evall(x); \ } \ } while(0) -/* Get two 32 bit ints from a double. */ +#define asuint(f) ((union{float _f; uint32_t _i;}){f})._i +#define asfloat(i) ((union{uint32_t _i; float _f;}){i})._f +#define asuint64(f) ((union{double _f; uint64_t _i;}){f})._i +#define asdouble(i) ((union{uint64_t _i; double _f;}){i})._f + #define EXTRACT_WORDS(hi,lo,d) \ do { \ - union {double f; uint64_t i;} __u; \ - __u.f = (d); \ - (hi) = __u.i >> 32; \ - (lo) = (uint32_t)__u.i; \ + uint64_t __u = asuint64(d); \ + (hi) = __u >> 32; \ + (lo) = (uint32_t)__u; \ } while (0) -/* Get the more significant 32 bit int from a double. */ #define GET_HIGH_WORD(hi,d) \ do { \ - union {double f; uint64_t i;} __u; \ - __u.f = (d); \ - (hi) = __u.i >> 32; \ + (hi) = asuint64(d) >> 32; \ } while (0) -/* Get the less significant 32 bit int from a double. */ #define GET_LOW_WORD(lo,d) \ do { \ - union {double f; uint64_t i;} __u; \ - __u.f = (d); \ - (lo) = (uint32_t)__u.i; \ + (lo) = (uint32_t)asuint64(d); \ } while (0) -/* Set a double from two 32 bit ints. */ #define INSERT_WORDS(d,hi,lo) \ do { \ - union {double f; uint64_t i;} __u; \ - __u.i = ((uint64_t)(hi)<<32) | (uint32_t)(lo); \ - (d) = __u.f; \ + (d) = asdouble(((uint64_t)(hi)<<32) | (uint32_t)(lo)); \ } while (0) -/* Set the more significant 32 bits of a double from an int. */ #define SET_HIGH_WORD(d,hi) \ -do { \ - union {double f; uint64_t i;} __u; \ - __u.f = (d); \ - __u.i &= 0xffffffff; \ - __u.i |= (uint64_t)(hi) << 32; \ - (d) = __u.f; \ -} while (0) + INSERT_WORDS(d, hi, (uint32_t)asuint64(d)) -/* Set the less significant 32 bits of a double from an int. */ #define SET_LOW_WORD(d,lo) \ -do { \ - union {double f; uint64_t i;} __u; \ - __u.f = (d); \ - __u.i &= 0xffffffff00000000ull; \ - __u.i |= (uint32_t)(lo); \ - (d) = __u.f; \ -} while (0) + INSERT_WORDS(d, asuint64(d)>>32, lo) -/* Get a 32 bit int from a float. */ #define GET_FLOAT_WORD(w,d) \ do { \ - union {float f; uint32_t i;} __u; \ - __u.f = (d); \ - (w) = __u.i; \ + (w) = asuint(d); \ } while (0) -/* Set a float from a 32 bit int. */ #define SET_FLOAT_WORD(d,w) \ do { \ - union {float f; uint32_t i;} __u; \ - __u.i = (w); \ - (d) = __u.f; \ + (d) = asfloat(w); \ } while (0) -#undef __CMPLX -#undef CMPLX -#undef CMPLXF -#undef CMPLXL - -#define __CMPLX(x, y, t) \ - ((union { _Complex t __z; t __xy[2]; }){.__xy = {(x),(y)}}.__z) - -#define CMPLX(x, y) __CMPLX(x, y, double) -#define CMPLXF(x, y) __CMPLX(x, y, float) -#define CMPLXL(x, y) __CMPLX(x, y, long double) - -/* fdlibm kernel functions */ - hidden int __rem_pio2_large(double*,double*,int,int,int); hidden int __rem_pio2(double,double*); hidden double __sin(double,double,int); hidden double __cos(double,double); hidden double __tan(double,double,int); -hidden double __expo2(double); -hidden double complex __ldexp_cexp(double complex,int); +hidden double __expo2(double,double); hidden int __rem_pio2f(float,double*); hidden float __sindf(double); hidden float __cosdf(double); hidden float __tandf(double,int); -hidden float __expo2f(float); -hidden float complex __ldexp_cexpf(float complex,int); +hidden float __expo2f(float,float); hidden int __rem_pio2l(long double, long double *); hidden long double __sinl(long double, long double, int); hidden long double __cosl(long double, long double); hidden long double __tanl(long double, long double, int); -/* polynomial evaluation */ hidden long double __polevll(long double, const long double *, int); hidden long double __p1evll(long double, const long double *, int); @@ -196,4 +256,19 @@ extern int __signgam; hidden double __lgamma_r(double, int *); hidden float __lgammaf_r(float, int *); +/* error handling functions */ +hidden float __math_xflowf(uint32_t, float); +hidden float __math_uflowf(uint32_t); +hidden float __math_oflowf(uint32_t); +hidden float __math_divzerof(uint32_t); +hidden float __math_invalidf(float); +hidden double __math_xflow(uint32_t, double); +hidden double __math_uflow(uint32_t); +hidden double __math_oflow(uint32_t); +hidden double __math_divzero(uint32_t); +hidden double __math_invalid(double); +#if LDBL_MANT_DIG != DBL_MANT_DIG +hidden long double __math_invalidl(long double); +#endif + #endif |