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