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 static void
acb_log_sin_pi_half(acb_t res,const acb_t z,slong prec,int upper)15 acb_log_sin_pi_half(acb_t res, const acb_t z, slong prec, int upper)
16 {
17 acb_t t, u, zmid;
18 arf_t n;
19 arb_t pi;
20
21 acb_init(t);
22 acb_init(u);
23 acb_init(zmid);
24 arf_init(n);
25 arb_init(pi);
26
27 arf_set(arb_midref(acb_realref(zmid)), arb_midref(acb_realref(z)));
28 arf_set(arb_midref(acb_imagref(zmid)), arb_midref(acb_imagref(z)));
29
30 arf_floor(n, arb_midref(acb_realref(zmid)));
31 arb_sub_arf(acb_realref(zmid), acb_realref(zmid), n, prec);
32
33 arb_const_pi(pi, prec);
34
35 if (arf_cmpabs_2exp_si(arb_midref(acb_imagref(zmid)), 2) < 1)
36 {
37 acb_sin_pi(t, zmid, prec);
38 acb_log(t, t, prec);
39 }
40 else /* i*pi*(z-0.5) + log((1-exp(-2i*pi*z))/2) */
41 {
42 acb_mul_2exp_si(t, zmid, 1);
43 acb_neg(t, t);
44
45 if (upper)
46 acb_conj(t, t);
47
48 acb_exp_pi_i(t, t, prec);
49 acb_sub_ui(t, t, 1, prec);
50 acb_neg(t, t);
51
52 acb_mul_2exp_si(t, t, -1);
53
54 acb_log(t, t, prec);
55 acb_one(u);
56 acb_mul_2exp_si(u, u, -1);
57 acb_sub(u, zmid, u, prec);
58 if (upper)
59 acb_conj(u, u);
60 acb_mul_onei(u, u);
61 acb_addmul_arb(t, u, pi, prec);
62 if (upper)
63 acb_conj(t, t);
64 }
65
66 if (upper)
67 arb_submul_arf(acb_imagref(t), pi, n, prec);
68 else
69 arb_addmul_arf(acb_imagref(t), pi, n, prec);
70
71 /* propagated error bound from the derivative pi cot(pi z) */
72 if (!acb_is_exact(z))
73 {
74 mag_t zm, um;
75
76 mag_init(zm);
77 mag_init(um);
78
79 acb_cot_pi(u, z, prec);
80 acb_mul_arb(u, u, pi, prec);
81
82 mag_hypot(zm, arb_radref(acb_realref(z)), arb_radref(acb_imagref(z)));
83 acb_get_mag(um, u);
84 mag_mul(um, um, zm);
85
86 acb_add_error_mag(t, um);
87
88 mag_clear(zm);
89 mag_clear(um);
90 }
91
92 acb_set(res, t);
93
94 acb_clear(t);
95 acb_clear(u);
96 acb_clear(zmid);
97 arf_clear(n);
98 arb_clear(pi);
99 }
100
101 void
acb_log_sin_pi(acb_t res,const acb_t z,slong prec)102 acb_log_sin_pi(acb_t res, const acb_t z, slong prec)
103 {
104 if (!acb_is_finite(z))
105 {
106 acb_indeterminate(res);
107 return;
108 }
109
110 if (arb_is_positive(acb_imagref(z)) ||
111 (arb_is_zero(acb_imagref(z)) && arb_is_negative(acb_realref(z))))
112 {
113 acb_log_sin_pi_half(res, z, prec, 1);
114 }
115 else if (arb_is_negative(acb_imagref(z)) ||
116 (arb_is_zero(acb_imagref(z)) && arb_is_positive(acb_realref(z))))
117 {
118 acb_log_sin_pi_half(res, z, prec, 0);
119 }
120 else
121 {
122 acb_t t;
123 acb_init(t);
124 acb_log_sin_pi_half(t, z, prec, 1);
125 acb_log_sin_pi_half(res, z, prec, 0);
126 arb_union(acb_realref(res), acb_realref(res), acb_realref(t), prec);
127 arb_union(acb_imagref(res), acb_imagref(res), acb_imagref(t), prec);
128 acb_clear(t);
129 }
130 }
131
132