1 /*
2 Copyright (C) 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 "acb.h"
13
14 void
acb_sqrt(acb_t y,const acb_t x,slong prec)15 acb_sqrt(acb_t y, const acb_t x, slong prec)
16 {
17 arb_t r, t, u;
18 slong wp;
19 int done;
20
21 #define a acb_realref(x)
22 #define b acb_imagref(x)
23 #define c acb_realref(y)
24 #define d acb_imagref(y)
25
26 if (arb_is_zero(b))
27 {
28 if (arb_is_nonnegative(a))
29 {
30 arb_sqrt(c, a, prec);
31 arb_zero(d);
32 return;
33 }
34 else if (arb_is_nonpositive(a))
35 {
36 arb_neg(d, a);
37 arb_sqrt(d, d, prec);
38 arb_zero(c);
39 return;
40 }
41 }
42
43 if (arb_is_zero(a))
44 {
45 if (arb_is_nonnegative(b))
46 {
47 arb_mul_2exp_si(c, b, -1);
48 arb_sqrt(c, c, prec);
49 arb_set(d, c);
50 return;
51 }
52 else if (arb_is_nonpositive(b))
53 {
54 arb_mul_2exp_si(c, b, -1);
55 arb_neg(c, c);
56 arb_sqrt(c, c, prec);
57 arb_neg(d, c);
58 return;
59 }
60 }
61
62 wp = prec + 4;
63
64 arb_init(r);
65 arb_init(t);
66 arb_init(u);
67
68 /* r = |a+bi| */
69 acb_abs(r, x, wp);
70
71 done = 0;
72
73 if (arf_sgn(arb_midref(a)) >= 0)
74 {
75 arb_add(t, r, a, wp);
76
77 if (arb_rel_accuracy_bits(t) > 8)
78 {
79 /* sqrt(a+bi) = sqrt((r+a)/2) + b/sqrt(2*(r+a))*i */
80 arb_mul_2exp_si(u, t, 1);
81 arb_sqrt(u, u, wp);
82 arb_div(d, b, u, prec);
83
84 arb_set_round(c, u, prec);
85 arb_mul_2exp_si(c, c, -1);
86 done = 1;
87 }
88 else
89 {
90 arb_sub(u, r, a, wp);
91 }
92 }
93 else if (!arb_contains_zero(b))
94 {
95 arb_sub(u, r, a, wp);
96
97 if (arb_rel_accuracy_bits(u) > 8)
98 {
99 /* sqrt(a+bi) = |b|/sqrt(2*(r-a)) + sgn(b)*sqrt((r-a)/2)*i */
100 int sgn = arf_sgn(arb_midref(b));
101
102 arb_mul_2exp_si(t, u, 1);
103 arb_sqrt(t, t, wp);
104 arb_div(c, b, t, prec);
105 arb_abs(c, c);
106
107 arb_set_round(d, t, prec);
108 arb_mul_2exp_si(d, d, -1);
109 if (sgn < 0)
110 arb_neg(d, d);
111 done = 1;
112 }
113 else
114 {
115 arb_add(t, r, a, wp);
116 }
117 }
118 else
119 {
120 arb_add(t, r, a, wp);
121 arb_sub(u, r, a, wp);
122 }
123
124 /* t = r+a, u = r-a */
125
126 if (!done)
127 {
128 /* sqrt(a+bi) = sqrt((r+a)/2) + (b/|b|)*sqrt((r-a)/2)*i
129 (sign) */
130 arb_mul_2exp_si(t, t, -1);
131 arb_mul_2exp_si(u, u, -1);
132 arb_sqrtpos(c, t, prec);
133
134 if (arb_is_nonnegative(b))
135 {
136 arb_sqrtpos(d, u, prec);
137 }
138 else if (arb_is_nonpositive(b))
139 {
140 arb_sqrtpos(d, u, prec);
141 arb_neg(d, d);
142 }
143 else
144 {
145 arb_sqrtpos(t, u, wp);
146 arb_neg(u, t);
147 arb_union(d, t, u, prec);
148 }
149 }
150
151 arb_clear(r);
152 arb_clear(t);
153 arb_clear(u);
154
155 #undef a
156 #undef b
157 #undef c
158 #undef d
159 }
160
161 void
acb_sqrt_analytic(acb_ptr res,const acb_t z,int analytic,slong prec)162 acb_sqrt_analytic(acb_ptr res, const acb_t z, int analytic, slong prec)
163 {
164 if (analytic && arb_contains_zero(acb_imagref(z)) &&
165 !arb_is_positive(acb_realref(z)))
166 {
167 acb_indeterminate(res);
168 }
169 else
170 {
171 acb_sqrt(res, z, prec);
172 }
173 }
174
175