add fmpcb_rsqrt (complex reciprocal square root)

doc/source/fmpcb.rst

@@ -404,6 +404,14 @@ Elementary functions
     `\sqrt{a+bi} = u/2 + ib/u, u = \sqrt{2(|a+bi|+a)}`,
     requiring two real square root extractions.
 
+.. function:: void fmpcb_rsqrt(fmpcb_t r, const fmpcb_t z, long prec)
+
+    Sets *r* to the reciprocal square root of *z*.
+    If either the real or imaginary part is exactly zero, only
+    a single real reciprocal square root is needed. Generally, we use the
+    formula `1/\sqrt{a+bi} = ((a+r) - bi)/v, r = |a+bi|, v = \sqrt{r |a+bi+r|^2}`,
+    requiring one real square root and one real reciprocal square root.
+
 .. function:: void fmpcb_invroot_newton(fmpcb_t r, const fmpcb_t a, ulong m, const fmpcb_t r0, long startprec, long prec)
 
     Given one inverse *m*-th root *r0* (with a valid error bound) of the complex

fmpcb.h

@@ -569,6 +569,7 @@ 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);
+void fmpcb_rsqrt(fmpcb_t y, const fmpcb_t x, long prec);
 
 void fmpcb_invroot_newton(fmpcb_t r, const fmpcb_t a, ulong m, const fmpcb_t r0, long startprec, long prec);
 void fmpcb_root_exp(fmpcb_t r, const fmpcb_t a, long m, long index, long prec);

fmpcb/rsqrt.c

@@ -0,0 +1,116 @@
+/*=============================================================================
+
+    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"
+#include "profiler.h"
+
+void
+fmpcb_rsqrt(fmpcb_t y, const fmpcb_t x, long prec)
+{
+    fmprb_t r, t, u, v;
+    long wp;
+
+#define a fmpcb_realref(x)
+#define b fmpcb_imagref(x)
+#define c fmpcb_realref(y)
+#define d fmpcb_imagref(y)
+
+    if (fmprb_is_zero(b))
+    {
+        if (fmprb_is_nonnegative(a))
+        {
+            fmprb_rsqrt(c, a, prec);
+            fmprb_zero(d);
+            return;
+        }
+        else if (fmprb_is_nonpositive(a))
+        {
+            fmprb_neg(d, a);
+            fmprb_rsqrt(d, d, prec);
+            fmprb_neg(d, d);
+            fmprb_zero(c);
+            return;
+        }
+    }
+
+    if (fmprb_is_zero(a))
+    {
+        if (fmprb_is_nonnegative(b))
+        {
+            fmprb_mul_2exp_si(c, b, 1);
+            fmprb_rsqrt(c, c, prec);
+            fmprb_neg(d, c);
+            return;
+        }
+        else if (fmprb_is_nonpositive(b))
+        {
+            fmprb_mul_2exp_si(c, b, 1);
+            fmprb_neg(c, c);
+            fmprb_rsqrt(c, c, prec);
+            fmprb_set(d, c);
+            return;
+        }
+    }
+
+    /* based on the identity sqrt(z) = sqrt(r) (z+r) / |z+r|: */
+    /* 1/sqrt(a+bi) = (1/v)((a+r) - b*i), r = |a+bi|, v = sqrt(r*(b^2+(a+r)^2)) */
+
+    wp = prec + 6;
+
+    fmprb_init(r);
+    fmprb_init(t);
+    fmprb_init(u);
+    fmprb_init(v);
+
+    /* u = b^2, r = |a+bi| */
+    fmprb_mul(t, a, a, wp);
+    fmprb_mul(u, b, b, wp);
+    fmprb_add(r, t, u, wp);
+    fmprb_sqrtpos(r, r, wp);
+
+    /* t = a+r, v = r*(b^2+(a+r)^2) */
+    fmprb_add(t, r, a, wp);
+    fmprb_mul(v, t, t, wp);
+    fmprb_add(v, v, u, wp);
+    fmprb_mul(v, v, r, wp);
+
+    /* v = 1/sqrt(v) */
+    fmprb_rsqrt(v, v, wp);
+
+    fmprb_mul(c, t, v, prec);
+    fmprb_mul(d, b, v, prec);
+    fmprb_neg(d, d);
+
+    fmprb_clear(r);
+    fmprb_clear(t);
+    fmprb_clear(u);
+    fmprb_clear(v);
+
+#undef a
+#undef b
+#undef c
+#undef d
+}
+

fmpcb/test/t-rsqrt.c

@@ -0,0 +1,85 @@
+/*=============================================================================
+
+    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"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("rsqrt....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    /* check (a^(-1/2))^(-2) = a */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        fmpcb_t a, b, c;
+        long prec;
+
+        fmpcb_init(a);
+        fmpcb_init(b);
+        fmpcb_init(c);
+
+        fmpcb_randtest(a, state, 1 + n_randint(state, 2000), 10);
+        fmpcb_randtest(b, state, 1 + n_randint(state, 2000), 10);
+
+        prec = 2 + n_randint(state, 2000);
+
+        fmpcb_rsqrt(b, a, prec);
+        fmpcb_inv(c, b, prec);
+        fmpcb_mul(c, c, c, prec);
+
+        if (!fmpcb_contains(c, a))
+        {
+            printf("FAIL: containment\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");
+            abort();
+        }
+
+        fmpcb_rsqrt(a, a, prec);
+        if (!fmpcb_equal(a, b))
+        {
+            printf("FAIL: aliasing\n\n");
+            printf("a = "); fmpcb_print(a); printf("\n\n");
+            printf("b = "); fmpcb_print(b); printf("\n\n");
+            abort();
+        }
+
+        fmpcb_clear(a);
+        fmpcb_clear(b);
+        fmpcb_clear(c);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

fmpcb_poly/rsqrt_series.c

@@ -31,9 +31,7 @@ _fmpcb_poly_rsqrt_series(fmpcb_ptr g,
 {
     hlen = FLINT_MIN(hlen, len);
 
-    /* TODO: write an fmpcb_rsqrt */
-    fmpcb_sqrt(g, h, prec);
-    fmpcb_inv(g, g, prec);
+    fmpcb_rsqrt(g, h, prec);
 
     if (hlen == 1)
     {