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