From c50370d897db158e87ed1f4d927299abcf49a558 Mon Sep 17 00:00:00 2001
From: Fredrik Johansson <fredrik.johansson@gmail.com>
Date: Mon, 10 Sep 2012 23:50:33 +0200
Subject: [PATCH] add test code for ball arithmetic

---
 fmprb/div.c           |   9 ++-
 fmprb/pow.c           |   8 +-
 fmprb/test/t-add.c    | 168 ++++++++++++++++++++++++++++++++++++++++
 fmprb/test/t-addmul.c | 171 +++++++++++++++++++++++++++++++++++++++++
 fmprb/test/t-div.c    | 174 ++++++++++++++++++++++++++++++++++++++++++
 fmprb/test/t-mul.c    | 168 ++++++++++++++++++++++++++++++++++++++++
 fmprb/test/t-sub.c    | 168 ++++++++++++++++++++++++++++++++++++++++
 fmprb/test/t-submul.c | 171 +++++++++++++++++++++++++++++++++++++++++
 fmprb/zeta_ui.c       |   4 +-
 9 files changed, 1034 insertions(+), 7 deletions(-)
 create mode 100644 fmprb/test/t-add.c
 create mode 100644 fmprb/test/t-addmul.c
 create mode 100644 fmprb/test/t-div.c
 create mode 100644 fmprb/test/t-mul.c
 create mode 100644 fmprb/test/t-sub.c
 create mode 100644 fmprb/test/t-submul.c

