diff --git a/double_interval.h b/double_interval.h index 6eac4d76..fed7e0ea 100644 --- a/double_interval.h +++ b/double_interval.h @@ -32,6 +32,13 @@ typedef struct } di_t; +#define DI_CHECK(__x) \ + if (!(__x.a <= __x.b)) \ + { \ + flint_printf("di_t endpoints %g, %g not ordered\n", __x.a, __x.b); \ + flint_abort(); \ + } \ + DOUBLE_INTERVAL_INLINE di_t di_interval(double a, double b) { @@ -53,12 +60,22 @@ double _di_below(double x) { double t; - t = x; - if (t < 0.0) - t = -t; + if (x <= 1e300) + { + t = x; + if (t < 0.0) + t = -t; - t += 1e-300; - return x - t * 4.440892098500626e-16; + t += 1e-300; + return x - t * 4.440892098500626e-16; + } + else + { + if (x != x) + return -D_INF; + + return 1e300; + } } DOUBLE_INTERVAL_INLINE @@ -66,12 +83,22 @@ double _di_above(double x) { double t; - t = x; - if (t < 0.0) - t = -t; + if (x >= -1e300) + { + t = x; + if (t < 0.0) + t = -t; - t += 1e-300; - return x + t * 4.440892098500626e-16; + t += 1e-300; + return x + t * 4.440892098500626e-16; + } + else + { + if (x != x) + return D_INF; + + return -1e300; + } } DOUBLE_INTERVAL_INLINE diff --git a/double_interval/fast_mul.c b/double_interval/fast_mul.c index 6da661f8..9abf63fe 100644 --- a/double_interval/fast_mul.c +++ b/double_interval/fast_mul.c @@ -15,22 +15,22 @@ di_t di_fast_mul(di_t x, di_t y) { di_t res; - if (x.a >= 0 && y.a >= 0) + if (x.a > 0 && y.a > 0) { res.a = x.a * y.a; res.b = x.b * y.b; } - else if (x.a >= 0 && y.b <= 0) + else if (x.a > 0 && y.b < 0) { res.a = x.b * y.a; res.b = x.a * y.b; } - else if (x.b <= 0 && y.a >= 0) + else if (x.b < 0 && y.a > 0) { res.a = x.a * y.b; res.b = x.b * y.a; } - else if (x.b <= 0 && y.b <= 0) + else if (x.b < 0 && y.b < 0) { res.a = x.b * y.b; res.b = x.a * y.a; diff --git a/double_interval/test/t-fast_add.c b/double_interval/test/t-fast_add.c index 51c318eb..d820eb8d 100644 --- a/double_interval/test/t-fast_add.c +++ b/double_interval/test/t-fast_add.c @@ -37,6 +37,8 @@ int main() z = di_fast_add(x, y); + DI_CHECK(z) + arf_set_d(a, x.a); arf_set_d(t, y.a); arf_add(a, a, t, ARF_PREC_EXACT, ARF_RND_DOWN); diff --git a/double_interval/test/t-fast_div.c b/double_interval/test/t-fast_div.c index 85d07bc7..6b85b9c6 100644 --- a/double_interval/test/t-fast_div.c +++ b/double_interval/test/t-fast_div.c @@ -58,6 +58,8 @@ int main() y = di_randtest(state); z = di_fast_div(x, y); + DI_CHECK(z) + arf_set_d(a, x.a); arf_set_d(b, x.b); arf_set_d(c, y.a); diff --git a/double_interval/test/t-fast_mul.c b/double_interval/test/t-fast_mul.c index d12d86b1..92596839 100644 --- a/double_interval/test/t-fast_mul.c +++ b/double_interval/test/t-fast_mul.c @@ -54,10 +54,23 @@ int main() arf_init(za); arf_init(zb); - x = di_randtest(state); - y = di_randtest(state); + if (iter == 0) + { + x.a = 3.5835915908219665e+103; + x.b = 4.0874812242073031e+295; + y.a = -8.3711609938763647e+298; + y.b = -3.414037107833399e+243; + } + else + { + x = di_randtest(state); + y = di_randtest(state); + } + z = di_fast_mul(x, y); + DI_CHECK(z) + arf_set_d(a, x.a); arf_set_d(b, x.b); arf_set_d(c, y.a);