1 /*
2 Copyright (C) 2013, 2018 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 /* r - |m| */
15 void
arb_get_mag_reverse(mag_t res,const arb_t x)16 arb_get_mag_reverse(mag_t res, const arb_t x)
17 {
18 mag_t t;
19 mag_init(t);
20 arf_get_mag_lower(t, arb_midref(x));
21 mag_sub(res, arb_radref(x), t);
22 mag_clear(t);
23 }
24
25 /* upper bound for re(rsqrt(x+yi)) / |rsqrt(x+yi)|,
26 given upper bound for x, lower bound for y */
27 void
mag_rsqrt_re_quadrant1_upper(mag_t res,const mag_t x,const mag_t y)28 mag_rsqrt_re_quadrant1_upper(mag_t res, const mag_t x, const mag_t y)
29 {
30 if (mag_is_zero(x))
31 {
32 mag_one(res);
33 mag_mul_2exp_si(res, res, -1);
34 }
35 else
36 {
37 mag_t t, u;
38 mag_init(t);
39 mag_init(u);
40
41 /* t = (y/x)^2 -- the result is a decreasing function of t */
42 mag_div_lower(t, y, x);
43 mag_mul_lower(t, t, t);
44
45 /* (rsqrt(t^2+1)+1)/2 */
46 mag_add_ui_lower(u, t, 1);
47 mag_rsqrt(u, u);
48 mag_add_ui(u, u, 1);
49 mag_mul_2exp_si(res, u, -1);
50
51 mag_clear(t);
52 mag_clear(u);
53 }
54
55 mag_sqrt(res, res);
56 }
57
58 /* lower bound for re(rsqrt(x+yi)) / |rsqrt(x+yi)|,
59 given lower bound for x, upper bound for y */
60 void
mag_rsqrt_re_quadrant1_lower(mag_t res,const mag_t x,const mag_t y)61 mag_rsqrt_re_quadrant1_lower(mag_t res, const mag_t x, const mag_t y)
62 {
63 if (mag_is_zero(x))
64 {
65 mag_one(res);
66 mag_mul_2exp_si(res, res, -1);
67 }
68 else
69 {
70 mag_t t, u;
71 mag_init(t);
72 mag_init(u);
73
74 /* t = (y/x)^2 -- the result is a decreasing function of t */
75 mag_div(t, y, x);
76 mag_mul(t, t, t);
77
78 /* (rsqrt(t^2+1)+1)/2 */
79 mag_add_ui(u, t, 1);
80 mag_rsqrt_lower(u, u);
81 mag_add_ui_lower(u, u, 1);
82 mag_mul_2exp_si(res, u, -1);
83
84 mag_clear(t);
85 mag_clear(u);
86 }
87
88 mag_sqrt_lower(res, res);
89 }
90
91 /* upper bound for re(rsqrt(-x+yi)) / |rsqrt(x+yi)|,
92 given lower bound for -x, upper bound for y */
93 void
mag_rsqrt_re_quadrant2_upper(mag_t res,const mag_t x,const mag_t y)94 mag_rsqrt_re_quadrant2_upper(mag_t res, const mag_t x, const mag_t y)
95 {
96 if (mag_is_zero(x))
97 {
98 mag_one(res);
99 mag_mul_2exp_si(res, res, -1);
100 }
101 else
102 {
103 mag_t t, u, v;
104 mag_init(t);
105 mag_init(u);
106 mag_init(v);
107
108 /* t = (y/x)^2 -- the result is an increasing function of t */
109 mag_div(t, y, x);
110 mag_mul(t, t, t);
111
112 /* t / (2*(t+1)*(rsqrt(t+1)+1)) */
113 mag_add_ui(u, t, 1);
114 mag_rsqrt_lower(v, u);
115 mag_add_ui_lower(v, v, 1);
116 mag_add_ui_lower(u, t, 1);
117 mag_mul_lower(v, v, u);
118 mag_mul_2exp_si(v, v, 1);
119 mag_div(res, t, v);
120
121 mag_clear(t);
122 mag_clear(u);
123 mag_clear(v);
124 }
125
126 mag_sqrt(res, res);
127 }
128
129 /* lower bound for re(rsqrt(-x+yi)) / |rsqrt(x+yi)|,
130 given upper bound for -x, lower bound for y */
131 void
mag_rsqrt_re_quadrant2_lower(mag_t res,const mag_t x,const mag_t y)132 mag_rsqrt_re_quadrant2_lower(mag_t res, const mag_t x, const mag_t y)
133 {
134 if (mag_is_zero(x))
135 {
136 mag_one(res);
137 mag_mul_2exp_si(res, res, -1);
138 }
139 else
140 {
141 mag_t t, u, v;
142 mag_init(t);
143 mag_init(u);
144 mag_init(v);
145
146 /* t = (y/x)^2 -- the result is an increasing function of t */
147 mag_div_lower(t, y, x);
148 mag_mul_lower(t, t, t);
149
150 /* t / (2*(t+1)*(rsqrt(t+1)+1)) */
151 mag_add_ui_lower(u, t, 1);
152 mag_rsqrt(v, u);
153 mag_add_ui(v, v, 1);
154 mag_add_ui(u, t, 1);
155 mag_mul(v, v, u);
156 mag_mul_2exp_si(v, v, 1);
157 mag_div_lower(res, t, v);
158
159 mag_clear(t);
160 mag_clear(u);
161 mag_clear(v);
162 }
163
164 mag_sqrt_lower(res, res);
165 }
166
167 void
acb_rsqrt_wide(acb_t res,const acb_t z,slong prec)168 acb_rsqrt_wide(acb_t res, const acb_t z, slong prec)
169 {
170 mag_t ax, ay, bx, by, cx, cy, dx, dy, am, bm;
171 mag_t one;
172
173 mag_init(ax); mag_init(ay); mag_init(bx); mag_init(by);
174 mag_init(cx); mag_init(cy); mag_init(dx); mag_init(dy);
175 mag_init(am); mag_init(bm);
176 mag_init(one);
177
178 mag_one(one);
179
180 /* magnitude */
181 acb_get_mag(am, z);
182 mag_rsqrt_lower(am, am);
183 acb_get_mag_lower(bm, z);
184 mag_rsqrt(bm, bm);
185
186 /* upper or lower half plane */
187 if (arb_is_nonnegative(acb_imagref(z)) || arb_is_negative(acb_imagref(z)))
188 {
189 if (arb_is_nonnegative(acb_realref(z)))
190 {
191 arb_get_mag_lower(ax, acb_realref(z));
192 arb_get_mag(ay, acb_imagref(z));
193 arb_get_mag(bx, acb_realref(z));
194 arb_get_mag_lower(by, acb_imagref(z));
195
196 mag_rsqrt_re_quadrant2_lower(cx, bx, by);
197 mag_rsqrt_re_quadrant2_upper(dx, ax, ay);
198
199 /* equivalent but more expensive than pythagoras
200 mag_rsqrt_re_quadrant1_lower(ax, ax, ay);
201 mag_rsqrt_re_quadrant1_upper(bx, bx, by);
202 */
203
204 mag_mul(ax, dx, dx);
205 mag_sub_lower(ax, one, ax);
206 mag_sqrt_lower(ax, ax);
207 mag_mul_lower(bx, cx, cx);
208 mag_sub(bx, one, bx);
209 mag_sqrt(bx, bx);
210 }
211 else
212 {
213 if (arb_is_nonpositive(acb_realref(z)))
214 {
215 arb_get_mag(ax, acb_realref(z));
216 arb_get_mag_lower(ay, acb_imagref(z));
217 arb_get_mag_lower(bx, acb_realref(z));
218 arb_get_mag(by, acb_imagref(z));
219
220 /* equivalent but more expensive than pythagoras
221 mag_rsqrt_re_quadrant1_lower(cx, bx, by);
222 mag_rsqrt_re_quadrant1_upper(dx, ax, ay);
223 */
224
225 mag_rsqrt_re_quadrant2_lower(ax, ax, ay);
226 mag_rsqrt_re_quadrant2_upper(bx, bx, by);
227 }
228 else if (arf_sgn(arb_midref(acb_realref(z))) >= 0)
229 {
230 arb_get_mag_reverse(ax, acb_realref(z));
231 arb_get_mag_lower(ay, acb_imagref(z));
232 arb_get_mag(bx, acb_realref(z));
233 arb_get_mag_lower(by, acb_imagref(z));
234
235 mag_rsqrt_re_quadrant2_lower(ax, ax, ay);
236 mag_rsqrt_re_quadrant1_upper(bx, bx, by);
237 }
238 else
239 {
240 arb_get_mag(ax, acb_realref(z));
241 arb_get_mag_lower(ay, acb_imagref(z));
242 arb_get_mag_reverse(bx, acb_realref(z));
243 arb_get_mag_lower(by, acb_imagref(z));
244
245 mag_rsqrt_re_quadrant2_lower(ax, ax, ay);
246 mag_rsqrt_re_quadrant1_upper(bx, bx, by);
247 }
248
249 /* pythagoras */
250 mag_mul(cx, bx, bx);
251 mag_sub_lower(cx, one, bx);
252 mag_sqrt_lower(cx, cx);
253 mag_mul_lower(dx, ax, ax);
254 mag_sub(dx, one, dx);
255 mag_sqrt(dx, dx);
256 }
257
258 mag_mul_lower(ax, ax, am);
259 mag_mul_lower(cx, cx, am);
260 mag_mul(bx, bx, bm);
261 mag_mul(dx, dx, bm);
262
263 if (arf_sgn(arb_midref(acb_imagref(z))) > 0)
264 {
265 arb_set_interval_mag(acb_realref(res), ax, bx, prec);
266 arb_set_interval_mag(acb_imagref(res), cx, dx, prec);
267 arf_neg(arb_midref(acb_imagref(res)), arb_midref(acb_imagref(res)));
268 }
269 else
270 {
271 arb_set_interval_mag(acb_realref(res), ax, bx, prec);
272 arb_set_interval_mag(acb_imagref(res), cx, dx, prec);
273 }
274 }
275 else if (arb_is_positive(acb_realref(z)))
276 {
277 /* right half plane, straddling real line */
278 int symmetric;
279
280 symmetric = arf_is_zero(arb_midref(acb_imagref(z)));
281
282 arb_get_mag_lower(ax, acb_realref(z));
283 arb_get_mag(dy, acb_imagref(z));
284 arb_get_mag_reverse(cy, acb_imagref(z));
285
286 if (!symmetric)
287 mag_rsqrt_re_quadrant2_lower(cx, ax, cy);
288 mag_rsqrt_re_quadrant2_upper(dx, ax, dy);
289
290 mag_one(bx);
291 /* mag_rsqrt_re_quadrant1_lower(ax, ax, dy); */
292 mag_mul(ax, dx, dx);
293 mag_sub_lower(ax, one, ax);
294 mag_sqrt_lower(ax, ax);
295
296 mag_mul_lower(ax, ax, am);
297 mag_mul(bx, bx, bm);
298 mag_mul(cx, cx, bm);
299 mag_mul(dx, dx, bm);
300
301 if (symmetric)
302 arb_set_interval_neg_pos_mag(acb_imagref(res), dx, dx, prec);
303 else if (arf_sgn(arb_midref(acb_imagref(z))) > 0)
304 arb_set_interval_neg_pos_mag(acb_imagref(res), dx, cx, prec);
305 else
306 arb_set_interval_neg_pos_mag(acb_imagref(res), cx, dx, prec);
307
308 arb_set_interval_mag(acb_realref(res), ax, bx, prec);
309 }
310 else /* left half plane, straddling branch cut */
311 {
312 mag_zero(ax);
313 arb_get_mag_lower(bx, acb_realref(z));
314 arb_get_mag(by, acb_imagref(z));
315 mag_rsqrt_re_quadrant2_upper(bx, bx, by);
316
317 mag_mul_lower(ax, ax, am);
318 mag_mul(bx, bx, bm);
319 arb_set_interval_mag(acb_realref(res), ax, bx, prec);
320
321 /* cx, dx = 1,1 */
322 arb_set_interval_neg_pos_mag(acb_imagref(res), bm, bm, prec);
323 }
324
325 mag_clear(ax); mag_clear(ay); mag_clear(bx); mag_clear(by);
326 mag_clear(cx); mag_clear(cy); mag_clear(dx); mag_clear(dy);
327 mag_clear(am); mag_clear(bm);
328 mag_clear(one);
329 }
330
331 void
acb_rsqrt_precise(acb_t y,const acb_t x,slong prec)332 acb_rsqrt_precise(acb_t y, const acb_t x, slong prec)
333 {
334 #define a acb_realref(x)
335 #define b acb_imagref(x)
336 #define c acb_realref(y)
337 #define d acb_imagref(y)
338
339 arb_t r, t, u, v;
340 slong wp;
341
342 /* based on the identity sqrt(z) = sqrt(r) (z+r) / |z+r|: */
343 /* 1/sqrt(a+bi) = (1/v)((a+r) - b*i), r = |a+bi|, v = sqrt(r*(b^2+(a+r)^2)) */
344
345 wp = prec + 6;
346
347 arb_init(r);
348 arb_init(t);
349 arb_init(u);
350 arb_init(v);
351
352 /* u = b^2, r = |a+bi| */
353 arb_mul(t, a, a, wp);
354 arb_mul(u, b, b, wp);
355 arb_add(r, t, u, wp);
356 arb_sqrtpos(r, r, wp);
357
358 /* t = a+r, v = r*(b^2+(a+r)^2) */
359 arb_add(t, r, a, wp);
360 arb_mul(v, t, t, wp);
361 arb_add(v, v, u, wp);
362 arb_mul(v, v, r, wp);
363
364 /* v = 1/sqrt(v) */
365 arb_rsqrt(v, v, wp);
366
367 arb_mul(c, t, v, prec);
368 arb_mul(d, b, v, prec);
369 arb_neg(d, d);
370
371 arb_clear(r);
372 arb_clear(t);
373 arb_clear(u);
374 arb_clear(v);
375
376 #undef a
377 #undef b
378 #undef c
379 #undef d
380 }
381
382 void
acb_rsqrt(acb_t y,const acb_t x,slong prec)383 acb_rsqrt(acb_t y, const acb_t x, slong prec)
384 {
385 slong acc;
386
387 #define a acb_realref(x)
388 #define b acb_imagref(x)
389 #define c acb_realref(y)
390 #define d acb_imagref(y)
391
392 if (acb_contains_zero(x))
393 {
394 acb_indeterminate(y);
395 return;
396 }
397
398 if (arb_is_zero(b))
399 {
400 if (arb_is_nonnegative(a))
401 {
402 arb_rsqrt(c, a, prec);
403 arb_zero(d);
404 return;
405 }
406 else if (arb_is_nonpositive(a))
407 {
408 arb_neg(d, a);
409 arb_rsqrt(d, d, prec);
410 arb_neg(d, d);
411 arb_zero(c);
412 return;
413 }
414 }
415
416 if (arb_is_zero(a))
417 {
418 if (arb_is_nonnegative(b))
419 {
420 arb_mul_2exp_si(c, b, 1);
421 arb_rsqrt(c, c, prec);
422 arb_neg(d, c);
423 return;
424 }
425 else if (arb_is_nonpositive(b))
426 {
427 arb_mul_2exp_si(c, b, 1);
428 arb_neg(c, c);
429 arb_rsqrt(c, c, prec);
430 arb_set(d, c);
431 return;
432 }
433 }
434
435 acc = acb_rel_accuracy_bits(x);
436
437 if (acc < 25)
438 {
439 acb_rsqrt_wide(y, x, prec);
440 }
441 else
442 {
443 if (arb_is_positive(a))
444 {
445 acb_rsqrt_precise(y, x, prec);
446 }
447 else if (arb_is_nonnegative(b))
448 {
449 acb_neg(y, x);
450 acb_rsqrt_precise(y, y, prec);
451 acb_div_onei(y, y);
452 }
453 else if (arb_is_negative(b))
454 {
455 acb_neg(y, x);
456 acb_rsqrt_precise(y, y, prec);
457 acb_mul_onei(y, y);
458 }
459 else
460 {
461 acb_rsqrt_wide(y, x, prec);
462 }
463 }
464
465 #undef a
466 #undef b
467 #undef c
468 #undef d
469 }
470
471 void
acb_rsqrt_analytic(acb_ptr res,const acb_t z,int analytic,slong prec)472 acb_rsqrt_analytic(acb_ptr res, const acb_t z, int analytic, slong prec)
473 {
474 if (analytic && arb_contains_zero(acb_imagref(z)) &&
475 !arb_is_positive(acb_realref(z)))
476 {
477 acb_indeterminate(res);
478 }
479 else
480 {
481 acb_rsqrt(res, z, prec);
482 }
483 }
484
485