synch with TL 37803
[luatex.git] / source / libs / mpfr / mpfr-3.1.3 / src / cbrt.c
blobd6858bc92faec05e36d4b86b17b0cab7f5b476d2
1 /* mpfr_cbrt -- cube root function.
3 Copyright 2002-2015 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
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.
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.
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 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
26 /* The computation of y = x^(1/3) is done as follows:
28 Let x = sign * m * 2^(3*e) where m is an integer
30 with 2^(3n-3) <= m < 2^(3n) where n = PREC(y)
32 and m = s^3 + r where 0 <= r and m < (s+1)^3
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
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
44 int
45 mpfr_cbrt (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
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);
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));
58 /* special values */
59 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
61 if (MPFR_IS_NAN (x))
63 MPFR_SET_NAN (y);
64 MPFR_RET_NAN;
66 else if (MPFR_IS_INF (x))
68 MPFR_SET_INF (y);
69 MPFR_SET_SAME_SIGN (y, x);
70 MPFR_RET (0);
72 /* case 0: cbrt(+/- 0) = +/- 0 */
73 else /* x is necessarily 0 */
75 MPFR_ASSERTD (MPFR_IS_ZERO (x));
76 MPFR_SET_ZERO (y);
77 MPFR_SET_SAME_SIGN (y, x);
78 MPFR_RET (0);
82 /* General case */
83 MPFR_SAVE_EXPO_MARK (expo);
84 mpz_init (m);
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) */
94 MPFR_MPZ_SIZEINBASE2 (size_m, m);
95 n = MPFR_PREC (y) + (rnd_mode == MPFR_RNDN);
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)
103 mpz_mul_2exp (m, m, sh);
104 e = e - sh;
106 else if (r > 0)
108 mpz_mul_2exp (m, m, r);
109 e = e - r;
112 /* invariant: x = m*2^e, with e divisible by 3 */
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;
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 */
122 inexact = inexact || ((mpfr_exp_t) mpz_scan1 (m, 0) < sh);
123 mpz_fdiv_q_2exp (m, m, sh);
124 e += 3 * sh;
127 if (inexact)
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;
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);
144 if (negative)
146 MPFR_CHANGE_SIGN (y);
147 inexact = -inexact;
150 mpz_clear (m);
151 MPFR_SAVE_EXPO_FREE (expo);
152 return mpfr_check_range (y, inexact, rnd_mode);