sysdeps/ieee754/dbl-64/math_config.h

   1 /* Configuration for double precision math routines.
   2    Copyright (C) 2018-2023 Free Software Foundation, Inc.
   3    This file is part of the GNU C Library.
   4
   5    The GNU C Library is free software; you can redistribute it and/or
   6    modify it under the terms of the GNU Lesser General Public
   7    License as published by the Free Software Foundation; either
   8    version 2.1 of the License, or (at your option) any later version.
   9
  10    The GNU C Library is distributed in the hope that it will be useful,
  11    but WITHOUT ANY WARRANTY; without even the implied warranty of
  12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13    Lesser General Public License for more details.
  14
  15    You should have received a copy of the GNU Lesser General Public
  16    License along with the GNU C Library; if not, see
  17    <https://www.gnu.org/licenses/>.  */
  18
  19 #ifndef _MATH_CONFIG_H
  20 #define _MATH_CONFIG_H
  21
  22 #include <math.h>
  23 #include <math_private.h>
  24 #include <nan-high-order-bit.h>
  25 #include <stdint.h>
  26
  27 #ifndef WANT_ROUNDING
  28 /* Correct special case results in non-nearest rounding modes.  */
  29 # define WANT_ROUNDING 1
  30 #endif
  31 #ifndef WANT_ERRNO
  32 /* Set errno according to ISO C with (math_errhandling & MATH_ERRNO) != 0.  */
  33 # define WANT_ERRNO 1
  34 #endif
  35 #ifndef WANT_ERRNO_UFLOW
  36 /* Set errno to ERANGE if result underflows to 0 (in all rounding modes).  */
  37 # define WANT_ERRNO_UFLOW (WANT_ROUNDING && WANT_ERRNO)
  38 #endif
  39
  40 #ifndef TOINT_INTRINSICS
  41 /* When set, the roundtoint and converttoint functions are provided with
  42    the semantics documented below.  */
  43 # define TOINT_INTRINSICS 0
  44 #endif
  45
  46 static inline int
  47 clz_uint64 (uint64_t x)
  48 {
  49   if (sizeof (uint64_t) == sizeof (unsigned long))
  50     return __builtin_clzl (x);
  51   else
  52     return __builtin_clzll (x);
  53 }
  54
  55 static inline int
  56 ctz_uint64 (uint64_t x)
  57 {
  58   if (sizeof (uint64_t) == sizeof (unsigned long))
  59     return __builtin_ctzl (x);
  60   else
  61     return __builtin_ctzll (x);
  62 }
  63
  64 #if TOINT_INTRINSICS
  65 /* Round x to nearest int in all rounding modes, ties have to be rounded
  66    consistently with converttoint so the results match.  If the result
  67    would be outside of [-2^31, 2^31-1] then the semantics is unspecified.  */
  68 static inline double_t
  69 roundtoint (double_t x);
  70
  71 /* Convert x to nearest int in all rounding modes, ties have to be rounded
  72    consistently with roundtoint.  If the result is not representible in an
  73    int32_t then the semantics is unspecified.  */
  74 static inline int32_t
  75 converttoint (double_t x);
  76 #endif
  77
  78 static inline uint64_t
  79 asuint64 (double f)
  80 {
  81   union
  82   {
  83     double f;
  84     uint64_t i;
  85   } u = {f};
  86   return u.i;
  87 }
  88
  89 static inline double
  90 asdouble (uint64_t i)
  91 {
  92   union
  93   {
  94     uint64_t i;
  95     double f;
  96   } u = {i};
  97   return u.f;
  98 }
  99
 100 static inline int
 101 issignaling_inline (double x)
 102 {
 103   uint64_t ix = asuint64 (x);
 104   if (HIGH_ORDER_BIT_IS_SET_FOR_SNAN)
 105     return (ix & 0x7ff8000000000000) == 0x7ff8000000000000;
 106   return 2 * (ix ^ 0x0008000000000000) > 2 * 0x7ff8000000000000ULL;
 107 }
 108
 109 #define BIT_WIDTH       64
 110 #define MANTISSA_WIDTH  52
 111 #define EXPONENT_WIDTH  11
 112 #define MANTISSA_MASK   UINT64_C(0x000fffffffffffff)
 113 #define EXPONENT_MASK   UINT64_C(0x7ff0000000000000)
 114 #define EXP_MANT_MASK   UINT64_C(0x7fffffffffffffff)
 115 #define QUIET_NAN_MASK  UINT64_C(0x0008000000000000)
 116 #define SIGN_MASK       UINT64_C(0x8000000000000000)
 117
 118 static inline bool
 119 is_nan (uint64_t x)
 120 {
 121   return (x & EXP_MANT_MASK) > EXPONENT_MASK;
 122 }
 123
 124 static inline uint64_t
 125 get_mantissa (uint64_t x)
 126 {
 127   return x & MANTISSA_MASK;
 128 }
 129
 130 /* Convert integer number X, unbiased exponent EP, and sign S to double:
 131
 132    result = X * 2^(EP+1 - exponent_bias)
 133
 134    NB: zero is not supported.  */
 135 static inline double
 136 make_double (uint64_t x, int64_t ep, uint64_t s)
 137 {
 138   int lz = clz_uint64 (x) - EXPONENT_WIDTH;
 139   x <<= lz;
 140   ep -= lz;
 141
 142   if (__glibc_unlikely (ep < 0 || x == 0))
 143     {
 144       x >>= -ep;
 145       ep = 0;
 146     }
 147
 148   return asdouble (s + x + (ep << MANTISSA_WIDTH));
 149 }
 150
 151
 152 #define NOINLINE __attribute__ ((noinline))
 153
 154 /* Error handling tail calls for special cases, with a sign argument.
 155    The sign of the return value is set if the argument is non-zero.  */
 156
 157 /* The result overflows.  */
 158 attribute_hidden double __math_oflow (uint32_t);
 159 /* The result underflows to 0 in nearest rounding mode.  */
 160 attribute_hidden double __math_uflow (uint32_t);
 161 /* The result underflows to 0 in some directed rounding mode only.  */
 162 attribute_hidden double __math_may_uflow (uint32_t);
 163 /* Division by zero.  */
 164 attribute_hidden double __math_divzero (uint32_t);
 165
 166 /* Error handling using input checking.  */
 167
 168 /* Invalid input unless it is a quiet NaN.  */
 169 attribute_hidden double __math_invalid (double);
 170
 171 /* Error handling using output checking, only for errno setting.  */
 172
 173 /* Check if the result overflowed to infinity.  */
 174 attribute_hidden double __math_check_oflow (double);
 175 /* Check if the result underflowed to 0.  */
 176 attribute_hidden double __math_check_uflow (double);
 177
 178 /* Check if the result overflowed to infinity.  */
 179 static inline double
 180 check_oflow (double x)
 181 {
 182   return WANT_ERRNO ? __math_check_oflow (x) : x;
 183 }
 184
 185 /* Check if the result underflowed to 0.  */
 186 static inline double
 187 check_uflow (double x)
 188 {
 189   return WANT_ERRNO ? __math_check_uflow (x) : x;
 190 }
 191
 192 #define EXP_TABLE_BITS 7
 193 #define EXP_POLY_ORDER 5
 194 #define EXP2_POLY_ORDER 5
 195 extern const struct exp_data
 196 {
 197   double invln2N;
 198   double shift;
 199   double negln2hiN;
 200   double negln2loN;
 201   double poly[4]; /* Last four coefficients.  */
 202   double exp2_shift;
 203   double exp2_poly[EXP2_POLY_ORDER];
 204   uint64_t tab[2*(1 << EXP_TABLE_BITS)];
 205 } __exp_data attribute_hidden;
 206
 207 #define LOG_TABLE_BITS 7
 208 #define LOG_POLY_ORDER 6
 209 #define LOG_POLY1_ORDER 12
 210 extern const struct log_data
 211 {
 212   double ln2hi;
 213   double ln2lo;
 214   double poly[LOG_POLY_ORDER - 1]; /* First coefficient is 1.  */
 215   double poly1[LOG_POLY1_ORDER - 1];
 216   /* See e_log_data.c for details.  */
 217   struct {double invc, logc;} tab[1 << LOG_TABLE_BITS];
 218 #ifndef __FP_FAST_FMA
 219   struct {double chi, clo;} tab2[1 << LOG_TABLE_BITS];
 220 #endif
 221 } __log_data attribute_hidden;
 222
 223 #define LOG2_TABLE_BITS 6
 224 #define LOG2_POLY_ORDER 7
 225 #define LOG2_POLY1_ORDER 11
 226 extern const struct log2_data
 227 {
 228   double invln2hi;
 229   double invln2lo;
 230   double poly[LOG2_POLY_ORDER - 1];
 231   double poly1[LOG2_POLY1_ORDER - 1];
 232   /* See e_log2_data.c for details.  */
 233   struct {double invc, logc;} tab[1 << LOG2_TABLE_BITS];
 234 #ifndef __FP_FAST_FMA
 235   struct {double chi, clo;} tab2[1 << LOG2_TABLE_BITS];
 236 #endif
 237 } __log2_data attribute_hidden;
 238
 239 #define POW_LOG_TABLE_BITS 7
 240 #define POW_LOG_POLY_ORDER 8
 241 extern const struct pow_log_data
 242 {
 243   double ln2hi;
 244   double ln2lo;
 245   double poly[POW_LOG_POLY_ORDER - 1]; /* First coefficient is 1.  */
 246   /* Note: the pad field is unused, but allows slightly faster indexing.  */
 247   /* See e_pow_log_data.c for details.  */
 248   struct {double invc, pad, logc, logctail;} tab[1 << POW_LOG_TABLE_BITS];
 249 } __pow_log_data attribute_hidden;
 250
 251 #endif