From d9c8a39a51d3ad189ab7053073cc9b18c033584e Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Thu, 18 Feb 2016 00:55:43 +0100 Subject: [PATCH] add complex sign function, acb_sgn --- acb.h | 2 + acb/sgn.c | 62 ++++++++++++++++++++++++++++++ acb/test/t-sgn.c | 95 ++++++++++++++++++++++++++++++++++++++++++++++ doc/source/acb.rst | 5 +++ 4 files changed, 164 insertions(+) create mode 100644 acb/sgn.c create mode 100644 acb/test/t-sgn.c diff --git a/acb.h b/acb.h index 2001a6b6..ff3b9e09 100644 --- a/acb.h +++ b/acb.h @@ -424,6 +424,8 @@ acb_get_rad_ubound_arf(arf_t u, const acb_t z, slong prec) void acb_arg(arb_t r, const acb_t z, slong prec); +void acb_sgn(acb_t res, const acb_t z, slong prec); + ACB_INLINE void acb_add(acb_t z, const acb_t x, const acb_t y, slong prec) { diff --git a/acb/sgn.c b/acb/sgn.c new file mode 100644 index 00000000..bd1585f4 --- /dev/null +++ b/acb/sgn.c @@ -0,0 +1,62 @@ +/*============================================================================= + + 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) 2016 Fredrik Johansson + +******************************************************************************/ + +#include "acb.h" + +void +acb_sgn(acb_t res, const acb_t z, slong prec) +{ + if (arb_is_zero(acb_imagref(z))) + { + arb_sgn(acb_realref(res), acb_realref(z)); + arb_zero(acb_imagref(res)); + } + else if (arb_is_zero(acb_realref(z))) + { + arb_sgn(acb_imagref(res), acb_imagref(z)); + arb_zero(acb_realref(res)); + } + else + { + arb_t t; + arb_init(t); + acb_abs(t, z, prec); + arb_inv(t, t, prec); + + if (arb_is_finite(t)) + { + acb_mul_arb(res, z, t, prec); + } + else + { + arf_zero(arb_midref(acb_realref(res))); + mag_one(arb_radref(acb_realref(res))); + arb_set(acb_imagref(res), acb_realref(res)); + } + + arb_clear(t); + } +} + diff --git a/acb/test/t-sgn.c b/acb/test/t-sgn.c new file mode 100644 index 00000000..6cb30b1c --- /dev/null +++ b/acb/test/t-sgn.c @@ -0,0 +1,95 @@ +/*============================================================================= + + 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) 2016 Fredrik Johansson + +******************************************************************************/ + +#include "acb.h" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("sgn...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + acb_t x, a, b; + slong prec; + + acb_init(x); + acb_init(a); + acb_init(b); + + acb_randtest_special(x, state, 1 + n_randint(state, 200), 2 + n_randint(state, 100)); + acb_randtest_special(a, state, 1 + n_randint(state, 200), 2 + n_randint(state, 100)); + + prec = 2 + n_randint(state, 200); + + acb_sgn(a, x, prec); + + if (acb_is_zero(x)) + { + acb_zero(b); + } + else + { + arb_zero(acb_realref(b)); + acb_arg(acb_imagref(b), x, prec); + acb_exp(b, b, prec); + } + + if (!acb_overlaps(a, b)) + { + flint_printf("FAIL: overlap\n\n"); + flint_printf("x = "); acb_printd(x, 15); flint_printf("\n\n"); + flint_printf("a = "); acb_printd(a, 15); flint_printf("\n\n"); + flint_printf("b = "); acb_printd(b, 15); flint_printf("\n\n"); + abort(); + } + + acb_sgn(x, x, prec); + + if (!acb_overlaps(a, x)) + { + flint_printf("FAIL: aliasing\n\n"); + flint_printf("x = "); acb_printd(x, 15); flint_printf("\n\n"); + flint_printf("a = "); acb_printd(a, 15); flint_printf("\n\n"); + flint_printf("b = "); acb_printd(b, 15); flint_printf("\n\n"); + abort(); + } + + acb_clear(x); + acb_clear(a); + acb_clear(b); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/doc/source/acb.rst b/doc/source/acb.rst index 0550205d..5849bfd0 100644 --- a/doc/source/acb.rst +++ b/doc/source/acb.rst @@ -331,6 +331,11 @@ Complex parts Sets *r* to the absolute value of *z*. +.. function:: void acb_sgn(arb_t r, const acb_t z, slong prec) + + Sets *r* to the complex sign of *z*, defined as 0 if *z* is exactly zero + and the projection onto the unit circle `z / |z| = \exp(i \arg(z))` otherwise. + Arithmetic -------------------------------------------------------------------------------