1 /* Test file for mpfr_ui_pow and mpfr_ui_pow_ui.
2
3 Copyright 2001-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 static void
26 test1 (void)
27 {
28 mpfr_t x, y, z, a;
29 int res1, res2;
30
31 mpfr_init2 (x, 32);
32 mpfr_init2 (y, 65);
33 mpfr_init2 (z, 17);
34 mpfr_init2 (a, 17);
35
36 mpfr_set_str_binary (x, "-0.101110001001011011011e-9");
37 mpfr_ui_pow (y, 7, x, MPFR_RNDN);
38
39 mpfr_set_prec (x, 40);
40 mpfr_set_str_binary (x, "-0.1100101100101111011001010010110011110110E-1");
41 mpfr_set_prec (y, 74);
42 mpfr_ui_pow (y, 8, x, MPFR_RNDN);
43 mpfr_set_prec (x, 74);
44 mpfr_set_str_binary (x, "0.11100000010100111101000011111011011010011000011000101011010011010101000011E-1");
45 if (mpfr_cmp (x, y))
46 {
47 printf ("Error for input of 40 bits, output of 74 bits\n");
48 exit (1);
49 }
50
51 /* Check for ui_pow_ui */
52 mpfr_ui_pow_ui (x, 0, 1, MPFR_RNDN);
53 MPFR_ASSERTN(mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_POS (x));
54 mpfr_ui_pow_ui (x, 0, 4, MPFR_RNDN);
55 MPFR_ASSERTN(mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_POS (x));
56 res1 = mpfr_ui_pow_ui (z, 17, 42, MPFR_RNDD);
57 mpfr_set_ui (x, 17, MPFR_RNDN);
58 mpfr_set_ui (y, 42, MPFR_RNDN);
59 res2 = mpfr_pow (a, x, y, MPFR_RNDD);
60 if (mpfr_cmp (z, a) || res1 != res2)
61 {
62 printf ("Error for ui_pow_ui for 17^42\n"
63 "Inexact1 = %d Inexact2 = %d\n", res1, res2);
64 mpfr_dump (z);
65 mpfr_dump (a);
66 exit (1);
67 }
68 mpfr_set_prec (x, 2);
69 mpfr_ui_pow_ui (x, 65537, 65535, MPFR_RNDN);
70 if (mpfr_cmp_str (x, "0.11E1048562", 2, MPFR_RNDN) != 0)
71 {
72 printf ("Error for ui_pow_ui for 65537 ^65535 with 2 bits of precision\n");
73 mpfr_dump (x);
74 exit (1);
75 }
76 mpfr_clear (x);
77 mpfr_clear (y);
78 mpfr_clear (z);
79 mpfr_clear (a);
80 }
81
82 static void
83 check1 (mpfr_ptr x, mpfr_prec_t prec, unsigned long nt, mpfr_rnd_t rnd)
84 {
85 mpfr_t y, z, t;
86 int inexact, compare, compare2;
87 mpfr_prec_t yprec;
88 mpfr_exp_t err;
89
90 yprec = prec + 10;
91
92 mpfr_init (y);
93 mpfr_init (z);
94 mpfr_init (t);
95 mpfr_set_prec (y, yprec);
96 mpfr_set_prec (z, prec);
97 mpfr_set_prec (t, prec);
98
99 compare = mpfr_ui_pow (y, nt, x, rnd);
100 err = (rnd == MPFR_RNDN) ? yprec + 1 : yprec;
101 if (mpfr_can_round (y, err, rnd, rnd, prec))
102 {
103 mpfr_set (t, y, rnd);
104 inexact = mpfr_ui_pow (z, nt, x, rnd);
105 if (mpfr_cmp (t, z))
106 {
107 printf ("results differ for x=");
108 mpfr_out_str (stdout, 2, prec, x, MPFR_RNDN);
109 printf (" n=%lu", nt);
110 printf (" prec=%u rnd_mode=%s\n", (unsigned) prec,
111 mpfr_print_rnd_mode (rnd));
112 printf ("got ");
113 mpfr_dump (z);
114 printf ("expected ");
115 mpfr_dump (t);
116 printf ("approx ");
117 mpfr_dump (y);
118 exit (1);
119 }
120 compare2 = mpfr_cmp (t, y);
121 /* if rounding to nearest, cannot know the sign of t - f(x)
122 because of composed rounding: y = o(f(x)) and t = o(y) */
123 if ((rnd != MPFR_RNDN) && (compare * compare2 >= 0))
124 compare = compare + compare2;
125 else
126 compare = inexact; /* cannot determine sign(t-f(x)) */
127 if (((inexact == 0) && (compare != 0)) ||
128 ((inexact > 0) && (compare <= 0)) ||
129 ((inexact < 0) && (compare >= 0)))
130 {
131 printf ("Wrong inexact flag for rnd=%s: expected %d, got %d\n",
132 mpfr_print_rnd_mode (rnd), compare, inexact);
133 printf ("x="); mpfr_dump (x);
134 printf ("y="); mpfr_dump (y);
135 printf ("t="); mpfr_dump (t);
136 exit (1);
137 }
138 }
139
140 mpfr_clear (y);
141 mpfr_clear (z);
142 mpfr_clear (t);
143 }
144
145 static void
146 huge (void)
147 {
148 mpfr_exp_t old_emin, old_emax;
149 mpfr_t x;
150
151 old_emin = mpfr_get_emin ();
152 old_emax = mpfr_get_emax ();
153
154 set_emin (MPFR_EMIN_MIN);
155 set_emax (MPFR_EMAX_MAX);
156
157 mpfr_init2 (x, 8);
158
159 /* The purpose of this test is more to check that mpfr_ui_pow_ui
160 terminates (without taking much memory) rather than checking
161 the value of x. On 2023-02-13, the +Inf case was not handled
162 in the Ziv iteration, yielding an infinite loop, affecting
163 mpfr_log10 in particular. See
164 commit 90de094f0d9c309daca707aa227470d810866616
165 */
166 mpfr_ui_pow_ui (x, 5, ULONG_MAX, MPFR_RNDN);
167 if (MPFR_EMAX_MAX <= ULONG_MAX) /* true with default _MPFR_EXP_FORMAT */
168 MPFR_ASSERTN (MPFR_IS_INF (x));
169
170 mpfr_clear (x);
171
172 set_emin (old_emin);
173 set_emax (old_emax);
174 }
175
176 static void
177 test2 (void)
178 {
179 mpfr_t x, y, z, t;
180 mpfr_prec_t prec;
181 mpfr_rnd_t rnd;
182 unsigned int n;
183 mpfr_prec_t p0 = 2, p1 = 100;
184 unsigned int N = 20;
185
186 mpfr_init2 (x, 2);
187 mpfr_init2 (y, 2);
188 mpfr_init2 (z, 38);
189 mpfr_init2 (t, 6);
190
191 /* check exact power */
192 mpfr_set_str_binary (t, "0.110000E5");
193 mpfr_ui_pow (z, 3, t, MPFR_RNDN);
194
195 mpfr_set_str (x, "-0.5", 10, MPFR_RNDZ);
196 mpfr_ui_pow (y, 4, x, MPFR_RNDD);
197 if (mpfr_cmp_ui_2exp (y, 1, -1))
198 {
199 printf ("Error for 4^(-0.5), prec=2, MPFR_RNDD\n");
200 printf ("expected 0.5, got ");
201 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN);
202 printf ("\n");
203 exit (1);
204 }
205
206 /* problem found by Kevin on spe175.testdrive.compaq.com
207 (03 Sep 2003), ia64 under HP-UX */
208 mpfr_set_str (x, "0.5", 10, MPFR_RNDN);
209 mpfr_ui_pow (y, 398441521, x, MPFR_RNDN);
210 if (mpfr_cmp_ui_2exp (y, 1, 14))
211 {
212 printf ("Error for 398441521^(0.5), prec=2, MPFR_RNDN\n");
213 printf ("expected 1.0e14, got ");
214 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN);
215 printf ("\n");
216 exit (1);
217 }
218
219 mpfr_set_str (x, "0.5", 10, MPFR_RNDN);
220 check1 (x, 2, 398441521, MPFR_RNDN); /* 398441521 = 19961^2 */
221
222 /* generic test */
223 for (prec = p0; prec <= p1; prec++)
224 {
225 mpfr_set_prec (x, prec);
226 for (n = 0; n < N; n++)
227 {
228 int nt;
229 nt = randlimb () & INT_MAX;
230 mpfr_urandomb (x, RANDS);
231 rnd = RND_RAND_NO_RNDF ();
232 check1 (x, prec, nt, rnd);
233 }
234 }
235
236 mpfr_clears (x, y, z, t, (mpfr_ptr) 0);
237 }
238
239 /* reverse the arguments n and x for tgeneric_ui.c */
240 static int
241 ui_pow_rev (mpfr_ptr y, mpfr_srcptr x, unsigned long n, mpfr_rnd_t rnd)
242 {
243 return mpfr_ui_pow (y, n, x, rnd);
244 }
245
246 #define TEST_FUNCTION ui_pow_rev
247 #define INTEGER_TYPE unsigned long
248 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), 1, RANDS)
249 #define INT_RAND_FUNCTION() \
250 (randlimb () % 16 == 0 ? randulong () : (unsigned long) (randlimb () % 32))
251 #include "tgeneric_ui.c"
252
253 int
254 main (int argc, char *argv[])
255 {
256 mpfr_t x, y;
257 unsigned long int n;
258
259 tests_start_mpfr ();
260
261 mpfr_init (x);
262 mpfr_init (y);
263
264 do n = randlimb (); while (n <= 1);
265
266 MPFR_SET_INF (x);
267 mpfr_ui_pow (y, n, x, MPFR_RNDN);
268 if (! MPFR_IS_INF (y))
269 {
270 printf ("evaluation of function at INF does not return INF\n");
271 exit (1);
272 }
273
274 MPFR_CHANGE_SIGN (x);
275 mpfr_ui_pow (y, n, x, MPFR_RNDN);
276 if (! MPFR_IS_ZERO (y))
277 {
278 printf ("evaluation of function in -INF does not return 0\n");
279 exit (1);
280 }
281
282 MPFR_SET_NAN (x);
283 mpfr_ui_pow (y, n, x, MPFR_RNDN);
284 if (! MPFR_IS_NAN (y))
285 {
286 printf ("evaluation of function in NAN does not return NAN\n");
287 exit (1);
288 }
289
290 mpfr_clear (x);
291 mpfr_clear (y);
292
293 test1 ();
294 test2 ();
295 huge ();
296
297 test_generic_ui (MPFR_PREC_MIN, 100, 100);
298
299 tests_end_mpfr ();
300 return 0;
301 }