Merge pull request #125 from argriffing/arb-min-max

ENH: arb min and max

arb.h

@@ -535,6 +535,9 @@ void arb_union(arb_t z, const arb_t x, const arb_t y, slong prec);
 int arb_intersection(arb_t z, const arb_t x, const arb_t y, slong prec);
 void arb_get_rand_fmpq(fmpq_t q, flint_rand_t state, const arb_t x, slong bits);
 
+void arb_min(arb_t z, const arb_t x, const arb_t y, slong prec);
+void arb_max(arb_t z, const arb_t x, const arb_t y, slong prec);
+
 int arb_can_round_arf(const arb_t x, slong prec, arf_rnd_t rnd);
 int arb_can_round_mpfr(const arb_t x, slong prec, mpfr_rnd_t rnd);
 

arb/max.c

@@ -0,0 +1,59 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2016 Arb authors
+
+******************************************************************************/
+
+#include "arb.h"
+
+void
+arb_max(arb_t z, const arb_t x, const arb_t y, slong prec)
+{
+    arf_t left, right, t, xr, yr;
+
+    if (arf_is_nan(arb_midref(x)) || arf_is_nan(arb_midref(y)))
+    {
+        arb_indeterminate(z);
+        return;
+    }
+
+    arf_init(left);
+    arf_init(right);
+    arf_init(t);
+
+    arf_init_set_mag_shallow(xr, arb_radref(x));
+    arf_init_set_mag_shallow(yr, arb_radref(y));
+
+    arf_sub(left, arb_midref(x), xr, prec, ARF_RND_FLOOR);
+    arf_sub(t, arb_midref(y), yr, prec, ARF_RND_FLOOR);
+    arf_max(left, left, t);
+
+    arf_add(right, arb_midref(x), xr, prec, ARF_RND_CEIL);
+    arf_add(t, arb_midref(y), yr, prec, ARF_RND_CEIL);
+    arf_max(right, right, t);
+
+    arb_set_interval_arf(z, left, right, prec);
+
+    arf_clear(left);
+    arf_clear(right);
+    arf_clear(t);
+}

arb/min.c

@@ -0,0 +1,59 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2016 Arb authors
+
+******************************************************************************/
+
+#include "arb.h"
+
+void
+arb_min(arb_t z, const arb_t x, const arb_t y, slong prec)
+{
+    arf_t left, right, t, xr, yr;
+
+    if (arf_is_nan(arb_midref(x)) || arf_is_nan(arb_midref(y)))
+    {
+        arb_indeterminate(z);
+        return;
+    }
+
+    arf_init(left);
+    arf_init(right);
+    arf_init(t);
+
+    arf_init_set_mag_shallow(xr, arb_radref(x));
+    arf_init_set_mag_shallow(yr, arb_radref(y));
+
+    arf_sub(left, arb_midref(x), xr, prec, ARF_RND_FLOOR);
+    arf_sub(t, arb_midref(y), yr, prec, ARF_RND_FLOOR);
+    arf_min(left, left, t);
+
+    arf_add(right, arb_midref(x), xr, prec, ARF_RND_CEIL);
+    arf_add(t, arb_midref(y), yr, prec, ARF_RND_CEIL);
+    arf_min(right, right, t);
+
+    arb_set_interval_arf(z, left, right, prec);
+
+    arf_clear(left);
+    arf_clear(right);
+    arf_clear(t);
+}

arb/test/t-max.c

