1/* mpfr_cbrt -- cube root function. 2 3Copyright 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4Contributed by the AriC and Caramel 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 20http://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 x = sign * m * 2^(3*e) where m is an integer 29 30 with 2^(3n-3) <= m < 2^(3n) where n = PREC(y) 31 32 and m = s^3 + r where 0 <= r and m < (s+1)^3 33 34 we want that s has n bits i.e. s >= 2^(n-1), or m >= 2^(3n-3) 35 i.e. m must have at least 3n-2 bits 36 37 then x^(1/3) = s * 2^e if r=0 38 x^(1/3) = (s+1) * 2^e if round up 39 x^(1/3) = (s-1) * 2^e if round down 40 x^(1/3) = s * 2^e if nearest and r < 3/2*s^2+3/4*s+1/8 41 (s+1) * 2^e otherwise 42 */ 43 44int 45mpfr_cbrt (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 46{ 47 mpz_t m; 48 mpfr_exp_t e, r, sh; 49 mpfr_prec_t n, size_m, tmp; 50 int inexact, negative; 51 MPFR_SAVE_EXPO_DECL (expo); 52 53 MPFR_LOG_FUNC ( 54 ("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode), 55 ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y, 56 inexact)); 57 58 /* special values */ 59 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 60 { 61 if (MPFR_IS_NAN (x)) 62 { 63 MPFR_SET_NAN (y); 64 MPFR_RET_NAN; 65 } 66 else if (MPFR_IS_INF (x)) 67 { 68 MPFR_SET_INF (y); 69 MPFR_SET_SAME_SIGN (y, x); 70 MPFR_RET (0); 71 } 72 /* case 0: cbrt(+/- 0) = +/- 0 */ 73 else /* x is necessarily 0 */ 74 { 75 MPFR_ASSERTD (MPFR_IS_ZERO (x)); 76 MPFR_SET_ZERO (y); 77 MPFR_SET_SAME_SIGN (y, x); 78 MPFR_RET (0); 79 } 80 } 81 82 /* General case */ 83 MPFR_SAVE_EXPO_MARK (expo); 84 mpz_init (m); 85 86 e = mpfr_get_z_2exp (m, x); /* x = m * 2^e */ 87 if ((negative = MPFR_IS_NEG(x))) 88 mpz_neg (m, m); 89 r = e % 3; 90 if (r < 0) 91 r += 3; 92 /* x = (m*2^r) * 2^(e-r) = (m*2^r) * 2^(3*q) */ 93 94 MPFR_MPZ_SIZEINBASE2 (size_m, m); 95 n = MPFR_PREC (y) + (rnd_mode == MPFR_RNDN); 96 97 /* we want 3*n-2 <= size_m + 3*sh + r <= 3*n 98 i.e. 3*sh + size_m + r <= 3*n */ 99 sh = (3 * (mpfr_exp_t) n - (mpfr_exp_t) size_m - r) / 3; 100 sh = 3 * sh + r; 101 if (sh >= 0) 102 { 103 mpz_mul_2exp (m, m, sh); 104 e = e - sh; 105 } 106 else if (r > 0) 107 { 108 mpz_mul_2exp (m, m, r); 109 e = e - r; 110 } 111 112 /* invariant: x = m*2^e, with e divisible by 3 */ 113 114 /* we reuse the variable m to store the cube root, since it is not needed 115 any more: we just need to know if the root is exact */ 116 inexact = mpz_root (m, m, 3) == 0; 117 118 MPFR_MPZ_SIZEINBASE2 (tmp, m); 119 sh = tmp - n; 120 if (sh > 0) /* we have to flush to 0 the last sh bits from m */ 121 { 122 inexact = inexact || ((mpfr_exp_t) mpz_scan1 (m, 0) < sh); 123 mpz_fdiv_q_2exp (m, m, sh); 124 e += 3 * sh; 125 } 126 127 if (inexact) 128 { 129 if (negative) 130 rnd_mode = MPFR_INVERT_RND (rnd_mode); 131 if (rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDA 132 || (rnd_mode == MPFR_RNDN && mpz_tstbit (m, 0))) 133 inexact = 1, mpz_add_ui (m, m, 1); 134 else 135 inexact = -1; 136 } 137 138 /* either inexact is not zero, and the conversion is exact, i.e. inexact 139 is not changed; or inexact=0, and inexact is set only when 140 rnd_mode=MPFR_RNDN and bit (n+1) from m is 1 */ 141 inexact += mpfr_set_z (y, m, MPFR_RNDN); 142 MPFR_SET_EXP (y, MPFR_GET_EXP (y) + e / 3); 143 144 if (negative) 145 { 146 MPFR_CHANGE_SIGN (y); 147 inexact = -inexact; 148 } 149 150 mpz_clear (m); 151 MPFR_SAVE_EXPO_FREE (expo); 152 return mpfr_check_range (y, inexact, rnd_mode); 153} 154