helpers for double interval arithmetic

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 = 

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*.

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

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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
+

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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;
+}

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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);
+}

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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;
+}

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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;
+}

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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;
+}
+

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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;
+}
+

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 <http://www.gnu.org/licenses/>.
+*/
+
+#define DOUBLE_INTERVAL_INLINES_C
+#include "double_interval.h"
+

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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;
+}
+

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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;
+}
+

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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;
+}
+