1 /* Test file for mpfr_cmp2. 2 3 Copyright 1999, 2000, 2001, 2002, 2003, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4 Contributed by the Arenaire and Cacao 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", 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", 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", l, expected); 117 exit (1); 118 } 119 set_bit (x, i + j + k + 2, 0); 120 set_bit (x, i + j + k + 3, 0); 121 set_bit (y, i + j + k + 3, 1); 122 l = 0; mpfr_cmp2 (x, y, &l); 123 expected = i + j + k + 2; 124 if (l != expected) 125 { 126 printf ("Error in mpfr_cmp2:\nx="); 127 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN); 128 printf ("\ny="); 129 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN); 130 printf ("\ngot %lu instead of %u\n", l, expected); 131 exit (1); 132 } 133 } 134 } 135 } 136 137 mpfr_clear (x); 138 mpfr_clear (y); 139 } 140 141 static void 142 tcmp2 (double x, double y, int i) 143 { 144 mpfr_t xx, yy; 145 mpfr_prec_t j; 146 147 if (i == -1) 148 { 149 if (x == y) 150 i = 53; 151 else 152 i = (int) (__gmpfr_floor_log2 (x) - __gmpfr_floor_log2 (x - y)); 153 } 154 mpfr_init2(xx, 53); mpfr_init2(yy, 53); 155 mpfr_set_d (xx, x, MPFR_RNDN); 156 mpfr_set_d (yy, y, MPFR_RNDN); 157 j = 0; 158 if (mpfr_cmp2 (xx, yy, &j) == 0) 159 { 160 if (x != y) 161 { 162 printf ("Error in mpfr_cmp2 for\nx="); 163 mpfr_out_str (stdout, 2, 0, xx, MPFR_RNDN); 164 printf ("\ny="); 165 mpfr_out_str (stdout, 2, 0, yy, MPFR_RNDN); 166 printf ("\ngot sign 0 for x != y\n"); 167 exit (1); 168 } 169 } 170 else if (j != (unsigned) i) 171 { 172 printf ("Error in mpfr_cmp2 for\nx="); 173 mpfr_out_str (stdout, 2, 0, xx, MPFR_RNDN); 174 printf ("\ny="); 175 mpfr_out_str (stdout, 2, 0, yy, MPFR_RNDN); 176 printf ("\ngot %lu instead of %d\n", j, i); 177 exit (1); 178 } 179 mpfr_clear(xx); mpfr_clear(yy); 180 } 181 182 static void 183 special (void) 184 { 185 mpfr_t x, y; 186 mpfr_prec_t j; 187 188 mpfr_init (x); mpfr_init (y); 189 190 /* bug found by Nathalie Revol, 21 March 2001 */ 191 mpfr_set_prec (x, 65); 192 mpfr_set_prec (y, 65); 193 mpfr_set_str_binary (x, "0.10000000000000000000000000000000000001110010010110100110011110000E1"); 194 mpfr_set_str_binary (y, "0.11100100101101001100111011111111110001101001000011101001001010010E-35"); 195 j = 0; 196 if (mpfr_cmp2 (x, y, &j) <= 0 || j != 1) 197 { 198 printf ("Error in mpfr_cmp2:\n"); 199 printf ("x="); 200 mpfr_print_binary (x); 201 puts (""); 202 printf ("y="); 203 mpfr_print_binary (y); 204 puts (""); 205 printf ("got %lu, expected 1\n", j); 206 exit (1); 207 } 208 209 mpfr_set_prec(x, 127); mpfr_set_prec(y, 127); 210 mpfr_set_str_binary(x, "0.1011010000110111111000000101011110110001000101101011011110010010011110010000101101000010011001100110010000000010110000101000101E6"); 211 mpfr_set_str_binary(y, "0.1011010000110111111000000101011011111100011101000011001111000010100010100110110100110010011001100110010000110010010110000010110E6"); 212 j = 0; 213 if (mpfr_cmp2(x, y, &j) <= 0 || j != 32) 214 { 215 printf ("Error in mpfr_cmp2:\n"); 216 printf ("x="); mpfr_print_binary(x); puts (""); 217 printf ("y="); mpfr_print_binary(y); puts (""); 218 printf ("got %lu, expected 32\n", j); 219 exit (1); 220 } 221 222 mpfr_set_prec (x, 128); mpfr_set_prec (y, 239); 223 mpfr_set_str_binary (x, "0.10001000110110000111011000101011111100110010010011001101000011111010010110001000000010100110100111111011011010101100100000000000E167"); 224 mpfr_set_str_binary (y, "0.10001000110110000111011000101011111100110010010011001101000011111010010110001000000010100110100111111011011010101100011111111111111111111111111111111111111111111111011111100101011100011001101000100111000010000000000101100110000111111000101E167"); 225 j = 0; 226 if (mpfr_cmp2(x, y, &j) <= 0 || j != 164) 227 { 228 printf ("Error in mpfr_cmp2:\n"); 229 printf ("x="); mpfr_print_binary(x); puts (""); 230 printf ("y="); mpfr_print_binary(y); puts (""); 231 printf ("got %lu, expected 164\n", j); 232 exit (1); 233 } 234 235 /* bug found by Nathalie Revol, 29 March 2001 */ 236 mpfr_set_prec (x, 130); mpfr_set_prec (y, 130); 237 mpfr_set_str_binary (x, "0.1100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E2"); 238 mpfr_set_str_binary (y, "0.1011111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111100E2"); 239 j = 0; 240 if (mpfr_cmp2(x, y, &j) <= 0 || j != 127) 241 { 242 printf ("Error in mpfr_cmp2:\n"); 243 printf ("x="); mpfr_print_binary(x); puts (""); 244 printf ("y="); mpfr_print_binary(y); puts (""); 245 printf ("got %lu, expected 127\n", j); 246 exit (1); 247 } 248 249 /* bug found by Nathalie Revol, 2 Apr 2001 */ 250 mpfr_set_prec (x, 65); mpfr_set_prec (y, 65); 251 mpfr_set_ui (x, 5, MPFR_RNDN); 252 mpfr_set_str_binary (y, "0.10011111111111111111111111111111111111111111111111111111111111101E3"); 253 j = 0; 254 if (mpfr_cmp2(x, y, &j) <= 0 || j != 63) 255 { 256 printf ("Error in mpfr_cmp2:\n"); 257 printf ("x="); mpfr_print_binary(x); puts (""); 258 printf ("y="); mpfr_print_binary(y); puts (""); 259 printf ("got %lu, expected 63\n", j); 260 exit (1); 261 } 262 263 /* bug found by Nathalie Revol, 2 Apr 2001 */ 264 mpfr_set_prec (x, 65); mpfr_set_prec (y, 65); 265 mpfr_set_str_binary (x, "0.10011011111000101001110000000000000000000000000000000000000000000E-69"); 266 mpfr_set_str_binary (y, "0.10011011111000101001101111111111111111111111111111111111111111101E-69"); 267 j = 0; 268 if (mpfr_cmp2(x, y, &j) <= 0 || j != 63) 269 { 270 printf ("Error in mpfr_cmp2:\n"); 271 printf ("x="); mpfr_print_binary(x); puts (""); 272 printf ("y="); mpfr_print_binary(y); puts (""); 273 printf ("got %lu, expected 63\n", j); 274 exit (1); 275 } 276 277 mpfr_set_prec (x, 65); 278 mpfr_set_prec (y, 32); 279 mpfr_set_str_binary (x, "1.1110111011110001110111011111111111101000011001011100101100101101"); 280 mpfr_set_str_binary (y, "0.11101110111100011101110111111111"); 281 if (mpfr_cmp2 (x, y, &j) <= 0 || j != 0) 282 { 283 printf ("Error in mpfr_cmp2 (1)\n"); 284 exit (1); 285 } 286 287 mpfr_set_prec (x, 2 * GMP_NUMB_BITS); 288 mpfr_set_prec (y, GMP_NUMB_BITS); 289 mpfr_set_ui (x, 1, MPFR_RNDN); /* x = 1 */ 290 mpfr_set_ui (y, 1, MPFR_RNDN); 291 mpfr_nextbelow (y); /* y = 1 - 2^(-GMP_NUMB_BITS) */ 292 mpfr_cmp2 (x, y, &j); 293 if (mpfr_cmp2 (x, y, &j) <= 0 || j != GMP_NUMB_BITS) 294 { 295 printf ("Error in mpfr_cmp2 (2)\n"); 296 exit (1); 297 } 298 299 mpfr_clear (x); 300 mpfr_clear (y); 301 } 302 303 int 304 main (void) 305 { 306 int i; 307 long j; 308 double x, y, z; 309 310 tests_start_mpfr (); 311 mpfr_test_init (); 312 313 worst_cases (); 314 special (); 315 tcmp2 (5.43885304644369510000e+185, -1.87427265794105340000e-57, 1); 316 tcmp2 (1.06022698059744327881e+71, 1.05824655795525779205e+71, -1); 317 tcmp2 (1.0, 1.0, 53); 318 /* warning: cmp2 does not allow 0 as input */ 319 for (x = 0.5, i = 1; i < 54; i++) 320 { 321 tcmp2 (1.0, 1.0-x, i); 322 x /= 2.0; 323 } 324 for (x = 0.5, i = 1; i < 100; i++) 325 { 326 tcmp2 (1.0, x, 1); 327 x /= 2.0; 328 } 329 for (j = 0; j < 100000; j++) 330 { 331 x = DBL_RAND (); 332 y = DBL_RAND (); 333 if (x < y) 334 { 335 z = x; 336 x = y; 337 y = z; 338 } 339 if (y != 0.0) 340 tcmp2 (x, y, -1); 341 } 342 343 tests_end_mpfr (); 344 345 return 0; 346 } 347