1 /*
2 Copyright (C) 2015 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 static void
arb_sqrt1pm1_tiny(arb_t r,const arb_t z,slong prec)15 arb_sqrt1pm1_tiny(arb_t r, const arb_t z, slong prec)
16 {
17 mag_t b, c;
18 arb_t t;
19
20 mag_init(b);
21 mag_init(c);
22 arb_init(t);
23
24 /* if |z| < 1, then |(sqrt(1+z)-1) - (z/2-z^2/8)| <= |z|^3/(1-|z|)/16 */
25 arb_get_mag(b, z);
26 mag_one(c);
27 mag_sub_lower(c, c, b);
28 mag_pow_ui(b, b, 3);
29 mag_div(b, b, c);
30 mag_mul_2exp_si(b, b, -4);
31
32 arb_mul(t, z, z, prec);
33 arb_mul_2exp_si(t, t, -2);
34 arb_sub(r, z, t, prec);
35 arb_mul_2exp_si(r, r, -1);
36
37 if (mag_is_finite(b))
38 arb_add_error_mag(r, b);
39 else
40 arb_indeterminate(r);
41
42 mag_clear(b);
43 mag_clear(c);
44 arb_clear(t);
45 }
46
47 void
arb_sqrt1pm1(arb_t r,const arb_t z,slong prec)48 arb_sqrt1pm1(arb_t r, const arb_t z, slong prec)
49 {
50 slong magz, wp;
51
52 if (arb_is_zero(z))
53 {
54 arb_zero(r);
55 return;
56 }
57
58 magz = arf_abs_bound_lt_2exp_si(arb_midref(z));
59
60 if (magz < -prec)
61 {
62 arb_sqrt1pm1_tiny(r, z, prec);
63 }
64 else
65 {
66 if (magz < 0)
67 wp = prec + (-magz) + 4;
68 else
69 wp = prec + 4;
70
71 arb_add_ui(r, z, 1, wp);
72 arb_sqrt(r, r, wp);
73 arb_sub_ui(r, r, 1, wp);
74 }
75 }
76
77