include/flint/double_interval.h

0001 /*
0002     Copyright (C) 2021 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_INTERVAL_H
0013 #define DOUBLE_INTERVAL_H
0014
0015 #ifdef DOUBLE_INTERVAL_INLINES_C
0016 #define DOUBLE_INTERVAL_INLINE
0017 #else
0018 #define DOUBLE_INTERVAL_INLINE static inline
0019 #endif
0020
0021 #include <math.h>
0022 #include "double_extras.h"
0023 #include "arb_types.h"
0024
0025 #ifdef __cplusplus
0026 extern "C" {
0027 #endif
0028
0029 typedef struct
0030 {
0031     double a;
0032     double b;
0033 }
0034 di_t;
0035
0036 #define DI_CHECK(__x) \
0037     if (!(__x.a <= __x.b)) \
0038     { \
0039         flint_throw(FLINT_ERROR, "di_t endpoints %g, %g not ordered\n", __x.a, __x.b); \
0040     } \
0041
0042 DOUBLE_INTERVAL_INLINE
0043 di_t di_interval(double a, double b)
0044 {
0045     di_t res;
0046
0047     if (!(a <= b))
0048     {
0049         flint_throw(FLINT_ERROR, "di_interval endpoints %g, %g not ordered\n", a, b);
0050     }
0051
0052     res.a = a;
0053     res.b = b;
0054     return res;
0055 }
0056
0057 DOUBLE_INTERVAL_INLINE
0058 double _di_below(double x)
0059 {
0060     double t;
0061
0062     if (x <= 1e300)
0063     {
0064         t = x;
0065         if (t < 0.0)
0066             t = -t;
0067
0068         t += 1e-300;
0069         return x - t * 4.440892098500626e-16;
0070     }
0071     else
0072     {
0073         if (x != x)
0074             return -D_INF;
0075
0076         return 1e300;
0077     }
0078 }
0079
0080 DOUBLE_INTERVAL_INLINE
0081 double _di_above(double x)
0082 {
0083     double t;
0084
0085     if (x >= -1e300)
0086     {
0087         t = x;
0088         if (t < 0.0)
0089             t = -t;
0090
0091         t += 1e-300;
0092         return x + t * 4.440892098500626e-16;
0093     }
0094     else
0095     {
0096         if (x != x)
0097             return D_INF;
0098
0099         return -1e300;
0100     }
0101 }
0102
0103 DOUBLE_INTERVAL_INLINE
0104 di_t di_neg(di_t x)
0105 {
0106     di_t res;
0107     res.a = -x.b;
0108     res.b = -x.a;
0109     return res;
0110 }
0111
0112 DOUBLE_INTERVAL_INLINE
0113 di_t di_fast_add(di_t x, di_t y)
0114 {
0115     di_t res;
0116     res.a = _di_below(x.a + y.a);
0117     res.b = _di_above(x.b + y.b);
0118     return res;
0119 }
0120
0121 DOUBLE_INTERVAL_INLINE
0122 di_t di_fast_sub(di_t x, di_t y)
0123 {
0124     di_t res;
0125     res.a = _di_below(x.a - y.b);
0126     res.b = _di_above(x.b - y.a);
0127     return res;
0128 }
0129
0130 di_t di_fast_mul(di_t x, di_t y);
0131 di_t di_fast_sqr(di_t x);
0132 di_t di_fast_div(di_t x, di_t y);
0133
0134 DOUBLE_INTERVAL_INLINE
0135 di_t di_fast_add_d(di_t x, double y)
0136 {
0137     return di_fast_add(x, di_interval(y, y));
0138 }
0139
0140 DOUBLE_INTERVAL_INLINE
0141 di_t di_fast_sub_d(di_t x, double y)
0142 {
0143     return di_fast_sub(x, di_interval(y, y));
0144 }
0145
0146 DOUBLE_INTERVAL_INLINE
0147 di_t di_fast_mul_d(di_t x, double y)
0148 {
0149     return di_fast_mul(x, di_interval(y, y));
0150 }
0151
0152 DOUBLE_INTERVAL_INLINE
0153 di_t di_fast_div_d(di_t x, double y)
0154 {
0155     return di_fast_div(x, di_interval(y, y));
0156 }
0157
0158 di_t di_fast_log_nonnegative(di_t x);
0159
0160 DOUBLE_INTERVAL_INLINE
0161 di_t di_fast_mid(di_t x)
0162 {
0163     di_t a, b;
0164     if (x.a == -D_INF || x.b == D_INF)
0165         return di_interval(-D_INF, D_INF);
0166     a = di_interval(x.a, x.a);
0167     b = di_interval(x.b, x.b);
0168     return di_fast_mul_d(di_fast_add(a, b), 0.5);
0169 }
0170
0171 DOUBLE_INTERVAL_INLINE
0172 double di_fast_ubound_radius(di_t x)
0173 {
0174     return _di_above((x.b - x.a) * 0.5);
0175 }
0176
0177 DOUBLE_INTERVAL_INLINE
0178 void di_print(di_t x)
0179 {
0180     flint_printf("[%.17g, %.17g]", x.a, x.b);
0181 }
0182
0183 di_t arb_get_di(const arb_t x);
0184 void arb_set_di(arb_t res, di_t x, slong prec);
0185
0186 DOUBLE_INTERVAL_INLINE
0187 double d_randtest2(flint_rand_t state)
0188 {
0189     double x;
0190
0191     x = d_randtest(state);
0192     if (n_randint(state, 2))
0193         x = -x;
0194
0195     return ldexp(x, n_randint(state, 2400) - 1200);
0196 }
0197
0198 DOUBLE_INTERVAL_INLINE
0199 di_t di_randtest(flint_rand_t state)
0200 {
0201     di_t res;
0202
0203     res.a = d_randtest2(state);
0204     res.b = d_randtest2(state);
0205
0206     if (res.a > res.b)
0207     {
0208         double t = res.a;
0209         res.a = res.b;
0210         res.b = t;
0211     }
0212
0213     return res;
0214 }
0215
0216 #ifdef __cplusplus
0217 }
0218 #endif
0219
0220 #endif