1 /* Test file for mpfr_pow_z -- power function x^z with z a MPZ
2
3 Copyright 2005-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 #include "mpfr-test.h"
24
25 #define PRINT_ERROR(str) \
26 do { printf ("Error for %s\n", str); exit (1); } while (0)
27
28 static void
29 check_special (void)
30 {
31 mpfr_t x, y;
32 mpz_t z;
33 int res;
34
35 mpfr_init (x);
36 mpfr_init (y);
37 mpz_init (z);
38
39 /* x^0 = 1 except for NAN */
40 mpfr_set_ui (x, 23, MPFR_RNDN);
41 mpz_set_ui (z, 0);
42 res = mpfr_pow_z (y, x, z, MPFR_RNDN);
43 if (res != 0 || mpfr_cmp_ui (y, 1) != 0)
44 PRINT_ERROR ("23^0");
45 mpfr_set_nan (x);
46 res = mpfr_pow_z (y, x, z, MPFR_RNDN);
47 if (res != 0 || mpfr_nan_p (y) || mpfr_cmp_si (y, 1) != 0)
48 PRINT_ERROR ("NAN^0");
49 mpfr_set_inf (x, 1);
50 res = mpfr_pow_z (y, x, z, MPFR_RNDN);
51 if (res != 0 || mpfr_cmp_ui (y, 1) != 0)
52 PRINT_ERROR ("INF^0");
53
54 /* sINF^N = INF if s==1 or n even if N > 0*/
55 mpz_set_ui (z, 42);
56 mpfr_set_inf (x, 1);
57 res = mpfr_pow_z (y, x, z, MPFR_RNDN);
58 if (res != 0 || mpfr_inf_p (y) == 0 || mpfr_sgn (y) <= 0)
59 PRINT_ERROR ("INF^42");
60 mpfr_set_inf (x, -1);
61 res = mpfr_pow_z (y, x, z, MPFR_RNDN);
62 if (res != 0 || mpfr_inf_p (y) == 0 || mpfr_sgn (y) <= 0)
63 PRINT_ERROR ("-INF^42");
64 mpz_set_ui (z, 17);
65 mpfr_set_inf (x, 1);
66 res = mpfr_pow_z (y, x, z, MPFR_RNDN);
67 if (res != 0 || mpfr_inf_p (y) == 0 || mpfr_sgn (y) <= 0)
68 PRINT_ERROR ("INF^17");
69 mpfr_set_inf (x, -1);
70 res = mpfr_pow_z (y, x, z, MPFR_RNDN);
71 if (res != 0 || mpfr_inf_p (y) == 0 || mpfr_sgn (y) >= 0)
72 PRINT_ERROR ("-INF^17");
73
74 mpz_set_si (z, -42);
75 mpfr_set_inf (x, 1);
76 res = mpfr_pow_z (y, x, z, MPFR_RNDN);
77 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_IS_NEG (y))
78 PRINT_ERROR ("INF^-42");
79 mpfr_set_inf (x, -1);
80 res = mpfr_pow_z (y, x, z, MPFR_RNDN);
81 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_IS_NEG (y))
82 PRINT_ERROR ("-INF^-42");
83 mpz_set_si (z, -17);
84 mpfr_set_inf (x, 1);
85 res = mpfr_pow_z (y, x, z, MPFR_RNDN);
86 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_IS_NEG (y))
87 PRINT_ERROR ("INF^-17");
88 mpfr_set_inf (x, -1);
89 res = mpfr_pow_z (y, x, z, MPFR_RNDN);
90 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_IS_POS (y))
91 PRINT_ERROR ("-INF^-17");
92
93 /* s0^N = +0 if s==+ or n even if N > 0*/
94 mpz_set_ui (z, 42);
95 MPFR_SET_ZERO (x); MPFR_SET_POS (x);
96 res = mpfr_pow_z (y, x, z, MPFR_RNDN);
97 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_IS_NEG (y))
98 PRINT_ERROR ("+0^42");
99 MPFR_SET_NEG (x);
100 res = mpfr_pow_z (y, x, z, MPFR_RNDN);
101 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_IS_NEG (y))
102 PRINT_ERROR ("-0^42");
103 mpz_set_ui (z, 17);
104 MPFR_SET_POS (x);
105 res = mpfr_pow_z (y, x, z, MPFR_RNDN);
106 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_IS_NEG (y))
107 PRINT_ERROR ("+0^17");
108 MPFR_SET_NEG (x);
109 res = mpfr_pow_z (y, x, z, MPFR_RNDN);
110 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_IS_POS (y))
111 PRINT_ERROR ("-0^17");
112
113 mpz_set_si (z, -42);
114 MPFR_SET_ZERO (x); MPFR_SET_POS (x);
115 res = mpfr_pow_z (y, x, z, MPFR_RNDN);
116 if (res != 0 || mpfr_inf_p (y) == 0 || MPFR_IS_NEG (y))
117 PRINT_ERROR ("+0^-42");
118 MPFR_SET_NEG (x);
119 res = mpfr_pow_z (y, x, z, MPFR_RNDN);
120 if (res != 0 || mpfr_inf_p (y) == 0 || MPFR_IS_NEG (y))
121 PRINT_ERROR ("-0^-42");
122 mpz_set_si (z, -17);
123 MPFR_SET_POS (x);
124 res = mpfr_pow_z (y, x, z, MPFR_RNDN);
125 if (res != 0 || mpfr_inf_p (y) == 0 || MPFR_IS_NEG (y))
126 PRINT_ERROR ("+0^-17");
127 MPFR_SET_NEG (x);
128 res = mpfr_pow_z (y, x, z, MPFR_RNDN);
129 if (res != 0 || mpfr_inf_p (y) == 0 || MPFR_IS_POS (y))
130 PRINT_ERROR ("-0^-17");
131
132 mpz_clear (z);
133 mpfr_clear (y);
134 mpfr_clear (x);
135 }
136
137 static void
138 check_integer (mpfr_prec_t begin, mpfr_prec_t end, unsigned long max)
139 {
140 mpfr_t x, y1, y2;
141 mpz_t z;
142 unsigned long i, n;
143 mpfr_prec_t p;
144 int res1, res2;
145 mpfr_rnd_t rnd;
146
147 mpfr_inits2 (begin, x, y1, y2, (mpfr_ptr) 0);
148 mpz_init (z);
149 for (p = begin ; p < end ; p+=4)
150 {
151 mpfr_set_prec (x, p);
152 mpfr_set_prec (y1, p);
153 mpfr_set_prec (y2, p);
154 for (i = 0 ; i < max ; i++)
155 {
156 mpz_urandomb (z, RANDS, GMP_NUMB_BITS);
157 if ((i & 1) != 0)
158 mpz_neg (z, z);
159 mpfr_urandomb (x, RANDS);
160 mpfr_mul_2ui (x, x, 1, MPFR_RNDN); /* 0 <= x < 2 */
161 rnd = RND_RAND ();
162 if (mpz_fits_slong_p (z))
163 {
164 n = mpz_get_si (z);
165 /* printf ("New test for x=%ld\nCheck Pow_si\n", n); */
166 res1 = mpfr_pow_si (y1, x, n, rnd);
167 /* printf ("Check pow_z\n"); */
168 res2 = mpfr_pow_z (y2, x, z, rnd);
169 if (! mpfr_equal_p (y1, y2))
170 {
171 printf ("Error for p = %lu, z = %lu, rnd = %s and x = ",
172 (unsigned long) p, n, mpfr_print_rnd_mode (rnd));
173 mpfr_dump (x);
174 printf ("Ypowsi = "); mpfr_dump (y1);
175 printf ("Ypowz = "); mpfr_dump (y2);
176 exit (1);
177 }
178 /* The ternary value is unspecified with MPFR_RNDF. */
179 if (rnd != MPFR_RNDF && ! SAME_SIGN (res1, res2))
180 {
181 printf ("Wrong ternary value for p = %lu, z = %lu, rnd = %s"
182 " and x = ", (unsigned long) p, n,
183 mpfr_print_rnd_mode (rnd));
184 mpfr_dump (x);
185 printf ("Ypowsi(inex = %2d) = ", res1); mpfr_dump (y1);
186 printf ("Ypowz (inex = %2d) = ", res2); mpfr_dump (y2);
187 exit (1);
188 }
189 }
190 } /* for i */
191 } /* for p */
192 mpfr_clears (x, y1, y2, (mpfr_ptr) 0);
193 mpz_clear (z);
194 }
195
196 static void
197 check_regression (void)
198 {
199 mpfr_t x, y;
200 mpz_t z;
201 int res1, res2;
202
203 mpz_init_set_ui (z, 2026876995);
204 mpfr_init2 (x, 122);
205 mpfr_init2 (y, 122);
206
207 mpfr_set_str_binary (x, "0.10000010010000111101001110100101101010011110011100001111000001001101000110011001001001001011001011010110110110101000111011E1");
208 res1 = mpfr_pow_z (y, x, z, MPFR_RNDU);
209 res2 = mpfr_pow_ui (x, x, 2026876995UL, MPFR_RNDU);
210 if (! mpfr_equal_p (x, y) || ! SAME_SIGN (res1, res2))
211 {
212 printf ("Regression (1) tested failed (%d=?%d)\n",res1, res2);
213 printf ("pow_ui: "); mpfr_dump (x);
214 printf ("pow_z: "); mpfr_dump (y);
215
216 exit (1);
217 }
218
219 mpfr_clear (x);
220 mpfr_clear (y);
221 mpz_clear (z);
222 }
223
224 /* Bug found by Kevin P. Rauch */
225 static void
226 bug20071104 (void)
227 {
228 mpfr_t x, y;
229 mpz_t z;
230 int inex;
231
232 mpz_init_set_si (z, -2);
233 mpfr_inits2 (20, x, y, (mpfr_ptr) 0);
234 mpfr_set_ui (x, 0, MPFR_RNDN);
235 mpfr_nextbelow (x); /* x = -2^(emin-1) */
236 mpfr_clear_flags ();
237 inex = mpfr_pow_z (y, x, z, MPFR_RNDN);
238 if (! mpfr_inf_p (y) || MPFR_IS_NEG (y))
239 {
240 printf ("Error in bug20071104: expected +Inf, got ");
241 mpfr_dump (y);
242 exit (1);
243 }
244 if (inex <= 0)
245 {
246 printf ("Error in bug20071104: bad ternary value (%d)\n", inex);
247 exit (1);
248 }
249 if (__gmpfr_flags != (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT))
250 {
251 printf ("Error in bug20071104: bad flags (%u)\n",
252 (unsigned int) __gmpfr_flags);
253 exit (1);
254 }
255 mpfr_clears (x, y, (mpfr_ptr) 0);
256 mpz_clear (z);
257 }
258
259 static void
260 check_overflow (void)
261 {
262 mpfr_t a;
263 mpz_t z;
264 unsigned long n;
265 int res;
266
267 mpfr_init2 (a, 53);
268
269 mpfr_set_str_binary (a, "1E10");
270 mpz_init_set_ui (z, ULONG_MAX);
271 res = mpfr_pow_z (a, a, z, MPFR_RNDN);
272 if (! MPFR_IS_INF (a) || MPFR_IS_NEG (a) || res <= 0)
273 {
274 printf ("Error for (1e10)^ULONG_MAX, expected +Inf,\ngot ");
275 mpfr_dump (a);
276 exit (1);
277 }
278
279 /* Bug in pow_z.c up to r5109: if x = y (same mpfr_t argument), the
280 input argument is negative and not a power of two, z is positive
281 and odd, an overflow or underflow occurs, and the temporary result
282 res is positive, then the result gets a wrong sign (positive
283 instead of negative). */
284 mpfr_set_str_binary (a, "-1.1E10");
285 n = (ULONG_MAX ^ (ULONG_MAX >> 1)) + 1;
286 mpz_set_ui (z, n);
287 res = mpfr_pow_z (a, a, z, MPFR_RNDN);
288 if (! MPFR_IS_INF (a) || MPFR_IS_POS (a) || res >= 0)
289 {
290 printf ("Error for (-1e10)^%lu, expected -Inf,\ngot ", n);
291 mpfr_dump (a);
292 exit (1);
293 }
294
295 mpfr_clear (a);
296 mpz_clear (z);
297 }
298
299 /* bug reported by Carl Witty (32-bit architecture) */
300 static void
301 bug20080223 (void)
302 {
303 mpfr_t a, exp, answer;
304
305 mpfr_init2 (a, 53);
306 mpfr_init2 (exp, 53);
307 mpfr_init2 (answer, 53);
308
309 mpfr_set_si (exp, -1073741824, MPFR_RNDN);
310
311 mpfr_set_str (a, "1.999999999", 10, MPFR_RNDN);
312 /* a = 562949953139837/2^48 */
313 mpfr_pow (answer, a, exp, MPFR_RNDN);
314 mpfr_set_str_binary (a, "0.110110101111011001110000111111100011101000111011101E-1073741823");
315 MPFR_ASSERTN(mpfr_cmp0 (answer, a) == 0);
316
317 mpfr_clear (a);
318 mpfr_clear (exp);
319 mpfr_clear (answer);
320 }
321
322 static void
323 bug20080904 (void)
324 {
325 mpz_t exp;
326 mpfr_t a, answer;
327 mpfr_exp_t emin_default;
328
329 mpz_init (exp);
330 mpfr_init2 (a, 70);
331 mpfr_init2 (answer, 70);
332
333 emin_default = mpfr_get_emin ();
334 set_emin (MPFR_EMIN_MIN);
335
336 mpz_set_str (exp, "-4eb92f8c7b7bf81e", 16);
337 mpfr_set_str_binary (a, "1.110000101110100110100011111000011110111101000011111001111001010011100");
338
339 mpfr_pow_z (answer, a, exp, MPFR_RNDN);
340 /* The correct result is near 2^(-2^62), so it underflows when
341 MPFR_EMIN_MIN > -2^62 (i.e. with 32 and 64 bits machines). */
342 mpfr_set_str (a, "AA500C0D7A69275DBp-4632850503556296886", 16, MPFR_RNDN);
343 if (! mpfr_equal_p (answer, a))
344 {
345 printf ("Error in bug20080904:\n");
346 printf ("Expected ");
347 mpfr_out_str (stdout, 16, 0, a, MPFR_RNDN);
348 putchar ('\n');
349 printf ("Got ");
350 mpfr_out_str (stdout, 16, 0, answer, MPFR_RNDN);
351 putchar ('\n');
352 exit (1);
353 }
354
355 set_emin (emin_default);
356
357 mpz_clear (exp);
358 mpfr_clear (a);
359 mpfr_clear (answer);
360 }
361
362 int
363 main (void)
364 {
365 tests_start_mpfr ();
366
367 check_special ();
368
369 check_integer (2, 163, 100);
370 check_regression ();
371 bug20071104 ();
372 bug20080223 ();
373 bug20080904 ();
374 check_overflow ();
375
376 tests_end_mpfr ();
377 return 0;
378 }