diff --git a/fmprb/div.c b/fmprb/div.c
index 8f752aac..1085b93e 100644
--- a/fmprb/div.c
+++ b/fmprb/div.c
@@ -30,6 +30,13 @@ fmprb_div(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec)
 {
     long r;
 
+    if (fmprb_contains_zero(y))
+    {
+        fmpr_zero(fmprb_midref(z));
+        fmpr_pos_inf(fmprb_radref(z));
+        return;
+    }
+
     if (fmprb_is_exact(y))
     {
         if (fmprb_is_exact(x))
@@ -52,7 +59,7 @@ fmprb_div(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec)
     else
     {
         /* x/y - (x+a)/(y-b)  =  a/(b-y) + x*b/(y*(b-y)) */
-        /* XXX: require y > b! */
+        /* where we assume y > b */
 
         fmpr_t t, u, by;
 
diff --git a/fmprb/pow.c b/fmprb/pow.c
index 77276daa..6dfbd836 100644
--- a/fmprb/pow.c
+++ b/fmprb/pow.c
@@ -30,7 +30,7 @@
 void
 fmprb_pow_ui(fmprb_t y, const fmprb_t b, ulong e, long prec)
 {
-    long i;
+    long i, wp;
 
     if (y == b)
     {
@@ -50,11 +50,13 @@ fmprb_pow_ui(fmprb_t y, const fmprb_t b, ulong e, long prec)
 
     fmprb_set(y, b);
 
+    wp = FMPR_PREC_ADD(prec, FLINT_BIT_COUNT(e));
+
     for (i = FLINT_BIT_COUNT(e) - 2; i >= 0; i--)
     {
-        fmprb_mul(y, y, y, prec);
+        fmprb_mul(y, y, y, wp);
         if (e & (1UL<<i))
-            fmprb_mul(y, y, b, prec);
+            fmprb_mul(y, y, b, wp);
     }
 }
 
diff --git a/fmprb/test/t-add.c b/fmprb/test/t-add.c
new file mode 100644
index 00000000..f84dd58f
--- /dev/null
+++ b/fmprb/test/t-add.c
@@ -0,0 +1,168 @@
+/*=============================================================================
+
+This file is part of fmprb.
+
+fmprb 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.
+
+fmprb 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 fmprb; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmprb.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("add....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 100000; iter++)
+    {
+        fmprb_t a, b, c;
+        fmpq_t x, y, z;
+
+        fmprb_init(a);
+        fmprb_init(b);
+        fmprb_init(c);
+
+        fmpq_init(x);
+        fmpq_init(y);
+        fmpq_init(z);
+
+        fmprb_randtest(a, state, 1 + n_randint(state, 200), 10);
+        fmprb_randtest(b, state, 1 + n_randint(state, 200), 10);
+        fmprb_randtest(c, state, 1 + n_randint(state, 200), 10);
+
+        fmprb_get_rand_fmpq(x, state, a, 1 + n_randint(state, 200));
+        fmprb_get_rand_fmpq(y, state, b, 1 + n_randint(state, 200));
+
+        fmprb_add(c, a, b, 2 + n_randint(state, 200));
+        fmpq_add(z, x, y);
+
+        if (!fmprb_contains_fmpq(c, z))
+        {
+            printf("FAIL: containment\n\n");
+            printf("a = "); fmprb_print(a); printf("\n\n");
+            printf("x = "); fmpq_print(x); printf("\n\n");
+            printf("b = "); fmprb_print(b); printf("\n\n");
+            printf("y = "); fmpq_print(y); printf("\n\n");
+            printf("c = "); fmprb_print(c); printf("\n\n");
+            printf("z = "); fmpq_print(z); printf("\n\n");
+            abort();
+        }
+
+        fmprb_clear(a);
+        fmprb_clear(b);
+        fmprb_clear(c);
+
+        fmpq_clear(x);
+        fmpq_clear(y);
+        fmpq_clear(z);
+    }
+
+    /* aliasing of c and a */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        fmprb_t a, b;
+        fmpq_t x, y, z;
+
+        fmprb_init(a);
+        fmprb_init(b);
+
+        fmpq_init(x);
+        fmpq_init(y);
+        fmpq_init(z);
+
+        fmprb_randtest(a, state, 1 + n_randint(state, 200), 10);
+        fmprb_randtest(b, state, 1 + n_randint(state, 200), 10);
+
+        fmprb_get_rand_fmpq(x, state, a, 1 + n_randint(state, 200));
+        fmprb_get_rand_fmpq(y, state, b, 1 + n_randint(state, 200));
+
+        fmprb_add(a, a, b, 2 + n_randint(state, 200));
+        fmpq_add(z, x, y);
+
+        if (!fmprb_contains_fmpq(a, z))
+        {
+            printf("FAIL: aliasing (c, a)\n\n");
+            printf("a = "); fmprb_print(a); printf("\n\n");
+            printf("x = "); fmpq_print(x); printf("\n\n");
+            printf("b = "); fmprb_print(b); printf("\n\n");
+            printf("y = "); fmpq_print(y); printf("\n\n");
+            printf("z = "); fmpq_print(z); printf("\n\n");
+            abort();
+        }
+
+        fmprb_clear(a);
+        fmprb_clear(b);
+
+        fmpq_clear(x);
+        fmpq_clear(y);
+        fmpq_clear(z);
+    }
+
+    /* aliasing of c and b */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        fmprb_t a, b;
+        fmpq_t x, y, z;
+
+        fmprb_init(a);
+        fmprb_init(b);
+
+        fmpq_init(x);
+        fmpq_init(y);
+        fmpq_init(z);
+
+        fmprb_randtest(a, state, 1 + n_randint(state, 200), 10);
+        fmprb_randtest(b, state, 1 + n_randint(state, 200), 10);
+
+        fmprb_get_rand_fmpq(x, state, a, 1 + n_randint(state, 200));
+        fmprb_get_rand_fmpq(y, state, b, 1 + n_randint(state, 200));
+
+        fmprb_add(b, a, b, 2 + n_randint(state, 200));
+        fmpq_add(z, x, y);
+
+        if (!fmprb_contains_fmpq(b, z))
+        {
+            printf("FAIL: aliasing (c, b)\n\n");
+            printf("a = "); fmprb_print(a); printf("\n\n");
+            printf("x = "); fmpq_print(x); printf("\n\n");
+            printf("b = "); fmprb_print(b); printf("\n\n");
+            printf("y = "); fmpq_print(y); printf("\n\n");
+            printf("z = "); fmpq_print(z); printf("\n\n");
+            abort();
+        }
+
+        fmprb_clear(a);
+        fmprb_clear(b);
+
+        fmpq_clear(x);
+        fmpq_clear(y);
+        fmpq_clear(z);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}
diff --git a/fmprb/test/t-addmul.c b/fmprb/test/t-addmul.c
new file mode 100644
index 00000000..7a6adccc
--- /dev/null
+++ b/fmprb/test/t-addmul.c
@@ -0,0 +1,171 @@
+/*=============================================================================
+
+This file is part of fmprb.
+
+fmprb 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.
+
+fmprb 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 fmprb; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmprb.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("addmul....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 100000; iter++)
+    {
+        fmprb_t a, b, c;
+        fmpq_t x, y, z;
+
+        fmprb_init(a);
+        fmprb_init(b);
+        fmprb_init(c);
+
+        fmpq_init(x);
+        fmpq_init(y);
+        fmpq_init(z);
+
+        fmprb_randtest(a, state, 1 + n_randint(state, 200), 10);
+        fmprb_randtest(b, state, 1 + n_randint(state, 200), 10);
+        fmprb_randtest(c, state, 1 + n_randint(state, 200), 10);
+
+        fmprb_get_rand_fmpq(x, state, a, 1 + n_randint(state, 200));
+        fmprb_get_rand_fmpq(y, state, b, 1 + n_randint(state, 200));
+        fmprb_get_rand_fmpq(z, state, c, 1 + n_randint(state, 200));
+
+        fmprb_addmul(c, a, b, 2 + n_randint(state, 200));
+        fmpq_addmul(z, x, y);
+
+        if (!fmprb_contains_fmpq(c, z))
+        {
+            printf("FAIL: containment\n\n");
+            printf("a = "); fmprb_print(a); printf("\n\n");
+            printf("x = "); fmpq_print(x); printf("\n\n");
+            printf("b = "); fmprb_print(b); printf("\n\n");
+            printf("y = "); fmpq_print(y); printf("\n\n");
+            printf("c = "); fmprb_print(c); printf("\n\n");
+            printf("z = "); fmpq_print(z); printf("\n\n");
+            abort();
+        }
+
+        fmprb_clear(a);
+        fmprb_clear(b);
+        fmprb_clear(c);
+
+        fmpq_clear(x);
+        fmpq_clear(y);
+        fmpq_clear(z);
+    }
+
+    /* aliasing of c and a */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        fmprb_t a, b;
+        fmpq_t x, y, z;
+
+        fmprb_init(a);
+        fmprb_init(b);
+
+        fmpq_init(x);
+        fmpq_init(y);
+        fmpq_init(z);
+
+        fmprb_randtest(a, state, 1 + n_randint(state, 200), 10);
+        fmprb_randtest(b, state, 1 + n_randint(state, 200), 10);
+
+        fmprb_get_rand_fmpq(x, state, a, 1 + n_randint(state, 200));
+        fmprb_get_rand_fmpq(y, state, b, 1 + n_randint(state, 200));
+        fmpq_set(z, x);
+
+        fmprb_addmul(a, a, b, 2 + n_randint(state, 200));
+        fmpq_addmul(z, x, y);
+
+        if (!fmprb_contains_fmpq(a, z))
+        {
+            printf("FAIL: aliasing (c, a)\n\n");
+            printf("a = "); fmprb_print(a); printf("\n\n");
+            printf("x = "); fmpq_print(x); printf("\n\n");
+            printf("b = "); fmprb_print(b); printf("\n\n");
+            printf("y = "); fmpq_print(y); printf("\n\n");
+            printf("z = "); fmpq_print(z); printf("\n\n");
+            abort();
+        }
+
+        fmprb_clear(a);
+        fmprb_clear(b);
+
+        fmpq_clear(x);
+        fmpq_clear(y);
+        fmpq_clear(z);
+    }
+
+    /* aliasing of c and b */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        fmprb_t a, b;
+        fmpq_t x, y, z;
+
+        fmprb_init(a);
+        fmprb_init(b);
+
+        fmpq_init(x);
+        fmpq_init(y);
+        fmpq_init(z);
+
+        fmprb_randtest(a, state, 1 + n_randint(state, 200), 10);
+        fmprb_randtest(b, state, 1 + n_randint(state, 200), 10);
+
+        fmprb_get_rand_fmpq(x, state, a, 1 + n_randint(state, 200));
+        fmprb_get_rand_fmpq(y, state, b, 1 + n_randint(state, 200));
+        fmpq_set(z, y);
+
+        fmprb_addmul(b, a, b, 2 + n_randint(state, 200));
+        fmpq_addmul(z, x, y);
+
+        if (!fmprb_contains_fmpq(b, z))
+        {
+            printf("FAIL: aliasing (c, b)\n\n");
+            printf("a = "); fmprb_print(a); printf("\n\n");
+            printf("x = "); fmpq_print(x); printf("\n\n");
+            printf("b = "); fmprb_print(b); printf("\n\n");
+            printf("y = "); fmpq_print(y); printf("\n\n");
+            printf("z = "); fmpq_print(z); printf("\n\n");
+            abort();
+        }
+
+        fmprb_clear(a);
+        fmprb_clear(b);
+
+        fmpq_clear(x);
+        fmpq_clear(y);
+        fmpq_clear(z);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}
diff --git a/fmprb/test/t-div.c b/fmprb/test/t-div.c
new file mode 100644
index 00000000..3ecdf48a
--- /dev/null
+++ b/fmprb/test/t-div.c
@@ -0,0 +1,174 @@
+/*=============================================================================
+
+This file is part of fmprb.
+
+fmprb 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.
+
+fmprb 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 fmprb; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmprb.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("div....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 100000; iter++)
+    {
+        fmprb_t a, b, c;
+        fmpq_t x, y, z;
+
+        fmprb_init(a);
+        fmprb_init(b);
+        fmprb_init(c);
+
+        fmpq_init(x);
+        fmpq_init(y);
+        fmpq_init(z);
+
+        do {
+            fmprb_randtest(a, state, 1 + n_randint(state, 200), 10);
+            fmprb_randtest(b, state, 1 + n_randint(state, 200), 10);
+            fmprb_randtest(c, state, 1 + n_randint(state, 200), 10);
+
+            fmprb_get_rand_fmpq(x, state, a, 1 + n_randint(state, 200));
+            fmprb_get_rand_fmpq(y, state, b, 1 + n_randint(state, 200));
+        } while (fmpq_is_zero(y));
+
+        fmprb_div(c, a, b, 2 + n_randint(state, 200));
+        fmpq_div(z, x, y);
+
+        if (!fmprb_contains_fmpq(c, z))
+        {
+            printf("FAIL: containment\n\n");
+            printf("a = "); fmprb_print(a); printf("\n\n");
+            printf("x = "); fmpq_print(x); printf("\n\n");
+            printf("b = "); fmprb_print(b); printf("\n\n");
+            printf("y = "); fmpq_print(y); printf("\n\n");
+            printf("c = "); fmprb_print(c); printf("\n\n");
+            printf("z = "); fmpq_print(z); printf("\n\n");
+            abort();
+        }
+
+        fmprb_clear(a);
+        fmprb_clear(b);
+        fmprb_clear(c);
+
+        fmpq_clear(x);
+        fmpq_clear(y);
+        fmpq_clear(z);
+    }
+
+    /* aliasing of c and a */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        fmprb_t a, b;
+        fmpq_t x, y, z;
+
+        fmprb_init(a);
+        fmprb_init(b);
+
+        fmpq_init(x);
+        fmpq_init(y);
+        fmpq_init(z);
+
+        do {
+            fmprb_randtest(a, state, 1 + n_randint(state, 200), 10);
+            fmprb_randtest(b, state, 1 + n_randint(state, 200), 10);
+
+            fmprb_get_rand_fmpq(x, state, a, 1 + n_randint(state, 200));
+            fmprb_get_rand_fmpq(y, state, b, 1 + n_randint(state, 200));
+        } while (fmpq_is_zero(y));
+
+        fmprb_div(a, a, b, 2 + n_randint(state, 200));
+        fmpq_div(z, x, y);
+
+        if (!fmprb_contains_fmpq(a, z))
+        {
+            printf("FAIL: aliasing (c, a)\n\n");
+            printf("a = "); fmprb_print(a); printf("\n\n");
+            printf("x = "); fmpq_print(x); printf("\n\n");
+            printf("b = "); fmprb_print(b); printf("\n\n");
+            printf("y = "); fmpq_print(y); printf("\n\n");
+            printf("z = "); fmpq_print(z); printf("\n\n");
+            abort();
+        }
+
+        fmprb_clear(a);
+        fmprb_clear(b);
+
+        fmpq_clear(x);
+        fmpq_clear(y);
+        fmpq_clear(z);
+    }
+
+    /* aliasing of c and b */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        fmprb_t a, b;
+        fmpq_t x, y, z;
+
+        fmprb_init(a);
+        fmprb_init(b);
+
+        fmpq_init(x);
+        fmpq_init(y);
+        fmpq_init(z);
+
+        do {
+            fmprb_randtest(a, state, 1 + n_randint(state, 200), 10);
+            fmprb_randtest(b, state, 1 + n_randint(state, 200), 10);
+
+            fmprb_get_rand_fmpq(x, state, a, 1 + n_randint(state, 200));
+            fmprb_get_rand_fmpq(y, state, b, 1 + n_randint(state, 200));
+        } while (fmpq_is_zero(y));
+
+        fmprb_div(b, a, b, 2 + n_randint(state, 200));
+        fmpq_div(z, x, y);
+
+        if (!fmprb_contains_fmpq(b, z))
+        {
+            printf("FAIL: aliasing (c, b)\n\n");
+            printf("a = "); fmprb_print(a); printf("\n\n");
+            printf("x = "); fmpq_print(x); printf("\n\n");
+            printf("b = "); fmprb_print(b); printf("\n\n");
+            printf("y = "); fmpq_print(y); printf("\n\n");
+            printf("z = "); fmpq_print(z); printf("\n\n");
+            abort();
+        }
+
+        fmprb_clear(a);
+        fmprb_clear(b);
+
+        fmpq_clear(x);
+        fmpq_clear(y);
+        fmpq_clear(z);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}
diff --git a/fmprb/test/t-mul.c b/fmprb/test/t-mul.c
new file mode 100644
index 00000000..8a25d9df
--- /dev/null
+++ b/fmprb/test/t-mul.c
@@ -0,0 +1,168 @@
+/*=============================================================================
+
+This file is part of fmprb.
+
+fmprb 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.
+
+fmprb 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 fmprb; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmprb.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("mul....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 100000; iter++)
+    {
+        fmprb_t a, b, c;
+        fmpq_t x, y, z;
+
+        fmprb_init(a);
+        fmprb_init(b);
+        fmprb_init(c);
+
+        fmpq_init(x);
+        fmpq_init(y);
+        fmpq_init(z);
+
+        fmprb_randtest(a, state, 1 + n_randint(state, 200), 10);
+        fmprb_randtest(b, state, 1 + n_randint(state, 200), 10);
+        fmprb_randtest(c, state, 1 + n_randint(state, 200), 10);
+
+        fmprb_get_rand_fmpq(x, state, a, 1 + n_randint(state, 200));
+        fmprb_get_rand_fmpq(y, state, b, 1 + n_randint(state, 200));
+
+        fmprb_mul(c, a, b, 2 + n_randint(state, 200));
+        fmpq_mul(z, x, y);
+
+        if (!fmprb_contains_fmpq(c, z))
+        {
+            printf("FAIL: containment\n\n");
+            printf("a = "); fmprb_print(a); printf("\n\n");
+            printf("x = "); fmpq_print(x); printf("\n\n");
+            printf("b = "); fmprb_print(b); printf("\n\n");
+            printf("y = "); fmpq_print(y); printf("\n\n");
+            printf("c = "); fmprb_print(c); printf("\n\n");
+            printf("z = "); fmpq_print(z); printf("\n\n");
+            abort();
+        }
+
+        fmprb_clear(a);
+        fmprb_clear(b);
+        fmprb_clear(c);
+
+        fmpq_clear(x);
+        fmpq_clear(y);
+        fmpq_clear(z);
+    }
+
+    /* aliasing of c and a */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        fmprb_t a, b;
+        fmpq_t x, y, z;
+
+        fmprb_init(a);
+        fmprb_init(b);
+
+        fmpq_init(x);
+        fmpq_init(y);
+        fmpq_init(z);
+
+        fmprb_randtest(a, state, 1 + n_randint(state, 200), 10);
+        fmprb_randtest(b, state, 1 + n_randint(state, 200), 10);
+
+        fmprb_get_rand_fmpq(x, state, a, 1 + n_randint(state, 200));
+        fmprb_get_rand_fmpq(y, state, b, 1 + n_randint(state, 200));
+
+        fmprb_mul(a, a, b, 2 + n_randint(state, 200));
+        fmpq_mul(z, x, y);
+
+        if (!fmprb_contains_fmpq(a, z))
+        {
+            printf("FAIL: aliasing (c, a)\n\n");
+            printf("a = "); fmprb_print(a); printf("\n\n");
+            printf("x = "); fmpq_print(x); printf("\n\n");
+            printf("b = "); fmprb_print(b); printf("\n\n");
+            printf("y = "); fmpq_print(y); printf("\n\n");
+            printf("z = "); fmpq_print(z); printf("\n\n");
+            abort();
+        }
+
+        fmprb_clear(a);
+        fmprb_clear(b);
+
+        fmpq_clear(x);
+        fmpq_clear(y);
+        fmpq_clear(z);
+    }
+
+    /* aliasing of c and b */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        fmprb_t a, b;
+        fmpq_t x, y, z;
+
+        fmprb_init(a);
+        fmprb_init(b);
+
+        fmpq_init(x);
+        fmpq_init(y);
+        fmpq_init(z);
+
+        fmprb_randtest(a, state, 1 + n_randint(state, 200), 10);
+        fmprb_randtest(b, state, 1 + n_randint(state, 200), 10);
+
+        fmprb_get_rand_fmpq(x, state, a, 1 + n_randint(state, 200));
+        fmprb_get_rand_fmpq(y, state, b, 1 + n_randint(state, 200));
+
+        fmprb_mul(b, a, b, 2 + n_randint(state, 200));
+        fmpq_mul(z, x, y);
+
+        if (!fmprb_contains_fmpq(b, z))
+        {
+            printf("FAIL: aliasing (c, b)\n\n");
+            printf("a = "); fmprb_print(a); printf("\n\n");
+            printf("x = "); fmpq_print(x); printf("\n\n");
+            printf("b = "); fmprb_print(b); printf("\n\n");
+            printf("y = "); fmpq_print(y); printf("\n\n");
+            printf("z = "); fmpq_print(z); printf("\n\n");
+            abort();
+        }
+
+        fmprb_clear(a);
+        fmprb_clear(b);
+
+        fmpq_clear(x);
+        fmpq_clear(y);
+        fmpq_clear(z);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}
diff --git a/fmprb/test/t-sub.c b/fmprb/test/t-sub.c
new file mode 100644
index 00000000..c12820be
--- /dev/null
+++ b/fmprb/test/t-sub.c
@@ -0,0 +1,168 @@
+/*=============================================================================
+
+This file is part of fmprb.
+
+fmprb 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.
+
+fmprb 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 fmprb; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmprb.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("sub....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 100000; iter++)
+    {
+        fmprb_t a, b, c;
+        fmpq_t x, y, z;
+
+        fmprb_init(a);
+        fmprb_init(b);
+        fmprb_init(c);
+
+        fmpq_init(x);
+        fmpq_init(y);
+        fmpq_init(z);
+
+        fmprb_randtest(a, state, 1 + n_randint(state, 200), 10);
+        fmprb_randtest(b, state, 1 + n_randint(state, 200), 10);
+        fmprb_randtest(c, state, 1 + n_randint(state, 200), 10);
+
+        fmprb_get_rand_fmpq(x, state, a, 1 + n_randint(state, 200));
+        fmprb_get_rand_fmpq(y, state, b, 1 + n_randint(state, 200));
+
+        fmprb_sub(c, a, b, 2 + n_randint(state, 200));
+        fmpq_sub(z, x, y);
+
+        if (!fmprb_contains_fmpq(c, z))
+        {
+            printf("FAIL: containment\n\n");
+            printf("a = "); fmprb_print(a); printf("\n\n");
+            printf("x = "); fmpq_print(x); printf("\n\n");
+            printf("b = "); fmprb_print(b); printf("\n\n");
+            printf("y = "); fmpq_print(y); printf("\n\n");
+            printf("c = "); fmprb_print(c); printf("\n\n");
+            printf("z = "); fmpq_print(z); printf("\n\n");
+            abort();
+        }
+
+        fmprb_clear(a);
+        fmprb_clear(b);
+        fmprb_clear(c);
+
+        fmpq_clear(x);
+        fmpq_clear(y);
+        fmpq_clear(z);
+    }
+
+    /* aliasing of c and a */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        fmprb_t a, b;
+        fmpq_t x, y, z;
+
+        fmprb_init(a);
+        fmprb_init(b);
+
+        fmpq_init(x);
+        fmpq_init(y);
+        fmpq_init(z);
+
+        fmprb_randtest(a, state, 1 + n_randint(state, 200), 10);
+        fmprb_randtest(b, state, 1 + n_randint(state, 200), 10);
+
+        fmprb_get_rand_fmpq(x, state, a, 1 + n_randint(state, 200));
+        fmprb_get_rand_fmpq(y, state, b, 1 + n_randint(state, 200));
+
+        fmprb_sub(a, a, b, 2 + n_randint(state, 200));
+        fmpq_sub(z, x, y);
+
+        if (!fmprb_contains_fmpq(a, z))
+        {
+            printf("FAIL: aliasing (c, a)\n\n");
+            printf("a = "); fmprb_print(a); printf("\n\n");
+            printf("x = "); fmpq_print(x); printf("\n\n");
+            printf("b = "); fmprb_print(b); printf("\n\n");
+            printf("y = "); fmpq_print(y); printf("\n\n");
+            printf("z = "); fmpq_print(z); printf("\n\n");
+            abort();
+        }
+
+        fmprb_clear(a);
+        fmprb_clear(b);
+
+        fmpq_clear(x);
+        fmpq_clear(y);
+        fmpq_clear(z);
+    }
+
+    /* aliasing of c and b */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        fmprb_t a, b;
+        fmpq_t x, y, z;
+
+        fmprb_init(a);
+        fmprb_init(b);
+
+        fmpq_init(x);
+        fmpq_init(y);
+        fmpq_init(z);
+
+        fmprb_randtest(a, state, 1 + n_randint(state, 200), 10);
+        fmprb_randtest(b, state, 1 + n_randint(state, 200), 10);
+
+        fmprb_get_rand_fmpq(x, state, a, 1 + n_randint(state, 200));
+        fmprb_get_rand_fmpq(y, state, b, 1 + n_randint(state, 200));
+
+        fmprb_sub(b, a, b, 2 + n_randint(state, 200));
+        fmpq_sub(z, x, y);
+
+        if (!fmprb_contains_fmpq(b, z))
+        {
+            printf("FAIL: aliasing (c, b)\n\n");
+            printf("a = "); fmprb_print(a); printf("\n\n");
+            printf("x = "); fmpq_print(x); printf("\n\n");
+            printf("b = "); fmprb_print(b); printf("\n\n");
+            printf("y = "); fmpq_print(y); printf("\n\n");
+            printf("z = "); fmpq_print(z); printf("\n\n");
+            abort();
+        }
+
+        fmprb_clear(a);
+        fmprb_clear(b);
+
+        fmpq_clear(x);
+        fmpq_clear(y);
+        fmpq_clear(z);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}
diff --git a/fmprb/test/t-submul.c b/fmprb/test/t-submul.c
new file mode 100644
index 00000000..aecc4a18
--- /dev/null
+++ b/fmprb/test/t-submul.c
@@ -0,0 +1,171 @@
+/*=============================================================================
+
+This file is part of fmprb.
+
+fmprb 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.
+
+fmprb 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 fmprb; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmprb.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("submul....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 100000; iter++)
+    {
+        fmprb_t a, b, c;
+        fmpq_t x, y, z;
+
+        fmprb_init(a);
+        fmprb_init(b);
+        fmprb_init(c);
+
+        fmpq_init(x);
+        fmpq_init(y);
+        fmpq_init(z);
+
+        fmprb_randtest(a, state, 1 + n_randint(state, 200), 10);
+        fmprb_randtest(b, state, 1 + n_randint(state, 200), 10);
+        fmprb_randtest(c, state, 1 + n_randint(state, 200), 10);
+
+        fmprb_get_rand_fmpq(x, state, a, 1 + n_randint(state, 200));
+        fmprb_get_rand_fmpq(y, state, b, 1 + n_randint(state, 200));
+        fmprb_get_rand_fmpq(z, state, c, 1 + n_randint(state, 200));
+
+        fmprb_submul(c, a, b, 2 + n_randint(state, 200));
+        fmpq_submul(z, x, y);
+
+        if (!fmprb_contains_fmpq(c, z))
+        {
+            printf("FAIL: containment\n\n");
+            printf("a = "); fmprb_print(a); printf("\n\n");
+            printf("x = "); fmpq_print(x); printf("\n\n");
+            printf("b = "); fmprb_print(b); printf("\n\n");
+            printf("y = "); fmpq_print(y); printf("\n\n");
+            printf("c = "); fmprb_print(c); printf("\n\n");
+            printf("z = "); fmpq_print(z); printf("\n\n");
+            abort();
+        }
+
+        fmprb_clear(a);
+        fmprb_clear(b);
+        fmprb_clear(c);
+
+        fmpq_clear(x);
+        fmpq_clear(y);
+        fmpq_clear(z);
+    }
+
+    /* aliasing of c and a */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        fmprb_t a, b;
+        fmpq_t x, y, z;
+
+        fmprb_init(a);
+        fmprb_init(b);
+
+        fmpq_init(x);
+        fmpq_init(y);
+        fmpq_init(z);
+
+        fmprb_randtest(a, state, 1 + n_randint(state, 200), 10);
+        fmprb_randtest(b, state, 1 + n_randint(state, 200), 10);
+
+        fmprb_get_rand_fmpq(x, state, a, 1 + n_randint(state, 200));
+        fmprb_get_rand_fmpq(y, state, b, 1 + n_randint(state, 200));
+        fmpq_set(z, x);
+
+        fmprb_submul(a, a, b, 2 + n_randint(state, 200));
+        fmpq_submul(z, x, y);
+
+        if (!fmprb_contains_fmpq(a, z))
+        {
+            printf("FAIL: aliasing (c, a)\n\n");
+            printf("a = "); fmprb_print(a); printf("\n\n");
+            printf("x = "); fmpq_print(x); printf("\n\n");
+            printf("b = "); fmprb_print(b); printf("\n\n");
+            printf("y = "); fmpq_print(y); printf("\n\n");
+            printf("z = "); fmpq_print(z); printf("\n\n");
+            abort();
+        }
+
+        fmprb_clear(a);
+        fmprb_clear(b);
+
+        fmpq_clear(x);
+        fmpq_clear(y);
+        fmpq_clear(z);
+    }
+
+    /* aliasing of c and b */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        fmprb_t a, b;
+        fmpq_t x, y, z;
+
+        fmprb_init(a);
+        fmprb_init(b);
+
+        fmpq_init(x);
+        fmpq_init(y);
+        fmpq_init(z);
+
+        fmprb_randtest(a, state, 1 + n_randint(state, 200), 10);
+        fmprb_randtest(b, state, 1 + n_randint(state, 200), 10);
+
+        fmprb_get_rand_fmpq(x, state, a, 1 + n_randint(state, 200));
+        fmprb_get_rand_fmpq(y, state, b, 1 + n_randint(state, 200));
+        fmpq_set(z, y);
+
+        fmprb_submul(b, a, b, 2 + n_randint(state, 200));
+        fmpq_submul(z, x, y);
+
+        if (!fmprb_contains_fmpq(b, z))
+        {
+            printf("FAIL: aliasing (c, b)\n\n");
+            printf("a = "); fmprb_print(a); printf("\n\n");
+            printf("x = "); fmpq_print(x); printf("\n\n");
+            printf("b = "); fmprb_print(b); printf("\n\n");
+            printf("y = "); fmpq_print(y); printf("\n\n");
+            printf("z = "); fmpq_print(z); printf("\n\n");
+            abort();
+        }
+
+        fmprb_clear(a);
+        fmprb_clear(b);
+
+        fmpq_clear(x);
+        fmpq_clear(y);
+        fmpq_clear(z);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}
diff --git a/fmprb/zeta_ui.c b/fmprb/zeta_ui.c
index 700cb512..cdfbe289 100644
--- a/fmprb/zeta_ui.c
+++ b/fmprb/zeta_ui.c
@@ -84,9 +84,7 @@ void
 fmprb_zeta_ui_vec_odd(fmprb_struct * x, ulong start, long num, long prec)
 {
     long i, num_borwein;
-    ulong end, cutoff;
-
-    end = start + num * 2;
+    ulong cutoff;
 
     cutoff = 40 + 0.3 * prec;