1 /*
2 Copyright (C) 2014, 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 "arb.h"
13
14 void
arb_div_arf(arb_t z,const arb_t x,const arf_t y,slong prec)15 arb_div_arf(arb_t z, const arb_t x, const arf_t y, slong prec)
16 {
17 mag_t zr, ym;
18 int inexact;
19
20 if (arf_is_special(y) || !arb_is_finite(x))
21 {
22 if (arf_is_inf(arb_midref(x)) && mag_is_finite(arb_radref(x)) &&
23 arf_is_finite(y) && !arf_is_zero(y))
24 {
25 /* inf / finite nonzero = inf */
26 arf_div(arb_midref(z), arb_midref(x), y, prec, ARF_RND_DOWN);
27 mag_zero(arb_radref(z));
28 }
29 else if (arb_is_finite(x) && arf_is_inf(y))
30 {
31 /* finite / inf = 0 */
32 arb_zero(z);
33 }
34 else if (!arf_is_nan(arb_midref(x)) && mag_is_inf(arb_radref(x)) &&
35 arf_is_finite(y) && !arf_is_zero(y))
36 {
37 /* [+/- inf] / finite nonzero = [+/- inf] */
38 arb_zero_pm_inf(z);
39 }
40 else
41 {
42 arb_indeterminate(z);
43 }
44 }
45 else if (mag_is_zero(arb_radref(x)))
46 {
47 inexact = arf_div(arb_midref(z), arb_midref(x), y, prec, ARB_RND);
48
49 if (inexact)
50 arf_mag_set_ulp(arb_radref(z), arb_midref(z), prec);
51 else
52 mag_zero(arb_radref(z));
53 }
54 else
55 {
56 mag_init(ym);
57 mag_init(zr);
58
59 arf_get_mag_lower(ym, y);
60 mag_div(zr, arb_radref(x), ym);
61
62 inexact = arf_div(arb_midref(z), arb_midref(x), y, prec, ARB_RND);
63
64 if (inexact)
65 arf_mag_add_ulp(arb_radref(z), zr, arb_midref(z), prec);
66 else
67 mag_swap(arb_radref(z), zr);
68
69 mag_clear(ym);
70 mag_clear(zr);
71 }
72 }
73
74 static void
arb_div_wide(arb_t z,const arb_t x,const arb_t y,slong prec)75 arb_div_wide(arb_t z, const arb_t x, const arb_t y, slong prec)
76 {
77 mag_t a, b, t, u;
78
79 mag_init(t);
80 arb_get_mag_lower(t, y);
81
82 if (mag_is_zero(t))
83 {
84 arb_indeterminate(z);
85 }
86 else if (arf_is_zero(arb_midref(x)))
87 {
88 mag_div(arb_radref(z), arb_radref(x), t);
89 arf_zero(arb_midref(z));
90 }
91 else
92 {
93 if (arf_cmpabs_mag(arb_midref(x), arb_radref(x)) >= 0)
94 {
95 /* [a,b] / [t,u] = [a/u, b/t]
96 [a,b] / -[t,u] = -[a/u, b/t]
97 -[a,b] / [t,u] = -[a/u, b/t]
98 -[a,b] / -[t,u] = [a/u, b/t] */
99 mag_init(a);
100 mag_init(b);
101 mag_init(u);
102
103 arb_get_mag_lower(a, x);
104 arb_get_mag(b, x);
105 arb_get_mag(u, y);
106
107 mag_div_lower(a, a, u);
108 mag_div(b, b, t);
109
110 if ((arf_sgn(arb_midref(x)) < 0) ^ (arf_sgn(arb_midref(y)) < 0))
111 {
112 arb_set_interval_mag(z, a, b, prec);
113 arb_neg(z, z);
114 }
115 else
116 {
117 arb_set_interval_mag(z, a, b, prec);
118 }
119
120 mag_clear(a);
121 mag_clear(b);
122 mag_clear(u);
123 }
124 else
125 {
126 /* [-a,b] / [t,u] = [-a/t, b/t]
127 [-a,b] / -[t,u] = [-b/t, a/t] */
128 mag_init(a);
129 mag_init(b);
130
131 arb_get_mag(b, x);
132 arf_get_mag_lower(a, arb_midref(x));
133 mag_sub(a, arb_radref(x), a);
134
135 if ((arf_sgn(arb_midref(x)) < 0) ^ (arf_sgn(arb_midref(y)) < 0))
136 mag_swap(a, b);
137
138 mag_div(a, a, t);
139 mag_div(b, b, t);
140 arb_set_interval_neg_pos_mag(z, a, b, prec);
141
142 mag_clear(a);
143 mag_clear(b);
144 }
145 }
146
147 mag_clear(t);
148 }
149
150 void
arb_div(arb_t z,const arb_t x,const arb_t y,slong prec)151 arb_div(arb_t z, const arb_t x, const arb_t y, slong prec)
152 {
153 int inexact;
154 slong acc, xacc, yacc;
155
156 if (mag_is_zero(arb_radref(y)))
157 {
158 arb_div_arf(z, x, arb_midref(y), prec);
159 }
160 else if (arf_is_zero(arb_midref(y))) /* anything / 0 = nan */
161 {
162 arb_indeterminate(z);
163 }
164 else if (ARB_IS_LAGOM(x) && ARB_IS_LAGOM(y)) /* fast case */
165 {
166 /* both finite, both midpoint and radius of y not special */
167 yacc = MAG_EXP(arb_midref(y)) - ARF_EXP(arb_radref(y));
168 xacc = arb_rel_accuracy_bits(x);
169 acc = FLINT_MIN(xacc, yacc);
170 acc = FLINT_MAX(acc, 0);
171 acc = FLINT_MIN(acc, prec);
172 prec = FLINT_MIN(prec, acc + MAG_BITS);
173 prec = FLINT_MAX(prec, 2);
174
175 if (acc <= 20)
176 {
177 arb_div_wide(z, x, y, prec);
178 }
179 else
180 {
181 mag_t t, u, v;
182 mag_init(t); /* no need to free */
183 mag_init(u);
184 mag_init(v);
185
186 /* (x*yrad + y*xrad)/(y*(y-yrad))
187 <= (1+eps) (x*yrad + y*xrad)/y^2
188 <= (1+eps) (x*yrad/y^2 + y*xrad/y^2)
189 <= (1+eps) ((x/y)*yrad/y + xrad/y)
190 <= (1+eps) ((x/y)*yrad + xrad)/y
191 */
192
193 /* t = y */
194 arf_get_mag_lower(t, arb_midref(y));
195
196 inexact = arf_div(arb_midref(z), arb_midref(x), arb_midref(y), prec, ARB_RND);
197
198 /* u = x/y */
199 arf_get_mag(u, arb_midref(z));
200
201 /* v = (x/y)*yrad + xrad */
202 MAG_MAN(v) = MAG_MAN(arb_radref(x));
203 MAG_EXP(v) = MAG_EXP(arb_radref(x));
204 mag_fast_addmul(v, arb_radref(y), u);
205 /* v = v / y */
206 mag_div(arb_radref(z), v, t);
207
208 /* correct for errors */
209 MAG_MAN(t) = MAG_ONE_HALF + (MAG_ONE_HALF >> 16);
210 MAG_EXP(t) = 1;
211 mag_fast_mul(arb_radref(z), arb_radref(z), t);
212
213 if (inexact)
214 arf_mag_fast_add_ulp(arb_radref(z), arb_radref(z), arb_midref(z), prec);
215 }
216 }
217 else if (!arb_is_finite(y))
218 {
219 /* finite / inf = 0 */
220 if (arf_is_inf(arb_midref(y)) && mag_is_finite(arb_radref(y)) &&
221 arb_is_finite(x))
222 arb_zero(z);
223 else
224 arb_indeterminate(z);
225 }
226 else if (!arb_is_finite(x))
227 {
228 if (arb_contains_zero(y) || arf_is_nan(arb_midref(x)))
229 {
230 arb_indeterminate(z);
231 }
232 else
233 {
234 /* +/- inf / finite nonzero = +/- inf */
235 if (arf_is_inf(arb_midref(x)) && mag_is_finite(arb_radref(x)))
236 {
237 arf_div(arb_midref(z), arb_midref(x), arb_midref(y), prec, ARF_RND_DOWN);
238 mag_zero(arb_radref(z));
239 }
240 else if (!arf_is_nan(arb_midref(x)) && mag_is_inf(arb_radref(x)))
241 {
242 arb_zero_pm_inf(z);
243 }
244 else
245 {
246 arb_indeterminate(z);
247 }
248 }
249 }
250 else
251 {
252 yacc = arb_rel_accuracy_bits(y);
253 xacc = arb_rel_accuracy_bits(x);
254 acc = FLINT_MIN(xacc, yacc);
255 acc = FLINT_MAX(acc, 0);
256 acc = FLINT_MIN(acc, prec);
257 prec = FLINT_MIN(prec, acc + MAG_BITS);
258 prec = FLINT_MAX(prec, 2);
259
260 if (acc <= 20)
261 {
262 arb_div_wide(z, x, y, prec);
263 }
264 else
265 {
266 mag_t zr, xm, ym, yl, yw;
267
268 mag_init_set_arf(xm, arb_midref(x));
269 mag_init_set_arf(ym, arb_midref(y));
270 mag_init(zr);
271 mag_init(yl);
272 mag_init(yw);
273
274 /* (|x|*yrad + |y|*xrad)/(|y|*(|y|-yrad)) */
275 mag_mul(zr, xm, arb_radref(y));
276 mag_addmul(zr, ym, arb_radref(x));
277 arb_get_mag_lower(yw, y);
278
279 arf_get_mag_lower(yl, arb_midref(y));
280 mag_mul_lower(yl, yl, yw);
281
282 mag_div(zr, zr, yl);
283
284 inexact = arf_div(arb_midref(z), arb_midref(x), arb_midref(y), prec, ARB_RND);
285
286 if (inexact)
287 arf_mag_add_ulp(arb_radref(z), zr, arb_midref(z), prec);
288 else
289 mag_swap(arb_radref(z), zr);
290
291 mag_clear(xm);
292 mag_clear(ym);
293 mag_clear(zr);
294 mag_clear(yl);
295 mag_clear(yw);
296 }
297 }
298 }
299
300 void
arb_div_si(arb_t z,const arb_t x,slong y,slong prec)301 arb_div_si(arb_t z, const arb_t x, slong y, slong prec)
302 {
303 arf_t t;
304 arf_init_set_si(t, y); /* no need to free */
305 arb_div_arf(z, x, t, prec);
306 }
307
308 void
arb_div_ui(arb_t z,const arb_t x,ulong y,slong prec)309 arb_div_ui(arb_t z, const arb_t x, ulong y, slong prec)
310 {
311 arf_t t;
312 arf_init_set_ui(t, y); /* no need to free */
313 arb_div_arf(z, x, t, prec);
314 }
315
316 void
arb_div_fmpz(arb_t z,const arb_t x,const fmpz_t y,slong prec)317 arb_div_fmpz(arb_t z, const arb_t x, const fmpz_t y, slong prec)
318 {
319 arf_t t;
320
321 if (!COEFF_IS_MPZ(*y))
322 {
323 arf_init_set_si(t, *y); /* no need to free */
324 arb_div_arf(z, x, t, prec);
325 }
326 else
327 {
328 arf_init(t);
329 arf_set_fmpz(t, y);
330 arb_div_arf(z, x, t, prec);
331 arf_clear(t);
332 }
333 }
334
335
336 void
arb_fmpz_div_fmpz(arb_t z,const fmpz_t x,const fmpz_t y,slong prec)337 arb_fmpz_div_fmpz(arb_t z, const fmpz_t x, const fmpz_t y, slong prec)
338 {
339 int inexact;
340
341 inexact = arf_fmpz_div_fmpz(arb_midref(z), x, y, prec, ARB_RND);
342
343 if (inexact)
344 arf_mag_set_ulp(arb_radref(z), arb_midref(z), prec);
345 else
346 mag_zero(arb_radref(z));
347 }
348
349 void
arb_ui_div(arb_t z,ulong x,const arb_t y,slong prec)350 arb_ui_div(arb_t z, ulong x, const arb_t y, slong prec)
351 {
352 arb_t t;
353 arb_init(t);
354 arb_set_ui(t, x);
355 arb_div(z, t, y, prec);
356 arb_clear(t);
357 }
358
359