Add a function to compute b/sqrt(a)

[fredrik-johansson]

An arb_div_sqrt function to compute b/sqrt(a) could be useful.

Dan Schultz proposed the following algorithm:

    x = 1/sqrt(a) + O(ep)
    y = b x + O(ep)
    b x + 1/2 y (1 - a x^2) = b/sqrt(a) + O(ep^2)

So something like this:

void
_arf_div_sqrt_approx(arf_t res, const arf_t b, const arf_t a, slong prec)
{
    slong wp = prec + GUARD_BITS;
    slong hp = prec / 2 + GUARD_BITS;

    arf_t t, bx, x, y;

    arf_init(t);
    arf_init(x);
    arf_init(bx);
    arf_init(y);

    /* x = 1/sqrt(a) + O(eps) */
    arf_set_round(x, a, hp, ARF_RND_DOWN);
    _arf_rsqrt_newton(x, x, hp);

    /* bx = bx + O(eps^2) */
    arf_mul(bx, b, x, wp, ARF_RND_DOWN);
    /* y = bx + O(eps) */
    arf_set_round(y, bx, hp, ARF_RND_DOWN);

    /* t = 1/2 y (1 - ax^2) */
    arf_mul(t, x, x, wp, ARF_RND_DOWN);
    arf_mul(t, t, a, wp, ARF_RND_DOWN);
    arf_sub_ui(t, t, 1, hp, ARF_RND_DOWN);
    arf_neg(t, t);
    arf_mul(t, t, y, hp, ARF_RND_DOWN);
    arf_mul_2exp_si(t, t, -1);

    arf_add(res, bx, t, prec, ARF_RND_DOWN);

    arf_clear(t);
    arf_clear(x);
    arf_clear(bx);
    arf_clear(y);
}

On my machine this is 5-8% faster than arb_rsqrt + arb_mul at high enough precision.