multiplication code

arf.h

@@ -21,6 +21,9 @@
 
     Copyright (C) 2014 Fredrik Johansson
 
+    2x2 mul code taken from MPFR 2.3.0
+    (Copyright (C) 1991-2007 Free Software Foundation, Inc.)
+
 ******************************************************************************/
 
 #ifndef ARF_H
@@ -38,6 +41,8 @@ extern "C" {
 #define LIMB_ONES (-(mp_limb_t) 1)
 #define LIMB_TOP (((mp_limb_t) 1) << (FLINT_BITS - 1))
 
+#define MASK_LIMB(n, c) ((n) & (LIMB_ONES << (c)))
+
 #define arf_rnd_t fmpr_rnd_t
 #define ARF_RND_DOWN FMPR_RND_DOWN
 #define ARF_RND_UP FMPR_RND_UP
@@ -45,6 +50,15 @@ extern "C" {
 #define ARF_RND_CEIL FMPR_RND_CEIL
 #define ARF_RND_NEAR FMPR_RND_NEAR
 
+static __inline__ int
+arf_rounds_down(arf_rnd_t rnd, int sgnbit)
+{
+    if (rnd == ARF_RND_DOWN) return 1;
+    if (rnd == ARF_RND_UP) return 0;
+    if (rnd == ARF_RND_FLOOR) return !sgnbit;
+    return sgnbit;
+}
+
 static __inline__ int
 arf_rounds_up(arf_rnd_t rnd, int sgnbit)
 {
@@ -678,7 +692,7 @@ arf_set_fmpz_2exp(arf_t x, const fmpz_t man, const fmpz_t exp)
 
 void arf_debug(const arf_t x);
 
-#define arf_print arf_debug
+void arf_print(const arf_t x);
 
 void arf_randtest(arf_t x, flint_rand_t state, long bits, long mag_bits);
 
@@ -686,6 +700,105 @@ void arf_randtest_not_zero(arf_t x, flint_rand_t state, long bits, long mag_bits
 
 void arf_randtest_special(arf_t x, flint_rand_t state, long bits, long mag_bits);
 
+#define ADD2_FAST_MAX (COEFF_MAX / 4)
+#define ADD2_FAST_MIN (-ADD2_FAST_MAX)
+
+static __inline__ void
+_fmpz_add2_fast(fmpz_t z, const fmpz_t x, const fmpz_t y, long c)
+{
+    fmpz ze, xe, ye;
+
+    ze = *z;
+    xe = *x;
+    ye = *y;
+
+    if (!COEFF_IS_MPZ(ze) && (xe > ADD2_FAST_MIN && xe < ADD2_FAST_MAX) &&
+                             (ye > ADD2_FAST_MIN && ye < ADD2_FAST_MAX))
+    {
+        *z = xe + ye + c;
+    }
+    else
+    {
+        fmpz_add(z, x, y);
+        fmpz_add_si(z, z, c);
+    }
+}
+
+#define MUL_MPFR_MIN_LIMBS 25
+#define MUL_MPFR_MAX_LIMBS 10000
+
+#define nn_mul_2x1(r2, r1, r0, a1, a0, b0)                  \
+    do {                                                    \
+        mp_limb_t t1;                                       \
+        umul_ppmm(r1, r0, a0, b0);                          \
+        umul_ppmm(r2, t1, a1, b0);                          \
+        add_ssaaaa(r2, r1, r2, r1, 0, t1);                  \
+    } while (0)
+
+#define nn_mul_2x2(r3, r2, r1, r0, a1, a0, b1, b0)          \
+    do {                                                    \
+        mp_limb_t t1, t2, t3;                               \
+        umul_ppmm(r1, r0, a0, b0);                          \
+        umul_ppmm(r2, t1, a1, b0);                          \
+        add_ssaaaa(r2, r1, r2, r1, 0, t1);                  \
+        umul_ppmm(t1, t2, a0, b1);                          \
+        umul_ppmm(r3, t3, a1, b1);                          \
+        add_ssaaaa(r3, t1, r3, t1, 0, t3);                  \
+        add_ssaaaa(r2, r1, r2, r1, t1, t2);                 \
+        r3 += r2 < t1;                                      \
+    } while (0)
+
+#define ARF_MUL_STACK_ALLOC 40
+#define ARF_MUL_TLS_ALLOC 1000
+
+extern TLS_PREFIX mp_ptr __arf_mul_tmp;
+extern TLS_PREFIX long __arf_mul_alloc;
+
+extern void _arf_mul_tmp_cleanup(void);
+
+#define ARF_MUL_TMP_DECL \
+    mp_limb_t tmp_stack[ARF_MUL_STACK_ALLOC]; \
+
+#define ARF_MUL_TMP_ALLOC(tmp, alloc) \
+    if (alloc <= ARF_MUL_STACK_ALLOC) \
+    { \
+        tmp = tmp_stack; \
+    } \
+    else if (alloc <= ARF_MUL_TLS_ALLOC) \
+    { \
+        if (__arf_mul_alloc < alloc) \
+        { \
+            if (__arf_mul_alloc == 0) \
+            { \
+                flint_register_cleanup_function(_arf_mul_tmp_cleanup); \
+            } \
+            __arf_mul_tmp = flint_realloc(__arf_mul_tmp, sizeof(mp_limb_t) * alloc); \
+            __arf_mul_alloc = alloc; \
+        } \
+        tmp = __arf_mul_tmp; \
+    } \
+    else \
+    { \
+        tmp = flint_malloc(sizeof(mp_limb_t) * alloc); \
+    }
+
+#define ARF_MUL_TMP_FREE(tmp, alloc) \
+    if (alloc > ARF_MUL_TLS_ALLOC) \
+        flint_free(tmp);
+
+void arf_mul_special(arf_t z, const arf_t x, const arf_t y);
+
+int arf_mul_via_mpfr(arf_t z, const arf_t x, const arf_t y, long prec, arf_rnd_t rnd);
+
+int arf_mul_rnd_any(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec, arf_rnd_t rnd);
+
+int arf_mul_rnd_down(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec);
+
+#define arf_mul(z, x, y, prec, rnd)              \
+    ((rnd == FMPR_RND_DOWN)                      \
+        ? arf_mul_rnd_down(z, x, y, prec)        \
+        : arf_mul_rnd_any(z, x, y, prec, rnd))
+
 #ifdef __cplusplus
 }
 #endif

arf/mul_rnd_any.c

@@ -0,0 +1,100 @@
+/*=============================================================================
+
+    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) 2014 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arf.h"
+
+TLS_PREFIX mp_ptr __arf_mul_tmp = NULL;
+TLS_PREFIX long __arf_mul_alloc = 0;
+
+void _arf_mul_tmp_cleanup(void)
+{
+    flint_free(__arf_mul_tmp);
+    __arf_mul_tmp = NULL;
+    __arf_mul_alloc = 0;
+}
+
+int
+arf_mul_rnd_any(arf_ptr z, arf_srcptr x, arf_srcptr y,
+        long prec, arf_rnd_t rnd)
+{
+    mp_size_t xn, yn;
+    int sgnbit, inexact;
+
+    xn = ARF_XSIZE(x);
+    yn = ARF_XSIZE(y);
+    sgnbit = (xn ^ yn) & 1;
+    xn >>= 1;
+    yn >>= 1;
+
+    if (yn > xn)
+    {
+        arf_srcptr __t; mp_size_t __u;
+        __t = x; x = y; y = __t;
+        __u = xn; xn = yn; yn = __u;
+    }
+
+    /* Either operand is a special value. */
+    if (yn == 0)
+    {
+        arf_mul_special(z, x, y);
+        return 0;
+    }
+    else
+    {
+        mp_size_t zn, alloc;
+        mp_srcptr xptr, yptr;
+        mp_ptr tmp;
+        ARF_MUL_TMP_DECL
+
+        ARF_GET_MPN_READONLY(xptr, xn, x);
+        ARF_GET_MPN_READONLY(yptr, yn, y);
+
+        alloc = zn = xn + yn;
+        ARF_MUL_TMP_ALLOC(tmp, alloc)
+
+        if (yn == 1)
+        {
+            tmp[zn - 1] = mpn_mul_1(tmp, xptr, xn, yptr[0]);
+        }
+        else if (xn == yn)
+        {
+            if (xptr == yptr)
+                mpn_sqr(tmp, xptr, xn);
+            else
+                mpn_mul_n(tmp, xptr, yptr, yn);
+        }
+        else 
+        {
+            mpn_mul(tmp, xptr, xn, yptr, yn);
+        }
+
+        _fmpz_add2_fast(ARF_EXPREF(z), ARF_EXPREF(x), ARF_EXPREF(y), 0);
+        inexact = arf_set_round_mpn(z, tmp, zn, sgnbit, ARF_EXPREF(z), prec, rnd);
+        ARF_MUL_TMP_FREE(tmp, alloc)
+
+        return inexact;
+    }
+}
+

arf/mul_rnd_down.c

@@ -0,0 +1,233 @@
+/*=============================================================================
+
+    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) 2014 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arf.h"
+
+int
+arf_mul_rnd_down(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec)
+{
+    mp_size_t xn, yn, zn;
+    mp_limb_t hi, lo;
+    int sgnbit, ret, fix;
+    mp_ptr zptr;
+
+    xn = ARF_XSIZE(x);
+    yn = ARF_XSIZE(y);
+    sgnbit = (xn ^ yn) & 1;
+    xn >>= 1;
+    yn >>= 1;
+
+    if (yn > xn)
+    {
+        arf_srcptr __t; mp_size_t __u;
+        __t = x; x = y; y = __t;
+        __u = xn; xn = yn; yn = __u;
+    }
+
+    /* Either operand is a special value. */
+    if (yn == 0)
+    {
+        arf_mul_special(z, x, y);
+        return 0;
+    }
+
+    /* xn == yn == 1 */
+    if (xn == 1)
+    {
+        lo = ARF_NOPTR_D(x)[0];
+        hi = ARF_NOPTR_D(y)[0];
+
+        umul_ppmm(hi, lo, hi, lo);
+
+        /* Shift so that the top bit is set (branch free version). */
+        fix = !(hi >> (FLINT_BITS - 1));
+        hi = (hi << fix) | ((lo >> (FLINT_BITS - 1)) & fix);
+        lo = (lo << fix);
+
+        ARF_DEMOTE(z);
+
+        if (lo == 0)
+        {
+            zn = 1;
+
+            if (prec >= FLINT_BITS)
+            {
+                lo = hi;
+                ret = 0;
+            }
+            else
+            {
+                lo = MASK_LIMB(hi, FLINT_BITS - prec);
+                ret = (lo != hi);
+            }
+        }
+        else
+        {
+            zn = 2;
+
+            if (prec <= FLINT_BITS)  /* Must be inexact. */
+            {
+                lo = MASK_LIMB(hi, FLINT_BITS - prec);
+                zn = ret = 1;
+            }
+            else if (prec >= 2 * FLINT_BITS) /* Must be exact. */
+            {
+                ret = 0;
+            }
+            else  /* prec < FLINT_BITS < 2 * FLINT_BITS */
+            {
+                ret = MASK_LIMB(lo, 2 * FLINT_BITS - prec) != lo;
+                lo = MASK_LIMB(lo, 2 * FLINT_BITS - prec);
+
+                if (lo == 0)
+                {
+                    zn = 1;
+                    lo = hi;
+                }
+            }
+        }
+
+        _fmpz_add2_fast(ARF_EXPREF(z), ARF_EXPREF(x), ARF_EXPREF(y), -fix);
+        ARF_XSIZE(z) = ARF_MAKE_XSIZE(zn, sgnbit);
+
+        zptr = ARF_NOPTR_D(z);
+        zptr[0] = lo;
+        zptr[1] = hi;
+
+        return ret;
+    }
+    else if (xn == 2)
+    {
+        mp_limb_t zz[4];
+        mp_limb_t x1, x0, y1, y0;
+
+        x0 = ARF_NOPTR_D(x)[0];
+        x1 = ARF_NOPTR_D(x)[1];
+
+        if (yn == 2)
+        {
+            y0 = ARF_NOPTR_D(y)[0];
+            y1 = ARF_NOPTR_D(y)[1];
+
+            nn_mul_2x2(zz[3], zz[2], zz[1], zz[0], x1, x0, y1, y0);
+
+            /* Likely case, must be inexact */
+            if (prec <= 2 * FLINT_BITS)
+            {
+                ARF_DEMOTE(z);
+
+                fix = !(zz[3] >> (FLINT_BITS - 1));
+                zz[3] = (zz[3] << fix) | ((zz[2] >> (FLINT_BITS - 1)) & fix);
+                zz[2] = (zz[2] << fix) | ((zz[1] >> (FLINT_BITS - 1)) & fix);
+
+                _fmpz_add2_fast(ARF_EXPREF(z), ARF_EXPREF(x), ARF_EXPREF(y), -fix);
+
+                /* Rounding */
+                if (prec != 2 * FLINT_BITS)
+                {
+                    if (prec > FLINT_BITS)
+                    {
+                        zz[2] &= (LIMB_ONES << (2 * FLINT_BITS - prec));
+                    }
+                    else if (prec == FLINT_BITS)
+                    {
+                        zz[2] = 0;
+                    }
+                    else
+                    {
+                        zz[3] &= (LIMB_ONES << (FLINT_BITS - prec));
+                        zz[2] = 0;
+                    }
+                }
+
+                if (zz[2] == 0)
+                {
+                    ARF_XSIZE(z) = ARF_MAKE_XSIZE(1, sgnbit);
+                    ARF_NOPTR_D(z)[0] = zz[3];
+                }
+                else
+                {
+                    ARF_XSIZE(z) = ARF_MAKE_XSIZE(2, sgnbit);
+                    ARF_NOPTR_D(z)[0] = zz[2];
+                    ARF_NOPTR_D(z)[1] = zz[3];
+                }
+
+                return 1;
+            }
+        }
+        else
+        {
+            y0 = ARF_NOPTR_D(y)[0];
+            nn_mul_2x1(zz[2], zz[1], zz[0], x1, x0, y0);
+        }
+
+        zn = xn + yn;
+        _fmpz_add2_fast(ARF_EXPREF(z), ARF_EXPREF(x), ARF_EXPREF(y), 0);
+        ret = arf_set_round_mpn(z, zz, zn, sgnbit, ARF_EXPREF(z), prec, ARF_RND_DOWN);
+        return ret;
+    }
+    else if (yn > MUL_MPFR_MIN_LIMBS && prec != ARF_PREC_EXACT
+                && xn + yn > 1.25 * prec / FLINT_BITS
+                && xn < MUL_MPFR_MAX_LIMBS)  /* FIXME: proper cutoffs */
+    {
+        return arf_mul_via_mpfr(z, x, y, prec, ARF_RND_DOWN);
+    }
+    else
+    {
+        mp_size_t zn, alloc;
+        mp_srcptr xptr, yptr;
+        mp_ptr tmp;
+        ARF_MUL_TMP_DECL
+
+        ARF_GET_MPN_READONLY(xptr, xn, x);
+        ARF_GET_MPN_READONLY(yptr, yn, y);
+
+        alloc = zn = xn + yn;
+        ARF_MUL_TMP_ALLOC(tmp, alloc)
+
+        if (xn == yn)
+        {
+            if (xptr == yptr)
+                mpn_sqr(tmp, xptr, xn);
+            else
+                mpn_mul_n(tmp, xptr, yptr, yn);
+        }
+        else if (yn == 1)
+        {
+            tmp[zn - 1] = mpn_mul_1(tmp, xptr, xn, yptr[0]);
+        }
+        else
+        {
+            mpn_mul(tmp, xptr, xn, yptr, yn);
+        }
+
+        _fmpz_add2_fast(ARF_EXPREF(z), ARF_EXPREF(x), ARF_EXPREF(y), 0);
+        ret = arf_set_round_mpn(z, tmp, zn, sgnbit, ARF_EXPREF(z), prec, ARF_RND_DOWN);
+        ARF_MUL_TMP_FREE(tmp, alloc)
+
+        return ret;
+    }
+}
+

arf/mul_special.c

@@ -0,0 +1,52 @@
+/*=============================================================================
+
+    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) 2014 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arf.h"
+
+void
+arf_mul_special(arf_t z, const arf_t x, const arf_t y)
+{
+    if (arf_is_zero(x))
+    {
+        if (arf_is_finite(y))
+            arf_zero(z);
+        else
+            arf_nan(z);
+    }
+    else if (arf_is_zero(y))
+    {
+        if (arf_is_finite(x))
+            arf_zero(z);
+        else
+            arf_nan(z);
+    }
+    else if (arf_is_nan(x) || arf_is_nan(y))
+        arf_nan(z);
+    else if (arf_sgn(x) == arf_sgn(y))
+        arf_pos_inf(z);
+    else
+        arf_neg_inf(z);
+}
+

arf/mul_via_mpfr.c

@@ -0,0 +1,94 @@
+/*=============================================================================
+
+    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) 2014 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arf.h"
+
+int
+arf_mul_via_mpfr(arf_ptr z, arf_srcptr x, arf_srcptr y,
+        long prec, arf_rnd_t rnd)
+{
+    mp_size_t xn, yn, zn, val;
+    mp_srcptr xptr, yptr;
+    mp_ptr tmp, zptr;
+    mpfr_t xf, yf, zf;
+    int ret;
+    ARF_MUL_TMP_DECL
+
+    if (arf_is_special(x) || arf_is_special(y))
+    {
+        arf_mul_special(z, x, y);
+        return 0;
+    }
+
+    ARF_GET_MPN_READONLY(xptr, xn, x);
+    ARF_GET_MPN_READONLY(yptr, yn, y);
+
+    prec = FLINT_MIN((xn + yn) * FLINT_BITS, prec);
+    zn = (prec + FLINT_BITS - 1) / FLINT_BITS;
+
+    ARF_MUL_TMP_ALLOC(tmp, zn)
+
+    zf->_mpfr_d = tmp;
+    zf->_mpfr_prec = prec;
+    zf->_mpfr_sign = 1;
+    zf->_mpfr_exp = 0;
+
+    xf->_mpfr_d = (mp_ptr) xptr;
+    xf->_mpfr_prec = xn * FLINT_BITS;
+    xf->_mpfr_sign = ARF_SGNBIT(x) ? -1 : 1;
+    xf->_mpfr_exp = 0;
+
+    if (x == y)
+    {
+        ret = mpfr_sqr(zf, xf, arf_rnd_to_mpfr(rnd));
+    }
+    else
+    {
+        yf->_mpfr_d = (mp_ptr) yptr;
+        yf->_mpfr_prec = yn * FLINT_BITS;
+        yf->_mpfr_sign = ARF_SGNBIT(y) ? -1 : 1;
+        yf->_mpfr_exp = 0;
+
+        ret = mpfr_mul(zf, xf, yf, arf_rnd_to_mpfr(rnd));
+    }
+
+    ret = (ret != 0);
+
+    fmpz_add2_fmpz_si_inline(ARF_EXPREF(z),
+        ARF_EXPREF(x), ARF_EXPREF(y), zf->_mpfr_exp);
+
+    val = 0;
+    while (tmp[val] == 0)
+        val++;
+
+    ARF_GET_MPN_WRITE(zptr, zn - val, z);
+    flint_mpn_copyi(zptr, tmp + val, zn - val);
+    ARF_XSIZE(z) |= (zf->_mpfr_sign < 0);
+
+    ARF_MUL_TMP_FREE(tmp, zn)
+
+    return ret;
+}
+

arf/print.c

@@ -0,0 +1,57 @@
+/*=============================================================================
+
+    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) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arf.h"
+
+void
+arf_print(const arf_t x)
+{
+    if (arf_is_normal(x))
+    {
+        fmpz_t man, exp;
+
+        fmpz_init(man);
+        fmpz_init(exp);
+
+        arf_get_fmpz_2exp(man, exp, x);
+
+        printf("(");
+        fmpz_print(man);
+        printf(" * 2^");
+        fmpz_print(exp);
+        printf(")");
+
+        fmpz_clear(man);
+        fmpz_clear(exp);
+    }
+    else
+    {
+        if (arf_is_zero(x)) printf("(0)");
+        else if (arf_is_pos_inf(x)) printf("(+inf)");
+        else if (arf_is_neg_inf(x)) printf("(-inf)");
+        else printf("(nan)");
+    }
+}
+

arf/test/t-mul.c

@@ -0,0 +1,178 @@
+/*=============================================================================
+
+    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) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arf.h"
+
+int
+arf_mul_naive(arf_t z, const arf_t x, const arf_t y, long prec, arf_rnd_t rnd)
+{
+    fmpr_t a, b;
+    long r;
+
+    fmpr_init(a);
+    fmpr_init(b);
+
+    arf_get_fmpr(a, x);
+    arf_get_fmpr(b, y);
+
+    r = fmpr_mul_naive(a, a, b, prec, rnd);
+
+    arf_set_fmpr(z, a);
+
+    fmpr_clear(a);
+    fmpr_clear(b);
+
+    return (r == FMPR_RESULT_EXACT) ? 0 : 1;
+}
+
+int main()
+{
+    long iter, iter2;
+    flint_rand_t state;
+
+    printf("mul....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 10000; iter++)
+    {
+        arf_t x, y, z, v;
+        long prec, r1, r2;
+        arf_rnd_t rnd;
+
+        arf_init(x);
+        arf_init(y);
+        arf_init(z);
+        arf_init(v);
+
+        for (iter2 = 0; iter2 < 100; iter2++)
+        {
+            arf_randtest_special(x, state, 2000, 100);
+            arf_randtest_special(y, state, 2000, 100);
+            prec = 2 + n_randint(state, 2000);
+
+            if (n_randint(state, 50) == 0)
+                prec = ARF_PREC_EXACT;
+
+            switch (n_randint(state, 4))
+            {
+                case 0:  rnd = ARF_RND_DOWN; break;
+                case 1:  rnd = ARF_RND_UP; break;
+                case 2:  rnd = ARF_RND_FLOOR; break;
+                default: rnd = ARF_RND_CEIL; break;
+            }
+
+            switch (n_randint(state, 5))
+            {
+            case 0:
+                r1 = arf_mul(z, x, y, prec, rnd);
+                r2 = arf_mul_naive(v, x, y, prec, rnd);
+                if (!arf_equal(z, v) || r1 != r2)
+                {
+                    printf("FAIL!\n");
+                    printf("prec = %ld, rnd = %d\n\n", prec, rnd);
+                    printf("x = "); arf_print(x); printf("\n\n");
+                    printf("y = "); arf_print(y); printf("\n\n");
+                    printf("z = "); arf_print(z); printf("\n\n");
+                    printf("v = "); arf_print(v); printf("\n\n");
+                    printf("r1 = %ld, r2 = %ld\n", r1, r2);
+                    abort();
+                }
+                break;
+
+            case 1:
+                r1 = arf_mul(z, x, x, prec, rnd);
+                r2 = arf_mul_naive(v, x, x, prec, rnd);
+                if (!arf_equal(z, v) || r1 != r2)
+                {
+                    printf("FAIL (aliasing 1)!\n");
+                    printf("prec = %ld, rnd = %d\n\n", prec, rnd);
+                    printf("x = "); arf_print(x); printf("\n\n");
+                    printf("z = "); arf_print(z); printf("\n\n");
+                    printf("v = "); arf_print(v); printf("\n\n");
+                    printf("r1 = %ld, r2 = %ld\n", r1, r2);
+                    abort();
+                }
+                break;
+
+            case 2:
+                r2 = arf_mul_naive(v, x, x, prec, rnd);
+                r1 = arf_mul(x, x, x, prec, rnd);
+                if (!arf_equal(v, x) || r1 != r2)
+                {
+                    printf("FAIL (aliasing 2)!\n");
+                    printf("prec = %ld, rnd = %d\n\n", prec, rnd);
+                    printf("x = "); arf_print(x); printf("\n\n");
+                    printf("z = "); arf_print(z); printf("\n\n");
+                    printf("v = "); arf_print(v); printf("\n\n");
+                    printf("r1 = %ld, r2 = %ld\n", r1, r2);
+                    abort();
+                }
+                break;
+
+            case 3:
+                r2 = arf_mul_naive(v, x, y, prec, rnd);
+                r1 = arf_mul(x, x, y, prec, rnd);
+                if (!arf_equal(x, v) || r1 != r2)
+                {
+                    printf("FAIL (aliasing 3)!\n");
+                    printf("prec = %ld, rnd = %d\n\n", prec, rnd);
+                    printf("x = "); arf_print(x); printf("\n\n");
+                    printf("y = "); arf_print(y); printf("\n\n");
+                    printf("v = "); arf_print(v); printf("\n\n");
+                    printf("r1 = %ld, r2 = %ld\n", r1, r2);
+                    abort();
+                }
+                break;
+
+            default:
+                r2 = arf_mul_naive(v, x, y, prec, rnd);
+                r1 = arf_mul(x, y, x, prec, rnd);
+                if (!arf_equal(x, v) || r1 != r2)
+                {
+                    printf("FAIL (aliasing 4)!\n");
+                    printf("prec = %ld, rnd = %d\n\n", prec, rnd);
+                    printf("x = "); arf_print(x); printf("\n\n");
+                    printf("y = "); arf_print(y); printf("\n\n");
+                    printf("v = "); arf_print(v); printf("\n\n");
+                    printf("r1 = %ld, r2 = %ld\n", r1, r2);
+                    abort();
+                }
+                break;
+            }
+        }
+
+        arf_clear(x);
+        arf_clear(y);
+        arf_clear(z);
+        arf_clear(v);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}

arf/test/t-mul_via_mpfr.c

@@ -0,0 +1,178 @@
+/*=============================================================================
+
+    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) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arf.h"
+
+int
+arf_mul_naive(arf_t z, const arf_t x, const arf_t y, long prec, arf_rnd_t rnd)
+{
+    fmpr_t a, b;
+    long r;
+
+    fmpr_init(a);
+    fmpr_init(b);
+
+    arf_get_fmpr(a, x);
+    arf_get_fmpr(b, y);
+
+    r = fmpr_mul_naive(a, a, b, prec, rnd);
+
+    arf_set_fmpr(z, a);
+
+    fmpr_clear(a);
+    fmpr_clear(b);
+
+    return (r == FMPR_RESULT_EXACT) ? 0 : 1;
+}
+
+int main()
+{
+    long iter, iter2;
+    flint_rand_t state;
+
+    printf("mul_via_mpfr....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 10000; iter++)
+    {
+        arf_t x, y, z, v;
+        long prec, r1, r2;
+        arf_rnd_t rnd;
+
+        arf_init(x);
+        arf_init(y);
+        arf_init(z);
+        arf_init(v);
+
+        for (iter2 = 0; iter2 < 100; iter2++)
+        {
+            arf_randtest_special(x, state, 2000, 100);
+            arf_randtest_special(y, state, 2000, 100);
+            prec = 2 + n_randint(state, 2000);
+
+            if (n_randint(state, 50) == 0)
+                prec = ARF_PREC_EXACT;
+
+            switch (n_randint(state, 4))
+            {
+                case 0:  rnd = ARF_RND_DOWN; break;
+                case 1:  rnd = ARF_RND_UP; break;
+                case 2:  rnd = ARF_RND_FLOOR; break;
+                default: rnd = ARF_RND_CEIL; break;
+            }
+
+            switch (n_randint(state, 5))
+            {
+            case 0:
+                r1 = arf_mul_via_mpfr(z, x, y, prec, rnd);
+                r2 = arf_mul_naive(v, x, y, prec, rnd);
+                if (!arf_equal(z, v) || r1 != r2)
+                {
+                    printf("FAIL!\n");
+                    printf("prec = %ld, rnd = %d\n\n", prec, rnd);
+                    printf("x = "); arf_print(x); printf("\n\n");
+                    printf("y = "); arf_print(y); printf("\n\n");
+                    printf("z = "); arf_print(z); printf("\n\n");
+                    printf("v = "); arf_print(v); printf("\n\n");
+                    printf("r1 = %ld, r2 = %ld\n", r1, r2);
+                    abort();
+                }
+                break;
+
+            case 1:
+                r1 = arf_mul_via_mpfr(z, x, x, prec, rnd);
+                r2 = arf_mul_naive(v, x, x, prec, rnd);
+                if (!arf_equal(z, v) || r1 != r2)
+                {
+                    printf("FAIL (aliasing 1)!\n");
+                    printf("prec = %ld, rnd = %d\n\n", prec, rnd);
+                    printf("x = "); arf_print(x); printf("\n\n");
+                    printf("z = "); arf_print(z); printf("\n\n");
+                    printf("v = "); arf_print(v); printf("\n\n");
+                    printf("r1 = %ld, r2 = %ld\n", r1, r2);
+                    abort();
+                }
+                break;
+
+            case 2:
+                r2 = arf_mul_naive(v, x, x, prec, rnd);
+                r1 = arf_mul_via_mpfr(x, x, x, prec, rnd);
+                if (!arf_equal(v, x) || r1 != r2)
+                {
+                    printf("FAIL (aliasing 2)!\n");
+                    printf("prec = %ld, rnd = %d\n\n", prec, rnd);
+                    printf("x = "); arf_print(x); printf("\n\n");
+                    printf("z = "); arf_print(z); printf("\n\n");
+                    printf("v = "); arf_print(v); printf("\n\n");
+                    printf("r1 = %ld, r2 = %ld\n", r1, r2);
+                    abort();
+                }
+                break;
+
+            case 3:
+                r2 = arf_mul_naive(v, x, y, prec, rnd);
+                r1 = arf_mul_via_mpfr(x, x, y, prec, rnd);
+                if (!arf_equal(x, v) || r1 != r2)
+                {
+                    printf("FAIL (aliasing 3)!\n");
+                    printf("prec = %ld, rnd = %d\n\n", prec, rnd);
+                    printf("x = "); arf_print(x); printf("\n\n");
+                    printf("y = "); arf_print(y); printf("\n\n");
+                    printf("v = "); arf_print(v); printf("\n\n");
+                    printf("r1 = %ld, r2 = %ld\n", r1, r2);
+                    abort();
+                }
+                break;
+
+            default:
+                r2 = arf_mul_naive(v, x, y, prec, rnd);
+                r1 = arf_mul_via_mpfr(x, y, x, prec, rnd);
+                if (!arf_equal(x, v) || r1 != r2)
+                {
+                    printf("FAIL (aliasing 4)!\n");
+                    printf("prec = %ld, rnd = %d\n\n", prec, rnd);
+                    printf("x = "); arf_print(x); printf("\n\n");
+                    printf("y = "); arf_print(y); printf("\n\n");
+                    printf("v = "); arf_print(v); printf("\n\n");
+                    printf("r1 = %ld, r2 = %ld\n", r1, r2);
+                    abort();
+                }
+                break;
+            }
+        }
+
+        arf_clear(x);
+        arf_clear(y);
+        arf_clear(z);
+        arf_clear(v);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}