more work on powers

doc/source/fmpcb.rst

@@ -168,7 +168,7 @@ Complex parts
     the special value `\operatorname{arg}(0) = 0`, and
     `\operatorname{arg}(a+0i) = \pi` for `a < 0`. Equivalently, if
     `z = a+bi`, the argument is given by `\operatorname{atan2}(b,a)`
-    (see *fmprb_atan2*).
+    (see :funct:`fmprb_atan2`).
 
 .. function:: void fmpcb_abs(fmprb_t r, const fmpcb_t z, long prec)
 
@@ -308,11 +308,19 @@ Elementary functions
 
 .. function:: void fmpcb_pow_ui(fmpcb_t y, const fmpcb_t b, ulong e, long prec)
 
-    Sets *y* to *b* raised to the power *e*, computed using binary exponentiation.
+    Sets `y = b^e` using binary exponentiation (with an initial division
+    if `e < 0`). Provided that *b* and *e*
+    are small enough and the exponent is positive, the exact power can be
+    computed by setting the precision to *FMPR_PREC_EXACT*.
 
-.. function:: void fmpcb_pow(fmpcb_t r, const fmpcb_t x, const fmpcb_t y, long prec)
+    Note that these functions can get slow if the exponent is
+    extremely large (in such cases :func:`fmpcb_pow` may be superior).
 
-    Sets *r* to *x* raised to the power *y*, computed as `x^y = \exp(y \log x)`.
+.. function:: void fmpcb_pow(fmpcb_t z, const fmpcb_t x, const fmpcb_t y, long prec)
+
+    Sets `z = x^y`, computed using binary exponentiation if `y` if
+    a small exact integer, as `z = (x^{1/2})^{2y}` if `y` is a small exact
+    half-integer, and generally as `z = \exp(y \log x)`.
 
 .. function:: void fmpcb_sqrt(fmpcb_t r, const fmpcb_t z, long prec)
 

fmpcb.h

@@ -461,28 +461,7 @@ fmpcb_submul_fmprb(fmpcb_t z, const fmpcb_t x, const fmprb_t y, long prec)
     fmprb_submul(fmpcb_imagref(z), fmpcb_imagref(x), y, prec);
 }
 
-static __inline__ void
+void fmpcb_inv(fmpcb_t z, const fmpcb_t x, long prec);
-fmpcb_inv(fmpcb_t z, const fmpcb_t x, long prec)
-{
-    fmprb_t t;
-    fmprb_init(t);
-
-#define a fmpcb_realref(x)
-#define b fmpcb_imagref(x)
-
-    fmprb_mul(t, a, a, prec);
-    fmprb_addmul(t, b, b, prec);
-
-    fmprb_div(fmpcb_realref(z), a, t, prec);
-    fmprb_div(fmpcb_imagref(z), b, t, prec);
-
-    fmprb_neg(fmpcb_imagref(z), fmpcb_imagref(z));
-
-#undef a
-#undef b
-
-    fmprb_clear(t);
-}
 
 static __inline__ void
 fmpcb_div(fmpcb_t z, const fmpcb_t x, const fmpcb_t y, long prec)
@@ -531,6 +510,7 @@ void fmpcb_sin_cos(fmpcb_t s, fmpcb_t c, const fmpcb_t z, long prec);
 void fmpcb_sin_pi(fmpcb_t r, const fmpcb_t z, long prec);
 void fmpcb_sin_cos_pi(fmpcb_t s, fmpcb_t c, const fmpcb_t z, long prec);
 
+void fmpcb_pow_fmprb(fmpcb_t z, const fmpcb_t x, const fmprb_t y, long prec);
 void fmpcb_pow(fmpcb_t r, const fmpcb_t x, const fmpcb_t y, long prec);
 
 void fmpcb_sqrt(fmpcb_t y, const fmpcb_t x, long prec);

fmpcb/inv.c

