1 /* mpfr_mpn_exp -- auxiliary function for mpfr_get_str and mpfr_set_str
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
24 #define MPFR_NEED_LONGLONG_H
25 #include "mpfr-impl.h"
26
27 /* count the number of significant bits of e, i.e.,
28 nbits(mpfr_exp_t) - count_leading_zeros (e) */
29 static int
nbits_mpfr_exp_t(mpfr_exp_t e)30 nbits_mpfr_exp_t (mpfr_exp_t e)
31 {
32 int nbits = 0;
33
34 MPFR_ASSERTD (e > 0);
35 while (e >= 0x10000)
36 {
37 e >>= 16;
38 nbits += 16;
39 }
40 MPFR_ASSERTD (e <= 0xffff);
41 if (e >= 0x100)
42 {
43 e >>= 8;
44 nbits += 8;
45 }
46 MPFR_ASSERTD (e <= 0xff);
47 if (e >= 0x10)
48 {
49 e >>= 4;
50 nbits += 4;
51 }
52 MPFR_ASSERTD (e <= 0xf);
53 if (e >= 4)
54 {
55 e >>= 2;
56 nbits += 2;
57 }
58 MPFR_ASSERTD (e <= 3);
59 /* now e = 1, 2, or 3 */
60 return nbits + 1 + (e >= 2);
61 }
62
63 /* this function computes an approximation to b^e in {a, n}, with exponent
64 stored in exp_r. The computed value is rounded toward zero (truncated).
65 It returns an integer f such that the final error is bounded by 2^f ulps,
66 that is:
67 a*2^exp_r <= b^e <= 2^exp_r (a + 2^f),
68 where a represents {a, n}, i.e. the integer
69 a[0] + a[1]*B + ... + a[n-1]*B^(n-1) where B=2^GMP_NUMB_BITS
70
71 Return -1 is the result is exact.
72 Return -2 if an overflow occurred in the computation of exp_r.
73 */
74
75 int
mpfr_mpn_exp(mp_limb_t * a,mpfr_exp_t * exp_r,int b,mpfr_exp_t e,size_t n)76 mpfr_mpn_exp (mp_limb_t *a, mpfr_exp_t *exp_r, int b, mpfr_exp_t e, size_t n)
77 {
78 mp_limb_t *c, B;
79 mpfr_exp_t f, h;
80 int i;
81 int t; /* number of bits in e */
82 size_t bits, n1;
83 unsigned int error; /* (number - 1) of loops a^2b inexact */
84 /* error == t means no error */
85 int err_s_a2 = 0;
86 int err_s_ab = 0; /* number of error when shift A^2, AB */
87 MPFR_TMP_DECL(marker);
88
89 MPFR_ASSERTN (n > 0 && n <= ((size_t) -1) / GMP_NUMB_BITS);
90 MPFR_ASSERTN (e > 0);
91 MPFR_ASSERTN (2 <= b && b <= 62);
92
93 MPFR_TMP_MARK(marker);
94
95 /* initialization of a, b, f, h */
96
97 /* normalize the base */
98 B = (mp_limb_t) b;
99 count_leading_zeros (h, B);
100
101 bits = GMP_NUMB_BITS - h;
102
103 B = B << h;
104 h = - h;
105
106 /* allocate space for A and set it to B */
107 c = MPFR_TMP_LIMBS_ALLOC (2 * n);
108 a [n - 1] = B;
109 MPN_ZERO (a, n - 1);
110 /* initial exponent for A: invariant is A = {a, n} * 2^f */
111 f = h - (n - 1) * GMP_NUMB_BITS;
112
113 /* determine number of bits in e */
114 t = nbits_mpfr_exp_t (e);
115
116 error = t;
117 /* t, error <= bitsize(mpfr_exp_t) */
118 MPFR_ASSERTD (error >= 0);
119
120 MPN_ZERO (c, 2 * n);
121
122 for (i = t - 2; i >= 0; i--)
123 {
124 /* determine precision needed */
125 bits = n * GMP_NUMB_BITS - mpn_scan1 (a, 0);
126 n1 = (n * GMP_NUMB_BITS - bits) / GMP_NUMB_BITS;
127
128 /* square of A : {c+2n1, 2(n-n1)} = {a+n1, n-n1}^2 */
129 /* TODO: we should use a short square here, but this needs to redo
130 the error analysis */
131 mpn_sqr (c + 2 * n1, a + n1, n - n1);
132
133 /* set {c+n, 2n1-n} to 0 : {c, n} = {a, n}^2*K^n */
134
135 /* check overflow on f */
136 if (MPFR_UNLIKELY (f < MPFR_EXP_MIN / 2 || f > MPFR_EXP_MAX / 2))
137 {
138 overflow:
139 MPFR_TMP_FREE(marker);
140 return -2;
141 }
142 /* FIXME: Could f = 2 * f + n * GMP_NUMB_BITS be used? */
143 f = 2 * f;
144 MPFR_SADD_OVERFLOW (f, f, n * GMP_NUMB_BITS,
145 mpfr_exp_t, mpfr_uexp_t,
146 MPFR_EXP_MIN, MPFR_EXP_MAX,
147 goto overflow, goto overflow);
148 if (MPFR_LIMB_MSB (c[2*n - 1]) == 0)
149 {
150 /* shift A by one bit to the left */
151 mpn_lshift (a, c + n, n, 1);
152 a[0] |= mpn_lshift (c + n - 1, c + n - 1, 1, 1);
153 f --;
154 if (error != t)
155 err_s_a2 ++;
156 }
157 else
158 MPN_COPY (a, c + n, n);
159
160 if (error == t && 2 * n1 <= n &&
161 mpn_scan1 (c + 2 * n1, 0) < (n - 2 * n1) * GMP_NUMB_BITS)
162 error = i;
163
164 if ((e >> i) & 1)
165 {
166 /* multiply A by B */
167 c[2 * n - 1] = mpn_mul_1 (c + n - 1, a, n, B);
168 f += h + GMP_NUMB_BITS;
169 if ((c[2 * n - 1] & MPFR_LIMB_HIGHBIT) == 0)
170 { /* shift A by one bit to the left */
171 mpn_lshift (a, c + n, n, 1);
172 a[0] |= mpn_lshift (c + n - 1, c + n - 1, 1, 1);
173 f --;
174 }
175 else
176 {
177 MPN_COPY (a, c + n, n);
178 if (error != t)
179 err_s_ab ++;
180 }
181 if (error == t && c[n - 1] != 0)
182 error = i;
183 }
184 }
185
186 MPFR_ASSERTD (error >= 0);
187 MPFR_ASSERTD (err_s_a2 >= 0);
188 MPFR_ASSERTD (err_s_ab >= 0);
189
190 MPFR_TMP_FREE(marker);
191
192 *exp_r = f;
193
194 if (error == t)
195 return -1; /* result is exact */
196 else /* error <= t-2 <= bitsize(mpfr_exp_t)-2
197 err_s_ab, err_s_a2 <= t-1 */
198 {
199 /* if there are p loops after the first inexact result, with
200 j shifts in a^2 and l shifts in a*b, then the final error is
201 at most 2^(p+ceil((j+1)/2)+l+1)*ulp(res).
202 This is bounded by 2^(5/2*t-1/2) where t is the number of bits of e.
203 */
204 return error + err_s_ab + err_s_a2 / 2 + 3; /* <= 5t/2-1/2 */
205 }
206 }
207