[WIP] Split arb_pow_fmpz_binexp into a ui and mpz version and implement proper arb_sqr

arb.h

@@ -413,6 +413,8 @@ void arb_mul_si(arb_t z, const arb_t x, slong y, slong prec);
 void arb_mul_ui(arb_t z, const arb_t x, ulong y, slong prec);
 void arb_mul_fmpz(arb_t z, const arb_t x, const fmpz_t y, slong prec);
 
+void arb_sqr(arb_t res, const arb_t x, slong prec);
+
 void arb_addmul(arb_t z, const arb_t x, const arb_t y, slong prec);
 void arb_addmul_arf(arb_t z, const arb_t x, const arf_t y, slong prec);
 void arb_addmul_si(arb_t z, const arb_t x, slong y, slong prec);
@@ -484,9 +486,27 @@ void arb_rsqrt(arb_t z, const arb_t x, slong prec);
 void arb_rsqrt_ui(arb_t z, ulong x, slong prec);
 void arb_sqrt1pm1(arb_t r, const arb_t z, slong prec);
 
-void arb_pow_fmpz_binexp(arb_t y, const arb_t b, const fmpz_t e, slong prec);
-void arb_pow_fmpz(arb_t y, const arb_t b, const fmpz_t e, slong prec);
-void arb_pow_ui(arb_t y, const arb_t b, ulong e, slong prec);
+void _arb_pow_mpz_binexp(arb_t y, const arb_t b, mpz_srcptr e, slong prec);
+void arb_pow_ui_binexp(arb_t y, const arb_t b, ulong e, slong prec);
+void arb_pow_si(arb_t y, const arb_t b, slong e, slong prec);
+ARB_INLINE void
+arb_pow_fmpz_binexp(arb_t y, const arb_t b, const fmpz_t e, slong prec)
+{
+    if (!COEFF_IS_MPZ(*e))
+        arb_pow_si(y, b, *e, prec);
+    else
+        _arb_pow_mpz_binexp(y, b, COEFF_TO_PTR(*e), prec);
+}
+ARB_INLINE void
+arb_pow_fmpz(arb_t y, const arb_t b, const fmpz_t e, slong prec)
+{
+    arb_pow_fmpz_binexp(y, b, e, prec);
+}
+ARB_INLINE void
+arb_pow_ui(arb_t y, const arb_t b, ulong e, slong prec)
+{
+    arb_pow_ui_binexp(y, b, e, prec);
+}
 void arb_ui_pow_ui(arb_t y, ulong b, ulong e, slong prec);
 void arb_si_pow_ui(arb_t y, slong b, ulong e, slong prec);
 void arb_pow_fmpq(arb_t y, const arb_t x, const fmpq_t a, slong prec);
@@ -620,12 +640,6 @@ void arb_primorial_ui(arb_t res, ulong n, slong prec);
 
 void arb_lambertw(arb_t res, const arb_t x, int flags, slong prec);
 
-ARB_INLINE void
-arb_sqr(arb_t res, const arb_t val, slong prec)
-{
-    arb_mul(res, val, val, prec);
-}
-
 #define ARB_DEF_CACHED_CONSTANT(name, comp_func) \
     TLS_PREFIX slong name ## _cached_prec = 0; \
     TLS_PREFIX arb_t name ## _cached_value; \

arb/pow_binexp.c