@@ -0,0 +1,68 @@
+/*=============================================================================
+
+    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) 2013 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpcb.h"
+
+void
+fmpcb_inv(fmpcb_t z, const fmpcb_t x, long prec)
+{
+#define a fmpcb_realref(x)
+#define b fmpcb_imagref(x)
+#define c fmpcb_realref(z)
+#define d fmpcb_imagref(z)
+
+    if (fmprb_is_zero(b))
+    {
+        fmprb_ui_div(c, 1, a, prec);
+        fmprb_zero(d);
+    }
+    else if (fmprb_is_zero(a))
+    {
+        fmprb_ui_div(d, 1, b, prec);
+        fmprb_neg(d, d);
+        fmprb_zero(c);
+    }
+    else
+    {
+        fmprb_t t;
+        fmprb_init(t);
+
+        fmprb_mul(t, a, a, prec);
+        fmprb_addmul(t, b, b, prec);
+
+        fmprb_div(c, a, t, prec);
+        fmprb_div(d, b, t, prec);
+
+        fmprb_neg(d, d);
+
+        fmprb_clear(t);
+    }
+
+#undef a
+#undef b
+#undef c
+#undef d
+}
+

fmpcb/pow.c

@@ -26,13 +26,25 @@
 #include "fmpcb.h"
 
 void
-fmpcb_pow_fmpz(fmpcb_t y, const fmpcb_t b, const fmpz_t e, long prec)
+fmpcb_pow_fmpz_binexp(fmpcb_t y, const fmpcb_t b, const fmpz_t e, long prec)
 {
     long i, wp, bits;
 
-    if (fmpz_is_zero(e))
+    if (-2L <= *e && *e <= 2L)
     {
-        fmpcb_one(y);
+        if (*e == 0L)
+            fmpcb_one(y);
+        else if (*e == 1L)
+            fmpcb_set_round(y, b, prec);
+        else if (*e == -1L)
+            fmpcb_inv(y, b, prec);
+        else if (*e == 2L)
+            fmpcb_mul(y, b, b, prec);
+        else
+        {
+            fmpcb_inv(y, b, prec);
+            fmpcb_mul(y, y, y, prec);
+        }
         return;
     }
 
@@ -41,7 +53,7 @@ fmpcb_pow_fmpz(fmpcb_t y, const fmpcb_t b, const fmpz_t e, long prec)
         fmpz_t f;
         fmpz_init(f);
         fmpz_neg(f, e);
-        fmpcb_pow_fmpz(y, b, f, prec + 2);
+        fmpcb_pow_fmpz_binexp(y, b, f, prec + 2);
         fmpcb_inv(y, y, prec);
         fmpz_clear(f);
         return;
@@ -52,7 +64,7 @@ fmpcb_pow_fmpz(fmpcb_t y, const fmpcb_t b, const fmpz_t e, long prec)
         fmpcb_t t;
         fmpcb_init(t);
         fmpcb_set(t, b);
-        fmpcb_pow_fmpz(y, t, e, prec);
+        fmpcb_pow_fmpz_binexp(y, t, e, prec);
         fmpcb_clear(t);
         return;
     }
@@ -70,6 +82,12 @@ fmpcb_pow_fmpz(fmpcb_t y, const fmpcb_t b, const fmpz_t e, long prec)
     }
 }
 
+void
+fmpcb_pow_fmpz(fmpcb_t y, const fmpcb_t b, const fmpz_t e, long prec)
+{
+    fmpcb_pow_fmpz_binexp(y, b, e, prec);
+}
+
 void
 fmpcb_pow_ui(fmpcb_t y, const fmpcb_t b, ulong e, long prec)
 {
@@ -90,13 +108,86 @@ fmpcb_pow_si(fmpcb_t y, const fmpcb_t b, long e, long prec)
 }
 
 void
-fmpcb_pow(fmpcb_t r, const fmpcb_t x, const fmpcb_t y, long prec)
+_fmpcb_pow_exp(fmpcb_t z, const fmpcb_t x, const fmpcb_t y, long prec)
 {
     fmpcb_t t;
     fmpcb_init(t);
     fmpcb_log(t, x, prec);
     fmpcb_mul(t, t, y, prec);
-    fmpcb_exp(r, t, prec);
+    fmpcb_exp(z, t, prec);
     fmpcb_clear(t);
 }
 
+void
+_fmpcb_pow_fmprb_exp(fmpcb_t z, const fmpcb_t x, const fmprb_t y, long prec)
+{
+    fmpcb_t t;
+    fmpcb_init(t);
+    fmpcb_log(t, x, prec);
+    fmpcb_mul_fmprb(t, t, y, prec);
+    fmpcb_exp(z, t, prec);
+    fmpcb_clear(t);
+}
+
+void
+fmpcb_pow_fmprb(fmpcb_t z, const fmpcb_t x, const fmprb_t y, long prec)
+{
+    if (fmprb_is_exact(y))
+    {
+        if (fmpr_is_zero(fmprb_midref(y)))
+        {
+            fmpcb_one(z);
+            return;
+        }
+
+        if (!fmpr_is_special(fmprb_midref(y)))
+        {
+            const fmpz * exp_exp = fmpr_expref(fmprb_midref(y));
+
+            /* smallish integer powers and square roots */
+            if (!COEFF_IS_MPZ(*exp_exp) && (*exp_exp >= -1L))
+            {
+                long exp_bits;
+
+                exp_bits = *exp_exp + fmpz_bits(fmpr_manref(fmprb_midref(y)));
+
+                if (exp_bits < 64)
+                {
+                    fmpz_t e;
+                    fmpz_init(e);
+
+                    if (*exp_exp == -1L)
+                    {
+                        fmpcb_sqrt(z, x, prec + exp_bits);
+                        fmpz_set(e, fmpr_manref(fmprb_midref(y)));
+                        fmpcb_pow_fmpz(z, z, e, prec);
+                    }
+                    else
+                    {
+                        fmpz_mul_2exp(e, fmpr_manref(fmprb_midref(y)), *exp_exp);
+                        fmpcb_pow_fmpz(z, x, e, prec);
+                    }
+
+                    fmpz_clear(e);
+                    return;
+                }
+            }
+        }
+    }
+
+    _fmpcb_pow_fmprb_exp(z, x, y, prec);
+}
+
+void
+fmpcb_pow(fmpcb_t z, const fmpcb_t x, const fmpcb_t y, long prec)
+{
+    if (fmprb_is_zero(fmpcb_imagref(y)))
+    {
+        fmpcb_pow_fmprb(z, x, fmpcb_realref(y), prec);
+    }
+    else
+    {
+        _fmpcb_pow_exp(z, x, y, prec);
+    }
+}
+