@@ -0,0 +1,142 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2016 Arb authors
+
+******************************************************************************/
+
+#include "arb.h"
+
+/* sample (x, y) so that x /in y */
+void _sample_arf_in_arb(arf_t x, arb_t y, flint_rand_t state)
+{
+    slong bits, prec, expbits;
+    arf_t a, b;
+    slong i, n;
+
+    arf_init(a);
+    arf_init(b);
+
+    bits = 2 + n_randint(state, 1000);
+    prec = 2 + n_randint(state, 1000);
+    expbits = n_randint(state, 14);
+    n = n_randint(state, 3);
+
+    arf_randtest(x, state, bits, expbits);
+    arf_set(a, x);
+    arf_set(b, x);
+    for (i = 0; i < n; i++)
+    {
+        arf_randtest(x, state, bits, expbits);
+        arf_min(a, a, x);
+        arf_max(b, b, x);
+    }
+
+    arb_set_interval_arf(y, a, b, prec);
+
+    arf_clear(a);
+    arf_clear(b);
+}
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("max....");
+    fflush(stdout);
+    flint_randinit(state);
+
+    for (iter = 0; iter < 100000; iter++)
+    {
+        arf_t a, b, c;
+        arb_t x, y, z;
+        slong prec;
+
+        arf_init(a);
+        arf_init(b);
+        arf_init(c);
+
+        arb_init(x);
+        arb_init(y);
+        arb_init(z);
+
+        _sample_arf_in_arb(a, x, state);
+        _sample_arf_in_arb(b, y, state);
+
+        prec = 2 + n_randint(state, 200);
+
+        arf_max(c, a, b);
+        arb_max(z, x, y, prec);
+
+        if (!arb_contains_arf(x, a) ||
+            !arb_contains_arf(y, b) ||
+            !arb_contains_arf(z, c))
+        {
+            flint_printf("FAIL: containment\n\n");
+            flint_printf("a = "); arf_print(a); flint_printf("\n\n");
+            flint_printf("b = "); arf_print(b); flint_printf("\n\n");
+            flint_printf("c = "); arf_print(c); flint_printf("\n\n");
+            flint_printf("x = "); arb_print(x); flint_printf("\n\n");
+            flint_printf("y = "); arb_print(y); flint_printf("\n\n");
+            flint_printf("z = "); arb_print(z); flint_printf("\n\n");
+            abort();
+        }
+
+        /* aliasing */
+        {
+            int alias;
+
+            if (n_randint(state, 2))
+            {
+                arb_max(x, x, y, prec);
+                alias = arb_equal(x, z);
+            }
+            else
+            {
+                arb_max(y, x, y, prec);
+                alias = arb_equal(y, z);
+            }
+
+            if (!alias)
+            {
+                flint_printf("FAIL: aliasing\n\n");
+                flint_printf("x = "); arb_print(x); flint_printf("\n\n");
+                flint_printf("y = "); arb_print(y); flint_printf("\n\n");
+                flint_printf("z = "); arb_print(z); flint_printf("\n\n");
+                abort();
+            }
+        }
+
+        arf_clear(a);
+        arf_clear(b);
+        arf_clear(c);
+
+        arb_clear(x);
+        arb_clear(y);
+        arb_clear(z);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

arb/test/t-min.c

@@ -0,0 +1,142 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2016 Arb authors
+
+******************************************************************************/
+
+#include "arb.h"
+
+/* sample (x, y) so that x /in y */
+void _sample_arf_in_arb(arf_t x, arb_t y, flint_rand_t state)
+{
+    slong bits, prec, expbits;
+    arf_t a, b;
+    slong i, n;
+
+    arf_init(a);
+    arf_init(b);
+
+    bits = 2 + n_randint(state, 1000);
+    prec = 2 + n_randint(state, 1000);
+    expbits = n_randint(state, 14);
+    n = n_randint(state, 3);
+
+    arf_randtest(x, state, bits, expbits);
+    arf_set(a, x);
+    arf_set(b, x);
+    for (i = 0; i < n; i++)
+    {
+        arf_randtest(x, state, bits, expbits);
+        arf_min(a, a, x);
+        arf_max(b, b, x);
+    }
+
+    arb_set_interval_arf(y, a, b, prec);
+
+    arf_clear(a);
+    arf_clear(b);
+}
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("min....");
+    fflush(stdout);
+    flint_randinit(state);
+
+    for (iter = 0; iter < 100000; iter++)
+    {
+        arf_t a, b, c;
+        arb_t x, y, z;
+        slong prec;
+
+        arf_init(a);
+        arf_init(b);
+        arf_init(c);
+
+        arb_init(x);
+        arb_init(y);
+        arb_init(z);
+
+        _sample_arf_in_arb(a, x, state);
+        _sample_arf_in_arb(b, y, state);
+
+        prec = 2 + n_randint(state, 200);
+
+        arf_min(c, a, b);
+        arb_min(z, x, y, prec);
+
+        if (!arb_contains_arf(x, a) ||
+            !arb_contains_arf(y, b) ||
+            !arb_contains_arf(z, c))
+        {
+            flint_printf("FAIL: containment\n\n");
+            flint_printf("a = "); arf_print(a); flint_printf("\n\n");
+            flint_printf("b = "); arf_print(b); flint_printf("\n\n");
+            flint_printf("c = "); arf_print(c); flint_printf("\n\n");
+            flint_printf("x = "); arb_print(x); flint_printf("\n\n");
+            flint_printf("y = "); arb_print(y); flint_printf("\n\n");
+            flint_printf("z = "); arb_print(z); flint_printf("\n\n");
+            abort();
+        }
+
+        /* aliasing */
+        {
+            int alias;
+
+            if (n_randint(state, 2))
+            {
+                arb_min(x, x, y, prec);
+                alias = arb_equal(x, z);
+            }
+            else
+            {
+                arb_min(y, x, y, prec);
+                alias = arb_equal(y, z);
+            }
+
+            if (!alias)
+            {
+                flint_printf("FAIL: aliasing\n\n");
+                flint_printf("x = "); arb_print(x); flint_printf("\n\n");
+                flint_printf("y = "); arb_print(y); flint_printf("\n\n");
+                flint_printf("z = "); arb_print(z); flint_printf("\n\n");
+                abort();
+            }
+        }
+
+        arf_clear(a);
+        arf_clear(b);
+        arf_clear(c);
+
+        arb_clear(x);
+        arb_clear(y);
+        arb_clear(z);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

doc/source/arb.rst

@@ -619,6 +619,12 @@ Arithmetic
     Sets *y* to the sign function of *x*. The result is `[0 \pm 1]` if
     *x* contains both zero and nonzero numbers.
 
+.. function:: void arb_min(arb_t z, const arb_t x, const arb_t y, slong prec)
+
+.. function:: void arb_max(arb_t z, const arb_t x, const arb_t y, slong prec)
+
+    Sets *z* respectively to the minimum and the maximum of *x* and *y*.
+
 .. function:: void arb_add(arb_t z, const arb_t x, const arb_t y, slong prec)
 
 .. function:: void arb_add_arf(arb_t z, const arb_t x, const arf_t y, slong prec)