1 /*
2 Copyright (C) 2017 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
_arb_arf_div_rounded_den(arb_t res,const arf_t x,const arf_t y,int yinexact,slong prec)15 _arb_arf_div_rounded_den(arb_t res, const arf_t x, const arf_t y, int yinexact, slong prec)
16 {
17 int inexact = arf_div(arb_midref(res), x, y, prec, ARB_RND);
18
19 if (yinexact && !arf_is_special(arb_midref(res)))
20 arf_mag_set_ulp(arb_radref(res), arb_midref(res), prec - 1);
21 else if (inexact)
22 arf_mag_set_ulp(arb_radref(res), arb_midref(res), prec);
23 else
24 mag_zero(arb_radref(res));
25 }
26
27 static void
_arb_arf_div_rounded_den_add_err(arb_t res,const arf_t x,const arf_t y,int yinexact,slong prec)28 _arb_arf_div_rounded_den_add_err(arb_t res, const arf_t x, const arf_t y, int yinexact, slong prec)
29 {
30 int inexact = arf_div(arb_midref(res), x, y, prec, ARB_RND);
31
32 if (yinexact && !arf_is_special(arb_midref(res)))
33 arf_mag_add_ulp(arb_radref(res), arb_radref(res), arb_midref(res), prec - 1);
34 else if (inexact)
35 arf_mag_add_ulp(arb_radref(res), arb_radref(res), arb_midref(res), prec);
36 }
37
38 void
acb_inv(acb_t res,const acb_t z,slong prec)39 acb_inv(acb_t res, const acb_t z, slong prec)
40 {
41 mag_t am, bm;
42 slong hprec;
43
44 #define a arb_midref(acb_realref(z))
45 #define b arb_midref(acb_imagref(z))
46 #define x arb_radref(acb_realref(z))
47 #define y arb_radref(acb_imagref(z))
48
49 /* choose precision for the floating-point approximation of a^2+b^2 so
50 that the double rounding result in less than
51 2 ulp error; also use at least MAG_BITS bits since the
52 value will be recycled for error bounds */
53 hprec = FLINT_MAX(prec + 3, MAG_BITS);
54
55 if (arb_is_zero(acb_imagref(z)))
56 {
57 arb_inv(acb_realref(res), acb_realref(z), prec);
58 arb_zero(acb_imagref(res));
59 return;
60 }
61
62 if (arb_is_zero(acb_realref(z)))
63 {
64 arb_inv(acb_imagref(res), acb_imagref(z), prec);
65 arb_neg(acb_imagref(res), acb_imagref(res));
66 arb_zero(acb_realref(res));
67 return;
68 }
69
70 if (!acb_is_finite(z))
71 {
72 acb_indeterminate(res);
73 return;
74 }
75
76 if (mag_is_zero(x) && mag_is_zero(y))
77 {
78 int inexact;
79
80 arf_t a2b2;
81 arf_init(a2b2);
82
83 inexact = arf_sosq(a2b2, a, b, hprec, ARF_RND_DOWN);
84
85 if (arf_is_special(a2b2))
86 {
87 acb_indeterminate(res);
88 }
89 else
90 {
91 _arb_arf_div_rounded_den(acb_realref(res), a, a2b2, inexact, prec);
92 _arb_arf_div_rounded_den(acb_imagref(res), b, a2b2, inexact, prec);
93 arf_neg(arb_midref(acb_imagref(res)), arb_midref(acb_imagref(res)));
94 }
95
96 arf_clear(a2b2);
97 return;
98 }
99
100 mag_init(am);
101 mag_init(bm);
102
103 /* first bound |a|-x, |b|-y */
104 arb_get_mag_lower(am, acb_realref(z));
105 arb_get_mag_lower(bm, acb_imagref(z));
106
107 if ((mag_is_zero(am) && mag_is_zero(bm)))
108 {
109 acb_indeterminate(res);
110 }
111 else
112 {
113 /*
114 The propagated error in the real part is given exactly by
115
116 (a+x')/((a+x')^2+(b+y'))^2 - a/(a^2+b^2) = P / Q,
117
118 P = [(b^2-a^2) x' - a (x'^2+y'^2 + 2y'b)]
119 Q = [(a^2+b^2)((a+x')^2+(b+y')^2)]
120
121 where |x'| <= x and |y'| <= y, and analogously for the imaginary part.
122 */
123 mag_t t, u, v, w;
124 arf_t a2b2;
125 int inexact;
126
127 mag_init(t);
128 mag_init(u);
129 mag_init(v);
130 mag_init(w);
131
132 arf_init(a2b2);
133
134 inexact = arf_sosq(a2b2, a, b, hprec, ARF_RND_DOWN);
135
136 /* compute denominator */
137 /* t = (|a|-x)^2 + (|b|-x)^2 (lower bound) */
138 mag_mul_lower(t, am, am);
139 mag_mul_lower(u, bm, bm);
140 mag_add_lower(t, t, u);
141 /* u = a^2 + b^2 (lower bound) */
142 arf_get_mag_lower(u, a2b2);
143 /* t = ((|a|-x)^2 + (|b|-x)^2)(a^2 + b^2) (lower bound) */
144 mag_mul_lower(t, t, u);
145
146 /* compute numerator */
147 /* real: |a^2-b^2| x + |a| ((x^2 + y^2) + 2 |b| y)) */
148 /* imag: |a^2-b^2| y + |b| ((x^2 + y^2) + 2 |a| x)) */
149 /* am, bm = upper bounds for a, b */
150 arf_get_mag(am, a);
151 arf_get_mag(bm, b);
152
153 /* v = x^2 + y^2 */
154 mag_mul(v, x, x);
155 mag_addmul(v, y, y);
156
157 /* u = |a| ((x^2 + y^2) + 2 |b| y) */
158 mag_mul_2exp_si(u, bm, 1);
159 mag_mul(u, u, y);
160 mag_add(u, u, v);
161 mag_mul(u, u, am);
162
163 /* v = |b| ((x^2 + y^2) + 2 |a| x) */
164 mag_mul_2exp_si(w, am, 1);
165 mag_addmul(v, w, x);
166 mag_mul(v, v, bm);
167
168 /* w = |b^2 - a^2| (upper bound) */
169 if (arf_cmpabs(a, b) >= 0)
170 mag_mul(w, am, am);
171 else
172 mag_mul(w, bm, bm);
173
174 mag_addmul(u, w, x);
175 mag_addmul(v, w, y);
176
177 mag_div(arb_radref(acb_realref(res)), u, t);
178 mag_div(arb_radref(acb_imagref(res)), v, t);
179
180 _arb_arf_div_rounded_den_add_err(acb_realref(res), a, a2b2, inexact, prec);
181 _arb_arf_div_rounded_den_add_err(acb_imagref(res), b, a2b2, inexact, prec);
182 arf_neg(arb_midref(acb_imagref(res)), arb_midref(acb_imagref(res)));
183
184 mag_clear(t);
185 mag_clear(u);
186 mag_clear(v);
187 mag_clear(w);
188
189 arf_clear(a2b2);
190 }
191
192 mag_clear(am);
193 mag_clear(bm);
194 #undef a
195 #undef b
196 #undef x
197 #undef y
198 }
199
200