From 9ded02956b38520a0a5ab0f9eb6a87b70dee46f8 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Mon, 21 Oct 2013 17:58:42 +0200 Subject: [PATCH] add fmpcb_rsqrt (complex reciprocal square root) --- doc/source/fmpcb.rst | 8 +++ fmpcb.h | 1 + fmpcb/rsqrt.c | 116 ++++++++++++++++++++++++++++++++++++++ fmpcb/test/t-rsqrt.c | 85 ++++++++++++++++++++++++++++ fmpcb_poly/rsqrt_series.c | 4 +- 5 files changed, 211 insertions(+), 3 deletions(-) create mode 100644 fmpcb/rsqrt.c create mode 100644 fmpcb/test/t-rsqrt.c diff --git a/doc/source/fmpcb.rst b/doc/source/fmpcb.rst index 99fa3c21..914471e9 100644 --- a/doc/source/fmpcb.rst +++ b/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 diff --git a/fmpcb.h b/fmpcb.h index fb62ba8d..aa79015f 100644 --- a/fmpcb.h +++ b/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); diff --git a/fmpcb/rsqrt.c b/fmpcb/rsqrt.c new file mode 100644 index 00000000..17d27b58 --- /dev/null +++ b/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 +} + diff --git a/fmpcb/test/t-rsqrt.c b/fmpcb/test/t-rsqrt.c new file mode 100644 index 00000000..20801edf --- /dev/null +++ b/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; +} + diff --git a/fmpcb_poly/rsqrt_series.c b/fmpcb_poly/rsqrt_series.c index d7e49963..dcd9ee6b 100644 --- a/fmpcb_poly/rsqrt_series.c +++ b/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) {