From ee9a5eb8fc8bce623e19cbe579e91b5f1b44f107 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Fri, 22 Mar 2013 13:32:48 +0100 Subject: [PATCH] add complex square roots --- doc/source/fmpcb.rst | 8 ++++ fmpcb.h | 2 + fmpcb/sqrt.c | 102 +++++++++++++++++++++++++++++++++++++++++++ fmpcb/test/t-sqrt.c | 84 +++++++++++++++++++++++++++++++++++ 4 files changed, 196 insertions(+) create mode 100644 fmpcb/sqrt.c create mode 100644 fmpcb/test/t-sqrt.c diff --git a/doc/source/fmpcb.rst b/doc/source/fmpcb.rst index c3d2a09b..dcf80a5a 100644 --- a/doc/source/fmpcb.rst +++ b/doc/source/fmpcb.rst @@ -314,6 +314,14 @@ Elementary functions Sets *r* to *x* raised to the power *y*, computed as `x^y = \exp(y \log x)`. +.. function:: void fmpcb_sqrt(fmpcb_t r, const fmpcb_t z, long prec) + + Sets *r* to the square root of *z*. + If either the real or imaginary part is exactly zero, only + a single real square root is needed. Generally, we use the formula + `\sqrt{a+bi} = u/2 + ib/u, u = \sqrt{2(|a+bi|+a)}`, + requiring two real square root extractions. + .. 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 cca141c6..10503528 100644 --- a/fmpcb.h +++ b/fmpcb.h @@ -533,6 +533,8 @@ void fmpcb_sin_cos_pi(fmpcb_t s, fmpcb_t c, const fmpcb_t z, 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_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); void fmpcb_root_newton(fmpcb_t r, const fmpcb_t a, long m, long index, long prec); diff --git a/fmpcb/sqrt.c b/fmpcb/sqrt.c new file mode 100644 index 00000000..eae65f49 --- /dev/null +++ b/fmpcb/sqrt.c @@ -0,0 +1,102 @@ +/*============================================================================= + + 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_sqrt(fmpcb_t y, const fmpcb_t x, long prec) +{ + fmprb_t r, t, u; + 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_sqrt(c, a, prec); + fmprb_zero(d); + return; + } + else if (fmprb_is_nonpositive(a)) + { + fmprb_neg(d, a); + fmprb_sqrt(d, d, prec); + fmprb_zero(c); + return; + } + } + + if (fmprb_is_zero(a)) + { + if (fmprb_is_nonnegative(b)) + { + fmprb_mul_2exp_si(c, b, -1); + fmprb_sqrt(c, c, prec); + fmprb_set(d, c); + return; + } + else if (fmprb_is_nonpositive(b)) + { + fmprb_mul_2exp_si(c, b, -1); + fmprb_neg(c, c); + fmprb_sqrt(c, c, prec); + fmprb_neg(d, c); + return; + } + } + + /* sqrt(a+bi) = sqrt((r+a)/2) + b/sqrt(2*(r+a))*i, r = |a+bi| */ + + wp = prec + 4; + + fmprb_init(r); + fmprb_init(t); + fmprb_init(u); + + fmpcb_abs(r, x, wp); + fmprb_add(t, r, a, wp); + + fmprb_mul_2exp_si(u, t, 1); + fmprb_sqrt(u, u, wp); + fmprb_div(d, b, u, prec); + + fmprb_set_round(c, u, prec); + fmprb_mul_2exp_si(c, c, -1); + + fmprb_clear(r); + fmprb_clear(t); + fmprb_clear(u); + +#undef a +#undef b +#undef c +#undef d +} + diff --git a/fmpcb/test/t-sqrt.c b/fmpcb/test/t-sqrt.c new file mode 100644 index 00000000..eeeb6f62 --- /dev/null +++ b/fmpcb/test/t-sqrt.c @@ -0,0 +1,84 @@ +/*============================================================================= + + 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("sqrt...."); + 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_sqrt(b, a, prec); + fmpcb_mul(c, b, b, 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_sqrt(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); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} +