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 } 178
Indexes created Fri Jun 05 00:26:10 UTC 2026