1/* Test file for mpfr_cmp2. 2 3Copyright 1999, 2000, 2001, 2002, 2003, 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 <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 */ 29static void 30set_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 */ 49static void 50worst_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 141static void 142tcmp2 (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 182static void 183special (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 303int 304main (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