fmpcb/test/t-pow.c

@@ -0,0 +1,103 @@
+/*=============================================================================
+
+    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, 2013 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpcb.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("pow....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    /* check large arguments */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        fmpcb_t a, b, c, d, e, f;
+        long prec1, prec2;
+
+        prec1 = 2 + n_randint(state, 1000);
+        prec2 = prec1 + 30;
+
+        fmpcb_init(a);
+        fmpcb_init(b);
+        fmpcb_init(c);
+        fmpcb_init(d);
+        fmpcb_init(e);
+        fmpcb_init(f);
+
+        fmpcb_randtest(a, state, 1 + n_randint(state, 1000), 200);
+        fmpcb_randtest(b, state, 1 + n_randint(state, 1000), 200);
+
+        fmpcb_pow(c, a, b, prec1);
+        fmpcb_pow(d, a, b, prec2);
+
+        if (!fmpcb_overlaps(c, d))
+        {
+            printf("FAIL: overlap\n\n");
+            printf("a = "); fmpcb_print(a); printf("\n\n");
+            printf("b = "); fmpcb_print(b); printf("\n\n");
+            printf("c = "); fmpcb_print(c); printf("\n\n");
+            printf("d = "); fmpcb_print(d); printf("\n\n");
+            abort();
+        }
+
+        fmpcb_randtest(c, state, 1 + n_randint(state, 1000), 200);
+
+        /* check a^(b+c) = a^b*a^c */
+        fmpcb_add(e, b, c, prec1);
+        fmpcb_pow(d, a, e, prec1);
+
+        fmpcb_pow(e, a, b, prec1);
+        fmpcb_pow(f, a, c, prec1);
+        fmpcb_mul(e, e, f, prec1);
+
+        if (!fmpcb_overlaps(d, e))
+        {
+            printf("FAIL: functional equation\n\n");
+            printf("a = "); fmpcb_print(a); printf("\n\n");
+            printf("b = "); fmpcb_print(b); printf("\n\n");
+            printf("c = "); fmpcb_print(c); printf("\n\n");
+            printf("d = "); fmpcb_print(d); printf("\n\n");
+            printf("e = "); fmpcb_print(e); printf("\n\n");
+            abort();
+        }
+
+        fmpcb_clear(a);
+        fmpcb_clear(b);
+        fmpcb_clear(c);
+        fmpcb_clear(d);
+        fmpcb_clear(e);
+        fmpcb_clear(f);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}

fmpcb/test/t-pow_fmpz.c

