1 1.1.1.2 mrg /* Test file for mpfr_add_[q,z], mpfr_sub_[q,z], mpfr_div_[q,z], 2 1.1.1.2 mrg mpfr_mul_[q,z], mpfr_cmp_[f,q,z] 3 1.1 mrg 4 1.1.1.6 mrg Copyright 2004-2023 Free Software Foundation, Inc. 5 1.1.1.3 mrg Contributed by the AriC and Caramba projects, INRIA. 6 1.1 mrg 7 1.1 mrg This file is part of the GNU MPFR Library. 8 1.1 mrg 9 1.1 mrg The GNU MPFR Library is free software; you can redistribute it and/or modify 10 1.1 mrg it under the terms of the GNU Lesser General Public License as published by 11 1.1 mrg the Free Software Foundation; either version 3 of the License, or (at your 12 1.1 mrg option) any later version. 13 1.1 mrg 14 1.1 mrg The GNU MPFR Library is distributed in the hope that it will be useful, but 15 1.1 mrg WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 16 1.1 mrg or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 17 1.1 mrg License for more details. 18 1.1 mrg 19 1.1 mrg You should have received a copy of the GNU Lesser General Public License 20 1.1 mrg along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 21 1.1.1.5 mrg https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 22 1.1 mrg 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 23 1.1 mrg 24 1.1 mrg #include "mpfr-test.h" 25 1.1 mrg 26 1.1.1.4 mrg #ifndef MPFR_USE_MINI_GMP 27 1.1.1.4 mrg 28 1.1.1.2 mrg #define CHECK_FOR(str, cond) \ 29 1.1.1.2 mrg if ((cond) == 0) { \ 30 1.1.1.2 mrg printf ("Special case error %s. Ternary value = %d, flags = %u\n", \ 31 1.1.1.2 mrg str, res, __gmpfr_flags); \ 32 1.1.1.2 mrg printf ("Got "); mpfr_dump (y); \ 33 1.1.1.2 mrg printf ("X = "); mpfr_dump (x); \ 34 1.1.1.4 mrg printf ("Q = "); mpz_out_str (stdout, 10, mpq_numref(q)); \ 35 1.1.1.4 mrg printf ("\n /"); mpz_out_str (stdout, 10, mpq_denref(q)); \ 36 1.1.1.4 mrg printf ("\n"); \ 37 1.1.1.2 mrg exit (1); \ 38 1.1.1.2 mrg } 39 1.1.1.2 mrg 40 1.1.1.2 mrg #define CHECK_FORZ(str, cond) \ 41 1.1.1.2 mrg if ((cond) == 0) { \ 42 1.1.1.2 mrg printf ("Special case error %s. Ternary value = %d, flags = %u\n", \ 43 1.1.1.2 mrg str, res, __gmpfr_flags); \ 44 1.1.1.2 mrg printf ("Got "); mpfr_dump (y); \ 45 1.1.1.2 mrg printf ("X = "); mpfr_dump (x); \ 46 1.1.1.4 mrg printf ("Z = "); mpz_out_str (stdout, 10, z); \ 47 1.1.1.4 mrg printf ("\n"); \ 48 1.1.1.2 mrg exit (1); \ 49 1.1.1.2 mrg } 50 1.1 mrg 51 1.1 mrg static void 52 1.1 mrg special (void) 53 1.1 mrg { 54 1.1 mrg mpfr_t x, y; 55 1.1.1.2 mrg mpq_t q; 56 1.1.1.2 mrg mpz_t z; 57 1.1 mrg int res = 0; 58 1.1 mrg 59 1.1 mrg mpfr_init (x); 60 1.1 mrg mpfr_init (y); 61 1.1.1.2 mrg mpq_init (q); 62 1.1.1.2 mrg mpz_init (z); 63 1.1 mrg 64 1.1 mrg /* cancellation in mpfr_add_q */ 65 1.1 mrg mpfr_set_prec (x, 60); 66 1.1 mrg mpfr_set_prec (y, 20); 67 1.1.1.2 mrg mpz_set_str (mpq_numref (q), "-187207494", 10); 68 1.1.1.2 mrg mpz_set_str (mpq_denref (q), "5721", 10); 69 1.1 mrg mpfr_set_str_binary (x, "11111111101001011011100101100011011110010011100010000100001E-44"); 70 1.1.1.2 mrg mpfr_add_q (y, x, q, MPFR_RNDN); 71 1.1.1.5 mrg CHECK_FOR ("cancellation in add_q", mpfr_cmp_ui_2exp (y, 256783, -64) == 0); 72 1.1 mrg 73 1.1 mrg mpfr_set_prec (x, 19); 74 1.1 mrg mpfr_set_str_binary (x, "0.1011110101110011100E0"); 75 1.1.1.2 mrg mpz_set_str (mpq_numref (q), "187207494", 10); 76 1.1.1.2 mrg mpz_set_str (mpq_denref (q), "5721", 10); 77 1.1 mrg mpfr_set_prec (y, 29); 78 1.1.1.2 mrg mpfr_add_q (y, x, q, MPFR_RNDD); 79 1.1 mrg mpfr_set_prec (x, 29); 80 1.1 mrg mpfr_set_str_binary (x, "11111111101001110011010001001E-14"); 81 1.1.1.5 mrg CHECK_FOR ("cancellation in add_q", mpfr_cmp (x,y) == 0); 82 1.1 mrg 83 1.1 mrg /* Inf */ 84 1.1 mrg mpfr_set_inf (x, 1); 85 1.1.1.2 mrg mpz_set_str (mpq_numref (q), "395877315", 10); 86 1.1.1.2 mrg mpz_set_str (mpq_denref (q), "3508975966", 10); 87 1.1 mrg mpfr_set_prec (y, 118); 88 1.1.1.2 mrg mpfr_add_q (y, x, q, MPFR_RNDU); 89 1.1 mrg CHECK_FOR ("inf", mpfr_inf_p (y) && mpfr_sgn (y) > 0); 90 1.1.1.2 mrg mpfr_sub_q (y, x, q, MPFR_RNDU); 91 1.1 mrg CHECK_FOR ("inf", mpfr_inf_p (y) && mpfr_sgn (y) > 0); 92 1.1 mrg 93 1.1 mrg /* Nan */ 94 1.1 mrg MPFR_SET_NAN (x); 95 1.1.1.2 mrg mpfr_add_q (y, x, q, MPFR_RNDU); 96 1.1 mrg CHECK_FOR ("nan", mpfr_nan_p (y)); 97 1.1.1.2 mrg mpfr_sub_q (y, x, q, MPFR_RNDU); 98 1.1 mrg CHECK_FOR ("nan", mpfr_nan_p (y)); 99 1.1 mrg 100 1.1 mrg /* Exact value */ 101 1.1 mrg mpfr_set_prec (x, 60); 102 1.1 mrg mpfr_set_prec (y, 60); 103 1.1 mrg mpfr_set_str1 (x, "0.5"); 104 1.1.1.2 mrg mpz_set_str (mpq_numref (q), "3", 10); 105 1.1.1.2 mrg mpz_set_str (mpq_denref (q), "2", 10); 106 1.1.1.2 mrg res = mpfr_add_q (y, x, q, MPFR_RNDU); 107 1.1 mrg CHECK_FOR ("0.5+3/2", mpfr_cmp_ui(y, 2)==0 && res==0); 108 1.1.1.2 mrg res = mpfr_sub_q (y, x, q, MPFR_RNDU); 109 1.1 mrg CHECK_FOR ("0.5-3/2", mpfr_cmp_si(y, -1)==0 && res==0); 110 1.1 mrg 111 1.1.1.6 mrg /* Inf rational */ 112 1.1.1.2 mrg mpq_set_ui (q, 1, 0); 113 1.1 mrg mpfr_set_str1 (x, "0.5"); 114 1.1.1.2 mrg res = mpfr_add_q (y, x, q, MPFR_RNDN); 115 1.1.1.4 mrg CHECK_FOR ("0.5+1/0", mpfr_inf_p (y) && MPFR_IS_POS (y) && res == 0); 116 1.1.1.2 mrg res = mpfr_sub_q (y, x, q, MPFR_RNDN); 117 1.1.1.4 mrg CHECK_FOR ("0.5-1/0", mpfr_inf_p (y) && MPFR_IS_NEG (y) && res == 0); 118 1.1.1.2 mrg mpq_set_si (q, -1, 0); 119 1.1.1.2 mrg res = mpfr_add_q (y, x, q, MPFR_RNDN); 120 1.1.1.4 mrg CHECK_FOR ("0.5+ -1/0", mpfr_inf_p (y) && MPFR_IS_NEG (y) && res == 0); 121 1.1.1.2 mrg res = mpfr_sub_q (y, x, q, MPFR_RNDN); 122 1.1.1.4 mrg CHECK_FOR ("0.5- -1/0", mpfr_inf_p (y) && MPFR_IS_POS (y) && res == 0); 123 1.1.1.2 mrg res = mpfr_div_q (y, x, q, MPFR_RNDN); 124 1.1.1.4 mrg CHECK_FOR ("0.5 / (-1/0)", mpfr_zero_p (y) && MPFR_IS_NEG (y) && res == 0); 125 1.1.1.2 mrg mpq_set_ui (q, 1, 0); 126 1.1.1.2 mrg mpfr_set_inf (x, 1); 127 1.1.1.2 mrg res = mpfr_add_q (y, x, q, MPFR_RNDN); 128 1.1.1.4 mrg CHECK_FOR ("+Inf + +Inf", mpfr_inf_p (y) && MPFR_IS_POS (y) && res == 0); 129 1.1.1.2 mrg res = mpfr_sub_q (y, x, q, MPFR_RNDN); 130 1.1.1.2 mrg CHECK_FOR ("+Inf - +Inf", MPFR_IS_NAN (y) && res == 0); 131 1.1.1.2 mrg mpfr_set_inf (x, -1); 132 1.1.1.2 mrg res = mpfr_add_q (y, x, q, MPFR_RNDN); 133 1.1.1.2 mrg CHECK_FOR ("-Inf + +Inf", MPFR_IS_NAN (y) && res == 0); 134 1.1.1.2 mrg res = mpfr_sub_q (y, x, q, MPFR_RNDN); 135 1.1.1.4 mrg CHECK_FOR ("-Inf - +Inf", mpfr_inf_p (y) && MPFR_IS_NEG (y) && res == 0); 136 1.1.1.2 mrg mpq_set_si (q, -1, 0); 137 1.1.1.2 mrg mpfr_set_inf (x, 1); 138 1.1.1.2 mrg res = mpfr_add_q (y, x, q, MPFR_RNDN); 139 1.1.1.2 mrg CHECK_FOR ("+Inf + -Inf", MPFR_IS_NAN (y) && res == 0); 140 1.1.1.2 mrg res = mpfr_sub_q (y, x, q, MPFR_RNDN); 141 1.1.1.4 mrg CHECK_FOR ("+Inf - -Inf", mpfr_inf_p (y) && MPFR_IS_POS (y) && res == 0); 142 1.1.1.2 mrg mpfr_set_inf (x, -1); 143 1.1.1.2 mrg res = mpfr_add_q (y, x, q, MPFR_RNDN); 144 1.1.1.4 mrg CHECK_FOR ("-Inf + -Inf", mpfr_inf_p (y) && MPFR_IS_NEG (y) && res == 0); 145 1.1.1.2 mrg res = mpfr_sub_q (y, x, q, MPFR_RNDN); 146 1.1.1.2 mrg CHECK_FOR ("-Inf - -Inf", MPFR_IS_NAN (y) && res == 0); 147 1.1 mrg 148 1.1 mrg /* 0 */ 149 1.1.1.2 mrg mpq_set_ui (q, 0, 1); 150 1.1 mrg mpfr_set_ui (x, 42, MPFR_RNDN); 151 1.1.1.2 mrg res = mpfr_add_q (y, x, q, MPFR_RNDN); 152 1.1 mrg CHECK_FOR ("42+0/1", mpfr_cmp_ui (y, 42) == 0 && res == 0); 153 1.1.1.2 mrg res = mpfr_sub_q (y, x, q, MPFR_RNDN); 154 1.1 mrg CHECK_FOR ("42-0/1", mpfr_cmp_ui (y, 42) == 0 && res == 0); 155 1.1.1.2 mrg res = mpfr_mul_q (y, x, q, MPFR_RNDN); 156 1.1.1.4 mrg CHECK_FOR ("42*0/1", mpfr_zero_p (y) && MPFR_IS_POS (y) && res == 0); 157 1.1.1.2 mrg mpfr_clear_flags (); 158 1.1.1.2 mrg res = mpfr_div_q (y, x, q, MPFR_RNDN); 159 1.1.1.4 mrg CHECK_FOR ("42/(0/1)", mpfr_inf_p (y) && MPFR_IS_POS (y) && res == 0 160 1.1.1.2 mrg && mpfr_divby0_p ()); 161 1.1.1.2 mrg mpz_set_ui (z, 0); 162 1.1.1.2 mrg mpfr_clear_flags (); 163 1.1.1.2 mrg res = mpfr_div_z (y, x, z, MPFR_RNDN); 164 1.1.1.4 mrg CHECK_FORZ ("42/0", mpfr_inf_p (y) && MPFR_IS_POS (y) && res == 0 165 1.1.1.2 mrg && mpfr_divby0_p ()); 166 1.1 mrg 167 1.1.1.2 mrg mpz_clear (z); 168 1.1.1.2 mrg mpq_clear (q); 169 1.1 mrg mpfr_clear (x); 170 1.1 mrg mpfr_clear (y); 171 1.1 mrg } 172 1.1 mrg 173 1.1 mrg static void 174 1.1 mrg check_for_zero (void) 175 1.1 mrg { 176 1.1 mrg /* Check that 0 is unsigned! */ 177 1.1 mrg mpq_t q; 178 1.1 mrg mpz_t z; 179 1.1 mrg mpfr_t x; 180 1.1 mrg int r; 181 1.1 mrg mpfr_sign_t i; 182 1.1 mrg 183 1.1 mrg mpfr_init (x); 184 1.1 mrg mpz_init (z); 185 1.1 mrg mpq_init (q); 186 1.1 mrg 187 1.1 mrg mpz_set_ui (z, 0); 188 1.1 mrg mpq_set_ui (q, 0, 1); 189 1.1 mrg 190 1.1 mrg MPFR_SET_ZERO (x); 191 1.1 mrg RND_LOOP (r) 192 1.1 mrg { 193 1.1 mrg for (i = MPFR_SIGN_NEG ; i <= MPFR_SIGN_POS ; 194 1.1 mrg i+=MPFR_SIGN_POS-MPFR_SIGN_NEG) 195 1.1 mrg { 196 1.1 mrg MPFR_SET_SIGN(x, i); 197 1.1 mrg mpfr_add_z (x, x, z, (mpfr_rnd_t) r); 198 1.1 mrg if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i) 199 1.1 mrg { 200 1.1 mrg printf("GMP Zero errors for add_z & rnd=%s & s=%d\n", 201 1.1 mrg mpfr_print_rnd_mode ((mpfr_rnd_t) r), i); 202 1.1 mrg mpfr_dump (x); 203 1.1 mrg exit (1); 204 1.1 mrg } 205 1.1 mrg mpfr_sub_z (x, x, z, (mpfr_rnd_t) r); 206 1.1 mrg if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i) 207 1.1 mrg { 208 1.1 mrg printf("GMP Zero errors for sub_z & rnd=%s & s=%d\n", 209 1.1 mrg mpfr_print_rnd_mode ((mpfr_rnd_t) r), i); 210 1.1 mrg mpfr_dump (x); 211 1.1 mrg exit (1); 212 1.1 mrg } 213 1.1 mrg mpfr_mul_z (x, x, z, (mpfr_rnd_t) r); 214 1.1 mrg if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i) 215 1.1 mrg { 216 1.1 mrg printf("GMP Zero errors for mul_z & rnd=%s & s=%d\n", 217 1.1 mrg mpfr_print_rnd_mode ((mpfr_rnd_t) r), i); 218 1.1 mrg mpfr_dump (x); 219 1.1 mrg exit (1); 220 1.1 mrg } 221 1.1 mrg mpfr_add_q (x, x, q, (mpfr_rnd_t) r); 222 1.1 mrg if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i) 223 1.1 mrg { 224 1.1 mrg printf("GMP Zero errors for add_q & rnd=%s & s=%d\n", 225 1.1 mrg mpfr_print_rnd_mode ((mpfr_rnd_t) r), i); 226 1.1 mrg mpfr_dump (x); 227 1.1 mrg exit (1); 228 1.1 mrg } 229 1.1 mrg mpfr_sub_q (x, x, q, (mpfr_rnd_t) r); 230 1.1 mrg if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i) 231 1.1 mrg { 232 1.1 mrg printf("GMP Zero errors for sub_q & rnd=%s & s=%d\n", 233 1.1 mrg mpfr_print_rnd_mode ((mpfr_rnd_t) r), i); 234 1.1 mrg mpfr_dump (x); 235 1.1 mrg exit (1); 236 1.1 mrg } 237 1.1 mrg } 238 1.1 mrg } 239 1.1 mrg 240 1.1 mrg mpq_clear (q); 241 1.1 mrg mpz_clear (z); 242 1.1 mrg mpfr_clear (x); 243 1.1 mrg } 244 1.1 mrg 245 1.1 mrg static void 246 1.1 mrg test_cmp_z (mpfr_prec_t pmin, mpfr_prec_t pmax, int nmax) 247 1.1 mrg { 248 1.1 mrg mpfr_t x, z; 249 1.1 mrg mpz_t y; 250 1.1 mrg mpfr_prec_t p; 251 1.1 mrg int res1, res2; 252 1.1 mrg int n; 253 1.1 mrg 254 1.1 mrg mpfr_init (x); 255 1.1 mrg mpfr_init2 (z, MPFR_PREC_MIN); 256 1.1 mrg mpz_init (y); 257 1.1.1.2 mrg 258 1.1.1.2 mrg /* check the erange flag when x is NaN */ 259 1.1.1.2 mrg mpfr_set_nan (x); 260 1.1.1.2 mrg mpz_set_ui (y, 17); 261 1.1.1.2 mrg mpfr_clear_erangeflag (); 262 1.1.1.2 mrg res1 = mpfr_cmp_z (x, y); 263 1.1.1.2 mrg if (res1 != 0 || mpfr_erangeflag_p () == 0) 264 1.1.1.2 mrg { 265 1.1.1.2 mrg printf ("Error for mpfr_cmp_z (NaN, 17)\n"); 266 1.1.1.2 mrg printf ("Return value: expected 0, got %d\n", res1); 267 1.1.1.2 mrg printf ("Erange flag: expected set, got %d\n", mpfr_erangeflag_p ()); 268 1.1.1.2 mrg exit (1); 269 1.1.1.2 mrg } 270 1.1.1.2 mrg 271 1.1.1.5 mrg for (p = pmin ; p < pmax ; p++) 272 1.1 mrg { 273 1.1 mrg mpfr_set_prec (x, p); 274 1.1.1.5 mrg for (n = 0 ; n < nmax ; n++) 275 1.1 mrg { 276 1.1 mrg mpfr_urandomb (x, RANDS); 277 1.1 mrg mpz_urandomb (y, RANDS, 1024); 278 1.1 mrg if (!MPFR_IS_SINGULAR (x)) 279 1.1 mrg { 280 1.1 mrg mpfr_sub_z (z, x, y, MPFR_RNDN); 281 1.1.1.4 mrg res1 = (mpfr_sgn) (z); 282 1.1 mrg res2 = mpfr_cmp_z (x, y); 283 1.1 mrg if (res1 != res2) 284 1.1 mrg { 285 1.1 mrg printf("Error for mpfr_cmp_z: res=%d sub_z gives %d\n", 286 1.1 mrg res2, res1); 287 1.1 mrg exit (1); 288 1.1 mrg } 289 1.1 mrg } 290 1.1 mrg } 291 1.1 mrg } 292 1.1 mrg mpz_clear (y); 293 1.1 mrg mpfr_clear (x); 294 1.1 mrg mpfr_clear (z); 295 1.1 mrg } 296 1.1 mrg 297 1.1 mrg static void 298 1.1 mrg test_cmp_q (mpfr_prec_t pmin, mpfr_prec_t pmax, int nmax) 299 1.1 mrg { 300 1.1 mrg mpfr_t x, z; 301 1.1 mrg mpq_t y; 302 1.1 mrg mpfr_prec_t p; 303 1.1 mrg int res1, res2; 304 1.1 mrg int n; 305 1.1 mrg 306 1.1 mrg mpfr_init (x); 307 1.1 mrg mpfr_init2 (z, MPFR_PREC_MIN); 308 1.1 mrg mpq_init (y); 309 1.1.1.2 mrg 310 1.1.1.5 mrg /* Check the flags when x is NaN: the erange flags must be set, and 311 1.1.1.5 mrg only this one. */ 312 1.1.1.2 mrg mpfr_set_nan (x); 313 1.1.1.2 mrg mpq_set_ui (y, 17, 1); 314 1.1.1.5 mrg mpfr_clear_flags (); 315 1.1.1.2 mrg res1 = mpfr_cmp_q (x, y); 316 1.1.1.5 mrg if (res1 != 0 || __gmpfr_flags != MPFR_FLAGS_ERANGE) 317 1.1.1.2 mrg { 318 1.1.1.2 mrg printf ("Error for mpfr_cmp_q (NaN, 17)\n"); 319 1.1.1.2 mrg printf ("Return value: expected 0, got %d\n", res1); 320 1.1.1.5 mrg printf ("Expected flags:"); 321 1.1.1.5 mrg flags_out (MPFR_FLAGS_ERANGE); 322 1.1.1.5 mrg printf ("Got flags: "); 323 1.1.1.5 mrg flags_out (__gmpfr_flags); 324 1.1.1.5 mrg exit (1); 325 1.1.1.5 mrg } 326 1.1.1.5 mrg 327 1.1.1.5 mrg /* Check the flags when y is NaN: the erange flags must be set, and 328 1.1.1.5 mrg only this one. */ 329 1.1.1.5 mrg mpfr_set_ui (x, 42, MPFR_RNDN); 330 1.1.1.5 mrg /* A NaN rational is represented by 0/0 (MPFR extension). */ 331 1.1.1.5 mrg mpz_set_ui (mpq_numref (y), 0); 332 1.1.1.5 mrg mpz_set_ui (mpq_denref (y), 0); 333 1.1.1.5 mrg mpfr_clear_flags (); 334 1.1.1.5 mrg res1 = mpfr_cmp_q (x, y); 335 1.1.1.5 mrg if (res1 != 0 || __gmpfr_flags != MPFR_FLAGS_ERANGE) 336 1.1.1.5 mrg { 337 1.1.1.5 mrg printf ("Error for mpfr_cmp_q (42, NaN)\n"); 338 1.1.1.5 mrg printf ("Return value: expected 0, got %d\n", res1); 339 1.1.1.5 mrg printf ("Expected flags:"); 340 1.1.1.5 mrg flags_out (MPFR_FLAGS_ERANGE); 341 1.1.1.5 mrg printf ("Got flags: "); 342 1.1.1.5 mrg flags_out (__gmpfr_flags); 343 1.1.1.2 mrg exit (1); 344 1.1.1.2 mrg } 345 1.1.1.2 mrg 346 1.1.1.5 mrg for (p = pmin ; p < pmax ; p++) 347 1.1 mrg { 348 1.1 mrg mpfr_set_prec (x, p); 349 1.1 mrg for (n = 0 ; n < nmax ; n++) 350 1.1 mrg { 351 1.1 mrg mpfr_urandomb (x, RANDS); 352 1.1 mrg mpq_set_ui (y, randlimb (), randlimb() ); 353 1.1 mrg if (!MPFR_IS_SINGULAR (x)) 354 1.1 mrg { 355 1.1 mrg mpfr_sub_q (z, x, y, MPFR_RNDN); 356 1.1.1.4 mrg res1 = (mpfr_sgn) (z); 357 1.1 mrg res2 = mpfr_cmp_q (x, y); 358 1.1 mrg if (res1 != res2) 359 1.1 mrg { 360 1.1 mrg printf("Error for mpfr_cmp_q: res=%d sub_z gives %d\n", 361 1.1 mrg res2, res1); 362 1.1 mrg exit (1); 363 1.1 mrg } 364 1.1 mrg } 365 1.1 mrg } 366 1.1 mrg } 367 1.1.1.5 mrg 368 1.1.1.5 mrg /* check for y = 1/0 */ 369 1.1.1.5 mrg mpz_set_ui (mpq_numref (y), 1); 370 1.1.1.5 mrg mpz_set_ui (mpq_denref (y), 0); 371 1.1.1.5 mrg mpfr_set_ui (x, 1, MPFR_RNDN); 372 1.1.1.5 mrg MPFR_ASSERTN(mpfr_cmp_q (x, y) < 0); 373 1.1.1.5 mrg mpfr_set_inf (x, -1); 374 1.1.1.5 mrg MPFR_ASSERTN(mpfr_cmp_q (x, y) < 0); 375 1.1.1.5 mrg mpfr_set_inf (x, +1); 376 1.1.1.5 mrg MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0); 377 1.1.1.5 mrg mpfr_set_nan (x); 378 1.1.1.5 mrg mpfr_clear_erangeflag (); 379 1.1.1.5 mrg MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0); 380 1.1.1.5 mrg MPFR_ASSERTN(mpfr_erangeflag_p ()); 381 1.1.1.5 mrg 382 1.1.1.5 mrg /* check for y = -1/0 */ 383 1.1.1.5 mrg mpz_set_si (mpq_numref (y), -1); 384 1.1.1.5 mrg mpz_set_ui (mpq_denref (y), 0); 385 1.1.1.5 mrg mpfr_set_ui (x, 1, MPFR_RNDN); 386 1.1.1.5 mrg MPFR_ASSERTN(mpfr_cmp_q (x, y) > 0); 387 1.1.1.5 mrg mpfr_set_inf (x, -1); 388 1.1.1.5 mrg MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0); 389 1.1.1.5 mrg mpfr_set_inf (x, +1); 390 1.1.1.5 mrg MPFR_ASSERTN(mpfr_cmp_q (x, y) > 0); 391 1.1.1.5 mrg mpfr_set_nan (x); 392 1.1.1.5 mrg mpfr_clear_erangeflag (); 393 1.1.1.5 mrg MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0); 394 1.1.1.5 mrg MPFR_ASSERTN(mpfr_erangeflag_p ()); 395 1.1.1.5 mrg 396 1.1.1.5 mrg /* check for y = 0/0 */ 397 1.1.1.5 mrg mpz_set_ui (mpq_numref (y), 0); 398 1.1.1.5 mrg mpz_set_ui (mpq_denref (y), 0); 399 1.1.1.5 mrg mpfr_set_ui (x, 1, MPFR_RNDN); 400 1.1.1.5 mrg mpfr_clear_erangeflag (); 401 1.1.1.5 mrg MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0); 402 1.1.1.5 mrg MPFR_ASSERTN(mpfr_erangeflag_p ()); 403 1.1.1.5 mrg mpfr_set_inf (x, -1); 404 1.1.1.5 mrg mpfr_clear_erangeflag (); 405 1.1.1.5 mrg MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0); 406 1.1.1.5 mrg MPFR_ASSERTN(mpfr_erangeflag_p ()); 407 1.1.1.5 mrg mpfr_set_inf (x, +1); 408 1.1.1.5 mrg mpfr_clear_erangeflag (); 409 1.1.1.5 mrg MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0); 410 1.1.1.5 mrg MPFR_ASSERTN(mpfr_erangeflag_p ()); 411 1.1.1.5 mrg mpfr_set_nan (x); 412 1.1.1.5 mrg mpfr_clear_erangeflag (); 413 1.1.1.5 mrg MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0); 414 1.1.1.5 mrg MPFR_ASSERTN(mpfr_erangeflag_p ()); 415 1.1.1.5 mrg 416 1.1 mrg mpq_clear (y); 417 1.1 mrg mpfr_clear (x); 418 1.1 mrg mpfr_clear (z); 419 1.1 mrg } 420 1.1 mrg 421 1.1 mrg static void 422 1.1 mrg test_cmp_f (mpfr_prec_t pmin, mpfr_prec_t pmax, int nmax) 423 1.1 mrg { 424 1.1 mrg mpfr_t x, z; 425 1.1 mrg mpf_t y; 426 1.1 mrg mpfr_prec_t p; 427 1.1 mrg int res1, res2; 428 1.1 mrg int n; 429 1.1 mrg 430 1.1 mrg mpfr_init (x); 431 1.1 mrg mpfr_init2 (z, pmax+GMP_NUMB_BITS); 432 1.1 mrg mpf_init2 (y, MPFR_PREC_MIN); 433 1.1.1.2 mrg 434 1.1.1.2 mrg /* check the erange flag when x is NaN */ 435 1.1.1.2 mrg mpfr_set_nan (x); 436 1.1.1.2 mrg mpf_set_ui (y, 17); 437 1.1.1.2 mrg mpfr_clear_erangeflag (); 438 1.1.1.2 mrg res1 = mpfr_cmp_f (x, y); 439 1.1.1.2 mrg if (res1 != 0 || mpfr_erangeflag_p () == 0) 440 1.1.1.2 mrg { 441 1.1.1.2 mrg printf ("Error for mpfr_cmp_f (NaN, 17)\n"); 442 1.1.1.2 mrg printf ("Return value: expected 0, got %d\n", res1); 443 1.1.1.2 mrg printf ("Erange flag: expected set, got %d\n", mpfr_erangeflag_p ()); 444 1.1.1.2 mrg exit (1); 445 1.1.1.2 mrg } 446 1.1.1.2 mrg 447 1.1.1.5 mrg for (p = pmin ; p < pmax ; p += 3) 448 1.1 mrg { 449 1.1 mrg mpfr_set_prec (x, p); 450 1.1 mrg mpf_set_prec (y, p); 451 1.1.1.5 mrg for (n = 0 ; n < nmax ; n++) 452 1.1 mrg { 453 1.1 mrg mpfr_urandomb (x, RANDS); 454 1.1 mrg mpf_urandomb (y, RANDS, p); 455 1.1 mrg if (!MPFR_IS_SINGULAR (x)) 456 1.1 mrg { 457 1.1 mrg mpfr_set_f (z, y, MPFR_RNDN); 458 1.1 mrg mpfr_sub (z, x, z, MPFR_RNDN); 459 1.1.1.4 mrg res1 = (mpfr_sgn) (z); 460 1.1 mrg res2 = mpfr_cmp_f (x, y); 461 1.1 mrg if (res1 != res2) 462 1.1 mrg { 463 1.1 mrg printf("Error for mpfr_cmp_f: res=%d sub gives %d\n", 464 1.1 mrg res2, res1); 465 1.1 mrg exit (1); 466 1.1 mrg } 467 1.1 mrg } 468 1.1 mrg } 469 1.1 mrg } 470 1.1 mrg mpf_clear (y); 471 1.1 mrg mpfr_clear (x); 472 1.1 mrg mpfr_clear (z); 473 1.1 mrg } 474 1.1 mrg 475 1.1 mrg static void 476 1.1 mrg test_specialz (int (*mpfr_func)(mpfr_ptr, mpfr_srcptr, mpz_srcptr, mpfr_rnd_t), 477 1.1 mrg void (*mpz_func)(mpz_ptr, mpz_srcptr, mpz_srcptr), 478 1.1 mrg const char *op) 479 1.1 mrg { 480 1.1 mrg mpfr_t x1, x2; 481 1.1 mrg mpz_t z1, z2; 482 1.1 mrg int res; 483 1.1 mrg 484 1.1 mrg mpfr_inits2 (128, x1, x2, (mpfr_ptr) 0); 485 1.1 mrg mpz_init (z1); mpz_init(z2); 486 1.1 mrg mpz_fac_ui (z1, 19); /* 19!+1 fits perfectly in a 128 bits mantissa */ 487 1.1 mrg mpz_add_ui (z1, z1, 1); 488 1.1 mrg mpz_fac_ui (z2, 20); /* 20!+1 fits perfectly in a 128 bits mantissa */ 489 1.1 mrg mpz_add_ui (z2, z2, 1); 490 1.1 mrg 491 1.1 mrg res = mpfr_set_z(x1, z1, MPFR_RNDN); 492 1.1 mrg if (res) 493 1.1 mrg { 494 1.1 mrg printf("Specialz %s: set_z1 error\n", op); 495 1.1 mrg exit(1); 496 1.1 mrg } 497 1.1 mrg mpfr_set_z (x2, z2, MPFR_RNDN); 498 1.1 mrg if (res) 499 1.1 mrg { 500 1.1 mrg printf("Specialz %s: set_z2 error\n", op); 501 1.1 mrg exit(1); 502 1.1 mrg } 503 1.1 mrg 504 1.1 mrg /* (19!+1) * (20!+1) fits in a 128 bits number */ 505 1.1 mrg res = mpfr_func(x1, x1, z2, MPFR_RNDN); 506 1.1 mrg if (res) 507 1.1 mrg { 508 1.1 mrg printf("Specialz %s: wrong inexact flag.\n", op); 509 1.1 mrg exit(1); 510 1.1 mrg } 511 1.1 mrg mpz_func(z1, z1, z2); 512 1.1 mrg res = mpfr_set_z (x2, z1, MPFR_RNDN); 513 1.1 mrg if (res) 514 1.1 mrg { 515 1.1 mrg printf("Specialz %s: set_z2 error\n", op); 516 1.1 mrg exit(1); 517 1.1 mrg } 518 1.1 mrg if (mpfr_cmp(x1, x2)) 519 1.1 mrg { 520 1.1 mrg printf("Specialz %s: results differ.\nx1=", op); 521 1.1.1.4 mrg mpfr_dump (x1); 522 1.1.1.4 mrg printf ("x2="); 523 1.1.1.4 mrg mpfr_dump (x2); 524 1.1.1.4 mrg printf ("Z2="); 525 1.1 mrg mpz_out_str (stdout, 2, z1); 526 1.1 mrg putchar('\n'); 527 1.1 mrg exit(1); 528 1.1 mrg } 529 1.1 mrg 530 1.1 mrg mpz_set_ui (z1, 1); 531 1.1 mrg mpz_set_ui (z2, 0); 532 1.1 mrg mpfr_set_ui (x1, 1, MPFR_RNDN); 533 1.1 mrg mpz_func (z1, z1, z2); 534 1.1 mrg res = mpfr_func(x1, x1, z2, MPFR_RNDN); 535 1.1 mrg mpfr_set_z (x2, z1, MPFR_RNDN); 536 1.1 mrg if (mpfr_cmp(x1, x2)) 537 1.1 mrg { 538 1.1 mrg printf("Specialz %s: results differ(2).\nx1=", op); 539 1.1.1.4 mrg mpfr_dump (x1); 540 1.1.1.4 mrg printf ("x2="); 541 1.1.1.4 mrg mpfr_dump (x2); 542 1.1 mrg exit(1); 543 1.1 mrg } 544 1.1 mrg 545 1.1 mrg mpz_clear (z1); mpz_clear(z2); 546 1.1 mrg mpfr_clears (x1, x2, (mpfr_ptr) 0); 547 1.1 mrg } 548 1.1 mrg 549 1.1 mrg static void 550 1.1.1.2 mrg test_special2z (int (*mpfr_func)(mpfr_ptr, mpz_srcptr, mpfr_srcptr, mpfr_rnd_t), 551 1.1.1.2 mrg void (*mpz_func)(mpz_ptr, mpz_srcptr, mpz_srcptr), 552 1.1.1.2 mrg const char *op) 553 1.1.1.2 mrg { 554 1.1.1.2 mrg mpfr_t x1, x2; 555 1.1.1.2 mrg mpz_t z1, z2; 556 1.1.1.2 mrg int res; 557 1.1.1.2 mrg 558 1.1.1.2 mrg mpfr_inits2 (128, x1, x2, (mpfr_ptr) 0); 559 1.1.1.2 mrg mpz_init (z1); mpz_init(z2); 560 1.1.1.2 mrg mpz_fac_ui (z1, 19); /* 19!+1 fits perfectly in a 128 bits mantissa */ 561 1.1.1.2 mrg mpz_add_ui (z1, z1, 1); 562 1.1.1.2 mrg mpz_fac_ui (z2, 20); /* 20!+1 fits perfectly in a 128 bits mantissa */ 563 1.1.1.2 mrg mpz_add_ui (z2, z2, 1); 564 1.1.1.2 mrg 565 1.1.1.2 mrg res = mpfr_set_z(x1, z1, MPFR_RNDN); 566 1.1.1.2 mrg if (res) 567 1.1.1.2 mrg { 568 1.1.1.2 mrg printf("Special2z %s: set_z1 error\n", op); 569 1.1.1.2 mrg exit(1); 570 1.1.1.2 mrg } 571 1.1.1.2 mrg mpfr_set_z (x2, z2, MPFR_RNDN); 572 1.1.1.2 mrg if (res) 573 1.1.1.2 mrg { 574 1.1.1.2 mrg printf("Special2z %s: set_z2 error\n", op); 575 1.1.1.2 mrg exit(1); 576 1.1.1.2 mrg } 577 1.1.1.2 mrg 578 1.1.1.2 mrg /* (19!+1) * (20!+1) fits in a 128 bits number */ 579 1.1.1.2 mrg res = mpfr_func(x1, z1, x2, MPFR_RNDN); 580 1.1.1.2 mrg if (res) 581 1.1.1.2 mrg { 582 1.1.1.2 mrg printf("Special2z %s: wrong inexact flag.\n", op); 583 1.1.1.2 mrg exit(1); 584 1.1.1.2 mrg } 585 1.1.1.2 mrg mpz_func(z1, z1, z2); 586 1.1.1.2 mrg res = mpfr_set_z (x2, z1, MPFR_RNDN); 587 1.1.1.2 mrg if (res) 588 1.1.1.2 mrg { 589 1.1.1.2 mrg printf("Special2z %s: set_z2 error\n", op); 590 1.1.1.2 mrg exit(1); 591 1.1.1.2 mrg } 592 1.1.1.2 mrg if (mpfr_cmp(x1, x2)) 593 1.1.1.2 mrg { 594 1.1.1.2 mrg printf("Special2z %s: results differ.\nx1=", op); 595 1.1.1.4 mrg mpfr_dump (x1); 596 1.1.1.4 mrg printf ("x2="); 597 1.1.1.4 mrg mpfr_dump (x2); 598 1.1.1.4 mrg printf ("Z2="); 599 1.1.1.2 mrg mpz_out_str (stdout, 2, z1); 600 1.1.1.2 mrg putchar('\n'); 601 1.1.1.2 mrg exit(1); 602 1.1.1.2 mrg } 603 1.1.1.2 mrg 604 1.1.1.2 mrg mpz_set_ui (z1, 0); 605 1.1.1.2 mrg mpz_set_ui (z2, 1); 606 1.1.1.2 mrg mpfr_set_ui (x2, 1, MPFR_RNDN); 607 1.1.1.2 mrg res = mpfr_func(x1, z1, x2, MPFR_RNDN); 608 1.1.1.2 mrg mpz_func (z1, z1, z2); 609 1.1.1.2 mrg mpfr_set_z (x2, z1, MPFR_RNDN); 610 1.1.1.2 mrg if (mpfr_cmp(x1, x2)) 611 1.1.1.2 mrg { 612 1.1.1.2 mrg printf("Special2z %s: results differ(2).\nx1=", op); 613 1.1.1.4 mrg mpfr_dump (x1); 614 1.1.1.4 mrg printf ("x2="); 615 1.1.1.4 mrg mpfr_dump (x2); 616 1.1.1.2 mrg exit(1); 617 1.1.1.2 mrg } 618 1.1.1.2 mrg 619 1.1.1.2 mrg mpz_clear (z1); mpz_clear(z2); 620 1.1.1.2 mrg mpfr_clears (x1, x2, (mpfr_ptr) 0); 621 1.1.1.2 mrg } 622 1.1.1.2 mrg 623 1.1.1.2 mrg static void 624 1.1 mrg test_genericz (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int N, 625 1.1 mrg int (*func)(mpfr_ptr, mpfr_srcptr, mpz_srcptr, mpfr_rnd_t), 626 1.1 mrg const char *op) 627 1.1 mrg { 628 1.1 mrg mpfr_prec_t prec; 629 1.1 mrg mpfr_t arg1, dst_big, dst_small, tmp; 630 1.1 mrg mpz_t arg2; 631 1.1 mrg mpfr_rnd_t rnd; 632 1.1 mrg int inexact, compare, compare2; 633 1.1 mrg unsigned int n; 634 1.1 mrg 635 1.1 mrg mpfr_inits (arg1, dst_big, dst_small, tmp, (mpfr_ptr) 0); 636 1.1 mrg mpz_init (arg2); 637 1.1 mrg 638 1.1 mrg for (prec = p0; prec <= p1; prec++) 639 1.1 mrg { 640 1.1 mrg mpfr_set_prec (arg1, prec); 641 1.1 mrg mpfr_set_prec (tmp, prec); 642 1.1 mrg mpfr_set_prec (dst_small, prec); 643 1.1 mrg 644 1.1 mrg for (n=0; n<N; n++) 645 1.1 mrg { 646 1.1 mrg mpfr_urandomb (arg1, RANDS); 647 1.1 mrg mpz_urandomb (arg2, RANDS, 1024); 648 1.1.1.4 mrg rnd = RND_RAND_NO_RNDF (); 649 1.1 mrg mpfr_set_prec (dst_big, 2*prec); 650 1.1.1.4 mrg compare = func (dst_big, arg1, arg2, rnd); 651 1.1 mrg if (mpfr_can_round (dst_big, 2*prec, rnd, rnd, prec)) 652 1.1 mrg { 653 1.1 mrg mpfr_set (tmp, dst_big, rnd); 654 1.1.1.4 mrg inexact = func (dst_small, arg1, arg2, rnd); 655 1.1 mrg if (mpfr_cmp (tmp, dst_small)) 656 1.1 mrg { 657 1.1 mrg printf ("Results differ for prec=%u rnd_mode=%s and %s_z:\n" 658 1.1 mrg "arg1=", 659 1.1 mrg (unsigned) prec, mpfr_print_rnd_mode (rnd), op); 660 1.1.1.4 mrg mpfr_dump (arg1); 661 1.1.1.4 mrg printf ("arg2="); 662 1.1.1.2 mrg mpz_out_str (stdout, 10, arg2); 663 1.1.1.2 mrg printf ("\ngot "); 664 1.1.1.2 mrg mpfr_dump (dst_small); 665 1.1.1.2 mrg printf ("expected "); 666 1.1.1.2 mrg mpfr_dump (tmp); 667 1.1.1.2 mrg printf ("approx "); 668 1.1.1.2 mrg mpfr_dump (dst_big); 669 1.1.1.2 mrg exit (1); 670 1.1.1.2 mrg } 671 1.1.1.2 mrg compare2 = mpfr_cmp (tmp, dst_big); 672 1.1.1.2 mrg /* if rounding to nearest, cannot know the sign of t - f(x) 673 1.1.1.2 mrg because of composed rounding: y = o(f(x)) and t = o(y) */ 674 1.1.1.2 mrg if (compare * compare2 >= 0) 675 1.1.1.2 mrg compare = compare + compare2; 676 1.1.1.2 mrg else 677 1.1.1.2 mrg compare = inexact; /* cannot determine sign(t-f(x)) */ 678 1.1.1.2 mrg if (((inexact == 0) && (compare != 0)) || 679 1.1.1.2 mrg ((inexact > 0) && (compare <= 0)) || 680 1.1.1.2 mrg ((inexact < 0) && (compare >= 0))) 681 1.1.1.2 mrg { 682 1.1.1.2 mrg printf ("Wrong inexact flag for rnd=%s and %s_z:\n" 683 1.1.1.2 mrg "expected %d, got %d\n", 684 1.1.1.2 mrg mpfr_print_rnd_mode (rnd), op, compare, inexact); 685 1.1.1.4 mrg printf ("arg1="); mpfr_dump (arg1); 686 1.1.1.4 mrg printf ("arg2="); mpz_out_str(stdout, 2, arg2); 687 1.1.1.4 mrg printf ("\ndstl="); mpfr_dump (dst_big); 688 1.1.1.4 mrg printf ("dsts="); mpfr_dump (dst_small); 689 1.1.1.4 mrg printf ("tmp ="); mpfr_dump (tmp); 690 1.1.1.2 mrg exit (1); 691 1.1.1.2 mrg } 692 1.1.1.2 mrg } 693 1.1.1.2 mrg } 694 1.1.1.2 mrg } 695 1.1.1.2 mrg 696 1.1.1.2 mrg mpz_clear (arg2); 697 1.1.1.2 mrg mpfr_clears (arg1, dst_big, dst_small, tmp, (mpfr_ptr) 0); 698 1.1.1.2 mrg } 699 1.1.1.2 mrg 700 1.1.1.2 mrg static void 701 1.1.1.2 mrg test_generic2z (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int N, 702 1.1.1.2 mrg int (*func)(mpfr_ptr, mpz_srcptr, mpfr_srcptr, mpfr_rnd_t), 703 1.1.1.2 mrg const char *op) 704 1.1.1.2 mrg { 705 1.1.1.2 mrg mpfr_prec_t prec; 706 1.1.1.2 mrg mpfr_t arg1, dst_big, dst_small, tmp; 707 1.1.1.2 mrg mpz_t arg2; 708 1.1.1.2 mrg mpfr_rnd_t rnd; 709 1.1.1.2 mrg int inexact, compare, compare2; 710 1.1.1.2 mrg unsigned int n; 711 1.1.1.2 mrg 712 1.1.1.2 mrg mpfr_inits (arg1, dst_big, dst_small, tmp, (mpfr_ptr) 0); 713 1.1.1.2 mrg mpz_init (arg2); 714 1.1.1.2 mrg 715 1.1.1.2 mrg for (prec = p0; prec <= p1; prec++) 716 1.1.1.2 mrg { 717 1.1.1.2 mrg mpfr_set_prec (arg1, prec); 718 1.1.1.2 mrg mpfr_set_prec (tmp, prec); 719 1.1.1.2 mrg mpfr_set_prec (dst_small, prec); 720 1.1.1.2 mrg 721 1.1.1.2 mrg for (n=0; n<N; n++) 722 1.1.1.2 mrg { 723 1.1.1.2 mrg mpfr_urandomb (arg1, RANDS); 724 1.1.1.2 mrg mpz_urandomb (arg2, RANDS, 1024); 725 1.1.1.4 mrg rnd = RND_RAND_NO_RNDF (); 726 1.1.1.2 mrg mpfr_set_prec (dst_big, 2*prec); 727 1.1.1.2 mrg compare = func(dst_big, arg2, arg1, rnd); 728 1.1.1.2 mrg if (mpfr_can_round (dst_big, 2*prec, rnd, rnd, prec)) 729 1.1.1.2 mrg { 730 1.1.1.2 mrg mpfr_set (tmp, dst_big, rnd); 731 1.1.1.2 mrg inexact = func(dst_small, arg2, arg1, rnd); 732 1.1.1.2 mrg if (mpfr_cmp (tmp, dst_small)) 733 1.1.1.2 mrg { 734 1.1.1.2 mrg printf ("Results differ for prec=%u rnd_mode=%s and %s_z:\n" 735 1.1.1.2 mrg "arg1=", 736 1.1.1.2 mrg (unsigned) prec, mpfr_print_rnd_mode (rnd), op); 737 1.1.1.4 mrg mpfr_dump (arg1); 738 1.1.1.4 mrg printf ("arg2="); 739 1.1.1.2 mrg mpz_out_str (stdout, 10, arg2); 740 1.1 mrg printf ("\ngot "); 741 1.1 mrg mpfr_dump (dst_small); 742 1.1 mrg printf ("expected "); 743 1.1 mrg mpfr_dump (tmp); 744 1.1 mrg printf ("approx "); 745 1.1 mrg mpfr_dump (dst_big); 746 1.1 mrg exit (1); 747 1.1 mrg } 748 1.1 mrg compare2 = mpfr_cmp (tmp, dst_big); 749 1.1 mrg /* if rounding to nearest, cannot know the sign of t - f(x) 750 1.1 mrg because of composed rounding: y = o(f(x)) and t = o(y) */ 751 1.1 mrg if (compare * compare2 >= 0) 752 1.1 mrg compare = compare + compare2; 753 1.1 mrg else 754 1.1 mrg compare = inexact; /* cannot determine sign(t-f(x)) */ 755 1.1 mrg if (((inexact == 0) && (compare != 0)) || 756 1.1 mrg ((inexact > 0) && (compare <= 0)) || 757 1.1 mrg ((inexact < 0) && (compare >= 0))) 758 1.1 mrg { 759 1.1 mrg printf ("Wrong inexact flag for rnd=%s and %s_z:\n" 760 1.1 mrg "expected %d, got %d\n", 761 1.1 mrg mpfr_print_rnd_mode (rnd), op, compare, inexact); 762 1.1.1.4 mrg printf ("arg1="); mpfr_dump (arg1); 763 1.1.1.4 mrg printf ("arg2="); mpz_out_str(stdout, 2, arg2); 764 1.1.1.4 mrg printf ("\ndstl="); mpfr_dump (dst_big); 765 1.1.1.4 mrg printf ("dsts="); mpfr_dump (dst_small); 766 1.1.1.4 mrg printf ("tmp ="); mpfr_dump (tmp); 767 1.1 mrg exit (1); 768 1.1 mrg } 769 1.1 mrg } 770 1.1 mrg } 771 1.1 mrg } 772 1.1 mrg 773 1.1 mrg mpz_clear (arg2); 774 1.1 mrg mpfr_clears (arg1, dst_big, dst_small, tmp, (mpfr_ptr) 0); 775 1.1 mrg } 776 1.1 mrg 777 1.1 mrg static void 778 1.1 mrg test_genericq (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int N, 779 1.1 mrg int (*func)(mpfr_ptr, mpfr_srcptr, mpq_srcptr, mpfr_rnd_t), 780 1.1 mrg const char *op) 781 1.1 mrg { 782 1.1 mrg mpfr_prec_t prec; 783 1.1 mrg mpfr_t arg1, dst_big, dst_small, tmp; 784 1.1 mrg mpq_t arg2; 785 1.1 mrg mpfr_rnd_t rnd; 786 1.1 mrg int inexact, compare, compare2; 787 1.1 mrg unsigned int n; 788 1.1 mrg 789 1.1 mrg mpfr_inits (arg1, dst_big, dst_small, tmp, (mpfr_ptr) 0); 790 1.1 mrg mpq_init (arg2); 791 1.1 mrg 792 1.1 mrg for (prec = p0; prec <= p1; prec++) 793 1.1 mrg { 794 1.1 mrg mpfr_set_prec (arg1, prec); 795 1.1 mrg mpfr_set_prec (tmp, prec); 796 1.1 mrg mpfr_set_prec (dst_small, prec); 797 1.1 mrg 798 1.1 mrg for (n=0; n<N; n++) 799 1.1 mrg { 800 1.1 mrg mpfr_urandomb (arg1, RANDS); 801 1.1 mrg mpq_set_ui (arg2, randlimb (), randlimb() ); 802 1.1 mrg mpq_canonicalize (arg2); 803 1.1.1.4 mrg rnd = RND_RAND_NO_RNDF (); 804 1.1 mrg mpfr_set_prec (dst_big, prec+10); 805 1.1 mrg compare = func(dst_big, arg1, arg2, rnd); 806 1.1 mrg if (mpfr_can_round (dst_big, prec+10, rnd, rnd, prec)) 807 1.1 mrg { 808 1.1 mrg mpfr_set (tmp, dst_big, rnd); 809 1.1 mrg inexact = func(dst_small, arg1, arg2, rnd); 810 1.1 mrg if (mpfr_cmp (tmp, dst_small)) 811 1.1 mrg { 812 1.1 mrg printf ("Results differ for prec=%u rnd_mode=%s and %s_q:\n" 813 1.1 mrg "arg1=", 814 1.1 mrg (unsigned) prec, mpfr_print_rnd_mode (rnd), op); 815 1.1.1.4 mrg mpfr_dump (arg1); 816 1.1.1.4 mrg printf ("arg2="); 817 1.1 mrg mpq_out_str(stdout, 2, arg2); 818 1.1 mrg printf ("\ngot "); 819 1.1.1.4 mrg mpfr_dump (dst_small); 820 1.1.1.4 mrg printf ("expected "); 821 1.1.1.4 mrg mpfr_dump (tmp); 822 1.1.1.4 mrg printf ("approx "); 823 1.1.1.4 mrg mpfr_dump (dst_big); 824 1.1 mrg exit (1); 825 1.1 mrg } 826 1.1 mrg compare2 = mpfr_cmp (tmp, dst_big); 827 1.1 mrg /* if rounding to nearest, cannot know the sign of t - f(x) 828 1.1 mrg because of composed rounding: y = o(f(x)) and t = o(y) */ 829 1.1 mrg if (compare * compare2 >= 0) 830 1.1 mrg compare = compare + compare2; 831 1.1 mrg else 832 1.1 mrg compare = inexact; /* cannot determine sign(t-f(x)) */ 833 1.1 mrg if (((inexact == 0) && (compare != 0)) || 834 1.1 mrg ((inexact > 0) && (compare <= 0)) || 835 1.1 mrg ((inexact < 0) && (compare >= 0))) 836 1.1 mrg { 837 1.1 mrg printf ("Wrong inexact flag for rnd=%s and %s_q:\n" 838 1.1 mrg "expected %d, got %d", 839 1.1 mrg mpfr_print_rnd_mode (rnd), op, compare, inexact); 840 1.1.1.4 mrg printf ("arg1="); mpfr_dump (arg1); 841 1.1.1.4 mrg printf ("arg2="); mpq_out_str(stdout, 2, arg2); 842 1.1.1.4 mrg printf ("\ndstl="); mpfr_dump (dst_big); 843 1.1.1.4 mrg printf ("dsts="); mpfr_dump (dst_small); 844 1.1.1.4 mrg printf ("tmp ="); mpfr_dump (tmp); 845 1.1 mrg exit (1); 846 1.1 mrg } 847 1.1 mrg } 848 1.1 mrg } 849 1.1 mrg } 850 1.1 mrg 851 1.1 mrg mpq_clear (arg2); 852 1.1 mrg mpfr_clears (arg1, dst_big, dst_small, tmp, (mpfr_ptr) 0); 853 1.1 mrg } 854 1.1 mrg 855 1.1 mrg static void 856 1.1 mrg test_specialq (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int N, 857 1.1 mrg int (*mpfr_func)(mpfr_ptr, mpfr_srcptr, mpq_srcptr, mpfr_rnd_t), 858 1.1 mrg void (*mpq_func)(mpq_ptr, mpq_srcptr, mpq_srcptr), 859 1.1 mrg const char *op) 860 1.1 mrg { 861 1.1 mrg mpfr_t fra, frb, frq; 862 1.1 mrg mpq_t q1, q2, qr; 863 1.1 mrg unsigned int n; 864 1.1 mrg mpfr_prec_t prec; 865 1.1 mrg 866 1.1 mrg for (prec = p0 ; prec < p1 ; prec++) 867 1.1 mrg { 868 1.1 mrg mpfr_inits2 (prec, fra, frb, frq, (mpfr_ptr) 0); 869 1.1 mrg mpq_init (q1); mpq_init(q2); mpq_init (qr); 870 1.1 mrg 871 1.1.1.5 mrg for (n = 0 ; n < N ; n++) 872 1.1 mrg { 873 1.1 mrg mpq_set_ui(q1, randlimb(), randlimb() ); 874 1.1 mrg mpq_set_ui(q2, randlimb(), randlimb() ); 875 1.1 mrg mpq_canonicalize (q1); 876 1.1 mrg mpq_canonicalize (q2); 877 1.1 mrg mpq_func (qr, q1, q2); 878 1.1 mrg mpfr_set_q (fra, q1, MPFR_RNDD); 879 1.1 mrg mpfr_func (fra, fra, q2, MPFR_RNDD); 880 1.1 mrg mpfr_set_q (frb, q1, MPFR_RNDU); 881 1.1 mrg mpfr_func (frb, frb, q2, MPFR_RNDU); 882 1.1 mrg mpfr_set_q (frq, qr, MPFR_RNDN); 883 1.1 mrg /* We should have fra <= qr <= frb */ 884 1.1 mrg if ( (mpfr_cmp(fra, frq) > 0) || (mpfr_cmp (frq, frb) > 0)) 885 1.1 mrg { 886 1.1.1.2 mrg printf("Range error for prec=%lu and %s", 887 1.1.1.2 mrg (unsigned long) prec, op); 888 1.1 mrg printf ("\nq1="); mpq_out_str(stdout, 2, q1); 889 1.1 mrg printf ("\nq2="); mpq_out_str(stdout, 2, q2); 890 1.1.1.4 mrg printf ("\nfr_dn="); mpfr_dump (fra); 891 1.1.1.4 mrg printf ("fr_q ="); mpfr_dump (frq); 892 1.1.1.4 mrg printf ("fr_up="); mpfr_dump (frb); 893 1.1 mrg exit (1); 894 1.1 mrg } 895 1.1 mrg } 896 1.1 mrg 897 1.1 mrg mpq_clear (q1); mpq_clear (q2); mpq_clear (qr); 898 1.1 mrg mpfr_clears (fra, frb, frq, (mpfr_ptr) 0); 899 1.1 mrg } 900 1.1 mrg } 901 1.1 mrg 902 1.1.1.2 mrg static void 903 1.1.1.2 mrg bug_mul_q_20100810 (void) 904 1.1.1.2 mrg { 905 1.1.1.2 mrg mpfr_t x; 906 1.1.1.2 mrg mpfr_t y; 907 1.1.1.2 mrg mpq_t q; 908 1.1.1.2 mrg int inexact; 909 1.1.1.2 mrg 910 1.1.1.2 mrg mpfr_init (x); 911 1.1.1.2 mrg mpfr_init (y); 912 1.1.1.2 mrg mpq_init (q); 913 1.1.1.2 mrg 914 1.1.1.2 mrg /* mpfr_mul_q: the inexact value must be set in case of overflow */ 915 1.1.1.2 mrg mpq_set_ui (q, 4096, 3); 916 1.1.1.2 mrg mpfr_set_inf (x, +1); 917 1.1.1.2 mrg mpfr_nextbelow (x); 918 1.1.1.2 mrg inexact = mpfr_mul_q (y, x, q, MPFR_RNDU); 919 1.1.1.2 mrg 920 1.1.1.2 mrg if (inexact <= 0) 921 1.1.1.2 mrg { 922 1.1.1.2 mrg printf ("Overflow error in mpfr_mul_q. "); 923 1.1.1.2 mrg printf ("Wrong inexact flag: got %d, should be positive.\n", inexact); 924 1.1.1.2 mrg 925 1.1.1.2 mrg exit (1); 926 1.1.1.2 mrg } 927 1.1.1.2 mrg if (!mpfr_inf_p (y)) 928 1.1.1.2 mrg { 929 1.1.1.2 mrg printf ("Overflow error in mpfr_mul_q (y, x, q, MPFR_RNDD). "); 930 1.1.1.2 mrg printf ("\nx = "); 931 1.1.1.2 mrg mpfr_out_str (stdout, 10, 0, x, MPFR_RNDD); 932 1.1.1.2 mrg printf ("\nq = "); 933 1.1.1.2 mrg mpq_out_str (stdout, 10, q); 934 1.1.1.2 mrg printf ("\ny = "); 935 1.1.1.2 mrg mpfr_out_str (stdout, 10, 0, y, MPFR_RNDD); 936 1.1.1.2 mrg printf (" (should be +infinity)\n"); 937 1.1.1.2 mrg 938 1.1.1.2 mrg exit (1); 939 1.1.1.2 mrg } 940 1.1.1.2 mrg 941 1.1.1.2 mrg mpq_clear (q); 942 1.1.1.2 mrg mpfr_clear (y); 943 1.1.1.2 mrg mpfr_clear (x); 944 1.1.1.2 mrg } 945 1.1.1.2 mrg 946 1.1.1.2 mrg static void 947 1.1.1.2 mrg bug_div_q_20100810 (void) 948 1.1.1.2 mrg { 949 1.1.1.2 mrg mpfr_t x; 950 1.1.1.2 mrg mpfr_t y; 951 1.1.1.2 mrg mpq_t q; 952 1.1.1.2 mrg int inexact; 953 1.1.1.2 mrg 954 1.1.1.2 mrg mpfr_init (x); 955 1.1.1.2 mrg mpfr_init (y); 956 1.1.1.2 mrg mpq_init (q); 957 1.1.1.2 mrg 958 1.1.1.2 mrg /* mpfr_div_q: the inexact value must be set in case of overflow */ 959 1.1.1.2 mrg mpq_set_ui (q, 3, 4096); 960 1.1.1.2 mrg mpfr_set_inf (x, +1); 961 1.1.1.2 mrg mpfr_nextbelow (x); 962 1.1.1.2 mrg inexact = mpfr_div_q (y, x, q, MPFR_RNDU); 963 1.1.1.2 mrg 964 1.1.1.2 mrg if (inexact <= 0) 965 1.1.1.2 mrg { 966 1.1.1.2 mrg printf ("Overflow error in mpfr_div_q. "); 967 1.1.1.2 mrg printf ("Wrong inexact flag: got %d, should be positive.\n", inexact); 968 1.1.1.2 mrg 969 1.1.1.2 mrg exit (1); 970 1.1.1.2 mrg } 971 1.1.1.2 mrg if (!mpfr_inf_p (y)) 972 1.1.1.2 mrg { 973 1.1.1.2 mrg printf ("Overflow error in mpfr_div_q (y, x, q, MPFR_RNDD). "); 974 1.1.1.2 mrg printf ("\nx = "); 975 1.1.1.2 mrg mpfr_out_str (stdout, 10, 0, x, MPFR_RNDD); 976 1.1.1.2 mrg printf ("\nq = "); 977 1.1.1.2 mrg mpq_out_str (stdout, 10, q); 978 1.1.1.2 mrg printf ("\ny = "); 979 1.1.1.2 mrg mpfr_out_str (stdout, 10, 0, y, MPFR_RNDD); 980 1.1.1.2 mrg printf (" (should be +infinity)\n"); 981 1.1.1.2 mrg 982 1.1.1.2 mrg exit (1); 983 1.1.1.2 mrg } 984 1.1.1.2 mrg 985 1.1.1.2 mrg mpq_clear (q); 986 1.1.1.2 mrg mpfr_clear (y); 987 1.1.1.2 mrg mpfr_clear (x); 988 1.1.1.2 mrg } 989 1.1.1.2 mrg 990 1.1.1.2 mrg static void 991 1.1.1.2 mrg bug_mul_div_q_20100818 (void) 992 1.1.1.2 mrg { 993 1.1.1.2 mrg mpq_t qa, qb; 994 1.1.1.2 mrg mpfr_t x1, x2, y1, y2, y3; 995 1.1.1.2 mrg mpfr_exp_t emin, emax, e; 996 1.1.1.2 mrg int inex; 997 1.1.1.2 mrg int rnd; 998 1.1.1.2 mrg 999 1.1.1.2 mrg emin = mpfr_get_emin (); 1000 1.1.1.2 mrg emax = mpfr_get_emax (); 1001 1.1.1.2 mrg set_emin (MPFR_EMIN_MIN); 1002 1.1.1.2 mrg set_emax (MPFR_EMAX_MAX); 1003 1.1.1.2 mrg 1004 1.1.1.2 mrg mpq_init (qa); 1005 1.1.1.2 mrg mpq_init (qb); 1006 1.1.1.2 mrg mpfr_inits2 (32, x1, x2, y1, y2, y3, (mpfr_ptr) 0); 1007 1.1.1.2 mrg 1008 1.1.1.2 mrg mpq_set_ui (qa, 3, 17); 1009 1.1.1.2 mrg mpq_set_ui (qb, 17, 3); 1010 1.1.1.2 mrg inex = mpfr_set_ui (x1, 7, MPFR_RNDN); 1011 1.1.1.2 mrg MPFR_ASSERTN (inex == 0); 1012 1.1.1.2 mrg 1013 1.1.1.2 mrg e = MPFR_EMAX_MAX - 3; 1014 1.1.1.2 mrg inex = mpfr_set_ui_2exp (x2, 7, e, MPFR_RNDN); /* x2 = x1 * 2^e */ 1015 1.1.1.2 mrg MPFR_ASSERTN (inex == 0); 1016 1.1.1.2 mrg 1017 1.1.1.2 mrg RND_LOOP(rnd) 1018 1.1.1.2 mrg { 1019 1.1.1.2 mrg mpfr_mul_q (y1, x1, qa, (mpfr_rnd_t) rnd); 1020 1.1.1.2 mrg mpfr_div_q (y3, x1, qb, (mpfr_rnd_t) rnd); 1021 1.1.1.2 mrg MPFR_ASSERTN (mpfr_equal_p (y1, y3)); 1022 1.1.1.2 mrg inex = mpfr_set_ui_2exp (y3, 1, e, MPFR_RNDN); 1023 1.1.1.2 mrg MPFR_ASSERTN (inex == 0); 1024 1.1.1.2 mrg inex = mpfr_mul (y3, y3, y1, MPFR_RNDN); /* y3 = y1 * 2^e */ 1025 1.1.1.2 mrg MPFR_ASSERTN (inex == 0); 1026 1.1.1.2 mrg mpfr_mul_q (y2, x2, qa, (mpfr_rnd_t) rnd); 1027 1.1.1.2 mrg if (! mpfr_equal_p (y2, y3)) 1028 1.1.1.2 mrg { 1029 1.1.1.2 mrg printf ("Error 1 in bug_mul_div_q_20100818 (rnd = %d)\n", rnd); 1030 1.1.1.2 mrg printf ("Expected "); mpfr_dump (y3); 1031 1.1.1.2 mrg printf ("Got "); mpfr_dump (y2); 1032 1.1.1.2 mrg exit (1); 1033 1.1.1.2 mrg } 1034 1.1.1.2 mrg mpfr_div_q (y2, x2, qb, (mpfr_rnd_t) rnd); 1035 1.1.1.2 mrg if (! mpfr_equal_p (y2, y3)) 1036 1.1.1.2 mrg { 1037 1.1.1.2 mrg printf ("Error 2 in bug_mul_div_q_20100818 (rnd = %d)\n", rnd); 1038 1.1.1.2 mrg printf ("Expected "); mpfr_dump (y3); 1039 1.1.1.2 mrg printf ("Got "); mpfr_dump (y2); 1040 1.1.1.2 mrg exit (1); 1041 1.1.1.2 mrg } 1042 1.1.1.2 mrg } 1043 1.1.1.2 mrg 1044 1.1.1.2 mrg e = MPFR_EMIN_MIN; 1045 1.1.1.2 mrg inex = mpfr_set_ui_2exp (x2, 7, e, MPFR_RNDN); /* x2 = x1 * 2^e */ 1046 1.1.1.2 mrg MPFR_ASSERTN (inex == 0); 1047 1.1.1.2 mrg 1048 1.1.1.2 mrg RND_LOOP(rnd) 1049 1.1.1.2 mrg { 1050 1.1.1.2 mrg mpfr_div_q (y1, x1, qa, (mpfr_rnd_t) rnd); 1051 1.1.1.2 mrg mpfr_mul_q (y3, x1, qb, (mpfr_rnd_t) rnd); 1052 1.1.1.2 mrg MPFR_ASSERTN (mpfr_equal_p (y1, y3)); 1053 1.1.1.2 mrg inex = mpfr_set_ui_2exp (y3, 1, e, MPFR_RNDN); 1054 1.1.1.2 mrg MPFR_ASSERTN (inex == 0); 1055 1.1.1.2 mrg inex = mpfr_mul (y3, y3, y1, MPFR_RNDN); /* y3 = y1 * 2^e */ 1056 1.1.1.2 mrg MPFR_ASSERTN (inex == 0); 1057 1.1.1.2 mrg mpfr_div_q (y2, x2, qa, (mpfr_rnd_t) rnd); 1058 1.1.1.2 mrg if (! mpfr_equal_p (y2, y3)) 1059 1.1.1.2 mrg { 1060 1.1.1.2 mrg printf ("Error 3 in bug_mul_div_q_20100818 (rnd = %d)\n", rnd); 1061 1.1.1.2 mrg printf ("Expected "); mpfr_dump (y3); 1062 1.1.1.2 mrg printf ("Got "); mpfr_dump (y2); 1063 1.1.1.2 mrg exit (1); 1064 1.1.1.2 mrg } 1065 1.1.1.2 mrg mpfr_mul_q (y2, x2, qb, (mpfr_rnd_t) rnd); 1066 1.1.1.2 mrg if (! mpfr_equal_p (y2, y3)) 1067 1.1.1.2 mrg { 1068 1.1.1.2 mrg printf ("Error 4 in bug_mul_div_q_20100818 (rnd = %d)\n", rnd); 1069 1.1.1.2 mrg printf ("Expected "); mpfr_dump (y3); 1070 1.1.1.2 mrg printf ("Got "); mpfr_dump (y2); 1071 1.1.1.2 mrg exit (1); 1072 1.1.1.2 mrg } 1073 1.1.1.2 mrg } 1074 1.1.1.2 mrg 1075 1.1.1.2 mrg mpq_clear (qa); 1076 1.1.1.2 mrg mpq_clear (qb); 1077 1.1.1.2 mrg mpfr_clears (x1, x2, y1, y2, y3, (mpfr_ptr) 0); 1078 1.1.1.2 mrg 1079 1.1.1.2 mrg set_emin (emin); 1080 1.1.1.2 mrg set_emax (emax); 1081 1.1.1.2 mrg } 1082 1.1.1.2 mrg 1083 1.1.1.2 mrg static void 1084 1.1.1.2 mrg reduced_expo_range (void) 1085 1.1.1.2 mrg { 1086 1.1.1.2 mrg mpfr_t x; 1087 1.1.1.2 mrg mpz_t z; 1088 1.1.1.2 mrg mpq_t q; 1089 1.1.1.2 mrg mpfr_exp_t emin; 1090 1.1.1.2 mrg int inex; 1091 1.1.1.2 mrg 1092 1.1.1.2 mrg emin = mpfr_get_emin (); 1093 1.1.1.2 mrg set_emin (4); 1094 1.1.1.2 mrg 1095 1.1.1.2 mrg mpfr_init2 (x, 32); 1096 1.1.1.2 mrg 1097 1.1.1.2 mrg mpz_init (z); 1098 1.1.1.2 mrg mpfr_clear_flags (); 1099 1.1.1.2 mrg inex = mpfr_set_ui (x, 17, MPFR_RNDN); 1100 1.1.1.2 mrg MPFR_ASSERTN (inex == 0); 1101 1.1.1.2 mrg mpz_set_ui (z, 3); 1102 1.1.1.2 mrg inex = mpfr_mul_z (x, x, z, MPFR_RNDN); 1103 1.1.1.2 mrg if (inex != 0 || MPFR_IS_NAN (x) || mpfr_cmp_ui (x, 51) != 0) 1104 1.1.1.2 mrg { 1105 1.1.1.2 mrg printf ("Error 1 in reduce_expo_range: expected 51 with inex = 0," 1106 1.1.1.2 mrg " got\n"); 1107 1.1.1.2 mrg mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN); 1108 1.1.1.2 mrg printf ("with inex = %d\n", inex); 1109 1.1.1.2 mrg exit (1); 1110 1.1.1.2 mrg } 1111 1.1.1.2 mrg inex = mpfr_div_z (x, x, z, MPFR_RNDN); 1112 1.1.1.2 mrg if (inex != 0 || MPFR_IS_NAN (x) || mpfr_cmp_ui (x, 17) != 0) 1113 1.1.1.2 mrg { 1114 1.1.1.2 mrg printf ("Error 2 in reduce_expo_range: expected 17 with inex = 0," 1115 1.1.1.2 mrg " got\n"); 1116 1.1.1.2 mrg mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN); 1117 1.1.1.2 mrg printf ("with inex = %d\n", inex); 1118 1.1.1.2 mrg exit (1); 1119 1.1.1.2 mrg } 1120 1.1.1.2 mrg inex = mpfr_add_z (x, x, z, MPFR_RNDN); 1121 1.1.1.2 mrg if (inex != 0 || MPFR_IS_NAN (x) || mpfr_cmp_ui (x, 20) != 0) 1122 1.1.1.2 mrg { 1123 1.1.1.2 mrg printf ("Error 3 in reduce_expo_range: expected 20 with inex = 0," 1124 1.1.1.2 mrg " got\n"); 1125 1.1.1.2 mrg mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN); 1126 1.1.1.2 mrg printf ("with inex = %d\n", inex); 1127 1.1.1.2 mrg exit (1); 1128 1.1.1.2 mrg } 1129 1.1.1.2 mrg inex = mpfr_sub_z (x, x, z, MPFR_RNDN); 1130 1.1.1.2 mrg if (inex != 0 || MPFR_IS_NAN (x) || mpfr_cmp_ui (x, 17) != 0) 1131 1.1.1.2 mrg { 1132 1.1.1.2 mrg printf ("Error 4 in reduce_expo_range: expected 17 with inex = 0," 1133 1.1.1.2 mrg " got\n"); 1134 1.1.1.2 mrg mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN); 1135 1.1.1.2 mrg printf ("with inex = %d\n", inex); 1136 1.1.1.2 mrg exit (1); 1137 1.1.1.2 mrg } 1138 1.1.1.2 mrg MPFR_ASSERTN (__gmpfr_flags == 0); 1139 1.1.1.2 mrg if (mpfr_cmp_z (x, z) <= 0) 1140 1.1.1.2 mrg { 1141 1.1.1.2 mrg printf ("Error 5 in reduce_expo_range: expected a positive value.\n"); 1142 1.1.1.2 mrg exit (1); 1143 1.1.1.2 mrg } 1144 1.1.1.2 mrg mpz_clear (z); 1145 1.1.1.2 mrg 1146 1.1.1.2 mrg mpq_init (q); 1147 1.1.1.2 mrg mpq_set_ui (q, 1, 1); 1148 1.1.1.2 mrg mpfr_set_ui (x, 16, MPFR_RNDN); 1149 1.1.1.2 mrg inex = mpfr_add_q (x, x, q, MPFR_RNDN); 1150 1.1.1.2 mrg if (inex != 0 || MPFR_IS_NAN (x) || mpfr_cmp_ui (x, 17) != 0) 1151 1.1.1.2 mrg { 1152 1.1.1.2 mrg printf ("Error in reduce_expo_range for 16 + 1/1," 1153 1.1.1.2 mrg " got inex = %d and\nx = ", inex); 1154 1.1.1.2 mrg mpfr_dump (x); 1155 1.1.1.2 mrg exit (1); 1156 1.1.1.2 mrg } 1157 1.1.1.2 mrg inex = mpfr_sub_q (x, x, q, MPFR_RNDN); 1158 1.1.1.2 mrg if (inex != 0 || MPFR_IS_NAN (x) || mpfr_cmp_ui (x, 16) != 0) 1159 1.1.1.2 mrg { 1160 1.1.1.2 mrg printf ("Error in reduce_expo_range for 17 - 1/1," 1161 1.1.1.2 mrg " got inex = %d and\nx = ", inex); 1162 1.1.1.2 mrg mpfr_dump (x); 1163 1.1.1.2 mrg exit (1); 1164 1.1.1.2 mrg } 1165 1.1.1.2 mrg mpq_clear (q); 1166 1.1.1.2 mrg 1167 1.1.1.2 mrg mpfr_clear (x); 1168 1.1.1.2 mrg 1169 1.1.1.2 mrg set_emin (emin); 1170 1.1.1.2 mrg } 1171 1.1.1.2 mrg 1172 1.1.1.2 mrg static void 1173 1.1.1.2 mrg addsubq_overflow_aux (mpfr_exp_t e) 1174 1.1.1.2 mrg { 1175 1.1.1.2 mrg mpfr_t x, y; 1176 1.1.1.2 mrg mpq_t q; 1177 1.1.1.2 mrg mpfr_exp_t emax; 1178 1.1.1.2 mrg int inex; 1179 1.1.1.2 mrg int rnd; 1180 1.1.1.2 mrg int sign, sub; 1181 1.1.1.2 mrg 1182 1.1.1.2 mrg MPFR_ASSERTN (e <= LONG_MAX); 1183 1.1.1.2 mrg emax = mpfr_get_emax (); 1184 1.1.1.2 mrg set_emax (e); 1185 1.1.1.2 mrg mpfr_inits2 (16, x, y, (mpfr_ptr) 0); 1186 1.1.1.2 mrg mpq_init (q); 1187 1.1.1.2 mrg 1188 1.1.1.2 mrg mpfr_set_inf (x, 1); 1189 1.1.1.2 mrg mpfr_nextbelow (x); 1190 1.1.1.2 mrg mpq_set_ui (q, 1, 1); 1191 1.1.1.2 mrg 1192 1.1.1.2 mrg for (sign = 0; sign <= 1; sign++) 1193 1.1.1.2 mrg { 1194 1.1.1.2 mrg for (sub = 0; sub <= 1; sub++) 1195 1.1.1.2 mrg { 1196 1.1.1.2 mrg RND_LOOP(rnd) 1197 1.1.1.2 mrg { 1198 1.1.1.2 mrg unsigned int flags, ex_flags; 1199 1.1.1.2 mrg int inf; 1200 1.1.1.2 mrg 1201 1.1.1.2 mrg inf = rnd == MPFR_RNDA || 1202 1.1.1.2 mrg rnd == (sign ? MPFR_RNDD : MPFR_RNDU); 1203 1.1.1.2 mrg ex_flags = MPFR_FLAGS_INEXACT | (inf ? MPFR_FLAGS_OVERFLOW : 0); 1204 1.1.1.2 mrg mpfr_clear_flags (); 1205 1.1.1.2 mrg inex = sub ? 1206 1.1.1.2 mrg mpfr_sub_q (y, x, q, (mpfr_rnd_t) rnd) : 1207 1.1.1.2 mrg mpfr_add_q (y, x, q, (mpfr_rnd_t) rnd); 1208 1.1.1.2 mrg flags = __gmpfr_flags; 1209 1.1.1.2 mrg if (inex == 0 || flags != ex_flags || 1210 1.1.1.2 mrg (inf ? ! mpfr_inf_p (y) : ! mpfr_equal_p (x, y))) 1211 1.1.1.2 mrg { 1212 1.1.1.2 mrg printf ("Error in addsubq_overflow_aux(%ld)," 1213 1.1.1.2 mrg " sign = %d, %s\n", (long) e, sign, 1214 1.1.1.2 mrg mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 1215 1.1.1.2 mrg printf ("Got inex = %d, y = ", inex); 1216 1.1.1.2 mrg mpfr_dump (y); 1217 1.1.1.2 mrg printf ("Expected flags:"); 1218 1.1.1.2 mrg flags_out (ex_flags); 1219 1.1.1.2 mrg printf ("Got flags: "); 1220 1.1.1.2 mrg flags_out (flags); 1221 1.1.1.2 mrg exit (1); 1222 1.1.1.2 mrg } 1223 1.1.1.2 mrg } 1224 1.1.1.2 mrg mpq_neg (q, q); 1225 1.1.1.2 mrg } 1226 1.1.1.2 mrg mpfr_neg (x, x, MPFR_RNDN); 1227 1.1.1.2 mrg mpq_neg (q, q); 1228 1.1.1.2 mrg } 1229 1.1.1.2 mrg 1230 1.1.1.2 mrg mpq_clear (q); 1231 1.1.1.2 mrg mpfr_clears (x, y, (mpfr_ptr) 0); 1232 1.1.1.2 mrg set_emax (emax); 1233 1.1.1.2 mrg } 1234 1.1.1.2 mrg 1235 1.1.1.2 mrg static void 1236 1.1.1.2 mrg addsubq_overflow (void) 1237 1.1.1.2 mrg { 1238 1.1.1.2 mrg addsubq_overflow_aux (4913); 1239 1.1.1.2 mrg addsubq_overflow_aux (MPFR_EMAX_MAX); 1240 1.1.1.2 mrg } 1241 1.1.1.2 mrg 1242 1.1.1.2 mrg static void 1243 1.1.1.2 mrg coverage_mpfr_mul_q_20110218 (void) 1244 1.1.1.2 mrg { 1245 1.1.1.2 mrg mpfr_t cmp, res, op1; 1246 1.1.1.2 mrg mpq_t op2; 1247 1.1.1.2 mrg int status; 1248 1.1.1.2 mrg 1249 1.1.1.2 mrg mpfr_init2 (cmp, MPFR_PREC_MIN); 1250 1.1.1.2 mrg mpfr_init2 (res, MPFR_PREC_MIN); 1251 1.1.1.2 mrg mpfr_init_set_si (op1, 1, MPFR_RNDN); 1252 1.1.1.2 mrg 1253 1.1.1.2 mrg mpq_init (op2); 1254 1.1.1.2 mrg mpq_set_si (op2, 0, 0); 1255 1.1.1.2 mrg mpz_set_si (mpq_denref (op2), 0); 1256 1.1.1.2 mrg 1257 1.1.1.2 mrg status = mpfr_mul_q (res, op1, op2, MPFR_RNDN); 1258 1.1.1.2 mrg 1259 1.1.1.2 mrg if ((status != 0) || (mpfr_cmp (cmp, res) != 0)) 1260 1.1.1.2 mrg { 1261 1.1.1.2 mrg printf ("Results differ %d.\nres=", status); 1262 1.1.1.4 mrg mpfr_dump (res); 1263 1.1.1.4 mrg printf ("cmp="); 1264 1.1.1.4 mrg mpfr_dump (cmp); 1265 1.1.1.2 mrg exit (1); 1266 1.1.1.2 mrg } 1267 1.1.1.2 mrg 1268 1.1.1.2 mrg mpfr_set_si (op1, 1, MPFR_RNDN); 1269 1.1.1.2 mrg mpq_set_si (op2, -1, 0); 1270 1.1.1.2 mrg 1271 1.1.1.2 mrg status = mpfr_mul_q (res, op1, op2, MPFR_RNDN); 1272 1.1.1.2 mrg 1273 1.1.1.2 mrg mpfr_set_inf (cmp, -1); 1274 1.1.1.2 mrg if ((status != 0) || (mpfr_cmp(res, cmp) != 0)) 1275 1.1.1.2 mrg { 1276 1.1.1.4 mrg printf ("mpfr_mul_q 1 * (-1/0) returned a wrong value:\n"); 1277 1.1.1.4 mrg printf (" expected "); 1278 1.1.1.4 mrg mpfr_dump (cmp); 1279 1.1.1.4 mrg printf (" got "); 1280 1.1.1.4 mrg mpfr_dump (res); 1281 1.1.1.4 mrg printf (" ternary value is %d\n", status); 1282 1.1.1.2 mrg exit (1); 1283 1.1.1.2 mrg } 1284 1.1.1.2 mrg 1285 1.1.1.2 mrg mpq_clear (op2); 1286 1.1.1.2 mrg mpfr_clear (op1); 1287 1.1.1.2 mrg mpfr_clear (res); 1288 1.1.1.2 mrg mpfr_clear (cmp); 1289 1.1.1.2 mrg } 1290 1.1.1.2 mrg 1291 1.1.1.5 mrg static void 1292 1.1.1.5 mrg coverage (void) 1293 1.1.1.5 mrg { 1294 1.1.1.5 mrg mpfr_exp_t emax, emin; 1295 1.1.1.5 mrg mpz_t z; 1296 1.1.1.5 mrg mpfr_t x; 1297 1.1.1.5 mrg int cmp; 1298 1.1.1.5 mrg 1299 1.1.1.5 mrg mpz_init (z); 1300 1.1.1.5 mrg mpfr_init2 (x, 5); 1301 1.1.1.5 mrg 1302 1.1.1.5 mrg /* coverage for mpfr_cmp_z in case of overflow */ 1303 1.1.1.5 mrg emax = mpfr_get_emax (); 1304 1.1.1.6 mrg set_emax (63); 1305 1.1.1.5 mrg mpz_set_str (z, "9223372036854775808", 10); /* 2^63 */ 1306 1.1.1.5 mrg mpfr_set_ui_2exp (x, 1, mpfr_get_emax (), MPFR_RNDZ); 1307 1.1.1.5 mrg /* x = (1-2^(-p))*2^emax */ 1308 1.1.1.5 mrg mpfr_clear_flags (); 1309 1.1.1.5 mrg cmp = mpfr_cmp_z (x, z); 1310 1.1.1.5 mrg MPFR_ASSERTN(cmp < 0); 1311 1.1.1.5 mrg MPFR_ASSERTN(!mpfr_overflow_p ()); 1312 1.1.1.6 mrg set_emax (emax); 1313 1.1.1.5 mrg 1314 1.1.1.5 mrg /* coverage for mpfr_cmp_z in case of underflow */ 1315 1.1.1.5 mrg mpz_set_str (z, "18446744073709551615", 10); /* 2^64-1 */ 1316 1.1.1.5 mrg emin = mpfr_get_emin (); 1317 1.1.1.6 mrg set_emin (65); /* xmin = 2^64 */ 1318 1.1.1.5 mrg mpfr_set_ui_2exp (x, 1, 64, MPFR_RNDN); 1319 1.1.1.5 mrg mpfr_clear_flags (); 1320 1.1.1.5 mrg cmp = mpfr_cmp_z (x, z); 1321 1.1.1.5 mrg MPFR_ASSERTN(cmp > 0); 1322 1.1.1.5 mrg MPFR_ASSERTN(!mpfr_underflow_p ()); 1323 1.1.1.6 mrg set_emin (emin); 1324 1.1.1.5 mrg 1325 1.1.1.5 mrg mpfr_clear (x); 1326 1.1.1.5 mrg mpz_clear (z); 1327 1.1.1.5 mrg } 1328 1.1.1.5 mrg 1329 1.1 mrg int 1330 1.1 mrg main (int argc, char *argv[]) 1331 1.1 mrg { 1332 1.1 mrg tests_start_mpfr (); 1333 1.1 mrg 1334 1.1.1.5 mrg coverage (); 1335 1.1 mrg special (); 1336 1.1 mrg 1337 1.1 mrg test_specialz (mpfr_add_z, mpz_add, "add"); 1338 1.1 mrg test_specialz (mpfr_sub_z, mpz_sub, "sub"); 1339 1.1 mrg test_specialz (mpfr_mul_z, mpz_mul, "mul"); 1340 1.1.1.4 mrg test_genericz (MPFR_PREC_MIN, 100, 100, mpfr_add_z, "add"); 1341 1.1.1.4 mrg test_genericz (MPFR_PREC_MIN, 100, 100, mpfr_sub_z, "sub"); 1342 1.1.1.4 mrg test_genericz (MPFR_PREC_MIN, 100, 100, mpfr_mul_z, "mul"); 1343 1.1.1.4 mrg test_genericz (MPFR_PREC_MIN, 100, 100, mpfr_div_z, "div"); 1344 1.1.1.2 mrg test_special2z (mpfr_z_sub, mpz_sub, "sub"); 1345 1.1.1.4 mrg test_generic2z (MPFR_PREC_MIN, 100, 100, mpfr_z_sub, "sub"); 1346 1.1 mrg 1347 1.1.1.4 mrg test_genericq (MPFR_PREC_MIN, 100, 100, mpfr_add_q, "add"); 1348 1.1.1.4 mrg test_genericq (MPFR_PREC_MIN, 100, 100, mpfr_sub_q, "sub"); 1349 1.1.1.4 mrg test_genericq (MPFR_PREC_MIN, 100, 100, mpfr_mul_q, "mul"); 1350 1.1.1.4 mrg test_genericq (MPFR_PREC_MIN, 100, 100, mpfr_div_q, "div"); 1351 1.1.1.4 mrg test_specialq (MPFR_PREC_MIN, 100, 100, mpfr_mul_q, mpq_mul, "mul"); 1352 1.1.1.4 mrg test_specialq (MPFR_PREC_MIN, 100, 100, mpfr_div_q, mpq_div, "div"); 1353 1.1.1.4 mrg test_specialq (MPFR_PREC_MIN, 100, 100, mpfr_add_q, mpq_add, "add"); 1354 1.1.1.4 mrg test_specialq (MPFR_PREC_MIN, 100, 100, mpfr_sub_q, mpq_sub, "sub"); 1355 1.1.1.4 mrg 1356 1.1.1.4 mrg test_cmp_z (MPFR_PREC_MIN, 100, 100); 1357 1.1.1.4 mrg test_cmp_q (MPFR_PREC_MIN, 100, 100); 1358 1.1.1.4 mrg test_cmp_f (MPFR_PREC_MIN, 100, 100); 1359 1.1 mrg 1360 1.1 mrg check_for_zero (); 1361 1.1 mrg 1362 1.1.1.2 mrg bug_mul_q_20100810 (); 1363 1.1.1.2 mrg bug_div_q_20100810 (); 1364 1.1.1.2 mrg bug_mul_div_q_20100818 (); 1365 1.1.1.2 mrg reduced_expo_range (); 1366 1.1.1.2 mrg addsubq_overflow (); 1367 1.1.1.2 mrg 1368 1.1.1.2 mrg coverage_mpfr_mul_q_20110218 (); 1369 1.1.1.2 mrg 1370 1.1 mrg tests_end_mpfr (); 1371 1.1 mrg return 0; 1372 1.1 mrg } 1373 1.1 mrg 1374 1.1.1.4 mrg #else 1375 1.1.1.4 mrg 1376 1.1.1.4 mrg int 1377 1.1.1.4 mrg main (void) 1378 1.1.1.4 mrg { 1379 1.1.1.4 mrg return 77; 1380 1.1.1.4 mrg } 1381 1.1.1.4 mrg 1382 1.1.1.4 mrg #endif 1383