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_root_arf(arb_t z,const arf_t x,ulong k,slong prec)15 arb_root_arf(arb_t z, const arf_t x, ulong k, slong prec)
16 {
17 int inexact = arf_root(arb_midref(z), x, k, prec, ARB_RND);
18
19 if (inexact)
20 arf_mag_set_ulp(arb_radref(z), arb_midref(z), prec);
21 else
22 mag_zero(arb_radref(z));
23 }
24
25 void
arb_root_ui_algebraic(arb_t res,const arb_t x,ulong k,slong prec)26 arb_root_ui_algebraic(arb_t res, const arb_t x, ulong k, slong prec)
27 {
28 mag_t r, msubr, m1k, t;
29
30 if (arb_is_exact(x))
31 {
32 arb_root_arf(res, arb_midref(x), k, prec);
33 return;
34 }
35
36 if (!arb_is_nonnegative(x))
37 {
38 arb_indeterminate(res);
39 return;
40 }
41
42 mag_init(r);
43 mag_init(msubr);
44 mag_init(m1k);
45 mag_init(t);
46
47 /* x = [m-r, m+r] */
48 mag_set(r, arb_radref(x));
49 /* m - r */
50 arb_get_mag_lower(msubr, x);
51
52 /* m^(1/k) */
53 arb_root_arf(res, arb_midref(x), k, prec);
54
55 /* bound for m^(1/k) */
56 arb_get_mag(m1k, res);
57
58 /* C = min(1, log(1+r/(m-r))/k) */
59 mag_div(t, r, msubr);
60 mag_log1p(t, t);
61 mag_div_ui(t, t, k);
62 if (mag_cmp_2exp_si(t, 0) > 0)
63 mag_one(t);
64
65 /* C m^(1/k) */
66 mag_mul(t, m1k, t);
67 mag_add(arb_radref(res), arb_radref(res), t);
68
69 mag_clear(r);
70 mag_clear(msubr);
71 mag_clear(m1k);
72 mag_clear(t);
73 }
74
75 void
arb_root_ui_exp(arb_t res,const arb_t x,ulong k,slong prec)76 arb_root_ui_exp(arb_t res, const arb_t x, ulong k, slong prec)
77 {
78 arb_log(res, x, prec + 4);
79 arb_div_ui(res, res, k, prec + 4);
80 arb_exp(res, res, prec);
81 }
82
83 void
arb_root_ui(arb_t res,const arb_t x,ulong k,slong prec)84 arb_root_ui(arb_t res, const arb_t x, ulong k, slong prec)
85 {
86 if (k == 0)
87 {
88 arb_indeterminate(res);
89 }
90 else if (k == 1)
91 {
92 arb_set_round(res, x, prec);
93 }
94 else if (k == 2)
95 {
96 arb_sqrt(res, x, prec);
97 }
98 else if (k == 4)
99 {
100 arb_sqrt(res, x, prec + 2);
101 arb_sqrt(res, res, prec);
102 }
103 else
104 {
105 if (k > 50 || prec < (WORD(1) << ((k / 8) + 8)))
106 arb_root_ui_exp(res, x, k, prec);
107 else
108 arb_root_ui_algebraic(res, x, k, prec);
109 }
110 }
111
112 /* backwards compatible alias */
113 void
arb_root(arb_t res,const arb_t x,ulong k,slong prec)114 arb_root(arb_t res, const arb_t x, ulong k, slong prec)
115 {
116 arb_root_ui(res, x, k, prec);
117 }
118
119