1/* mpfr_cbrt -- cube root function. 2 3Copyright 2002-2023 Free Software Foundation, Inc. 4Contributed by the AriC and Caramba projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 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 48int 49mpfr_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[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode), 59 ("y[%Pu]=%.*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} 178