OLD | NEW |
| (Empty) |
1 /* mpfr_cbrt -- cube root function. | |
2 | |
3 Copyright 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundatio
n, Inc. | |
4 Contributed by the Arenaire and Cacao 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 2.1 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.LIB. If not, write to | |
20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, | |
21 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 | |
44 int | |
45 mpfr_cbrt (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) | |
46 { | |
47 mpz_t m; | |
48 mp_exp_t e, r, sh; | |
49 mp_prec_t n, size_m, tmp; | |
50 int inexact, negative; | |
51 MPFR_SAVE_EXPO_DECL (expo); | |
52 | |
53 /* special values */ | |
54 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) | |
55 { | |
56 if (MPFR_IS_NAN (x)) | |
57 { | |
58 MPFR_SET_NAN (y); | |
59 MPFR_RET_NAN; | |
60 } | |
61 else if (MPFR_IS_INF (x)) | |
62 { | |
63 MPFR_SET_INF (y); | |
64 MPFR_SET_SAME_SIGN (y, x); | |
65 MPFR_RET (0); | |
66 } | |
67 /* case 0: cbrt(+/- 0) = +/- 0 */ | |
68 else /* x is necessarily 0 */ | |
69 { | |
70 MPFR_ASSERTD (MPFR_IS_ZERO (x)); | |
71 MPFR_SET_ZERO (y); | |
72 MPFR_SET_SAME_SIGN (y, x); | |
73 MPFR_RET (0); | |
74 } | |
75 } | |
76 | |
77 /* General case */ | |
78 MPFR_SAVE_EXPO_MARK (expo); | |
79 mpz_init (m); | |
80 | |
81 e = mpfr_get_z_exp (m, x); /* x = m * 2^e */ | |
82 if ((negative = MPFR_IS_NEG(x))) | |
83 mpz_neg (m, m); | |
84 r = e % 3; | |
85 if (r < 0) | |
86 r += 3; | |
87 /* x = (m*2^r) * 2^(e-r) = (m*2^r) * 2^(3*q) */ | |
88 | |
89 MPFR_MPZ_SIZEINBASE2 (size_m, m); | |
90 n = MPFR_PREC (y) + (rnd_mode == GMP_RNDN); | |
91 | |
92 /* we want 3*n-2 <= size_m + 3*sh + r <= 3*n | |
93 i.e. 3*sh + size_m + r <= 3*n */ | |
94 sh = (3 * (mp_exp_t) n - (mp_exp_t) size_m - r) / 3; | |
95 sh = 3 * sh + r; | |
96 if (sh >= 0) | |
97 { | |
98 mpz_mul_2exp (m, m, sh); | |
99 e = e - sh; | |
100 } | |
101 else if (r > 0) | |
102 { | |
103 mpz_mul_2exp (m, m, r); | |
104 e = e - r; | |
105 } | |
106 | |
107 /* invariant: x = m*2^e, with e divisible by 3 */ | |
108 | |
109 /* we reuse the variable m to store the cube root, since it is not needed | |
110 any more: we just need to know if the root is exact */ | |
111 inexact = mpz_root (m, m, 3) == 0; | |
112 | |
113 MPFR_MPZ_SIZEINBASE2 (tmp, m); | |
114 sh = tmp - n; | |
115 if (sh > 0) /* we have to flush to 0 the last sh bits from m */ | |
116 { | |
117 inexact = inexact || ((mp_exp_t) mpz_scan1 (m, 0) < sh); | |
118 mpz_div_2exp (m, m, sh); | |
119 e += 3 * sh; | |
120 } | |
121 | |
122 if (inexact) | |
123 { | |
124 if (negative) | |
125 rnd_mode = MPFR_INVERT_RND (rnd_mode); | |
126 if (rnd_mode == GMP_RNDU | |
127 || (rnd_mode == GMP_RNDN && mpz_tstbit (m, 0))) | |
128 inexact = 1, mpz_add_ui (m, m, 1); | |
129 else | |
130 inexact = -1; | |
131 } | |
132 | |
133 /* either inexact is not zero, and the conversion is exact, i.e. inexact | |
134 is not changed; or inexact=0, and inexact is set only when | |
135 rnd_mode=GMP_RNDN and bit (n+1) from m is 1 */ | |
136 inexact += mpfr_set_z (y, m, GMP_RNDN); | |
137 MPFR_SET_EXP (y, MPFR_GET_EXP (y) + e / 3); | |
138 | |
139 if (negative) | |
140 { | |
141 MPFR_CHANGE_SIGN (y); | |
142 inexact = -inexact; | |
143 } | |
144 | |
145 mpz_clear (m); | |
146 MPFR_SAVE_EXPO_FREE (expo); | |
147 return mpfr_check_range (y, inexact, rnd_mode); | |
148 } | |
OLD | NEW |