1 /*
2 Copyright (C) 2020 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 mag_agm(mag_t res, const mag_t x, const mag_t y);
15
16 static void
agm_helper(acb_t res,const acb_t a,const acb_t b,slong prec)17 agm_helper(acb_t res, const acb_t a, const acb_t b, slong prec)
18 {
19 if (acb_rel_accuracy_bits(b) >= acb_rel_accuracy_bits(a))
20 {
21 acb_div(res, a, b, prec);
22 acb_agm1(res, res, prec);
23 acb_mul(res, res, b, prec);
24 }
25 else
26 {
27 acb_div(res, b, a, prec);
28 acb_agm1(res, res, prec);
29 acb_mul(res, res, a, prec);
30 }
31 }
32
33 void
acb_agm(acb_t res,const acb_t a,const acb_t b,slong prec)34 acb_agm(acb_t res, const acb_t a, const acb_t b, slong prec)
35 {
36 acb_t t, u, v;
37
38 if (!acb_is_finite(a) || !acb_is_finite(b))
39 {
40 acb_indeterminate(res);
41 return;
42 }
43
44 if (acb_is_zero(a) || acb_is_zero(b))
45 {
46 acb_zero(res);
47 return;
48 }
49
50 if (arb_is_zero(acb_imagref(a)) && arb_is_zero(acb_imagref(b)))
51 {
52 if (arb_is_nonnegative(acb_realref(a)) && arb_is_nonnegative(acb_realref(b)))
53 {
54 arb_agm(acb_realref(res), acb_realref(a), acb_realref(b), prec);
55 arb_zero(acb_imagref(res));
56 return;
57 }
58 }
59
60 if (acb_contains_zero(a) || acb_contains_zero(b))
61 {
62 mag_t ra, rb;
63
64 mag_init(ra);
65 mag_init(rb);
66
67 acb_get_mag(ra, a);
68 acb_get_mag(rb, b);
69 mag_agm(ra, ra, rb);
70 acb_zero(res);
71 acb_add_error_mag(res, ra);
72
73 mag_clear(ra);
74 mag_clear(rb);
75
76 return;
77 }
78
79 acb_init(t);
80
81 acb_add(t, a, b, prec);
82 acb_mul_2exp_si(t, t, -1);
83
84 /* a ~= -b; bound magnitude */
85 if (acb_contains_zero(t))
86 {
87 mag_t ra, rb;
88
89 mag_init(ra);
90 mag_init(rb);
91
92 acb_get_mag(ra, a);
93 acb_get_mag(rb, b);
94 mag_mul(rb, ra, rb);
95 mag_sqrt(rb, rb);
96
97 acb_get_mag(ra, t);
98 mag_agm(ra, ra, rb);
99
100 acb_zero(res);
101 acb_add_error_mag(res, ra);
102
103 mag_clear(ra);
104 mag_clear(rb);
105
106 acb_clear(t);
107 return;
108 }
109
110 /* Do the initial step with the optimal square root, reducing to agm1 */
111
112 acb_init(u);
113 acb_init(v);
114
115 acb_mul(u, a, b, prec);
116
117 /* we can compute either square root here; avoid the branch cut */
118 if (arf_sgn(arb_midref(acb_realref(u))) >= 0)
119 {
120 acb_sqrt(u, u, prec);
121 }
122 else
123 {
124 acb_neg(u, u);
125 acb_sqrt(u, u, prec);
126 acb_mul_onei(u, u);
127 }
128
129 acb_div(v, t, u, prec);
130
131 if (arb_is_nonnegative(acb_realref(v)))
132 {
133 agm_helper(res, t, u, prec);
134 }
135 else if (arb_is_negative(acb_realref(v)))
136 {
137 acb_neg(u, u);
138 agm_helper(res, t, u, prec);
139 }
140 else
141 {
142 agm_helper(v, t, u, prec);
143 acb_neg(u, u);
144 agm_helper(res, t, u, prec);
145 acb_union(res, res, v, prec);
146 }
147
148 acb_clear(t);
149 acb_clear(u);
150 acb_clear(v);
151 }
152
153