Home | History | Annotate | Line # | Download | only in libquadmath
      1 /* GCC Quad-Precision Math Library
      2    Copyright (C) 2010, 2011 Free Software Foundation, Inc.
      3    Written by Francois-Xavier Coudert  <fxcoudert (at) gcc.gnu.org>
      4 
      5 This file is part of the libquadmath library.
      6 Libquadmath is free software; you can redistribute it and/or
      7 modify it under the terms of the GNU Library General Public
      8 License as published by the Free Software Foundation; either
      9 version 2 of the License, or (at your option) any later version.
     10 
     11 Libquadmath is distributed in the hope that it will be useful,
     12 but WITHOUT ANY WARRANTY; without even the implied warranty of
     13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     14 Library General Public License for more details.
     15 
     16 You should have received a copy of the GNU Library General Public
     17 License along with libquadmath; see the file COPYING.LIB.  If
     18 not, write to the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
     19 Boston, MA 02110-1301, USA.  */
     20 
     21 #ifndef QUADMATH_IMP_H
     22 #define QUADMATH_IMP_H
     23 
     24 #include <errno.h>
     25 #include <limits.h>
     26 #include <stdbool.h>
     27 #include <stdint.h>
     28 #include <stdlib.h>
     29 #include "quadmath.h"
     30 #include "config.h"
     31 #ifdef HAVE_FENV_H
     32 # include <fenv.h>
     33 #endif
     34 
     35 
     36 /* Under IEEE 754, an architecture may determine tininess of
     37    floating-point results either "before rounding" or "after
     38    rounding", but must do so in the same way for all operations
     39    returning binary results.  Define TININESS_AFTER_ROUNDING to 1 for
     40    "after rounding" architectures, 0 for "before rounding"
     41    architectures.  */
     42 
     43 #define TININESS_AFTER_ROUNDING   1
     44 
     45 #define HIGH_ORDER_BIT_IS_SET_FOR_SNAN 0
     46 
     47 #define FIX_FLT128_LONG_CONVERT_OVERFLOW 0
     48 #define FIX_FLT128_LLONG_CONVERT_OVERFLOW 0
     49 
     50 /* Prototypes for internal functions.  */
     51 extern int32_t __quadmath_rem_pio2q (__float128, __float128 *);
     52 extern void __quadmath_kernel_sincosq (__float128, __float128, __float128 *,
     53 				       __float128 *, int);
     54 extern __float128 __quadmath_kernel_sinq (__float128, __float128, int);
     55 extern __float128 __quadmath_kernel_cosq (__float128, __float128);
     56 extern __float128 __quadmath_kernel_tanq (__float128, __float128, int);
     57 extern __float128 __quadmath_gamma_productq (__float128, __float128, int,
     58 					     __float128 *);
     59 extern __float128 __quadmath_gammaq_r (__float128, int *);
     60 extern __float128 __quadmath_lgamma_negq (__float128, int *);
     61 extern __float128 __quadmath_lgamma_productq (__float128, __float128,
     62 					      __float128, int);
     63 extern __float128 __quadmath_lgammaq_r (__float128, int *);
     64 extern __float128 __quadmath_x2y2m1q (__float128 x, __float128 y);
     65 extern __complex128 __quadmath_kernel_casinhq (__complex128, int);
     66 
     67 static inline void
     68 mul_splitq (__float128 *hi, __float128 *lo, __float128 x, __float128 y)
     69 {
     70   /* Fast built-in fused multiply-add.  */
     71   *hi = x * y;
     72   *lo = fmaq (x, y, -*hi);
     73 }
     74 
     75 
     76 
     77 
     78 /* Frankly, if you have __float128, you have 64-bit integers, right?  */
     79 #ifndef UINT64_C
     80 # error "No way!"
     81 #endif
     82 
     83 
     84 /* Main union type we use to manipulate the floating-point type.  */
     85 typedef union
     86 {
     87   __float128 value;
     88 
     89   struct
     90 #ifdef __MINGW32__
     91   /* On mingw targets the ms-bitfields option is active by default.
     92      Therefore enforce gnu-bitfield style.  */
     93   __attribute__ ((gcc_struct))
     94 #endif
     95   {
     96 #if __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__
     97     unsigned negative:1;
     98     unsigned exponent:15;
     99     unsigned mantissa0:16;
    100     unsigned mantissa1:32;
    101     unsigned mantissa2:32;
    102     unsigned mantissa3:32;
    103 #else
    104     unsigned mantissa3:32;
    105     unsigned mantissa2:32;
    106     unsigned mantissa1:32;
    107     unsigned mantissa0:16;
    108     unsigned exponent:15;
    109     unsigned negative:1;
    110 #endif
    111   } ieee;
    112 
    113   struct
    114   {
    115 #if __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__
    116     uint64_t high;
    117     uint64_t low;
    118 #else
    119     uint64_t low;
    120     uint64_t high;
    121 #endif
    122   } words64;
    123 
    124   struct
    125   {
    126 #if __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__
    127     uint32_t w0;
    128     uint32_t w1;
    129     uint32_t w2;
    130     uint32_t w3;
    131 #else
    132     uint32_t w3;
    133     uint32_t w2;
    134     uint32_t w1;
    135     uint32_t w0;
    136 #endif
    137   } words32;
    138 
    139   struct
    140 #ifdef __MINGW32__
    141   /* Make sure we are using gnu-style bitfield handling.  */
    142   __attribute__ ((gcc_struct))
    143 #endif
    144   {
    145 #if __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__
    146     unsigned negative:1;
    147     unsigned exponent:15;
    148     unsigned quiet_nan:1;
    149     unsigned mantissa0:15;
    150     unsigned mantissa1:32;
    151     unsigned mantissa2:32;
    152     unsigned mantissa3:32;
    153 #else
    154     unsigned mantissa3:32;
    155     unsigned mantissa2:32;
    156     unsigned mantissa1:32;
    157     unsigned mantissa0:15;
    158     unsigned quiet_nan:1;
    159     unsigned exponent:15;
    160     unsigned negative:1;
    161 #endif
    162   } ieee_nan;
    163 
    164 } ieee854_float128;
    165 
    166 
    167 /* Get two 64 bit ints from a long double.  */
    168 #define GET_FLT128_WORDS64(ix0,ix1,d)  \
    169 do {                                   \
    170   ieee854_float128 u;                  \
    171   u.value = (d);                       \
    172   (ix0) = u.words64.high;              \
    173   (ix1) = u.words64.low;               \
    174 } while (0)
    175 
    176 /* Set a long double from two 64 bit ints.  */
    177 #define SET_FLT128_WORDS64(d,ix0,ix1)  \
    178 do {                                   \
    179   ieee854_float128 u;                  \
    180   u.words64.high = (ix0);              \
    181   u.words64.low = (ix1);               \
    182   (d) = u.value;                       \
    183 } while (0)
    184 
    185 /* Get the more significant 64 bits of a long double mantissa.  */
    186 #define GET_FLT128_MSW64(v,d)          \
    187 do {                                   \
    188   ieee854_float128 u;                  \
    189   u.value = (d);                       \
    190   (v) = u.words64.high;                \
    191 } while (0)
    192 
    193 /* Set the more significant 64 bits of a long double mantissa from an int.  */
    194 #define SET_FLT128_MSW64(d,v)          \
    195 do {                                   \
    196   ieee854_float128 u;                  \
    197   u.value = (d);                       \
    198   u.words64.high = (v);                \
    199   (d) = u.value;                       \
    200 } while (0)
    201 
    202 /* Get the least significant 64 bits of a long double mantissa.  */
    203 #define GET_FLT128_LSW64(v,d)          \
    204 do {                                   \
    205   ieee854_float128 u;                  \
    206   u.value = (d);                       \
    207   (v) = u.words64.low;                 \
    208 } while (0)
    209 
    210 
    211 #define IEEE854_FLOAT128_BIAS 0x3fff
    212 
    213 #define QUADFP_NAN		0
    214 #define QUADFP_INFINITE		1
    215 #define QUADFP_ZERO		2
    216 #define QUADFP_SUBNORMAL	3
    217 #define QUADFP_NORMAL		4
    218 #define fpclassifyq(x) \
    219   __builtin_fpclassify (QUADFP_NAN, QUADFP_INFINITE, QUADFP_NORMAL, \
    220 			QUADFP_SUBNORMAL, QUADFP_ZERO, x)
    221 
    222 #ifndef math_opt_barrier
    223 # define math_opt_barrier(x) \
    224 ({ __typeof (x) __x = (x); __asm ("" : "+m" (__x)); __x; })
    225 # define math_force_eval(x) \
    226 ({ __typeof (x) __x = (x); __asm __volatile__ ("" : : "m" (__x)); })
    227 #endif
    228 
    229 /* math_narrow_eval reduces its floating-point argument to the range
    230    and precision of its semantic type.  (The original evaluation may
    231    still occur with excess range and precision, so the result may be
    232    affected by double rounding.)  */
    233 #define math_narrow_eval(x) (x)
    234 
    235 /* If X (which is not a NaN) is subnormal, force an underflow
    236    exception.  */
    237 #define math_check_force_underflow(x)				\
    238   do								\
    239     {								\
    240       __float128 force_underflow_tmp = (x);			\
    241       if (fabsq (force_underflow_tmp) < FLT128_MIN)		\
    242 	{							\
    243 	  __float128 force_underflow_tmp2			\
    244 	    = force_underflow_tmp * force_underflow_tmp;	\
    245 	  math_force_eval (force_underflow_tmp2);		\
    246 	}							\
    247     }								\
    248   while (0)
    249 /* Likewise, but X is also known to be nonnegative.  */
    250 #define math_check_force_underflow_nonneg(x)			\
    251   do								\
    252     {								\
    253       __float128 force_underflow_tmp = (x);			\
    254       if (force_underflow_tmp < FLT128_MIN)			\
    255 	{							\
    256 	  __float128 force_underflow_tmp2			\
    257 	    = force_underflow_tmp * force_underflow_tmp;	\
    258 	  math_force_eval (force_underflow_tmp2);		\
    259 	}							\
    260     }								\
    261   while (0)
    262 
    263 /* Likewise, for both real and imaginary parts of a complex
    264    result.  */
    265 #define math_check_force_underflow_complex(x)				\
    266   do									\
    267     {									\
    268       __typeof (x) force_underflow_complex_tmp = (x);			\
    269       math_check_force_underflow (__real__ force_underflow_complex_tmp); \
    270       math_check_force_underflow (__imag__ force_underflow_complex_tmp); \
    271     }									\
    272   while (0)
    273 
    274 #ifndef HAVE_FENV_H
    275 # define feraiseexcept(arg) ((void) 0)
    276 typedef int fenv_t;
    277 # define feholdexcept(arg) ((void) 0)
    278 # define fesetround(arg) ((void) 0)
    279 # define feupdateenv(arg) ((void) (arg))
    280 # define fesetenv(arg) ((void) (arg))
    281 # define fetestexcept(arg) 0
    282 # define feclearexcept(arg) ((void) 0)
    283 #else
    284 # ifndef HAVE_FEHOLDEXCEPT
    285 #  define feholdexcept(arg) ((void) 0)
    286 # endif
    287 # ifndef HAVE_FESETROUND
    288 #  define fesetround(arg) ((void) 0)
    289 # endif
    290 # ifndef HAVE_FEUPDATEENV
    291 #  define feupdateenv(arg) ((void) (arg))
    292 # endif
    293 # ifndef HAVE_FESETENV
    294 #  define fesetenv(arg) ((void) (arg))
    295 # endif
    296 # ifndef HAVE_FETESTEXCEPT
    297 #  define fetestexcept(arg) 0
    298 # endif
    299 #endif
    300 
    301 #ifndef __glibc_likely
    302 # define __glibc_likely(cond)	__builtin_expect ((cond), 1)
    303 #endif
    304 
    305 #ifndef __glibc_unlikely
    306 # define __glibc_unlikely(cond)	__builtin_expect ((cond), 0)
    307 #endif
    308 
    309 #if defined HAVE_FENV_H && defined HAVE_FESETROUND && defined HAVE_FEUPDATEENV
    310 struct rm_ctx
    311 {
    312   fenv_t env;
    313   bool updated_status;
    314 };
    315 
    316 # define SET_RESTORE_ROUNDF128(RM)					\
    317   struct rm_ctx ctx __attribute__((cleanup (libc_feresetround_ctx)));	\
    318   libc_feholdsetround_ctx (&ctx, (RM))
    319 
    320 static inline __attribute__ ((always_inline)) void
    321 libc_feholdsetround_ctx (struct rm_ctx *ctx, int round)
    322 {
    323   ctx->updated_status = false;
    324 
    325   /* Update rounding mode only if different.  */
    326   if (__glibc_unlikely (round != fegetround ()))
    327     {
    328       ctx->updated_status = true;
    329       fegetenv (&ctx->env);
    330       fesetround (round);
    331     }
    332 }
    333 
    334 static inline __attribute__ ((always_inline)) void
    335 libc_feresetround_ctx (struct rm_ctx *ctx)
    336 {
    337   /* Restore the rounding mode if updated.  */
    338   if (__glibc_unlikely (ctx->updated_status))
    339     feupdateenv (&ctx->env);
    340 }
    341 #else
    342 # define SET_RESTORE_ROUNDF128(RM) ((void) 0)
    343 #endif
    344 
    345 #endif
    346