include/flint/double_extras.h

0001 /*
0002     Copyright (C) 2012 Fredrik Johansson
0003
0004     This file is part of FLINT.
0005
0006     FLINT is free software: you can redistribute it and/or modify it under
0007     the terms of the GNU Lesser General Public License (LGPL) as published
0008     by the Free Software Foundation; either version 3 of the License, or
0009     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
0010 */
0011
0012 #ifndef DOUBLE_EXTRAS_H
0013 #define DOUBLE_EXTRAS_H
0014
0015 #ifdef DOUBLE_EXTRAS_INLINES_C
0016 #define DOUBLE_EXTRAS_INLINE
0017 #else
0018 #define DOUBLE_EXTRAS_INLINE static inline
0019 #endif
0020
0021 #include <stdint.h>
0022 #include <math.h>
0023 #include "flint.h"
0024
0025 #ifdef __cplusplus
0026  extern "C" {
0027 #endif
0028
0029 #define D_BITS 53
0030 #define D_EPS 2.2204460492503130808e-16
0031 #define D_INF HUGE_VAL
0032 #define D_NAN (HUGE_VAL - HUGE_VAL)
0033
0034 double d_randtest(flint_rand_t state);
0035
0036 double d_randtest_signed(flint_rand_t state, slong minexp, slong maxexp);
0037
0038 double d_randtest_special(flint_rand_t state, slong minexp, slong maxexp);
0039
0040 DOUBLE_EXTRAS_INLINE
0041 double d_polyval(const double * poly, int len, double x)
0042 {
0043     double t;
0044     int i;
0045
0046     for (t = poly[len-1], i = len-2; i >= 0; i--)
0047         t = poly[i] + x * t;
0048
0049     return t;
0050 }
0051
0052 double d_lambertw(double x);
0053
0054 DOUBLE_EXTRAS_INLINE
0055 int d_is_nan(double x)
0056 {
0057     return x != x;
0058 }
0059
0060 double d_log2(double x);
0061
0062 typedef union
0063 {
0064     double f;
0065     uint64_t i;
0066 }
0067 double_uint64_u;
0068
0069 #define D_MIN_NORMAL_EXPONENT -1022
0070 #define D_MAX_NORMAL_EXPONENT 1023
0071 #define D_EXPONENT_BIAS 1023
0072 #define D_EXPONENT_SHIFT 52
0073
0074 /* Assumes that 2^i is in the normal exponent range. */
0075 FLINT_FORCE_INLINE double d_mul_2exp_inrange(double x, int i)
0076 {
0077     FLINT_ASSERT(i >= D_MIN_NORMAL_EXPONENT && i <= D_MAX_NORMAL_EXPONENT);
0078     double_uint64_u u;
0079     u.i = ((uint64_t) (i + D_EXPONENT_BIAS)) << D_EXPONENT_SHIFT;
0080     return x * u.f;
0081 }
0082
0083 /* Assumes that 2^i, x and x*2^i are all in the normal exponent range. */
0084 /* In particular, also assumes x != 0. */
0085 FLINT_FORCE_INLINE double d_mul_2exp_inrange2(double x, int i)
0086 {
0087     FLINT_ASSERT(i >= D_MIN_NORMAL_EXPONENT && i <= D_MAX_NORMAL_EXPONENT);
0088     FLINT_ASSERT(x != 0);
0089
0090     double_uint64_u u;
0091     u.f = x;
0092     u.i += ((uint64_t) i) << D_EXPONENT_SHIFT;
0093     return u.f;
0094 }
0095
0096 FLINT_FORCE_INLINE double d_mul_2exp(double x, int i)
0097 {
0098     if (i >= D_MIN_NORMAL_EXPONENT && i <= D_MAX_NORMAL_EXPONENT)
0099         return d_mul_2exp_inrange(x, i);
0100     else
0101         return ldexp(x, i);
0102 }
0103
0104 #ifdef __cplusplus
0105 }
0106 #endif
0107
0108 #endif