1 /*
2 Copyright (C) 2014 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 #include "acb_poly.h"
14
15 void mag_agm(mag_t res, const mag_t x, const mag_t y);
16
17 /* Checks that |arg(z)| <= 3 pi / 4 */
18 static int
acb_check_arg(const acb_t z)19 acb_check_arg(const acb_t z)
20 {
21 mag_t re, im;
22 int res;
23
24 if (!arb_contains_negative(acb_realref(z)))
25 return 1;
26
27 mag_init(re);
28 mag_init(im);
29
30 arb_get_mag(re, acb_realref(z));
31 arb_get_mag_lower(im, acb_imagref(z));
32
33 res = mag_cmp(re, im) < 0;
34
35 mag_clear(re);
36 mag_clear(im);
37
38 return res;
39 }
40
41 static void
sqrtmul(acb_t c,const acb_t a,const acb_t b,slong prec)42 sqrtmul(acb_t c, const acb_t a, const acb_t b, slong prec)
43 {
44 if (arb_is_positive(acb_realref(a)) &&
45 arb_is_positive(acb_realref(b)))
46 {
47 acb_mul(c, a, b, prec);
48 acb_sqrt(c, c, prec);
49 }
50 else if (arb_is_nonnegative(acb_imagref(a)) &&
51 arb_is_nonnegative(acb_imagref(b)))
52 {
53 acb_mul(c, a, b, prec);
54 acb_neg(c, c);
55 acb_sqrt(c, c, prec);
56 acb_mul_onei(c, c);
57 }
58 else if (arb_is_nonpositive(acb_imagref(a)) &&
59 arb_is_nonpositive(acb_imagref(b)))
60 {
61 acb_mul(c, a, b, prec);
62 acb_neg(c, c);
63 acb_sqrt(c, c, prec);
64 acb_mul_onei(c, c);
65 acb_neg(c, c);
66 }
67 else
68 {
69 acb_t d;
70 acb_init(d);
71 acb_sqrt(c, a, prec);
72 acb_sqrt(d, b, prec);
73 acb_mul(c, c, d, prec);
74 acb_clear(d);
75 }
76 }
77
78 static void
acb_agm_close_taylor(acb_t res,acb_t z,acb_t z2,const acb_t aplusb,const acb_t aminusb,const mag_t err,slong prec)79 acb_agm_close_taylor(acb_t res, acb_t z, acb_t z2,
80 const acb_t aplusb, const acb_t aminusb,
81 const mag_t err, slong prec)
82 {
83 acb_div(z, aminusb, aplusb, prec);
84 acb_sqr(z, z, prec);
85 acb_sqr(z2, z, prec);
86
87 acb_mul_si(res, z2, -469, prec);
88 acb_addmul_si(res, z, -704, prec);
89 acb_mul(res, res, z2, prec);
90 acb_addmul_si(res, z2, -1280, prec);
91 acb_mul_2exp_si(z, z, 12);
92 acb_sub(res, res, z, prec);
93 acb_add_ui(res, res, 16384, prec);
94 acb_mul_2exp_si(res, res, -15);
95
96 acb_add_error_mag(res, err);
97
98 acb_mul(res, res, aplusb, prec);
99 }
100
101 static void
acb_agm1_around_zero(acb_t res,const acb_t z,slong prec)102 acb_agm1_around_zero(acb_t res, const acb_t z, slong prec)
103 {
104 mag_t a, b;
105 mag_init(a);
106 mag_init(b);
107 mag_one(a);
108 acb_get_mag(b, z);
109 mag_agm(a, a, b);
110 acb_zero(res);
111 acb_add_error_mag(res, a);
112 mag_clear(a);
113 mag_clear(b);
114 }
115
116 void
acb_agm1_basecase(acb_t res,const acb_t z,slong prec)117 acb_agm1_basecase(acb_t res, const acb_t z, slong prec)
118 {
119 acb_t a, b, t, u;
120 mag_t err, err2;
121 int isreal;
122
123 isreal = acb_is_real(z) && arb_is_nonnegative(acb_realref(z));
124
125 if (isreal)
126 {
127 acb_init(a);
128 acb_one(a);
129 arb_agm(acb_realref(res), acb_realref(a), acb_realref(z), prec);
130 arb_zero(acb_imagref(res));
131 acb_clear(a);
132 return;
133 }
134
135 if (acb_is_zero(z))
136 {
137 acb_zero(res);
138 return;
139 }
140
141 if (acb_is_one(z))
142 {
143 acb_one(res);
144 return;
145 }
146
147 if (!acb_check_arg(z))
148 {
149 acb_agm1_around_zero(res, z, prec);
150 return;
151 }
152
153 acb_init(a);
154 acb_init(b);
155 acb_init(t);
156 acb_init(u);
157 mag_init(err);
158 mag_init(err2);
159
160 acb_one(a);
161 acb_set_round(b, z, prec);
162
163 while (1)
164 {
165 acb_sub(u, a, b, prec);
166
167 if (acb_contains_zero(u))
168 {
169 /* Dupont's thesis, p. 87: |M(z) - a_n| <= |a_n - b_n| */
170 acb_set(res, a);
171 acb_get_mag(err, u);
172 acb_add_error_mag(res, err);
173 break;
174 }
175
176 acb_add(t, a, b, prec);
177
178 acb_get_mag(err, u);
179 acb_get_mag_lower(err2, t);
180 mag_div(err, err, err2);
181 mag_geom_series(err, err, 10);
182 mag_mul_2exp_si(err, err, -6);
183
184 if (mag_cmp_2exp_si(err, -prec) < 0)
185 {
186 acb_agm_close_taylor(res, a, b, t, u, err, prec);
187 break;
188 }
189
190 acb_mul_2exp_si(t, t, -1);
191
192 sqrtmul(u, a, b, prec);
193
194 acb_swap(t, a);
195 acb_swap(u, b);
196 }
197
198 acb_clear(a);
199 acb_clear(b);
200 acb_clear(t);
201 acb_clear(u);
202 mag_clear(err);
203 mag_clear(err2);
204 }
205
206 /*
207 Computes (M(z), M'(z)) using a finite difference.
208 Assumes z exact, |arg(z)| <= 3 pi / 4.
209 */
210 void
acb_agm1_deriv_diff(acb_t Mz,acb_t Mzp,const acb_t z,slong prec)211 acb_agm1_deriv_diff(acb_t Mz, acb_t Mzp, const acb_t z, slong prec)
212 {
213 mag_t err, t, C;
214 fmpz_t rexp, hexp;
215 acb_t u, v;
216 slong wp, qexp;
217 int isreal;
218
219 if (!acb_is_exact(z) || !acb_is_finite(z) ||
220 acb_is_zero(z) || !acb_check_arg(z))
221 {
222 acb_indeterminate(Mz);
223 acb_indeterminate(Mzp);
224 return;
225 }
226
227 isreal = acb_is_real(z) && arb_is_nonnegative(acb_realref(z));
228
229 /*
230 |M^(k)(z) / k!| <= C * D^k where
231 C = max(1, |z| + r),
232 D = 1/r,
233 and 0 < r < |z|
234
235 M(z+h) - M(z-h)
236 |--------------- - M'(z)| <= D^3 h^2 / (1 - D h)
237 2h
238
239 M(z+h) + M(z-h)
240 |--------------- - M(z)| <= D^2 h^2 / (1 - D h)
241 2
242
243 h D < 1.
244 */
245
246 fmpz_init(hexp);
247 fmpz_init(rexp);
248 mag_init(err);
249 mag_init(t);
250 mag_init(C);
251 acb_init(u);
252 acb_init(v);
253
254 /* choose r = 2^rexp such that r < |z| */
255 acb_get_mag_lower(t, z);
256 fmpz_sub_ui(rexp, MAG_EXPREF(t), 2);
257
258 /* Choose h = r/q = 2^hexp = 2^(rexp-qexp)
259 with qexp = floor(prec/2) + 5
260
261 D = 1/r = 2^-rexp
262
263 f(z) error <= C D^2 h^2 / (1-Dh)
264 f'(z) error <= C D^3 h^2 / (1-Dh)
265
266 1/(1-Dh) < 2, hence:
267
268 f(z) error < 2 C D^2 h^2 = C 2^(1-2*qexp)
269 f'(z) error < 2 C D^3 h^2 = C 2^(1-rexp-2*qexp)
270 */
271
272 /* C = max(1, |z| + r) */
273 acb_get_mag(C, z);
274 mag_one(t);
275 mag_mul_2exp_fmpz(t, t, rexp);
276 mag_add(C, C, t);
277 mag_one(t);
278 mag_max(C, C, t);
279
280 qexp = prec / 2 + 5;
281 /*
282 if (fmpz_sgn(rexp) < 0)
283 qexp += fmpz_bits(rexp);
284 */
285
286 /* compute h = 2^hexp */
287 fmpz_sub_ui(hexp, rexp, qexp);
288
289 /* compute finite differences */
290 wp = prec + qexp + 5;
291
292 acb_one(u);
293 acb_mul_2exp_fmpz(u, u, hexp);
294 acb_add(u, z, u, wp);
295 acb_agm1_basecase(u, u, wp);
296
297 acb_one(v);
298 acb_mul_2exp_fmpz(v, v, hexp);
299 acb_sub(v, z, v, wp);
300 acb_agm1_basecase(v, v, wp);
301
302 acb_add(Mz, u, v, prec);
303 acb_sub(Mzp, u, v, prec);
304 acb_mul_2exp_si(Mz, Mz, -1);
305 acb_mul_2exp_si(Mzp, Mzp, -1);
306 fmpz_neg(hexp, hexp);
307 acb_mul_2exp_fmpz(Mzp, Mzp, hexp);
308
309 /* add error */
310 mag_mul_2exp_si(err, C, 1 - 2 * qexp);
311
312 if (isreal)
313 arb_add_error_mag(acb_realref(Mz), err);
314 else
315 acb_add_error_mag(Mz, err);
316
317 fmpz_neg(rexp, rexp);
318 mag_mul_2exp_fmpz(err, err, rexp);
319
320 if (isreal)
321 arb_add_error_mag(acb_realref(Mzp), err);
322 else
323 acb_add_error_mag(Mzp, err);
324
325 fmpz_clear(hexp);
326 fmpz_clear(rexp);
327 mag_clear(err);
328 mag_clear(t);
329 mag_clear(C);
330 acb_clear(u);
331 acb_clear(v);
332 }
333
334 /*
335 For input z + eps
336
337 First derivative bound: max(1, |z|+|eps|+r) / r
338 Second derivative bound: 2 max(1, |z|+|eps|+r) / r^2
339
340 This is assuming that the circle at z with radius |eps| + r
341 does not cross the negative half axis, which we check.
342 */
343
344 void
acb_agm1_deriv_right(acb_t Mz,acb_t Mzp,const acb_t z,slong prec)345 acb_agm1_deriv_right(acb_t Mz, acb_t Mzp, const acb_t z, slong prec)
346 {
347 if (acb_is_exact(z))
348 {
349 acb_agm1_deriv_diff(Mz, Mzp, z, prec);
350 }
351 else
352 {
353 if (!acb_is_finite(z) || !acb_check_arg(z))
354 {
355 acb_indeterminate(Mz);
356 acb_indeterminate(Mzp);
357 }
358 else
359 {
360 acb_t t;
361 mag_t r, eps, err, one;
362 int isreal;
363
364 acb_init(t);
365 mag_init(r);
366 mag_init(err);
367 mag_init(one);
368 mag_init(eps);
369
370 isreal = acb_is_real(z) && arb_is_nonnegative(acb_realref(z));
371
372 mag_hypot(eps, arb_radref(acb_realref(z)), arb_radref(acb_imagref(z)));
373
374 /* choose r avoiding overlap with negative half axis */
375 if (arf_sgn(arb_midref(acb_realref(z))) < 0)
376 arb_get_mag_lower(r, acb_imagref(z));
377 else
378 acb_get_mag_lower(r, z);
379
380 mag_mul_2exp_si(r, r, -1);
381
382 if (mag_is_zero(r))
383 {
384 acb_indeterminate(Mz);
385 acb_indeterminate(Mzp);
386 }
387 else
388 {
389 acb_set(t, z);
390 mag_zero(arb_radref(acb_realref(t)));
391 mag_zero(arb_radref(acb_imagref(t)));
392
393 acb_get_mag(err, z);
394 mag_add(err, err, r);
395 mag_add(err, err, eps);
396 mag_one(one);
397 mag_max(err, err, one);
398 mag_mul(err, err, eps);
399
400 acb_agm1_deriv_diff(Mz, Mzp, t, prec);
401
402 mag_div(err, err, r);
403
404 if (isreal)
405 arb_add_error_mag(acb_realref(Mz), err);
406 else
407 acb_add_error_mag(Mz, err);
408
409 mag_div(err, err, r);
410 mag_mul_2exp_si(err, err, 1);
411
412 if (isreal)
413 arb_add_error_mag(acb_realref(Mzp), err);
414 else
415 acb_add_error_mag(Mzp, err);
416 }
417
418 acb_clear(t);
419 mag_clear(r);
420 mag_clear(err);
421 mag_clear(one);
422 mag_clear(eps);
423 }
424 }
425 }
426
427 void
acb_agm1(acb_t res,const acb_t z,slong prec)428 acb_agm1(acb_t res, const acb_t z, slong prec)
429 {
430 if (acb_is_zero(z))
431 {
432 acb_zero(res);
433 }
434 else if (!acb_is_finite(z))
435 {
436 acb_indeterminate(res);
437 }
438 else if (acb_contains_zero(z))
439 {
440 acb_agm1_around_zero(res, z, prec);
441 }
442 else if (arf_sgn(arb_midref(acb_realref(z))) >= 0)
443 {
444 acb_agm1_basecase(res, z, prec);
445 }
446 else if (acb_equal_si(z, -1))
447 {
448 acb_zero(res);
449 }
450 else
451 {
452 /* use M(1,z) = M((z+1)/2, sqrt(z))
453 = (z+1)/2 * M(1, 2 sqrt(z) / (z+1))
454 = sqrt(z) * M(1, (z+1) / (2 sqrt(z)) */
455 acb_t t;
456
457 acb_init(t);
458 acb_add_ui(t, z, 1, prec);
459 acb_mul_2exp_si(t, t, -1);
460
461 if (acb_contains_zero(t))
462 {
463 mag_t ra, rb;
464
465 mag_init(ra);
466 mag_init(rb);
467
468 acb_get_mag(ra, t);
469 acb_get_mag(rb, z);
470 mag_sqrt(rb, rb);
471
472 mag_agm(ra, ra, rb);
473
474 acb_zero(res);
475 acb_add_error_mag(res, ra);
476
477 mag_clear(ra);
478 mag_clear(rb);
479 }
480 else if (acb_rel_accuracy_bits(t) > acb_rel_accuracy_bits(z))
481 {
482 acb_sqrt(res, z, prec);
483 acb_div(res, res, t, prec);
484 acb_agm1_basecase(res, res, prec);
485 acb_mul(res, res, t, prec);
486 }
487 else
488 {
489 acb_sqrt(res, z, prec);
490 acb_div(t, t, res, prec);
491 acb_agm1_basecase(t, t, prec);
492 acb_mul(res, res, t, prec);
493 }
494
495 acb_clear(t);
496 }
497 }
498
499 void
acb_agm1_deriv(acb_t Mz,acb_t Mzp,const acb_t z,slong prec)500 acb_agm1_deriv(acb_t Mz, acb_t Mzp, const acb_t z, slong prec)
501 {
502 /*
503 u = 2 sqrt(z) / (1+z)
504
505 Mz = (1+z) M(u) / 2
506 Mzp = [M(u) - (z-1) M'(u) / ((1+z) sqrt(z))] / 2
507 */
508
509 if (arf_sgn(arb_midref(acb_realref(z))) >= 0)
510 {
511 if (acb_is_one(z))
512 {
513 acb_one(Mz);
514 acb_mul_2exp_si(Mzp, Mz, -1);
515 }
516 else
517 acb_agm1_deriv_right(Mz, Mzp, z, prec);
518 }
519 else
520 {
521 acb_t t, u, zp1, zm1;
522
523 acb_init(t);
524 acb_init(u);
525 acb_init(zp1);
526 acb_init(zm1);
527
528 acb_sqrt(t, z, prec);
529 acb_add_ui(zp1, z, 1, prec);
530 acb_sub_ui(zm1, z, 1, prec);
531
532 acb_div(u, t, zp1, prec);
533 acb_mul_2exp_si(u, u, 1);
534
535 acb_agm1_deriv_right(Mz, Mzp, u, prec);
536
537 acb_mul(Mzp, Mzp, zm1, prec);
538 acb_mul(t, t, zp1, prec);
539 acb_div(Mzp, Mzp, t, prec);
540 acb_sub(Mzp, Mz, Mzp, prec);
541 acb_mul_2exp_si(Mzp, Mzp, -1);
542
543 acb_mul(Mz, Mz, zp1, prec);
544 acb_mul_2exp_si(Mz, Mz, -1);
545
546 acb_clear(t);
547 acb_clear(u);
548 acb_clear(zp1);
549 acb_clear(zm1);
550 }
551 }
552
553 void
acb_agm1_cpx(acb_ptr m,const acb_t z,slong len,slong prec)554 acb_agm1_cpx(acb_ptr m, const acb_t z, slong len, slong prec)
555 {
556 if (len < 1)
557 return;
558
559 if (len == 1)
560 {
561 acb_agm1(m, z, prec);
562 return;
563 }
564
565 if (len == 2)
566 {
567 acb_agm1_deriv(m, m + 1, z, prec);
568 return;
569 }
570
571 if (len >= 3)
572 {
573 acb_t t, u, v;
574 acb_ptr w;
575 slong k, n;
576
577 acb_init(t);
578 acb_init(u);
579 acb_init(v);
580 w = _acb_vec_init(len);
581
582 acb_agm1_deriv(w, w + 1, z, prec);
583
584 /* invert series */
585 acb_inv(w, w, prec);
586 acb_mul(t, w, w, prec);
587 acb_mul(w + 1, w + 1, t, prec);
588 acb_neg(w + 1, w + 1);
589
590 if (acb_is_one(z))
591 {
592 for (k = 2; k < len; k++)
593 {
594 n = k - 2;
595
596 acb_mul_ui(w + k, w + n + 0, (n+1)*(n+1), prec);
597 acb_addmul_ui(w + k, w + n + 1, 7+3*n*(3+n), prec);
598 acb_div_ui(w + k, w + k, 2*(n+2)*(n+2), prec);
599 acb_neg(w + k, w + k);
600 }
601 }
602 else
603 {
604 /* t = 3z^2 - 1 */
605 /* u = -1 / (z^3 - z) */
606 acb_mul(t, z, z, prec);
607 acb_mul(u, t, z, prec);
608 acb_mul_ui(t, t, 3, prec);
609 acb_sub_ui(t, t, 1, prec);
610 acb_sub(u, u, z, prec);
611 acb_inv(u, u, prec);
612 acb_neg(u, u);
613
614 /* use differential equation for second derivative */
615 acb_mul(w + 2, z, w + 0, prec);
616 acb_addmul(w + 2, t, w + 1, prec);
617 acb_mul(w + 2, w + 2, u, prec);
618 acb_mul_2exp_si(w + 2, w + 2, -1);
619
620 /* recurrence */
621 for (k = 3; k < len; k++)
622 {
623 n = k - 3;
624 acb_mul_ui(w + k, w + n + 0, (n+1)*(n+1), prec);
625 acb_mul(v, w + n + 1, z, prec);
626 acb_addmul_ui(w + k, v, 7+3*n*(3+n), prec);
627 acb_mul(v, w + n + 2, t, prec);
628 acb_addmul_ui(w + k, v, (n+2)*(n+2), prec);
629 acb_mul(w + k, w + k, u, prec);
630 acb_div_ui(w + k, w + k, (n+2)*(n+3), prec);
631 }
632 }
633
634 /* invert series */
635 _acb_poly_inv_series(m, w, len, len, prec);
636
637 acb_clear(t);
638 acb_clear(u);
639 acb_clear(v);
640 _acb_vec_clear(w, len);
641 }
642 }
643
644