diff --git a/doc/source/fmpcb.rst b/doc/source/fmpcb.rst index a83e255b..23483b2e 100644 --- a/doc/source/fmpcb.rst +++ b/doc/source/fmpcb.rst @@ -152,10 +152,12 @@ Complex parts .. function:: void fmpcb_arg(fmprb_t r, const fmpcb_t z, long prec) - Sets *r* to a real interval containing the complex argument of *z*. We - define the complex argument have a discontinuity on `(-\infty,0]`, with + Sets *r* to a real interval containing the complex argument (phase) of *z*. + We define the complex argument have a discontinuity on `(-\infty,0]`, with the special value `\operatorname{arg}(0) = 0`, and - `\operatorname{arg}(x+0i) = \pi` for `x < 0`. + `\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*). .. function:: void fmpcb_abs(fmprb_t r, const fmpcb_t z, long prec) diff --git a/doc/source/fmprb.rst b/doc/source/fmprb.rst index 240d0352..e396b915 100644 --- a/doc/source/fmprb.rst +++ b/doc/source/fmprb.rst @@ -478,6 +478,13 @@ Elementary functions the propagated error is bounded by `r / (1 + d^2)` (this could be tightened). +.. function:: void fmprb_atan2(fmprb_t r, const fmprb_t b, const fmprb_t a, long prec) + + Sets *r* to an the argument (phase) of the complex number + `a + bi`, with the branch cut discontinuity on `(-\infty,0]`. + We define `\operatorname{atan2}(0,0) = 0`, and for `a < 0`, + `\operatorname{atan2}(0,a) = \pi`. + .. function:: void fmprb_sinh(fmprb_t s, const fmprb_t x, long prec) .. function:: void fmprb_cosh(fmprb_t c, const fmprb_t x, long prec) diff --git a/fmpcb/arg.c b/fmpcb/arg.c index 8eb85c0b..a3d858f1 100644 --- a/fmpcb/arg.c +++ b/fmpcb/arg.c @@ -28,112 +28,6 @@ void fmpcb_arg(fmprb_t r, const fmpcb_t z, long prec) { -#define a fmpcb_realref(z) -#define b fmpcb_imagref(z) -#define am fmprb_midref(a) -#define ar fmprb_radref(a) -#define bm fmprb_midref(b) -#define br fmprb_radref(b) - - /* a real number */ - if (fmprb_is_zero(b)) - { - /* exactly zero */ - if (fmprb_is_zero(a)) - { - /* define arg(0) = 0 by convention */ - fmprb_zero(r); - } - /* interval contains only nonnegative numbers */ - else if (fmpr_sgn(am) > 0 && fmpr_cmpabs(am, ar) >= 0) - { - fmprb_zero(r); - } - /* interval contains only positive numbers */ - else if (fmpr_sgn(am) < 0 && fmpr_cmpabs(am, ar) > 0) - { - fmprb_const_pi(r, prec); - } - else - { - /* both positive and negative -- argument will be in [0, pi] */ - fmprb_t t; - fmprb_init(t); - fmprb_const_pi(t, prec); - fmprb_mul_2exp_si(t, t, -1); - fmprb_zero(r); - fmprb_add_error(r, t); - fmprb_clear(t); - } - } - /* an imaginary number */ - else if (fmprb_is_zero(a)) - { - /* interval contains only positive numbers */ - if (fmpr_sgn(bm) > 0 && fmpr_cmpabs(bm, br) > 0) - { - fmprb_const_pi(r, prec); - fmprb_mul_2exp_si(r, r, -1); - } - /* interval contains only negative numbers */ - else if (fmpr_sgn(bm) < 0 && fmpr_cmpabs(bm, br) > 0) - { - fmprb_const_pi(r, prec); - fmprb_neg(r, r); - fmprb_mul_2exp_si(r, r, -1); - } - else - { - /* both positive and negative -- argument will be in 0 +/- pi/2 */ - fmprb_t t; - fmprb_init(t); - fmprb_const_pi(t, FMPRB_RAD_PREC); - fmprb_mul_2exp_si(t, t, -1); - fmprb_zero(r); - fmprb_add_error(r, t); - fmprb_clear(t); - } - } - /* strictly in the right half-plane -- atan(b/a) */ - else if (fmpr_sgn(am) > 0 && fmpr_cmpabs(am, ar) > 0) - { - fmprb_div(r, b, a, prec); - fmprb_atan(r, r, prec); - } - /* strictly in the upper half-plane -- pi/2 - atan(a/b) */ - else if (fmpr_sgn(bm) > 0 && fmpr_cmpabs(bm, br) > 0) - { - fmprb_t t; - fmprb_init(t); - fmprb_div(r, a, b, prec); - fmprb_atan(r, r, prec); - fmprb_const_pi(t, prec); - fmprb_mul_2exp_si(t, t, -1); - fmprb_sub(r, t, r, prec); - fmprb_clear(t); - } - /* strictly in the lower half-plane -- -pi/2 - atan(a/b) */ - else if (fmpr_sgn(bm) < 0 && fmpr_cmpabs(bm, br) > 0) - { - fmprb_t t; - fmprb_init(t); - fmprb_div(r, a, b, prec); - fmprb_atan(r, r, prec); - fmprb_const_pi(t, prec); - fmprb_mul_2exp_si(t, t, -1); - fmprb_add(r, t, r, prec); - fmprb_neg(r, r); - fmprb_clear(t); - } - /* overlaps the nonpositive half-axis -- [-pi, pi] */ - else - { - fmprb_t t; - fmprb_init(t); - fmprb_const_pi(t, FMPRB_RAD_PREC); - fmprb_zero(r); - fmprb_add_error(r, t); - fmprb_clear(t); - } + fmprb_atan2(r, fmpcb_imagref(z), fmpcb_realref(z), prec); } diff --git a/fmprb.h b/fmprb.h index f6ce6698..5330f120 100644 --- a/fmprb.h +++ b/fmprb.h @@ -276,6 +276,7 @@ void fmprb_cosh(fmprb_t z, const fmprb_t x, long prec); void fmprb_sinh_cosh(fmprb_t s, fmprb_t c, const fmprb_t x, long prec); void fmprb_atan(fmprb_t z, const fmprb_t x, long prec); +void fmprb_atan2(fmprb_t z, const fmprb_t b, const fmprb_t a, long prec); void fmprb_fac_ui(fmprb_t x, ulong n, long prec); diff --git a/fmprb/atan2.c b/fmprb/atan2.c new file mode 100644 index 00000000..86ddc26c --- /dev/null +++ b/fmprb/atan2.c @@ -0,0 +1,137 @@ +/*============================================================================= + + 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" + +void +fmprb_atan2(fmprb_t r, const fmprb_t b, const fmprb_t a, long prec) +{ +#define am fmprb_midref(a) +#define ar fmprb_radref(a) +#define bm fmprb_midref(b) +#define br fmprb_radref(b) + + /* a + bi is a real number */ + if (fmprb_is_zero(b)) + { + /* exactly zero */ + if (fmprb_is_zero(a)) + { + /* define arg(0) = 0 by convention */ + fmprb_zero(r); + } + /* interval contains only nonnegative numbers */ + else if (fmpr_sgn(am) > 0 && fmpr_cmpabs(am, ar) >= 0) + { + fmprb_zero(r); + } + /* interval contains only negative numbers */ + else if (fmpr_sgn(am) < 0 && fmpr_cmpabs(am, ar) > 0) + { + fmprb_const_pi(r, prec); + } + else + { + /* both positive and negative -- argument will be in [0, pi] */ + fmprb_t t; + fmprb_init(t); + fmprb_const_pi(t, prec); + fmprb_mul_2exp_si(t, t, -1); + fmprb_set(r, t); + fmprb_add_error(r, t); + fmprb_clear(t); + } + } + /* an imaginary number */ + else if (fmprb_is_zero(a)) + { + /* interval contains only positive numbers */ + if (fmpr_sgn(bm) > 0 && fmpr_cmpabs(bm, br) > 0) + { + fmprb_const_pi(r, prec); + fmprb_mul_2exp_si(r, r, -1); + } + /* interval contains only negative numbers */ + else if (fmpr_sgn(bm) < 0 && fmpr_cmpabs(bm, br) > 0) + { + fmprb_const_pi(r, prec); + fmprb_neg(r, r); + fmprb_mul_2exp_si(r, r, -1); + } + else + { + /* both positive and negative -- argument will be in 0 +/- pi/2 */ + fmprb_t t; + fmprb_init(t); + fmprb_const_pi(t, FMPRB_RAD_PREC); + fmprb_mul_2exp_si(t, t, -1); + fmprb_zero(r); + fmprb_add_error(r, t); + fmprb_clear(t); + } + } + /* strictly in the right half-plane -- atan(b/a) */ + else if (fmpr_sgn(am) > 0 && fmpr_cmpabs(am, ar) > 0) + { + fmprb_div(r, b, a, prec); + fmprb_atan(r, r, prec); + } + /* strictly in the upper half-plane -- pi/2 - atan(a/b) */ + else if (fmpr_sgn(bm) > 0 && fmpr_cmpabs(bm, br) > 0) + { + fmprb_t t; + fmprb_init(t); + fmprb_div(r, a, b, prec); + fmprb_atan(r, r, prec); + fmprb_const_pi(t, prec); + fmprb_mul_2exp_si(t, t, -1); + fmprb_sub(r, t, r, prec); + fmprb_clear(t); + } + /* strictly in the lower half-plane -- -pi/2 - atan(a/b) */ + else if (fmpr_sgn(bm) < 0 && fmpr_cmpabs(bm, br) > 0) + { + fmprb_t t; + fmprb_init(t); + fmprb_div(r, a, b, prec); + fmprb_atan(r, r, prec); + fmprb_const_pi(t, prec); + fmprb_mul_2exp_si(t, t, -1); + fmprb_add(r, t, r, prec); + fmprb_neg(r, r); + fmprb_clear(t); + } + /* overlaps the nonpositive half-axis -- [-pi, pi] */ + else + { + fmprb_t t; + fmprb_init(t); + fmprb_const_pi(t, FMPRB_RAD_PREC); + fmprb_zero(r); + fmprb_add_error(r, t); + fmprb_clear(t); + } +} + diff --git a/fmprb/test/t-atan2.c b/fmprb/test/t-atan2.c new file mode 100644 index 00000000..b2fd2f5a --- /dev/null +++ b/fmprb/test/t-atan2.c @@ -0,0 +1,88 @@ +/*============================================================================= + + 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 "fmprb.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("atan2...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000; iter++) + { + fmprb_t a, b, c; + fmpq_t q, r; + mpfr_t t, u; + long prec = 2 + n_randint(state, 200); + + fmprb_init(a); + fmprb_init(b); + fmprb_init(c); + fmpq_init(q); + fmpq_init(r); + mpfr_init2(t, prec + 100); + mpfr_init2(u, prec + 100); + + fmprb_randtest(a, state, 1 + n_randint(state, 200), 3); + fmprb_randtest(b, state, 1 + n_randint(state, 200), 3); + fmprb_randtest(c, state, 1 + n_randint(state, 200), 3); + + fmprb_get_rand_fmpq(q, state, a, 1 + n_randint(state, 200)); + fmprb_get_rand_fmpq(r, state, b, 1 + n_randint(state, 200)); + + fmpq_get_mpfr(t, q, MPFR_RNDN); + fmpq_get_mpfr(u, r, MPFR_RNDN); + mpfr_atan2(t, u, t, MPFR_RNDN); + + fmprb_atan2(c, b, a, prec); + + if (!fmprb_contains_mpfr(c, t)) + { + printf("FAIL: containment\n\n"); + printf("a = "); fmprb_print(a); printf("\n\n"); + printf("b = "); fmprb_print(b); printf("\n\n"); + printf("c = "); fmprb_print(c); printf("\n\n"); + abort(); + } + + fmprb_clear(a); + fmprb_clear(b); + fmprb_clear(c); + fmpq_clear(q); + fmpq_clear(r); + mpfr_clear(t); + mpfr_clear(u); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +}