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