1 /*
2 Copyright (C) 2012, 2013 Fredrik Johansson
3
4 This file is part of Arb.
5
6 Arb is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include "arb.h"
13
14 void
arb_atan2(arb_t r,const arb_t b,const arb_t a,slong prec)15 arb_atan2(arb_t r, const arb_t b, const arb_t a, slong prec)
16 {
17 #define am arb_midref(a)
18 #define ar arb_radref(a)
19 #define bm arb_midref(b)
20 #define br arb_radref(b)
21
22 /* a + bi is a real number */
23 if (arb_is_zero(b))
24 {
25 /* exactly zero */
26 if (arb_is_zero(a))
27 {
28 /* define arg(0) = 0 by convention */
29 arb_zero(r);
30 }
31 /* interval contains only nonnegative numbers */
32 else if (arf_sgn(am) > 0 && arf_cmpabs_mag(am, ar) >= 0)
33 {
34 arb_zero(r);
35 }
36 /* interval contains only negative numbers */
37 else if (arf_sgn(am) < 0 && arf_cmpabs_mag(am, ar) > 0)
38 {
39 arb_const_pi(r, prec);
40 }
41 else
42 {
43 /* both positive and negative -- argument will be in [0, pi] */
44 mag_const_pi(arb_radref(r));
45 arf_set_mag(arb_midref(r), arb_radref(r));
46 arb_set_round(r, r, prec);
47 arb_mul_2exp_si(r, r, -1);
48 }
49 }
50 /* an imaginary number */
51 else if (arb_is_zero(a))
52 {
53 /* interval contains only positive numbers */
54 if (arf_sgn(bm) > 0 && arf_cmpabs_mag(bm, br) > 0)
55 {
56 arb_const_pi(r, prec);
57 arb_mul_2exp_si(r, r, -1);
58 }
59 /* interval contains only negative numbers */
60 else if (arf_sgn(bm) < 0 && arf_cmpabs_mag(bm, br) > 0)
61 {
62 arb_const_pi(r, prec);
63 arb_neg(r, r);
64 arb_mul_2exp_si(r, r, -1);
65 }
66 else
67 {
68 /* both positive and negative -- argument will be in 0 +/- pi/2 */
69 arf_zero(arb_midref(r));
70 mag_const_pi(arb_radref(r));
71 mag_mul_2exp_si(arb_radref(r), arb_radref(r), -1);
72 }
73 }
74 /* strictly in the right half-plane -- atan(b/a) */
75 else if (arf_sgn(am) > 0 && arf_cmpabs_mag(am, ar) > 0)
76 {
77 arb_div(r, b, a, prec);
78 arb_atan(r, r, prec);
79 }
80 /* strictly in the upper half-plane -- pi/2 - atan(a/b) */
81 else if (arf_sgn(bm) > 0 && arf_cmpabs_mag(bm, br) > 0)
82 {
83 arb_t t;
84 arb_init(t);
85 arb_div(r, a, b, prec);
86 arb_atan(r, r, prec);
87 arb_const_pi(t, prec);
88 arb_mul_2exp_si(t, t, -1);
89 arb_sub(r, t, r, prec);
90 arb_clear(t);
91 }
92 /* strictly in the lower half-plane -- -pi/2 - atan(a/b) */
93 else if (arf_sgn(bm) < 0 && arf_cmpabs_mag(bm, br) > 0)
94 {
95 arb_t t;
96 arb_init(t);
97 arb_div(r, a, b, prec);
98 arb_atan(r, r, prec);
99 arb_const_pi(t, prec);
100 arb_mul_2exp_si(t, t, -1);
101 arb_add(r, t, r, prec);
102 arb_neg(r, r);
103 arb_clear(t);
104 }
105 /* overlaps the nonpositive half-axis -- [-pi, pi] */
106 else
107 {
108 arf_zero(arb_midref(r));
109 mag_const_pi(arb_radref(r));
110 }
111 }
112
113