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