1 /* mpfr_cbrt -- cube root function.
2
3 Copyright 2002-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba projects, INRIA.
5
6 This file is part of the GNU MPFR Library.
7
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
17
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
25
26 /* The computation of y = x^(1/3) is done as follows.
27
28 Let n = PREC(y), or PREC(y) + 1 if the rounding mode is MPFR_RNDN.
29 We seek to compute an integer cube root in precision n and the
30 associated inexact bit (non-zero iff the remainder is non-zero).
31
32 Let us write x, possibly truncated, under the form sign * m * 2^(3*e)
33 where m is an integer such that 2^(3n-3) <= m < 2^(3n), i.e. m has
34 between 3n-2 and 3n bits.
35
36 Let s be the integer cube root of m, i.e. the maximum integer such that
37 m = s^3 + t with t >= 0. Thus 2^(n-1) <= s < 2^n, i.e. s has n bits.
38
39 Then |x|^(1/3) = s * 2^e or (s+1) * 2^e depending on the rounding mode,
40 the sign, and whether s is "inexact" (i.e. t > 0 or the truncation of x
41 was not equal to x).
42
43 Note: The truncation of x was allowed because any breakpoint has n bits
44 and its cube has at most 3n bits. Thus the truncation of x cannot yield
45 a cube root below RNDZ(x^(1/3)) in precision n. [TODO: add details.]
46 */
47
48 int
49 mpfr_cbrt (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
50 {
51 mpz_t m;
52 mpfr_exp_t e, d, sh;
53 mpfr_prec_t n, size_m;
54 int inexact, inexact2, negative, r;
55 MPFR_SAVE_EXPO_DECL (expo);
56
57 MPFR_LOG_FUNC (
58 ("x[%Pd]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode),
59 ("y[%Pd]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y,
60 inexact));
61
62 /* special values */
63 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
64 {
65 if (MPFR_IS_NAN (x))
66 {
67 MPFR_SET_NAN (y);
68 MPFR_RET_NAN;
69 }
70 else if (MPFR_IS_INF (x))
71 {
72 MPFR_SET_INF (y);
73 MPFR_SET_SAME_SIGN (y, x);
74 MPFR_RET (0);
75 }
76 /* case 0: cbrt(+/- 0) = +/- 0 */
77 else /* x is necessarily 0 */
78 {
79 MPFR_ASSERTD (MPFR_IS_ZERO (x));
80 MPFR_SET_ZERO (y);
81 MPFR_SET_SAME_SIGN (y, x);
82 MPFR_RET (0);
83 }
84 }
85
86 /* General case */
87 MPFR_SAVE_EXPO_MARK (expo);
88 mpz_init (m);
89
90 e = mpfr_get_z_2exp (m, x); /* x = m * 2^e */
91 if ((negative = MPFR_IS_NEG(x)))
92 mpz_neg (m, m);
93 r = e % 3;
94 if (r < 0)
95 r += 3;
96 MPFR_ASSERTD (r >= 0 && r < 3 && (e - r) % 3 == 0);
97
98 /* x = (m*2^r) * 2^(e-r) = (m*2^r) * 2^(3*q) */
99
100 MPFR_LOG_MSG (("e=%" MPFR_EXP_FSPEC "d r=%d\n", (mpfr_eexp_t) e, r));
101
102 MPFR_MPZ_SIZEINBASE2 (size_m, m);
103 n = MPFR_PREC (y) + (rnd_mode == MPFR_RNDN);
104
105 /* We will need to multiply m by 2^(r'), truncated if r' < 0, and
106 subtract r' from e, so that m has between 3n-2 and 3n bits and
107 e becomes a multiple of 3.
108 Since r = e % 3, we write r' = 3 * sh + r.
109 We want 3 * n - 2 <= size_m + 3 * sh + r <= 3 * n.
110 Let d = 3 * n - size_m - r. Thus we want 0 <= d - 3 * sh <= 2,
111 i.e. sh = floor(d/3). */
112 d = 3 * (mpfr_exp_t) n - (mpfr_exp_t) size_m - r;
113 sh = d >= 0 ? d / 3 : - ((2 - d) / 3); /* floor(d/3) */
114 r += 3 * sh; /* denoted r' above */
115
116 e -= r;
117 MPFR_ASSERTD (e % 3 == 0);
118 e /= 3;
119
120 inexact = 0;
121
122 if (r > 0)
123 {
124 mpz_mul_2exp (m, m, r);
125 }
126 else if (r < 0)
127 {
128 r = -r;
129 inexact = mpz_scan1 (m, 0) < r;
130 mpz_fdiv_q_2exp (m, m, r);
131 }
132
133 /* we reuse the variable m to store the cube root, since it is not needed
134 any more: we just need to know if the root is exact */
135 inexact = ! mpz_root (m, m, 3) || inexact;
136
137 #if MPFR_WANT_ASSERT > 0
138 {
139 mpfr_prec_t tmp;
140
141 MPFR_MPZ_SIZEINBASE2 (tmp, m);
142 MPFR_ASSERTN (tmp == n);
143 }
144 #endif
145
146 if (inexact)
147 {
148 if (negative)
149 rnd_mode = MPFR_INVERT_RND (rnd_mode);
150 if (rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDA
151 || (rnd_mode == MPFR_RNDN && mpz_tstbit (m, 0)))
152 {
153 inexact = 1;
154 mpz_add_ui (m, m, 1);
155 }
156 else
157 inexact = -1;
158 }
159
160 /* either inexact is not zero, and the conversion is exact, i.e. inexact
161 is not changed; or inexact=0, and inexact is set only when
162 rnd_mode=MPFR_RNDN and bit (n+1) from m is 1 */
163 inexact2 = mpfr_set_z (y, m, MPFR_RNDN);
164 MPFR_ASSERTD (inexact == 0 || inexact2 == 0);
165 inexact += inexact2;
166 MPFR_SET_EXP (y, MPFR_GET_EXP (y) + e);
167
168 if (negative)
169 {
170 MPFR_CHANGE_SIGN (y);
171 inexact = -inexact;
172 }
173
174 mpz_clear (m);
175 MPFR_SAVE_EXPO_FREE (expo);
176 return mpfr_check_range (y, inexact, rnd_mode);
177 }