1 /* mpfr_pow_ui -- compute the power of a floating-point by a machine integer
2
3 Copyright 1999-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba 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 3 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.LESSER. If not, see
20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
25
26 #ifndef POW_U
27 #define POW_U mpfr_pow_ui
28 #define MPZ_SET_U mpz_set_ui
29 #define UTYPE unsigned long int
30 #define FSPEC "l"
31 #endif
32
33 /* sets y to x^n, and return 0 if exact, non-zero otherwise */
34 int
35 POW_U (mpfr_ptr y, mpfr_srcptr x, UTYPE n, mpfr_rnd_t rnd)
36 {
37 UTYPE m;
38 mpfr_t res;
39 mpfr_prec_t prec, err, nlen;
40 int inexact;
41 mpfr_rnd_t rnd1;
42 MPFR_SAVE_EXPO_DECL (expo);
43 MPFR_ZIV_DECL (loop);
44 MPFR_BLOCK_DECL (flags);
45
46 MPFR_LOG_FUNC
47 (("x[%Pd]=%.*Rg n=%" FSPEC "u rnd=%d",
48 mpfr_get_prec (x), mpfr_log_prec, x, n, rnd),
49 ("y[%Pd]=%.*Rg inexact=%d",
50 mpfr_get_prec (y), mpfr_log_prec, y, inexact));
51
52 /* x^0 = 1 for any x, even a NaN */
53 if (MPFR_UNLIKELY (n == 0))
54 return mpfr_set_ui (y, 1, rnd);
55
56 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
57 {
58 if (MPFR_IS_NAN (x))
59 {
60 MPFR_SET_NAN (y);
61 MPFR_RET_NAN;
62 }
63 else if (MPFR_IS_INF (x))
64 {
65 /* Inf^n = Inf, (-Inf)^n = Inf for n even, -Inf for n odd */
66 if (MPFR_IS_NEG (x) && (n & 1) == 1)
67 MPFR_SET_NEG (y);
68 else
69 MPFR_SET_POS (y);
70 MPFR_SET_INF (y);
71 MPFR_RET (0);
72 }
73 else /* x is zero */
74 {
75 MPFR_ASSERTD (MPFR_IS_ZERO (x));
76 /* 0^n = 0 for any n */
77 MPFR_SET_ZERO (y);
78 if (MPFR_IS_POS (x) || (n & 1) == 0)
79 MPFR_SET_POS (y);
80 else
81 MPFR_SET_NEG (y);
82 MPFR_RET (0);
83 }
84 }
85 else if (MPFR_UNLIKELY (n <= 2))
86 {
87 if (n < 2)
88 /* x^1 = x */
89 return mpfr_set (y, x, rnd);
90 else
91 /* x^2 = sqr(x) */
92 return mpfr_sqr (y, x, rnd);
93 }
94
95 /* Augment exponent range */
96 MPFR_SAVE_EXPO_MARK (expo);
97
98 for (m = n, nlen = 0; m != 0; nlen++, m >>= 1)
99 ;
100 /* 2^(nlen-1) <= n < 2^nlen */
101
102 /* set up initial precision */
103 prec = MPFR_PREC (y) + 3 + GMP_NUMB_BITS
104 + MPFR_INT_CEIL_LOG2 (MPFR_PREC (y));
105 if (prec <= nlen)
106 prec = nlen + 1;
107 mpfr_init2 (res, prec);
108
109 rnd1 = MPFR_IS_POS (x) ? MPFR_RNDU : MPFR_RNDD; /* away */
110
111 MPFR_ZIV_INIT (loop, prec);
112 for (;;)
113 {
114 int i;
115
116 MPFR_ASSERTD (prec > nlen);
117 err = prec - 1 - nlen;
118 /* First step: compute square from x */
119 MPFR_BLOCK (flags,
120 inexact = mpfr_sqr (res, x, MPFR_RNDU);
121 MPFR_ASSERTD (nlen >= 2 && nlen <= INT_MAX);
122 i = nlen;
123 if (n & ((UTYPE) 1 << (i-2)))
124 inexact |= mpfr_mul (res, res, x, rnd1);
125 for (i -= 3; i >= 0 && !MPFR_BLOCK_EXCEP; i--)
126 {
127 inexact |= mpfr_sqr (res, res, MPFR_RNDU);
128 if (n & ((UTYPE) 1 << i))
129 inexact |= mpfr_mul (res, res, x, rnd1);
130 });
131 /* let r(n) be the number of roundings: we have r(2)=1, r(3)=2,
132 and r(2n)=2r(n)+1, r(2n+1)=2r(n)+2, thus r(n)=n-1.
133 Using Higham's method, to each rounding corresponds a factor
134 (1-theta) with 0 <= theta <= 2^(1-p), thus at the end the
135 absolute error is bounded by (n-1)*2^(1-p)*res <= 2*(n-1)*ulp(res)
136 since 2^(-p)*x <= ulp(x). Since n < 2^i, this gives a maximal
137 error of 2^(1+i)*ulp(res).
138 */
139 if (MPFR_LIKELY (inexact == 0
140 || MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags)
141 || MPFR_CAN_ROUND (res, err, MPFR_PREC (y), rnd)))
142 break;
143 /* Actualisation of the precision */
144 MPFR_ZIV_NEXT (loop, prec);
145 mpfr_set_prec (res, prec);
146 }
147 MPFR_ZIV_FREE (loop);
148
149 if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags)))
150 {
151 mpz_t z;
152
153 /* Internal overflow or underflow. However, the approximation error has
154 * not been taken into account. So, let's solve this problem by using
155 * mpfr_pow_z, which can handle it. This case could be improved in the
156 * future, without having to use mpfr_pow_z.
157 */
158 MPFR_LOG_MSG (("Internal overflow or underflow,"
159 " let's use mpfr_pow_z.\n", 0));
160 mpfr_clear (res);
161 MPFR_SAVE_EXPO_FREE (expo);
162 mpz_init (z);
163 MPZ_SET_U (z, n);
164 inexact = mpfr_pow_z (y, x, z, rnd);
165 mpz_clear (z);
166 return inexact;
167 }
168
169 inexact = mpfr_set (y, res, rnd);
170 mpfr_clear (res);
171
172 MPFR_SAVE_EXPO_FREE (expo);
173 return mpfr_check_range (y, inexact, rnd);
174 }