@@ -0,0 +1,100 @@
+/*
+    Copyright (C) 2012, 2013 Fredrik Johansson
+    Copyright (C) 2021 Albin Ahlbäck
+
+    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/>.
+*/
+
+#include "arb.h"
+
+void
+_arb_pow_mpz_binexp(arb_t y, const arb_t b, mpz_srcptr e, slong prec)
+{
+    slong i, wp, bits;
+
+    if (e->_mp_size < 0)
+    {
+        __mpz_struct f = { 1, 0, NULL };
+        f._mp_size = -(e->_mp_size);
+        f._mp_d = e->_mp_d;
+
+        if (arb_is_exact(b))
+        {
+            _arb_pow_mpz_binexp(y, b, &f, prec + 2);
+            arb_inv(y, y, prec);
+        }
+        else
+        {
+            arb_inv(y, b, prec + mpz_sizeinbase(e, 2) + 2);
+            _arb_pow_mpz_binexp(y, y, &f, prec);
+        }
+
+        return;
+    }
+
+    if (y == b)
+    {
+        arb_t t;
+        arb_init(t);
+        arb_set(t, b);
+        _arb_pow_mpz_binexp(y, t, e, prec);
+        arb_clear(t);
+        return;
+    }
+
+    arb_set(y, b);
+
+    bits = mpz_sizeinbase(e, 2);
+    wp = ARF_PREC_ADD(prec, bits);
+
+    for (i = bits - 2; i >= 0; i--)
+    {
+        arb_sqr(y, y, wp);
+        if (mpz_tstbit(e, i))
+            arb_mul(y, y, b, wp);
+    }
+}
+
+void
+arb_pow_ui_binexp(arb_t y, const arb_t b, ulong e, slong prec)
+{
+    slong i, wp, bits;
+
+    if (e <= 2)
+    {
+        if (e == 0)
+            arb_one(y);
+        else if (e == 1)
+            arb_set_round(y, b, prec);
+        else
+            arb_sqr(y, b, prec);
+        return;
+    }
+
+    if (y == b)
+    {
+        arb_t t;
+        arb_init(t);
+        arb_set(t, b);
+        arb_pow_ui_binexp(y, t, e, prec);
+        arb_clear(t);
+        return;
+    }
+
+    arb_set(y, b);
+
+    bits = FLINT_BIT_COUNT(e);
+    wp = ARF_PREC_ADD(prec, bits);
+
+    for (i = bits - 2; i >= 0; i--)
+    {
+        arb_sqr(y, y, wp);
+        if (((WORD(1) << i) & e) != 0)
+            arb_mul(y, y, b, wp);
+    }
+}

arb/pow_si.c

@@ -0,0 +1,37 @@
+/*
+    Copyright (C) 2021 Albin Ahlbäck
+
+    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/>.
+*/
+
+#include "arb.h"
+
+void
+arb_pow_si(arb_t y, const arb_t b, slong e, slong prec)
+{
+    if (e >= 0)
+        arb_pow_ui_binexp(y, b, e, prec);
+    else if (e == -1)
+        arb_inv(y, b, prec);
+    else if (e == -2)
+    {
+        arb_inv(y, b, prec);
+        arb_sqr(y, y, prec);
+    }
+    else if (arb_is_exact(b))
+    {
+        arb_pow_ui_binexp(y, b, -e, prec + 2);
+        arb_inv(y, y, prec);
+    }
+    else
+    {
+        e = -e;
+        arb_inv(y, b, prec + FLINT_BIT_COUNT(e) + 2);
+        arb_pow_ui_binexp(y, y, e, prec);
+    }
+}

arb/sqr.c

@@ -0,0 +1,67 @@
+/*
+    Copyright (C) 2021 Albin Ahlbäck
+
+    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/>.
+*/
+
+#include "arb.h"
+
+void
+arb_sqr(arb_t res, const arb_t x, slong prec)
+{
+    mag_t resr, xm;
+    int inexact;
+
+    if (arb_is_exact(x))
+    {
+        inexact = arf_sqr_rnd_down(arb_midref(res), arb_midref(x), prec);
+
+        if (inexact)
+            arf_mag_set_ulp(arb_radref(res), arb_midref(res), prec);
+        else
+            mag_zero(arb_radref(res));
+    }
+    else if (ARB_IS_LAGOM(x) && ARB_IS_LAGOM(res))
+    {
+        mag_fast_init_set_arf(xm, arb_midref(x));
+
+        mag_init(resr);
+        mag_fast_mul(resr, xm, arb_radref(x));
+        if (resr->exp < COEFF_MAX && resr->exp >= COEFF_MIN)
+            (resr->exp)++;
+        else
+            fmpz_add_ui(&(resr->exp), &(resr->exp), 1);
+        mag_fast_addmul(resr, arb_radref(x), arb_radref(x));
+
+        inexact = arf_sqr_rnd_down(arb_midref(res), arb_midref(x), prec);
+
+        if (inexact)
+            arf_mag_fast_add_ulp(resr, resr, arb_midref(res), prec);
+
+        *arb_radref(res) = *resr;
+    }
+    else
+    {
+        mag_init_set_arf(xm, arb_midref(x));
+
+        mag_init(resr);
+        mag_mul(resr, xm, arb_radref(x));
+        fmpz_add_ui(&(resr->exp), &(resr->exp), 1);
+        mag_addmul(resr, arb_radref(x), arb_radref(x));
+
+        inexact = arf_sqr_rnd_down(arb_midref(res), arb_midref(x), prec);
+
+        if (inexact)
+            arf_mag_add_ulp(arb_radref(res), resr, arb_midref(res), prec);
+        else
+            mag_swap(arb_radref(res), resr);
+
+        mag_clear(xm);
+        mag_clear(resr);
+    }
+}