mul_via_mpfr.c 1.81 KB
Newer Older
1
/*
2 3
    Copyright (C) 2014 Fredrik Johansson

4 5 6 7 8 9 10
    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/>.
*/
11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79

#include "arf.h"

int
arf_mul_via_mpfr(arf_ptr z, arf_srcptr x, arf_srcptr y,
        slong 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_fast(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;
}