hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-05 09:21:38 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

add fmprb_atan2 with test code; use for and fix bug in fmpcb_arg

Browse source
This commit is contained in:
Fredrik Johansson 2013-02-19 11:17:21 +01:00
parent 800663b950
commit 6a6136fa99
6 changed files with 239 additions and 110 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

8
doc/source/fmpcb.rst
View file

@ -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)

7
doc/source/fmprb.rst
View file

@ -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)

108
fmpcb/arg.c
View file

@ -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);
}

1
fmprb.h
View file

@ -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);

137
fmprb/atan2.c Normal file
View file

@ -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);
}
}

88
fmprb/test/t-atan2.c Normal file
View file

@ -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;
}