1 /* Test file for mpfr_cmp2.
2
3 Copyright 1999-2003, 2006-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 /* set bit n of x to b, where bit 0 is the most significant one */
26 static void
set_bit(mpfr_ptr x,unsigned int n,int b)27 set_bit (mpfr_ptr x, unsigned int n, int b)
28 {
29 unsigned l;
30 mp_size_t xn;
31
32 xn = (MPFR_PREC(x) - 1) / mp_bits_per_limb;
33 l = n / mp_bits_per_limb;
34 n %= mp_bits_per_limb;
35 n = mp_bits_per_limb - 1 - n;
36 if (b)
37 MPFR_MANT(x)[xn - l] |= MPFR_LIMB_ONE << n;
38 else
39 MPFR_MANT(x)[xn - l] &= ~(MPFR_LIMB_ONE << n);
40 }
41
42 /* check that for x = 1.u 1 v 0^k low(x)
43 y = 1.u 0 v 1^k low(y)
44 mpfr_cmp2 (x, y) returns 1 + |u| + |v| + k for low(x) >= low(y),
45 and 1 + |u| + |v| + k + 1 otherwise */
46 static void
worst_cases(void)47 worst_cases (void)
48 {
49 mpfr_t x, y;
50 unsigned int i, j, k, b, expected;
51 mpfr_prec_t l;
52
53 mpfr_init2 (x, 200);
54 mpfr_init2 (y, 200);
55
56 mpfr_set_ui (y, 1, MPFR_RNDN);
57 for (i = 1; i < MPFR_PREC(x); i++)
58 {
59 mpfr_set_ui (x, 1, MPFR_RNDN);
60 mpfr_div_2ui (y, y, 1, MPFR_RNDN); /* y = 1/2^i */
61
62 l = 0;
63 if (mpfr_cmp2 (x, y, &l) <= 0 || l != 1)
64 {
65 printf ("Error in mpfr_cmp2:\nx=");
66 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
67 printf ("\ny=");
68 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN);
69 printf ("\ngot %lu instead of 1\n", (unsigned long) l);
70 exit (1);
71 }
72
73 mpfr_add (x, x, y, MPFR_RNDN); /* x = 1 + 1/2^i */
74 l = 0;
75 if (mpfr_cmp2 (x, y, &l) <= 0 || l != 0)
76 {
77 printf ("Error in mpfr_cmp2:\nx=");
78 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
79 printf ("\ny=");
80 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN);
81 printf ("\ngot %lu instead of 0\n", (unsigned long) l);
82 exit (1);
83 }
84 }
85
86 for (i = 0; i < 64; i++) /* |u| = i */
87 {
88 mpfr_urandomb (x, RANDS);
89 mpfr_set (y, x, MPFR_RNDN);
90 set_bit (x, i + 1, 1);
91 set_bit (y, i + 1, 0);
92 for (j = 0; j < 64; j++) /* |v| = j */
93 {
94 b = RAND_BOOL ();
95 set_bit (x, i + j + 2, b);
96 set_bit (y, i + j + 2, b);
97 for (k = 0; k < 64; k++)
98 {
99 if (k)
100 set_bit (x, i + j + k + 1, 0);
101 if (k)
102 set_bit (y, i + j + k + 1, 1);
103 set_bit (x, i + j + k + 2, 1);
104 set_bit (y, i + j + k + 2, 0);
105 l = 0; mpfr_cmp2 (x, y, &l);
106 expected = i + j + k + 1;
107 if (l != expected)
108 {
109 printf ("Error in mpfr_cmp2:\nx=");
110 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
111 printf ("\ny=");
112 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN);
113 printf ("\ngot %lu instead of %u\n",
114 (unsigned long) l, expected);
115 exit (1);
116 }
117 set_bit (x, i + j + k + 2, 0);
118 set_bit (x, i + j + k + 3, 0);
119 set_bit (y, i + j + k + 3, 1);
120 l = 0; mpfr_cmp2 (x, y, &l);
121 expected = i + j + k + 2;
122 if (l != expected)
123 {
124 printf ("Error in mpfr_cmp2:\nx=");
125 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
126 printf ("\ny=");
127 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN);
128 printf ("\ngot %lu instead of %u\n",
129 (unsigned long) l, expected);
130 exit (1);
131 }
132 }
133 }
134 }
135
136 mpfr_clear (x);
137 mpfr_clear (y);
138 }
139
140 static void
tcmp2(double x,double y,int i)141 tcmp2 (double x, double y, int i)
142 {
143 mpfr_t xx, yy;
144 mpfr_prec_t j;
145
146 if (i == -1)
147 {
148 if (x == y)
149 i = 53;
150 else
151 i = (int) (__gmpfr_floor_log2 (x) - __gmpfr_floor_log2 (x - y));
152 }
153 mpfr_init2(xx, 53); mpfr_init2(yy, 53);
154 mpfr_set_d (xx, x, MPFR_RNDN);
155 mpfr_set_d (yy, y, MPFR_RNDN);
156 j = 0;
157 if (mpfr_cmp2 (xx, yy, &j) == 0)
158 {
159 if (x != y)
160 {
161 printf ("Error in mpfr_cmp2 for\nx=");
162 mpfr_out_str (stdout, 2, 0, xx, MPFR_RNDN);
163 printf ("\ny=");
164 mpfr_out_str (stdout, 2, 0, yy, MPFR_RNDN);
165 printf ("\ngot sign 0 for x != y\n");
166 exit (1);
167 }
168 }
169 else if (j != (unsigned) i)
170 {
171 printf ("Error in mpfr_cmp2 for\nx=");
172 mpfr_out_str (stdout, 2, 0, xx, MPFR_RNDN);
173 printf ("\ny=");
174 mpfr_out_str (stdout, 2, 0, yy, MPFR_RNDN);
175 printf ("\ngot %lu instead of %d\n", (unsigned long) j, i);
176 exit (1);
177 }
178 mpfr_clear(xx); mpfr_clear(yy);
179 }
180
181 static void
special(void)182 special (void)
183 {
184 mpfr_t x, y;
185 mpfr_prec_t j;
186
187 mpfr_init (x); mpfr_init (y);
188
189 /* bug found by Nathalie Revol, 21 March 2001 */
190 mpfr_set_prec (x, 65);
191 mpfr_set_prec (y, 65);
192 mpfr_set_str_binary (x, "0.10000000000000000000000000000000000001110010010110100110011110000E1");
193 mpfr_set_str_binary (y, "0.11100100101101001100111011111111110001101001000011101001001010010E-35");
194 j = 0;
195 if (mpfr_cmp2 (x, y, &j) <= 0 || j != 1)
196 {
197 printf ("Error in mpfr_cmp2:\n");
198 printf ("x=");
199 mpfr_dump (x);
200 printf ("y=");
201 mpfr_dump (y);
202 printf ("got %lu, expected 1\n", (unsigned long) j);
203 exit (1);
204 }
205
206 mpfr_set_prec(x, 127); mpfr_set_prec(y, 127);
207 mpfr_set_str_binary(x, "0.1011010000110111111000000101011110110001000101101011011110010010011110010000101101000010011001100110010000000010110000101000101E6");
208 mpfr_set_str_binary(y, "0.1011010000110111111000000101011011111100011101000011001111000010100010100110110100110010011001100110010000110010010110000010110E6");
209 j = 0;
210 if (mpfr_cmp2(x, y, &j) <= 0 || j != 32)
211 {
212 printf ("Error in mpfr_cmp2:\n");
213 printf ("x="); mpfr_dump (x);
214 printf ("y="); mpfr_dump (y);
215 printf ("got %lu, expected 32\n", (unsigned long) j);
216 exit (1);
217 }
218
219 mpfr_set_prec (x, 128); mpfr_set_prec (y, 239);
220 mpfr_set_str_binary (x, "0.10001000110110000111011000101011111100110010010011001101000011111010010110001000000010100110100111111011011010101100100000000000E167");
221 mpfr_set_str_binary (y, "0.10001000110110000111011000101011111100110010010011001101000011111010010110001000000010100110100111111011011010101100011111111111111111111111111111111111111111111111011111100101011100011001101000100111000010000000000101100110000111111000101E167");
222 j = 0;
223 if (mpfr_cmp2(x, y, &j) <= 0 || j != 164)
224 {
225 printf ("Error in mpfr_cmp2:\n");
226 printf ("x="); mpfr_dump (x);
227 printf ("y="); mpfr_dump (y);
228 printf ("got %lu, expected 164\n", (unsigned long) j);
229 exit (1);
230 }
231
232 /* bug found by Nathalie Revol, 29 March 2001 */
233 mpfr_set_prec (x, 130); mpfr_set_prec (y, 130);
234 mpfr_set_str_binary (x, "0.1100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E2");
235 mpfr_set_str_binary (y, "0.1011111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111100E2");
236 j = 0;
237 if (mpfr_cmp2(x, y, &j) <= 0 || j != 127)
238 {
239 printf ("Error in mpfr_cmp2:\n");
240 printf ("x="); mpfr_dump (x);
241 printf ("y="); mpfr_dump (y);
242 printf ("got %lu, expected 127\n", (unsigned long) j);
243 exit (1);
244 }
245
246 /* bug found by Nathalie Revol, 2 Apr 2001 */
247 mpfr_set_prec (x, 65); mpfr_set_prec (y, 65);
248 mpfr_set_ui (x, 5, MPFR_RNDN);
249 mpfr_set_str_binary (y, "0.10011111111111111111111111111111111111111111111111111111111111101E3");
250 j = 0;
251 if (mpfr_cmp2(x, y, &j) <= 0 || j != 63)
252 {
253 printf ("Error in mpfr_cmp2:\n");
254 printf ("x="); mpfr_dump (x);
255 printf ("y="); mpfr_dump (y);
256 printf ("got %lu, expected 63\n", (unsigned long) j);
257 exit (1);
258 }
259
260 /* bug found by Nathalie Revol, 2 Apr 2001 */
261 mpfr_set_prec (x, 65); mpfr_set_prec (y, 65);
262 mpfr_set_str_binary (x, "0.10011011111000101001110000000000000000000000000000000000000000000E-69");
263 mpfr_set_str_binary (y, "0.10011011111000101001101111111111111111111111111111111111111111101E-69");
264 j = 0;
265 if (mpfr_cmp2(x, y, &j) <= 0 || j != 63)
266 {
267 printf ("Error in mpfr_cmp2:\n");
268 printf ("x="); mpfr_dump (x);
269 printf ("y="); mpfr_dump (y);
270 printf ("got %lu, expected 63\n", (unsigned long) j);
271 exit (1);
272 }
273
274 mpfr_set_prec (x, 65);
275 mpfr_set_prec (y, 32);
276 mpfr_set_str_binary (x, "1.1110111011110001110111011111111111101000011001011100101100101101");
277 mpfr_set_str_binary (y, "0.11101110111100011101110111111111");
278 if (mpfr_cmp2 (x, y, &j) <= 0 || j != 0)
279 {
280 printf ("Error in mpfr_cmp2 (1)\n");
281 exit (1);
282 }
283
284 mpfr_set_prec (x, 2 * GMP_NUMB_BITS);
285 mpfr_set_prec (y, GMP_NUMB_BITS);
286 mpfr_set_ui (x, 1, MPFR_RNDN); /* x = 1 */
287 mpfr_set_ui (y, 1, MPFR_RNDN);
288 mpfr_nextbelow (y); /* y = 1 - 2^(-GMP_NUMB_BITS) */
289 mpfr_cmp2 (x, y, &j);
290 if (mpfr_cmp2 (x, y, &j) <= 0 || j != GMP_NUMB_BITS)
291 {
292 printf ("Error in mpfr_cmp2 (2)\n");
293 exit (1);
294 }
295
296 mpfr_clear (x);
297 mpfr_clear (y);
298 }
299
300 /* Compare (m,kx) and (m,ky), where (m,k) means m fixed limbs followed by
301 k zero limbs. */
302 static void
test_equal(void)303 test_equal (void)
304 {
305 mpfr_t w, x, y;
306 int m, kx, ky, inex;
307 mpfr_prec_t j;
308
309 for (m = 1; m <= 4; m++)
310 {
311 mpfr_init2 (w, m * GMP_NUMB_BITS);
312 for (kx = 0; kx <= 4; kx++)
313 for (ky = 0; ky <= 4; ky++)
314 {
315 do mpfr_urandomb (w, RANDS); while (mpfr_zero_p (w));
316 mpfr_init2 (x, (m + kx) * GMP_NUMB_BITS
317 - (kx == 0 ? 0 : randlimb () % GMP_NUMB_BITS));
318 mpfr_init2 (y, (m + ky) * GMP_NUMB_BITS
319 - (ky == 0 ? 0 : randlimb () % GMP_NUMB_BITS));
320 inex = mpfr_set (x, w, MPFR_RNDN);
321 MPFR_ASSERTN (inex == 0);
322 inex = mpfr_set (y, w, MPFR_RNDN);
323 MPFR_ASSERTN (inex == 0);
324 MPFR_ASSERTN (mpfr_equal_p (x, y));
325 if (RAND_BOOL ())
326 mpfr_neg (x, x, MPFR_RNDN);
327 if (RAND_BOOL ())
328 mpfr_neg (y, y, MPFR_RNDN);
329 if (mpfr_cmp2 (x, y, &j) != 0)
330 {
331 printf ("Error in test_equal for m = %d, kx = %d, ky = %d\n",
332 m, kx, ky);
333 printf (" x = ");
334 mpfr_dump (x);
335 printf (" y = ");
336 mpfr_dump (y);
337 exit (1);
338 }
339 mpfr_clears (x, y, (mpfr_ptr) 0);
340 }
341 mpfr_clear (w);
342 }
343 }
344
345 int
main(void)346 main (void)
347 {
348 int i;
349 long j;
350 double x, y, z;
351
352 tests_start_mpfr ();
353 mpfr_test_init ();
354
355 worst_cases ();
356 special ();
357 tcmp2 (5.43885304644369510000e+185, -1.87427265794105340000e-57, 1);
358 tcmp2 (1.06022698059744327881e+71, 1.05824655795525779205e+71, -1);
359 tcmp2 (1.0, 1.0, 53);
360 /* warning: cmp2 does not allow 0 as input */
361 for (x = 0.5, i = 1; i < 54; i++)
362 {
363 tcmp2 (1.0, 1.0-x, i);
364 x /= 2.0;
365 }
366 for (x = 0.5, i = 1; i < 100; i++)
367 {
368 tcmp2 (1.0, x, 1);
369 x /= 2.0;
370 }
371 for (j = 0; j < 100000; j++)
372 {
373 x = DBL_RAND ();
374 y = DBL_RAND ();
375 if (x < y)
376 {
377 z = x;
378 x = y;
379 y = z;
380 }
381 if (y != 0.0)
382 tcmp2 (x, y, -1);
383 }
384
385 test_equal ();
386
387 tests_end_mpfr ();
388
389 return 0;
390 }
391