1 /*
2 Copyright (C) 2015 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 "acb.h"
13
14 /* branch cuts are on (+/-i)*[1,inf] */
15 int
acb_atan_on_branch_cut(const acb_t z)16 acb_atan_on_branch_cut(const acb_t z)
17 {
18 arb_t unit;
19 int result;
20
21 if (!acb_is_finite(z))
22 return 1;
23
24 if (arb_is_nonzero(acb_realref(z)))
25 return 0;
26
27 if (arb_contains_si(acb_imagref(z), 1) || arb_contains_si(acb_imagref(z), -1))
28 return 1;
29
30 arb_init(unit);
31 mag_one(arb_radref(unit));
32 result = !arb_contains(unit, acb_imagref(z));
33 arb_clear(unit);
34
35 return result;
36 }
37
38 void
acb_atan(acb_t r,const acb_t z,slong prec)39 acb_atan(acb_t r, const acb_t z, slong prec)
40 {
41 if (arb_is_zero(acb_imagref(z)))
42 {
43 arb_atan(acb_realref(r), acb_realref(z), prec);
44 arb_zero(acb_imagref(r));
45 }
46 else if (!acb_is_finite(z))
47 {
48 acb_indeterminate(r);
49 }
50 else
51 {
52 acb_t t, u;
53
54 acb_init(t);
55 acb_init(u);
56
57 if (acb_atan_on_branch_cut(z))
58 {
59 acb_mul_onei(u, z);
60 acb_neg(t, u);
61 acb_log1p(t, t, prec);
62 acb_log1p(u, u, prec);
63 acb_sub(t, t, u, prec);
64 acb_mul_onei(t, t);
65 acb_mul_2exp_si(r, t, -1);
66 }
67 else if (acb_is_exact(z))
68 {
69 acb_onei(t);
70 acb_sub(t, t, z, prec);
71 acb_div(t, z, t, prec);
72 acb_mul_2exp_si(t, t, 1);
73 acb_log1p(t, t, prec);
74 acb_mul_onei(t, t);
75 acb_mul_2exp_si(r, t, -1);
76 }
77 else
78 {
79 mag_t err, err2;
80 mag_init(err);
81 mag_init(err2);
82
83 /* atan'(z) = 1/(1+z^2) = 1/((z+i)(z-i)) */
84 acb_onei(t);
85 acb_add(t, z, t, prec);
86 acb_get_mag_lower(err, t);
87
88 acb_onei(t);
89 acb_sub(t, z, t, prec);
90 acb_get_mag_lower(err2, t);
91
92 mag_mul_lower(err, err, err2);
93
94 mag_hypot(err2, arb_radref(acb_realref(z)), arb_radref(acb_imagref(z)));
95 mag_div(err, err2, err);
96
97 arf_set(arb_midref(acb_realref(t)), arb_midref(acb_realref(z)));
98 arf_set(arb_midref(acb_imagref(t)), arb_midref(acb_imagref(z)));
99 mag_zero(arb_radref(acb_realref(t)));
100 mag_zero(arb_radref(acb_imagref(t)));
101 acb_atan(r, t, prec);
102 acb_add_error_mag(r, err);
103
104 mag_clear(err);
105 mag_clear(err2);
106 }
107
108 acb_clear(t);
109 acb_clear(u);
110 }
111 }
112
113