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 "arf.h"
13
14 int
arf_root(arf_ptr z,arf_srcptr x,ulong k,slong prec,arf_rnd_t rnd)15 arf_root(arf_ptr z, arf_srcptr x, ulong k, slong prec, arf_rnd_t rnd)
16 {
17 mp_size_t xn, zn, val;
18 mp_srcptr xptr;
19 mp_ptr tmp, zptr;
20 mpfr_t xf, zf;
21 fmpz_t q, r;
22 int inexact;
23
24 if (k == 0)
25 {
26 arf_nan(z);
27 return 0;
28 }
29
30 if (k == 1)
31 return arf_set_round(z, x, prec, rnd);
32
33 if (k == 2)
34 return arf_sqrt(z, x, prec, rnd);
35
36 if (arf_is_special(x))
37 {
38 if (arf_is_neg_inf(x))
39 arf_nan(z);
40 else
41 arf_set(z, x);
42 return 0;
43 }
44
45 if (ARF_SGNBIT(x))
46 {
47 arf_nan(z);
48 return 0;
49 }
50
51 fmpz_init(q);
52 fmpz_init(r);
53
54 /* x = m * 2^e where e = qk + r */
55 /* x^(1/k) = (m * 2^(qk+r))^(1/k) */
56 /* x^(1/k) = (m * 2^r)^(1/k) * 2^q */
57 fmpz_set_ui(r, k);
58 fmpz_fdiv_qr(q, r, ARF_EXPREF(x), r);
59
60 ARF_GET_MPN_READONLY(xptr, xn, x);
61 zn = (prec + FLINT_BITS - 1) / FLINT_BITS;
62
63 zf->_mpfr_d = tmp = flint_malloc(zn * sizeof(mp_limb_t));
64 zf->_mpfr_prec = prec;
65 zf->_mpfr_sign = 1;
66 zf->_mpfr_exp = 0;
67
68 xf->_mpfr_d = (mp_ptr) xptr;
69 xf->_mpfr_prec = xn * FLINT_BITS;
70 xf->_mpfr_sign = 1;
71 xf->_mpfr_exp = fmpz_get_ui(r);
72
73 #if MPFR_VERSION_MAJOR >= 4
74 inexact = mpfr_rootn_ui(zf, xf, k, arf_rnd_to_mpfr(rnd));
75 #else
76 inexact = mpfr_root(zf, xf, k, arf_rnd_to_mpfr(rnd));
77 #endif
78 inexact = (inexact != 0);
79
80 val = 0;
81 while (tmp[val] == 0)
82 val++;
83
84 ARF_GET_MPN_WRITE(zptr, zn - val, z);
85 flint_mpn_copyi(zptr, tmp + val, zn - val);
86
87 fmpz_add_si(ARF_EXPREF(z), q, zf->_mpfr_exp);
88
89 flint_free(tmp);
90 fmpz_clear(q);
91 fmpz_clear(r);
92
93 return inexact;
94 }
95
96