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