@@ -0,0 +1,104 @@
+/*=============================================================================
+
+    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, 2013 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpcb.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("pow_fmpz....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 5000; iter++)
+    {
+        fmpcb_t a, b, c, d;
+        fmpz_t e1, e2, e3;
+        long prec1, prec2;
+
+        prec1 = 2 + n_randint(state, 1000);
+        prec2 = prec1 + 30;
+
+        fmpcb_init(a);
+        fmpcb_init(b);
+        fmpcb_init(c);
+        fmpcb_init(d);
+        fmpz_init(e1);
+        fmpz_init(e2);
+        fmpz_init(e3);
+
+        fmpcb_randtest(a, state, 1 + n_randint(state, 1000), 200);
+        fmpcb_randtest(b, state, 1 + n_randint(state, 1000), 200);
+        fmpz_randtest(e1, state, 200);
+        fmpz_randtest(e2, state, 200);
+
+        fmpcb_pow_fmpz(b, a, e1, prec1);
+        fmpcb_pow_fmpz(c, a, e1, prec2);
+
+        if (!fmpcb_overlaps(b, c))
+        {
+            printf("FAIL: overlap\n\n");
+            printf("a = "); fmpcb_print(a); printf("\n\n");
+            printf("b = "); fmpcb_print(b); printf("\n\n");
+            printf("c = "); fmpcb_print(c); printf("\n\n");
+            printf("e1 = "); fmpz_print(e1); printf("\n\n");
+            abort();
+        }
+
+        /* check a^(e1+e2) = a^e1*a^e2 */
+        fmpcb_pow_fmpz(c, a, e2, prec1);
+        fmpcb_mul(d, b, c, prec1);
+        fmpz_add(e3, e1, e2);
+        fmpcb_pow_fmpz(c, a, e3, prec1);
+
+        if (!fmpcb_overlaps(c, d))
+        {
+            printf("FAIL: functional equation\n\n");
+            printf("a = "); fmpcb_print(a); printf("\n\n");
+            printf("b = "); fmpcb_print(b); printf("\n\n");
+            printf("c = "); fmpcb_print(c); printf("\n\n");
+            printf("d = "); fmpcb_print(d); printf("\n\n");
+            printf("e1 = "); fmpz_print(e1); printf("\n\n");
+            printf("e2 = "); fmpz_print(e2); printf("\n\n");
+            abort();
+        }
+
+        fmpcb_clear(a);
+        fmpcb_clear(b);
+        fmpcb_clear(c);
+        fmpcb_clear(d);
+        fmpz_clear(e1);
+        fmpz_clear(e2);
+        fmpz_clear(e3);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}

fmprb/pow_fmpz_binexp.c

@@ -53,7 +53,7 @@ fmprb_pow_fmpz_binexp(fmprb_t y, const fmprb_t b, const fmpz_t e, long prec)
         fmpz_t f;
         fmpz_init(f);
         fmpz_neg(f, e);
-        fmprb_pow_fmpz(y, b, f, prec + 2);
+        fmprb_pow_fmpz_binexp(y, b, f, prec + 2);
         fmprb_ui_div(y, 1UL, y, prec);
         fmpz_clear(f);
         return;
@@ -64,7 +64,7 @@ fmprb_pow_fmpz_binexp(fmprb_t y, const fmprb_t b, const fmpz_t e, long prec)
         fmprb_t t;
         fmprb_init(t);
         fmprb_set(t, b);
-        fmprb_pow_fmpz(y, t, e, prec);
+        fmprb_pow_fmpz_binexp(y, t, e, prec);
         fmprb_clear(t);
         return;
     }

fmprb/test/t-pow.c

@@ -51,8 +51,8 @@ int main()
         fmprb_init(e);
         fmprb_init(f);
 
-        fmprb_randtest_precise(a, state, 1 + n_randint(state, 1000), 200);
+        fmprb_randtest(a, state, 1 + n_randint(state, 1000), 200);
-        fmprb_randtest_precise(b, state, 1 + n_randint(state, 1000), 200);
+        fmprb_randtest(b, state, 1 + n_randint(state, 1000), 200);
 
         fmprb_pow(c, a, b, prec1);
         fmprb_pow(d, a, b, prec2);
@@ -67,7 +67,7 @@ int main()
             abort();
         }
 
-        fmprb_randtest_precise(c, state, 1 + n_randint(state, 1000), 200);
+        fmprb_randtest(c, state, 1 + n_randint(state, 1000), 200);
 
         /* check a^(b+c) = a^b*a^c */
         fmprb_add(e, b, c, prec1);