From 3345c7839cbec4d53b87be4e506733bb0c5b71c2 Mon Sep 17 00:00:00 2001 From: fredrik Date: Wed, 8 Dec 2021 11:12:52 +0100 Subject: [PATCH] helpers for double interval arithmetic --- Makefile.in | 2 +- doc/source/double_interval.rst | 98 +++++++++++++ doc/source/index.rst | 1 + double_interval.h | 193 +++++++++++++++++++++++++ double_interval/arb_get_di.c | 34 +++++ double_interval/arb_set_di.c | 25 ++++ double_interval/fast_div.c | 64 ++++++++ double_interval/fast_log_nonnegative.c | 30 ++++ double_interval/fast_mul.c | 64 ++++++++ double_interval/fast_sqr.c | 41 ++++++ double_interval/inlines.c | 14 ++ double_interval/test/t-fast_add.c | 75 ++++++++++ double_interval/test/t-fast_div.c | 114 +++++++++++++++ double_interval/test/t-fast_mul.c | 114 +++++++++++++++ 14 files changed, 868 insertions(+), 1 deletion(-) create mode 100644 doc/source/double_interval.rst create mode 100644 double_interval.h create mode 100644 double_interval/arb_get_di.c create mode 100644 double_interval/arb_set_di.c create mode 100644 double_interval/fast_div.c create mode 100644 double_interval/fast_log_nonnegative.c create mode 100644 double_interval/fast_mul.c create mode 100644 double_interval/fast_sqr.c create mode 100644 double_interval/inlines.c create mode 100644 double_interval/test/t-fast_add.c create mode 100644 double_interval/test/t-fast_div.c create mode 100644 double_interval/test/t-fast_mul.c diff --git a/Makefile.in b/Makefile.in index 5bcc72e2..89d52a08 100644 --- a/Makefile.in +++ b/Makefile.in @@ -15,7 +15,7 @@ AT=@ BUILD_DIRS = fmpr arf mag arb arb_mat arb_poly arb_calc acb acb_mat acb_poly \ acb_dft acb_calc acb_hypgeom acb_elliptic acb_modular dirichlet acb_dirichlet \ arb_hypgeom bernoulli hypgeom fmpz_extras bool_mat partitions dlog \ - arb_fmpz_poly arb_fpwrap \ + double_interval arb_fmpz_poly arb_fpwrap \ $(EXTRA_BUILD_DIRS) TEMPLATE_DIRS = diff --git a/doc/source/double_interval.rst b/doc/source/double_interval.rst new file mode 100644 index 00000000..c53e5724 --- /dev/null +++ b/doc/source/double_interval.rst @@ -0,0 +1,98 @@ +.. _double_interval: + +**double_interval.h** -- double-precision interval arithmetic and helpers +=============================================================================== + +This module provides helper functions for computing fast enclosures +using ``double`` arithmetic. + + +Types, macros and constants +------------------------------------------------------------------------------- + +.. type:: di_t + + Holds two ``double`` endpoints ``a`` and ``b`` representing + the extended real interval `[a, b]`. We generally assume that + `a \le b` and that neither endpoint is NaN. + +Basic manipulation +------------------------------------------------------------------------------- + +.. function:: di_t di_interval(double a, double b) + + Returns the interval `[a, b]`. We require that the endpoints + are ordered and not NaN. + +.. function:: di_t arb_get_di(const arb_t x) + + Returns the ball *x* converted to a double-precision interval. + +.. function:: void arb_set_di(arb_t res, di_t x, slong prec) + + Sets the ball *res* to the double-precision interval *x*, + rounded to *prec* bits. + +.. function:: void di_print(di_t x) + + Prints *x* to standard output. This simply prints decimal + representations of the floating-point endpoints; the + decimals are not guaranteed to be rounded outward. + +.. function:: double d_randtest2(flint_rand_t state) + + Returns a random non-NaN ``double`` with any exponent. + The value can be infinite or subnormal. + +.. function:: di_t di_randtest(flint_rand_t state) + + Returns an interval with random endpoints. + +Arithmetic +------------------------------------------------------------------------------- + +.. function:: di_t di_neg(di_t x) + + Returns the exact negation of *x*. + +Fast arithmetic +------------------------------------------------------------------------------- + +The following methods perform fast but sloppy interval arithmetic: +we manipulate the endpoints with default rounding and then add +or subtract generic perturbations regardless of whether the +operations were exact. +It is currently assumed that the CPU rounding mode is to nearest. + +.. function:: di_t di_fast_add(di_t x, di_t y) + di_t di_fast_sub(di_t x, di_t y) + di_t di_fast_mul(di_t x, di_t y) + di_t di_fast_div(di_t x, di_t y) + + Returns the sum, difference, product or quotient of *x* and *y*. + Division by zero is currently defined to return `[-\infty, +\infty]`. + +.. function:: di_t di_fast_sqr(di_t x) + + Returns the square of *x*. The output is clamped to + be nonnegative. + +.. function:: di_t di_fast_add_d(di_t x, double y) + di_t di_fast_sub_d(di_t x, double y) + di_t di_fast_mul_d(di_t x, double y) + di_t di_fast_div_d(di_t x, double y) + + Arithmetic with an exact ``double`` operand. + +.. function:: di_t di_fast_log_nonnegative(di_t x) + + Returns an enclosure of `\log(x)`. The lower endpoint of *x* + is rounded up to 0 if it is negative. + +.. function:: di_t di_fast_mid(di_t x) + + Returns an enclosure of the midpoint of *x*. + +.. function:: double di_fast_ubound_radius(di_t x) + + Returns an upper bound for the radius of *x*. diff --git a/doc/source/index.rst b/doc/source/index.rst index 5829c9de..3b94791a 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -217,6 +217,7 @@ Mainly for internal use. .. toctree:: :maxdepth: 1 + double_interval.rst fmpz_extras.rst bool_mat.rst dlog.rst diff --git a/double_interval.h b/double_interval.h new file mode 100644 index 00000000..032cea39 --- /dev/null +++ b/double_interval.h @@ -0,0 +1,193 @@ +/* + Copyright (C) 2021 Fredrik Johansson + + This file is part of Arb. + + Arb is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. See . +*/ + +#ifndef DOUBLE_INTERVAL_H +#define DOUBLE_INTERVAL_H + +#ifdef DOUBLE_INTERVAL_INLINES_C +#define DOUBLE_INTERVAL_INLINE +#else +#define DOUBLE_INTERVAL_INLINE static __inline__ +#endif + +#include "flint/double_extras.h" +#include "arb.h" + +#ifdef __cplusplus +extern "C" { +#endif + +typedef struct +{ + double a; + double b; +} +di_t; + +DOUBLE_INTERVAL_INLINE +di_t di_interval(double a, double b) +{ + di_t res; + + if (!(a < b)) + { + flint_printf("di_interval endpoints not ordered\n"); + flint_abort(); + } + + res.a = a; + res.b = b; + return res; +} + +DOUBLE_INTERVAL_INLINE +double _di_below(double x) +{ + double t; + + t = x; + if (t < 0.0) + t = -t; + + t += 1e-300; + return x - t * 4.440892098500626e-16; +} + +DOUBLE_INTERVAL_INLINE +double _di_above(double x) +{ + double t; + + t = x; + if (t < 0.0) + t = -t; + + t += 1e-300; + return x + t * 4.440892098500626e-16; +} + +DOUBLE_INTERVAL_INLINE +di_t di_neg(di_t x) +{ + di_t res; + res.a = -x.b; + res.b = -x.a; + return res; +} + +DOUBLE_INTERVAL_INLINE +di_t di_fast_add(di_t x, di_t y) +{ + di_t res; + res.a = _di_below(x.a + y.a); + res.b = _di_above(x.b + y.b); + return res; +} + +DOUBLE_INTERVAL_INLINE +di_t di_fast_sub(di_t x, di_t y) +{ + di_t res; + res.a = _di_below(x.a - y.b); + res.b = _di_above(x.b - y.a); + return res; +} + +di_t di_fast_mul(di_t x, di_t y); +di_t di_fast_sqr(di_t x); +di_t di_fast_div(di_t x, di_t y); + +DOUBLE_INTERVAL_INLINE +di_t di_fast_add_d(di_t x, double y) +{ + return di_fast_add(x, di_interval(y, y)); +} + +DOUBLE_INTERVAL_INLINE +di_t di_fast_sub_d(di_t x, double y) +{ + return di_fast_sub(x, di_interval(y, y)); +} + +DOUBLE_INTERVAL_INLINE +di_t di_fast_mul_d(di_t x, double y) +{ + return di_fast_mul(x, di_interval(y, y)); +} + +DOUBLE_INTERVAL_INLINE +di_t di_fast_div_d(di_t x, double y) +{ + return di_fast_div(x, di_interval(y, y)); +} + +di_t di_fast_log_nonnegative(di_t x); + +DOUBLE_INTERVAL_INLINE +di_t di_fast_mid(di_t x) +{ + di_t a, b; + a = di_interval(x.a, x.a); + b = di_interval(x.b, x.b); + return di_fast_mul_d(di_fast_add(a, b), 0.5); +} + +DOUBLE_INTERVAL_INLINE +double di_fast_ubound_radius(di_t x) +{ + return _di_above((x.b - x.a) * 0.5); +} + +DOUBLE_INTERVAL_INLINE +void di_print(di_t x) +{ + flint_printf("[%.17g, %.17g]", x.a, x.b); +} + +di_t arb_get_di(const arb_t x); +void arb_set_di(arb_t res, di_t x, slong prec); + +DOUBLE_INTERVAL_INLINE +double d_randtest2(flint_rand_t state) +{ + double x; + + x = d_randtest(state); + if (n_randint(state, 2)) + x = -x; + + return ldexp(x, n_randint(state, 2400) - 1200); +} + +DOUBLE_INTERVAL_INLINE +di_t di_randtest(flint_rand_t state) +{ + di_t res; + + res.a = d_randtest2(state); + res.b = d_randtest2(state); + + if (res.a > res.b) + { + double t = res.a; + res.a = res.b; + res.b = t; + } + + return res; +} + +#ifdef __cplusplus +} +#endif + +#endif + diff --git a/double_interval/arb_get_di.c b/double_interval/arb_get_di.c new file mode 100644 index 00000000..eaa81de0 --- /dev/null +++ b/double_interval/arb_get_di.c @@ -0,0 +1,34 @@ +/* + Copyright (C) 2021 Fredrik Johansson + + This file is part of Arb. + + Arb is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "double_interval.h" + +di_t +arb_get_di(const arb_t x) +{ + di_t res; + if (arf_is_nan(arb_midref(x))) + { + res.a = -D_INF; + res.b = D_INF; + } + else + { + arf_t t; + arf_init(t); + arb_get_lbound_arf(t, x, 53); + res.a = arf_get_d(t, ARF_RND_FLOOR); + arb_get_ubound_arf(t, x, 53); + res.b = arf_get_d(t, ARF_RND_CEIL); + arf_clear(t); + } + return res; +} diff --git a/double_interval/arb_set_di.c b/double_interval/arb_set_di.c new file mode 100644 index 00000000..0b91c0cb --- /dev/null +++ b/double_interval/arb_set_di.c @@ -0,0 +1,25 @@ +/* + Copyright (C) 2021 Fredrik Johansson + + This file is part of Arb. + + Arb is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "double_interval.h" + +void +arb_set_di(arb_t res, di_t x, slong prec) +{ + arf_t t, u; + arf_init(t); + arf_init(u); + arf_set_d(t, x.a); + arf_set_d(u, x.b); + arb_set_interval_arf(res, t, u, prec); + arf_clear(t); + arf_clear(u); +} diff --git a/double_interval/fast_div.c b/double_interval/fast_div.c new file mode 100644 index 00000000..e1257ba4 --- /dev/null +++ b/double_interval/fast_div.c @@ -0,0 +1,64 @@ +/* + Copyright (C) 2021 Fredrik Johansson + + This file is part of Arb. + + Arb is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "double_interval.h" + +di_t di_fast_div(di_t x, di_t y) +{ + di_t res; + + if (y.a > 0) + { + if (x.a >= 0) + { + res.a = x.a / y.b; + res.b = x.b / y.a; + } + else if (x.b <= 0) + { + res.a = x.a / y.a; + res.b = x.b / y.b; + } + else + { + res.a = x.a / y.a; + res.b = x.b / y.a; + } + } + else if (y.b < 0) + { + if (x.a >= 0) + { + res.a = x.b / y.b; + res.b = x.a / y.a; + } + else if (x.b <= 0) + { + res.a = x.b / y.a; + res.b = x.a / y.b; + } + else + { + res.a = x.b / y.b; + res.b = x.a / y.b; + } + } + else + { + res.a = -D_INF; + res.b = D_INF; + } + + res.a = _di_below(res.a); + res.b = _di_above(res.b); + + return res; +} diff --git a/double_interval/fast_log_nonnegative.c b/double_interval/fast_log_nonnegative.c new file mode 100644 index 00000000..a459c67b --- /dev/null +++ b/double_interval/fast_log_nonnegative.c @@ -0,0 +1,30 @@ +/* + Copyright (C) 2021 Fredrik Johansson + + This file is part of Arb. + + Arb is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "double_interval.h" +#include "mag.h" + +double mag_d_log_lower_bound(double x); +double mag_d_log_upper_bound(double x); + +di_t di_fast_log_nonnegative(di_t x) +{ + di_t res; + + if (x.a <= 0.0) + res.a = -D_INF; + else + res.a = mag_d_log_lower_bound(x.a); + + res.b = mag_d_log_upper_bound(x.b); + + return res; +} diff --git a/double_interval/fast_mul.c b/double_interval/fast_mul.c new file mode 100644 index 00000000..6da661f8 --- /dev/null +++ b/double_interval/fast_mul.c @@ -0,0 +1,64 @@ +/* + Copyright (C) 2021 Fredrik Johansson + + This file is part of Arb. + + Arb is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "double_interval.h" + +di_t di_fast_mul(di_t x, di_t y) +{ + di_t res; + + 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) + { + res.a = x.b * y.a; + res.b = x.a * y.b; + } + 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) + { + res.a = x.b * y.b; + res.b = x.a * y.a; + } + else + { + double a, b, c, d; + + a = x.a * y.a; + b = x.a * y.b; + c = x.b * y.a; + d = x.b * y.b; + + if (a != a || b != b || c != c || d != d) + { + res.a = -D_INF; + res.b = D_INF; + } + else + { + res.a = FLINT_MIN(FLINT_MIN(a, b), FLINT_MIN(c, d)); + res.b = FLINT_MAX(FLINT_MAX(a, b), FLINT_MAX(c, d)); + } + } + + res.a = _di_below(res.a); + res.b = _di_above(res.b); + + return res; +} + diff --git a/double_interval/fast_sqr.c b/double_interval/fast_sqr.c new file mode 100644 index 00000000..a5aa5fb6 --- /dev/null +++ b/double_interval/fast_sqr.c @@ -0,0 +1,41 @@ +/* + Copyright (C) 2021 Fredrik Johansson + + This file is part of Arb. + + Arb is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "double_interval.h" + +di_t di_fast_sqr(di_t x) +{ + di_t res; + + if (x.a >= 0) + { + res.a = x.a * x.a; + res.b = x.b * x.b; + } + else if (x.b <= 0) + { + res.a = x.b * x.b; + res.b = x.a * x.a; + } + else + { + res.a = 0.0; + res.b = FLINT_MAX(x.a * x.a, x.b * x.b); + } + + if (res.a != 0.0) + res.a = _di_below(res.a); + + res.b = _di_above(res.b); + + return res; +} + diff --git a/double_interval/inlines.c b/double_interval/inlines.c new file mode 100644 index 00000000..21f6c667 --- /dev/null +++ b/double_interval/inlines.c @@ -0,0 +1,14 @@ +/* + Copyright (C) 2021 Fredrik Johansson + + This file is part of Arb. + + Arb is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. See . +*/ + +#define DOUBLE_INTERVAL_INLINES_C +#include "double_interval.h" + diff --git a/double_interval/test/t-fast_add.c b/double_interval/test/t-fast_add.c new file mode 100644 index 00000000..51c318eb --- /dev/null +++ b/double_interval/test/t-fast_add.c @@ -0,0 +1,75 @@ +/* + Copyright (C) 2021 Fredrik Johansson + + This file is part of Arb. + + Arb is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "double_interval.h" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("fast_add...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000000 * arb_test_multiplier(); iter++) + { + di_t x, y, z; + arf_t a, b, t, za, zb; + + arf_init(a); + arf_init(b); + arf_init(t); + arf_init(za); + arf_init(zb); + + x = di_randtest(state); + y = di_randtest(state); + + z = di_fast_add(x, y); + + arf_set_d(a, x.a); + arf_set_d(t, y.a); + arf_add(a, a, t, ARF_PREC_EXACT, ARF_RND_DOWN); + arf_set_d(b, x.b); + arf_set_d(t, y.b); + arf_add(b, b, t, ARF_PREC_EXACT, ARF_RND_DOWN); + + arf_set_d(za, z.a); + arf_set_d(zb, z.b); + + if (arf_cmp(a, za) < 0 || arf_cmp(b, zb) > 0) + { + flint_printf("FAIL\n"); + flint_printf("x = "); di_print(x); printf("\n"); + flint_printf("y = "); di_print(y); printf("\n"); + flint_printf("z = "); di_print(z); printf("\n"); + flint_printf("a = "); arf_printd(a, 20); printf("\n"); + flint_printf("b = "); arf_printd(b, 20); printf("\n"); + flint_printf("za = "); arf_printd(za, 20); printf("\n"); + flint_printf("zb = "); arf_printd(zb, 20); printf("\n"); + flint_abort(); + } + + arf_clear(a); + arf_clear(b); + arf_clear(t); + arf_clear(za); + arf_clear(zb); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/double_interval/test/t-fast_div.c b/double_interval/test/t-fast_div.c new file mode 100644 index 00000000..85d07bc7 --- /dev/null +++ b/double_interval/test/t-fast_div.c @@ -0,0 +1,114 @@ +/* + Copyright (C) 2021 Fredrik Johansson + + This file is part of Arb. + + Arb is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "double_interval.h" + +void +arf_min2(arf_t res, const arf_t x, const arf_t y) +{ + if (arf_is_nan(x)) + arf_set(res, y); + else if (arf_is_nan(y)) + arf_set(res, x); + else + arf_min(res, x, y); +} + +void +arf_max2(arf_t res, const arf_t x, const arf_t y) +{ + if (arf_is_nan(x)) + arf_set(res, y); + else if (arf_is_nan(y)) + arf_set(res, x); + else + arf_max(res, x, y); +} + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("fast_div...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000000 * arb_test_multiplier(); iter++) + { + di_t x, y, z; + arf_t a, b, c, d, ra, rb, rc, rd, va, vb, za, zb; + + arf_init(a); arf_init(b); arf_init(c); arf_init(d); + arf_init(ra); arf_init(rb); arf_init(rc); arf_init(rd); + arf_init(va); arf_init(vb); + arf_init(za); + arf_init(zb); + + x = di_randtest(state); + y = di_randtest(state); + z = di_fast_div(x, y); + + arf_set_d(a, x.a); + arf_set_d(b, x.b); + arf_set_d(c, y.a); + arf_set_d(d, y.b); + + arf_div(ra, a, c, 106, ARF_RND_DOWN); + arf_div(rb, a, d, 106, ARF_RND_DOWN); + arf_div(rc, b, c, 106, ARF_RND_DOWN); + arf_div(rd, b, d, 106, ARF_RND_DOWN); + + arf_min2(va, ra, rb); + arf_min2(va, va, rc); + arf_min2(va, va, rd); + arf_max2(vb, ra, rb); + arf_max2(vb, vb, rc); + arf_max2(vb, vb, rd); + + arf_set_d(za, z.a); + arf_set_d(zb, z.b); + + if (arf_cmp(va, za) < 0 || arf_cmp(vb, zb) > 0) + { + flint_printf("FAIL\n"); + flint_printf("x = "); di_print(x); printf("\n"); + flint_printf("y = "); di_print(y); printf("\n"); + flint_printf("z = "); di_print(z); printf("\n"); + flint_printf("a = "); arf_printd(a, 20); printf("\n"); + flint_printf("b = "); arf_printd(b, 20); printf("\n"); + flint_printf("c = "); arf_printd(c, 20); printf("\n"); + flint_printf("d = "); arf_printd(d, 20); printf("\n"); + flint_printf("ra = "); arf_printd(ra, 20); printf("\n"); + flint_printf("rb = "); arf_printd(rb, 20); printf("\n"); + flint_printf("rc = "); arf_printd(rc, 20); printf("\n"); + flint_printf("rd = "); arf_printd(rd, 20); printf("\n"); + flint_printf("va = "); arf_printd(va, 20); printf("\n"); + flint_printf("vb = "); arf_printd(vb, 20); printf("\n"); + flint_printf("za = "); arf_printd(za, 20); printf("\n"); + flint_printf("zb = "); arf_printd(zb, 20); printf("\n"); + flint_abort(); + } + + arf_clear(a); arf_clear(b); arf_clear(c); arf_clear(d); + arf_clear(ra); arf_clear(rb); arf_clear(rc); arf_clear(rd); + arf_clear(va); arf_clear(vb); + arf_clear(za); + arf_clear(zb); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/double_interval/test/t-fast_mul.c b/double_interval/test/t-fast_mul.c new file mode 100644 index 00000000..d12d86b1 --- /dev/null +++ b/double_interval/test/t-fast_mul.c @@ -0,0 +1,114 @@ +/* + Copyright (C) 2021 Fredrik Johansson + + This file is part of Arb. + + Arb is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "double_interval.h" + +void +arf_min2(arf_t res, const arf_t x, const arf_t y) +{ + if (arf_is_nan(x)) + arf_set(res, y); + else if (arf_is_nan(y)) + arf_set(res, x); + else + arf_min(res, x, y); +} + +void +arf_max2(arf_t res, const arf_t x, const arf_t y) +{ + if (arf_is_nan(x)) + arf_set(res, y); + else if (arf_is_nan(y)) + arf_set(res, x); + else + arf_max(res, x, y); +} + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("fast_mul...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000000 * arb_test_multiplier(); iter++) + { + di_t x, y, z; + arf_t a, b, c, d, ra, rb, rc, rd, va, vb, za, zb; + + arf_init(a); arf_init(b); arf_init(c); arf_init(d); + arf_init(ra); arf_init(rb); arf_init(rc); arf_init(rd); + arf_init(va); arf_init(vb); + arf_init(za); + arf_init(zb); + + x = di_randtest(state); + y = di_randtest(state); + z = di_fast_mul(x, y); + + arf_set_d(a, x.a); + arf_set_d(b, x.b); + arf_set_d(c, y.a); + arf_set_d(d, y.b); + + arf_mul(ra, a, c, ARF_PREC_EXACT, ARF_RND_DOWN); + arf_mul(rb, a, d, ARF_PREC_EXACT, ARF_RND_DOWN); + arf_mul(rc, b, c, ARF_PREC_EXACT, ARF_RND_DOWN); + arf_mul(rd, b, d, ARF_PREC_EXACT, ARF_RND_DOWN); + + arf_min2(va, ra, rb); + arf_min2(va, va, rc); + arf_min2(va, va, rd); + arf_max2(vb, ra, rb); + arf_max2(vb, vb, rc); + arf_max2(vb, vb, rd); + + arf_set_d(za, z.a); + arf_set_d(zb, z.b); + + if (arf_cmp(va, za) < 0 || arf_cmp(vb, zb) > 0) + { + flint_printf("FAIL\n"); + flint_printf("x = "); di_print(x); printf("\n"); + flint_printf("y = "); di_print(y); printf("\n"); + flint_printf("z = "); di_print(z); printf("\n"); + flint_printf("a = "); arf_printd(a, 20); printf("\n"); + flint_printf("b = "); arf_printd(b, 20); printf("\n"); + flint_printf("c = "); arf_printd(c, 20); printf("\n"); + flint_printf("d = "); arf_printd(d, 20); printf("\n"); + flint_printf("ra = "); arf_printd(ra, 20); printf("\n"); + flint_printf("rb = "); arf_printd(rb, 20); printf("\n"); + flint_printf("rc = "); arf_printd(rc, 20); printf("\n"); + flint_printf("rd = "); arf_printd(rd, 20); printf("\n"); + flint_printf("va = "); arf_printd(va, 20); printf("\n"); + flint_printf("vb = "); arf_printd(vb, 20); printf("\n"); + flint_printf("za = "); arf_printd(za, 20); printf("\n"); + flint_printf("zb = "); arf_printd(zb, 20); printf("\n"); + flint_abort(); + } + + arf_clear(a); arf_clear(b); arf_clear(c); arf_clear(d); + arf_clear(ra); arf_clear(rb); arf_clear(rc); arf_clear(rd); + arf_clear(va); arf_clear(vb); + arf_clear(za); + arf_clear(zb); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} +