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 "arb.h"
13
14 void
arb_sqrt_ui(arb_t z,ulong x,slong prec)15 arb_sqrt_ui(arb_t z, ulong x, slong prec)
16 {
17 arf_t t;
18 arf_init_set_ui(t, x); /* no need to free */
19 arb_sqrt_arf(z, t, prec);
20 }
21
22 void
arb_sqrt_fmpz(arb_t z,const fmpz_t x,slong prec)23 arb_sqrt_fmpz(arb_t z, const fmpz_t x, slong prec)
24 {
25 arf_t t;
26 arf_init(t);
27 arf_set_fmpz(t, x);
28 arb_sqrt_arf(z, t, prec);
29 arf_clear(t);
30 }
31
32 void
arb_sqrt_arf(arb_t z,const arf_t x,slong prec)33 arb_sqrt_arf(arb_t z, const arf_t x, slong prec)
34 {
35 if (arf_sgn(x) < 0 || arf_is_nan(x))
36 {
37 arb_indeterminate(z);
38 }
39 else
40 {
41 int inexact;
42
43 inexact = arf_sqrt(arb_midref(z), x, prec, ARB_RND);
44
45 if (inexact)
46 arf_mag_set_ulp(arb_radref(z), arb_midref(z), prec);
47 else
48 mag_zero(arb_radref(z));
49 }
50 }
51
52 void
arb_sqrt(arb_t z,const arb_t x,slong prec)53 arb_sqrt(arb_t z, const arb_t x, slong prec)
54 {
55 mag_t rx, zr;
56 int inexact;
57
58 if (mag_is_zero(arb_radref(x)))
59 {
60 arb_sqrt_arf(z, arb_midref(x), prec);
61 }
62 else if (arf_is_special(arb_midref(x)) ||
63 arf_sgn(arb_midref(x)) < 0 || mag_is_inf(arb_radref(x)))
64 {
65 if (arf_is_pos_inf(arb_midref(x)) && mag_is_finite(arb_radref(x)))
66 arb_sqrt_arf(z, arb_midref(x), prec);
67 else
68 arb_indeterminate(z);
69 }
70 else /* now both mid and rad are non-special values, mid > 0 */
71 {
72 slong acc;
73
74 acc = _fmpz_sub_small(ARF_EXPREF(arb_midref(x)), MAG_EXPREF(arb_radref(x)));
75 acc = FLINT_MIN(acc, prec);
76 prec = FLINT_MIN(prec, acc + MAG_BITS);
77 prec = FLINT_MAX(prec, 2);
78
79 if (acc < 0)
80 {
81 arb_indeterminate(z);
82 }
83 else if (acc <= 20)
84 {
85 mag_t t, u;
86
87 mag_init(t);
88 mag_init(u);
89
90 arb_get_mag_lower(t, x);
91
92 if (mag_is_zero(t) && arb_contains_negative(x))
93 {
94 arb_indeterminate(z);
95 }
96 else
97 {
98 arb_get_mag(u, x);
99 mag_sqrt_lower(t, t);
100 mag_sqrt(u, u);
101 arb_set_interval_mag(z, t, u, prec);
102 }
103
104 mag_clear(t);
105 mag_clear(u);
106 }
107 else if (ARB_IS_LAGOM(x)) /* small exponents, acc *and* prec >= 20 */
108 {
109 mag_t t;
110 mag_init(t); /* no need to free */
111
112 inexact = arf_sqrt(arb_midref(z), arb_midref(x), prec, ARB_RND);
113
114 /* sqrt(x) - sqrt(x-r) <= 0.5 * r * rsqrt(x-r) */
115 /* we have rsqrt(x-r) ~= 1/sqrt(x) */
116 arf_get_mag_lower(t, arb_midref(z));
117
118 /* note: we need to write rad(z) first to use fast_mul later */
119 mag_div(arb_radref(z), arb_radref(x), t);
120
121 /* We are guaranteed to have acc and prec >= 20. */
122 /* 0.5 + eps corrects for errors */
123 MAG_MAN(t) = MAG_ONE_HALF + (MAG_ONE_HALF >> 16);
124 MAG_EXP(t) = 0;
125 mag_fast_mul(arb_radref(z), arb_radref(z), t);
126
127 if (inexact)
128 arf_mag_fast_add_ulp(arb_radref(z), arb_radref(z), arb_midref(z), prec);
129 }
130 else
131 {
132 mag_init(zr);
133 mag_init(rx);
134
135 /* rx = upper bound for r / x */
136 arf_get_mag_lower(rx, arb_midref(x));
137 mag_div(rx, arb_radref(x), rx);
138
139 inexact = arf_sqrt(arb_midref(z), arb_midref(x), prec, ARB_RND);
140
141 /* zr = upper bound for sqrt(x) */
142 arf_get_mag(zr, arb_midref(z));
143 if (inexact)
144 arf_mag_add_ulp(zr, zr, arb_midref(z), prec);
145
146 /* propagated error: sqrt(x) - sqrt(x-r)
147 = sqrt(x) * [1 - sqrt(1 - r/x)]
148 <= sqrt(x) * 0.5 * (rx + rx^2) */
149 mag_addmul(rx, rx, rx);
150 mag_mul(zr, zr, rx);
151 mag_mul_2exp_si(zr, zr, -1);
152
153 /* add the rounding error */
154 if (inexact)
155 arf_mag_add_ulp(arb_radref(z), zr, arb_midref(z), prec);
156 else
157 mag_swap(arb_radref(z), zr);
158
159 mag_clear(zr);
160 mag_clear(rx);
161 }
162 }
163 }
164
165