1 /* Test file for mpfr_div (and some mpfr_div_ui, etc. tests).
2
3 Copyright 1999, 2001-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 static void
check_equal(mpfr_srcptr a,mpfr_srcptr a2,const char * s,mpfr_srcptr b,mpfr_srcptr c,mpfr_rnd_t r)26 check_equal (mpfr_srcptr a, mpfr_srcptr a2, const char *s,
27 mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t r)
28 {
29 if (SAME_VAL (a, a2))
30 return;
31 if (r == MPFR_RNDF) /* RNDF might return different values */
32 return;
33 printf ("Error in %s\n", mpfr_print_rnd_mode (r));
34 printf ("b = ");
35 mpfr_dump (b);
36 printf ("c = ");
37 mpfr_dump (c);
38 printf ("mpfr_div result: ");
39 mpfr_dump (a);
40 printf ("%s result: ", s);
41 mpfr_dump (a2);
42 exit (1);
43 }
44
45 static int
mpfr_all_div(mpfr_ptr a,mpfr_srcptr b,mpfr_srcptr c,mpfr_rnd_t r)46 mpfr_all_div (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t r)
47 {
48 mpfr_t a2;
49 unsigned int oldflags, newflags;
50 int inex, inex2;
51
52 oldflags = __gmpfr_flags;
53 inex = mpfr_div (a, b, c, r);
54
55 /* this test makes no sense for RNDF, since it compares the ternary value
56 and the flags */
57 if (a == b || a == c || r == MPFR_RNDF)
58 return inex;
59
60 newflags = __gmpfr_flags;
61
62 mpfr_init2 (a2, MPFR_PREC (a));
63
64 if (mpfr_integer_p (b) && ! (MPFR_IS_ZERO (b) && MPFR_IS_NEG (b)))
65 {
66 /* b is an integer, but not -0 (-0 is rejected as
67 it becomes +0 when converted to an integer). */
68 if (mpfr_fits_ulong_p (b, MPFR_RNDA))
69 {
70 __gmpfr_flags = oldflags;
71 inex2 = mpfr_ui_div (a2, mpfr_get_ui (b, MPFR_RNDN), c, r);
72 if (!SAME_SIGN (inex2, inex))
73 {
74 printf ("Error for ternary value (rnd=%s), mpfr_div %d, mpfr_ui_div %d\n",
75 mpfr_print_rnd_mode (r), inex, inex2);
76 exit (1);
77 }
78 if (__gmpfr_flags != newflags)
79 {
80 printf ("Error for flags, mpfr_div %d, mpfr_ui_div %d\n",
81 newflags, __gmpfr_flags);
82 exit (1);
83 }
84 check_equal (a, a2, "mpfr_ui_div", b, c, r);
85 }
86 if (mpfr_fits_slong_p (b, MPFR_RNDA))
87 {
88 __gmpfr_flags = oldflags;
89 inex2 = mpfr_si_div (a2, mpfr_get_si (b, MPFR_RNDN), c, r);
90 MPFR_ASSERTN (SAME_SIGN (inex2, inex));
91 MPFR_ASSERTN (__gmpfr_flags == newflags);
92 check_equal (a, a2, "mpfr_si_div", b, c, r);
93 }
94 }
95
96 if (mpfr_integer_p (c) && ! (MPFR_IS_ZERO (c) && MPFR_IS_NEG (c)))
97 {
98 /* c is an integer, but not -0 (-0 is rejected as
99 it becomes +0 when converted to an integer). */
100 if (mpfr_fits_ulong_p (c, MPFR_RNDA))
101 {
102 __gmpfr_flags = oldflags;
103 inex2 = mpfr_div_ui (a2, b, mpfr_get_ui (c, MPFR_RNDN), r);
104 MPFR_ASSERTN (SAME_SIGN (inex2, inex));
105 MPFR_ASSERTN (__gmpfr_flags == newflags);
106 check_equal (a, a2, "mpfr_div_ui", b, c, r);
107 }
108 if (mpfr_fits_slong_p (c, MPFR_RNDA))
109 {
110 __gmpfr_flags = oldflags;
111 inex2 = mpfr_div_si (a2, b, mpfr_get_si (c, MPFR_RNDN), r);
112 MPFR_ASSERTN (SAME_SIGN (inex2, inex));
113 MPFR_ASSERTN (__gmpfr_flags == newflags);
114 check_equal (a, a2, "mpfr_div_si", b, c, r);
115 }
116 }
117
118 mpfr_clear (a2);
119
120 return inex;
121 }
122
123 #ifdef CHECK_EXTERNAL
124 static int
test_div(mpfr_ptr a,mpfr_srcptr b,mpfr_srcptr c,mpfr_rnd_t rnd_mode)125 test_div (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
126 {
127 int res;
128 int ok = rnd_mode == MPFR_RNDN && mpfr_number_p (b) && mpfr_number_p (c);
129 if (ok)
130 {
131 mpfr_print_raw (b);
132 printf (" ");
133 mpfr_print_raw (c);
134 }
135 res = mpfr_all_div (a, b, c, rnd_mode);
136 if (ok)
137 {
138 printf (" ");
139 mpfr_print_raw (a);
140 printf ("\n");
141 }
142 return res;
143 }
144 #else
145 #define test_div mpfr_all_div
146 #endif
147
148 #define check53(n, d, rnd, res) check4(n, d, rnd, 53, res)
149
150 /* return 0 iff a and b are of the same sign */
151 static int
inex_cmp(int a,int b)152 inex_cmp (int a, int b)
153 {
154 if (a > 0)
155 return (b > 0) ? 0 : 1;
156 else if (a == 0)
157 return (b == 0) ? 0 : 1;
158 else
159 return (b < 0) ? 0 : 1;
160 }
161
162 static void
check4(const char * Ns,const char * Ds,mpfr_rnd_t rnd_mode,int p,const char * Qs)163 check4 (const char *Ns, const char *Ds, mpfr_rnd_t rnd_mode, int p,
164 const char *Qs)
165 {
166 mpfr_t q, n, d;
167
168 mpfr_inits2 (p, q, n, d, (mpfr_ptr) 0);
169 mpfr_set_str1 (n, Ns);
170 mpfr_set_str1 (d, Ds);
171 test_div(q, n, d, rnd_mode);
172 if (mpfr_cmp_str (q, Qs, ((p==53) ? 10 : 2), MPFR_RNDN) )
173 {
174 printf ("mpfr_div failed for n=%s, d=%s, p=%d, rnd_mode=%s\n",
175 Ns, Ds, p, mpfr_print_rnd_mode (rnd_mode));
176 printf ("got ");
177 mpfr_dump (q);
178 mpfr_set_str (q, Qs, ((p==53) ? 10 : 2), MPFR_RNDN);
179 printf ("expected ");
180 mpfr_dump (q);
181 exit (1);
182 }
183 mpfr_clears (q, n, d, (mpfr_ptr) 0);
184 }
185
186 static void
check24(const char * Ns,const char * Ds,mpfr_rnd_t rnd_mode,const char * Qs)187 check24 (const char *Ns, const char *Ds, mpfr_rnd_t rnd_mode, const char *Qs)
188 {
189 mpfr_t q, n, d;
190
191 mpfr_inits2 (24, q, n, d, (mpfr_ptr) 0);
192
193 mpfr_set_str1 (n, Ns);
194 mpfr_set_str1 (d, Ds);
195 test_div(q, n, d, rnd_mode);
196 if (mpfr_cmp_str1 (q, Qs) )
197 {
198 printf ("mpfr_div failed for n=%s, d=%s, prec=24, rnd_mode=%s\n",
199 Ns, Ds, mpfr_print_rnd_mode(rnd_mode));
200 printf ("expected quotient is %s, got ", Qs);
201 mpfr_out_str(stdout,10,0,q, MPFR_RNDN); putchar('\n');
202 exit (1);
203 }
204 mpfr_clears (q, n, d, (mpfr_ptr) 0);
205 }
206
207 /* the following examples come from the paper "Number-theoretic Test
208 Generation for Directed Rounding" from Michael Parks, Table 2 */
209 static void
check_float(void)210 check_float(void)
211 {
212 check24("70368760954880.0", "8388609.0", MPFR_RNDN, "8.388609e6");
213 check24("140737479966720.0", "16777213.0", MPFR_RNDN, "8.388609e6");
214 check24("70368777732096.0", "8388611.0", MPFR_RNDN, "8.388609e6");
215 check24("105553133043712.0", "12582911.0", MPFR_RNDN, "8.38861e6");
216 /* the exponent for the following example was forgotten in
217 the Arith'14 version of Parks' paper */
218 check24 ("12582913.0", "12582910.0", MPFR_RNDN, "1.000000238");
219 check24 ("105553124655104.0", "12582910.0", MPFR_RNDN, "8388610.0");
220 check24("140737479966720.0", "8388609.0", MPFR_RNDN, "1.6777213e7");
221 check24("70368777732096.0", "8388609.0", MPFR_RNDN, "8.388611e6");
222 check24("105553133043712.0", "8388610.0", MPFR_RNDN, "1.2582911e7");
223 check24("105553124655104.0", "8388610.0", MPFR_RNDN, "1.258291e7");
224
225 check24("70368760954880.0", "8388609.0", MPFR_RNDZ, "8.388608e6");
226 check24("140737479966720.0", "16777213.0", MPFR_RNDZ, "8.388609e6");
227 check24("70368777732096.0", "8388611.0", MPFR_RNDZ, "8.388608e6");
228 check24("105553133043712.0", "12582911.0", MPFR_RNDZ, "8.38861e6");
229 check24("12582913.0", "12582910.0", MPFR_RNDZ, "1.000000238");
230 check24 ("105553124655104.0", "12582910.0", MPFR_RNDZ, "8388610.0");
231 check24("140737479966720.0", "8388609.0", MPFR_RNDZ, "1.6777213e7");
232 check24("70368777732096.0", "8388609.0", MPFR_RNDZ, "8.38861e6");
233 check24("105553133043712.0", "8388610.0", MPFR_RNDZ, "1.2582911e7");
234 check24("105553124655104.0", "8388610.0", MPFR_RNDZ, "1.258291e7");
235
236 check24("70368760954880.0", "8388609.0", MPFR_RNDU, "8.388609e6");
237 check24("140737479966720.0", "16777213.0", MPFR_RNDU, "8.38861e6");
238 check24("70368777732096.0", "8388611.0", MPFR_RNDU, "8.388609e6");
239 check24("105553133043712.0", "12582911.0", MPFR_RNDU, "8.388611e6");
240 check24("12582913.0", "12582910.0", MPFR_RNDU, "1.000000357");
241 check24 ("105553124655104.0", "12582910.0", MPFR_RNDU, "8388611.0");
242 check24("140737479966720.0", "8388609.0", MPFR_RNDU, "1.6777214e7");
243 check24("70368777732096.0", "8388609.0", MPFR_RNDU, "8.388611e6");
244 check24("105553133043712.0", "8388610.0", MPFR_RNDU, "1.2582912e7");
245 check24("105553124655104.0", "8388610.0", MPFR_RNDU, "1.2582911e7");
246
247 check24("70368760954880.0", "8388609.0", MPFR_RNDD, "8.388608e6");
248 check24("140737479966720.0", "16777213.0", MPFR_RNDD, "8.388609e6");
249 check24("70368777732096.0", "8388611.0", MPFR_RNDD, "8.388608e6");
250 check24("105553133043712.0", "12582911.0", MPFR_RNDD, "8.38861e6");
251 check24("12582913.0", "12582910.0", MPFR_RNDD, "1.000000238");
252 check24 ("105553124655104.0", "12582910.0", MPFR_RNDD, "8388610.0");
253 check24("140737479966720.0", "8388609.0", MPFR_RNDD, "1.6777213e7");
254 check24("70368777732096.0", "8388609.0", MPFR_RNDD, "8.38861e6");
255 check24("105553133043712.0", "8388610.0", MPFR_RNDD, "1.2582911e7");
256 check24("105553124655104.0", "8388610.0", MPFR_RNDD, "1.258291e7");
257
258 check24("70368760954880.0", "8388609.0", MPFR_RNDA, "8.388609e6");
259 }
260
261 static void
check_double(void)262 check_double(void)
263 {
264 check53("0.0", "1.0", MPFR_RNDZ, "0.0");
265 check53("-7.4988969224688591e63", "4.8816866450288732e306", MPFR_RNDD,
266 "-1.5361282826510687291e-243");
267 check53("-1.33225773037748601769e+199", "3.63449540676937123913e+79",
268 MPFR_RNDZ, "-3.6655920045905428978e119");
269 check53("9.89438396044940256501e-134", "5.93472984109987421717e-67",MPFR_RNDU,
270 "1.6672003992376663654e-67");
271 check53("9.89438396044940256501e-134", "5.93472984109987421717e-67",MPFR_RNDA,
272 "1.6672003992376663654e-67");
273 check53("9.89438396044940256501e-134", "-5.93472984109987421717e-67",
274 MPFR_RNDU, "-1.6672003992376663654e-67");
275 check53("-4.53063926135729747564e-308", "7.02293374921793516813e-84",
276 MPFR_RNDD, "-6.4512060388748850857e-225");
277 check53("6.25089225176473806123e-01","-2.35527154824420243364e-230",
278 MPFR_RNDD, "-2.6540006635008291192e229");
279 check53("6.25089225176473806123e-01","-2.35527154824420243364e-230",
280 MPFR_RNDA, "-2.6540006635008291192e229");
281 check53("6.52308934689126e15", "-1.62063546601505417497e273", MPFR_RNDN,
282 "-4.0250194961676020848e-258");
283 check53("1.04636807108079349236e-189", "3.72295730823253012954e-292",
284 MPFR_RNDZ, "2.810583051186143125e102");
285 /* problems found by Kevin under HP-PA */
286 check53 ("2.861044553323177e-136", "-1.1120354257068143e+45", MPFR_RNDZ,
287 "-2.5727998292003016e-181");
288 check53 ("-4.0559157245809205e-127", "-1.1237723844524865e+77", MPFR_RNDN,
289 "3.6091968273068081e-204");
290 check53 ("-1.8177943561493235e-93", "-8.51233984260364e-104", MPFR_RNDU,
291 "2.1354814184595821e+10");
292 }
293
294 static void
check_64(void)295 check_64(void)
296 {
297 mpfr_t x,y,z;
298
299 mpfr_inits2 (64, x, y, z, (mpfr_ptr) 0);
300
301 mpfr_set_str_binary(x, "1.00100100110110101001010010101111000001011100100101010000000000E54");
302 mpfr_set_str_binary(y, "1.00000000000000000000000000000000000000000000000000000000000000E584");
303 test_div(z, x, y, MPFR_RNDU);
304 if (mpfr_cmp_str (z, "0.1001001001101101010010100101011110000010111001001010100000000000E-529", 2, MPFR_RNDN))
305 {
306 printf ("Error for tdiv for MPFR_RNDU and p=64\nx=");
307 mpfr_dump (x);
308 printf ("y=");
309 mpfr_dump (y);
310 printf ("got ");
311 mpfr_dump (z);
312 printf ("expected 0.1001001001101101010010100101011110000010111001001010100000000000E-529\n");
313 exit (1);
314 }
315
316 mpfr_clears (x, y, z, (mpfr_ptr) 0);
317 }
318
319 static void
check_convergence(void)320 check_convergence (void)
321 {
322 mpfr_t x, y; int i, j;
323
324 mpfr_init2(x, 130);
325 mpfr_set_str_binary(x, "0.1011111101011010101000001010011111101000011100011101010011111011000011001010000000111100100111110011001010110100100001001000111001E6944");
326 mpfr_init2(y, 130);
327 mpfr_set_ui(y, 5, MPFR_RNDN);
328 test_div(x, x, y, MPFR_RNDD); /* exact division */
329
330 mpfr_set_prec(x, 64);
331 mpfr_set_prec(y, 64);
332 mpfr_set_str_binary(x, "0.10010010011011010100101001010111100000101110010010101E55");
333 mpfr_set_str_binary(y, "0.1E585");
334 test_div(x, x, y, MPFR_RNDN);
335 mpfr_set_str_binary(y, "0.10010010011011010100101001010111100000101110010010101E-529");
336 if (mpfr_cmp (x, y))
337 {
338 printf ("Error in mpfr_div for prec=64, rnd=MPFR_RNDN\n");
339 printf ("got "); mpfr_dump (x);
340 printf ("instead of "); mpfr_dump (y);
341 exit(1);
342 }
343
344 for (i=32; i<=64; i+=32)
345 {
346 mpfr_set_prec(x, i);
347 mpfr_set_prec(y, i);
348 mpfr_set_ui(x, 1, MPFR_RNDN);
349 RND_LOOP(j)
350 {
351 mpfr_set_ui (y, 1, MPFR_RNDN);
352 test_div (y, x, y, (mpfr_rnd_t) j);
353 if (mpfr_cmp_ui (y, 1))
354 {
355 printf ("mpfr_div failed for x=1.0, y=1.0, prec=%d rnd=%s\n",
356 i, mpfr_print_rnd_mode ((mpfr_rnd_t) j));
357 printf ("got "); mpfr_dump (y);
358 exit (1);
359 }
360 }
361 }
362
363 mpfr_clear (x);
364 mpfr_clear (y);
365 }
366
367 #define KMAX 10000
368
369 /* given y = o(x/u), x, u, find the inexact flag by
370 multiplying y by u */
371 static int
get_inexact(mpfr_ptr y,mpfr_ptr x,mpfr_ptr u)372 get_inexact (mpfr_ptr y, mpfr_ptr x, mpfr_ptr u)
373 {
374 mpfr_t xx;
375 int inex;
376 mpfr_init2 (xx, mpfr_get_prec (y) + mpfr_get_prec (u));
377 mpfr_mul (xx, y, u, MPFR_RNDN); /* exact */
378 inex = mpfr_cmp (xx, x);
379 mpfr_clear (xx);
380 return inex;
381 }
382
383 static void
check_hard(void)384 check_hard (void)
385 {
386 mpfr_t u, v, q, q2;
387 mpfr_prec_t precu, precv, precq;
388 int rnd;
389 int inex, inex2, i, j;
390
391 mpfr_init (q);
392 mpfr_init (q2);
393 mpfr_init (u);
394 mpfr_init (v);
395
396 for (precq = MPFR_PREC_MIN; precq <= 64; precq ++)
397 {
398 mpfr_set_prec (q, precq);
399 mpfr_set_prec (q2, precq + 1);
400 for (j = 0; j < 2; j++)
401 {
402 if (j == 0)
403 {
404 do
405 {
406 mpfr_urandomb (q2, RANDS);
407 }
408 while (mpfr_cmp_ui (q2, 0) == 0);
409 }
410 else /* use q2=1 */
411 mpfr_set_ui (q2, 1, MPFR_RNDN);
412 for (precv = precq; precv <= 10 * precq; precv += precq)
413 {
414 mpfr_set_prec (v, precv);
415 do
416 {
417 mpfr_urandomb (v, RANDS);
418 }
419 while (mpfr_cmp_ui (v, 0) == 0);
420 for (precu = precq; precu <= 10 * precq; precu += precq)
421 {
422 mpfr_set_prec (u, precu);
423 mpfr_mul (u, v, q2, MPFR_RNDN);
424 mpfr_nextbelow (u);
425 for (i = 0; i <= 2; i++)
426 {
427 RND_LOOP_NO_RNDF (rnd)
428 {
429 inex = test_div (q, u, v, (mpfr_rnd_t) rnd);
430 inex2 = get_inexact (q, u, v);
431 if (inex_cmp (inex, inex2))
432 {
433 printf ("Wrong inexact flag for rnd=%s: expected %d, got %d\n",
434 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), inex2, inex);
435 printf ("u= "); mpfr_dump (u);
436 printf ("v= "); mpfr_dump (v);
437 printf ("q= "); mpfr_dump (q);
438 mpfr_set_prec (q2, precq + precv);
439 mpfr_mul (q2, q, v, MPFR_RNDN);
440 printf ("q*v="); mpfr_dump (q2);
441 exit (1);
442 }
443 }
444 mpfr_nextabove (u);
445 }
446 }
447 }
448 }
449 }
450
451 mpfr_clear (q);
452 mpfr_clear (q2);
453 mpfr_clear (u);
454 mpfr_clear (v);
455 }
456
457 static void
check_lowr(void)458 check_lowr (void)
459 {
460 mpfr_t x, y, z, z2, z3, tmp;
461 int k, c, c2;
462
463
464 mpfr_init2 (x, 1000);
465 mpfr_init2 (y, 100);
466 mpfr_init2 (tmp, 850);
467 mpfr_init2 (z, 10);
468 mpfr_init2 (z2, 10);
469 mpfr_init2 (z3, 50);
470
471 for (k = 1; k < KMAX; k++)
472 {
473 do
474 {
475 mpfr_urandomb (z, RANDS);
476 }
477 while (mpfr_cmp_ui (z, 0) == 0);
478 do
479 {
480 mpfr_urandomb (tmp, RANDS);
481 }
482 while (mpfr_cmp_ui (tmp, 0) == 0);
483 mpfr_mul (x, z, tmp, MPFR_RNDN); /* exact */
484 c = test_div (z2, x, tmp, MPFR_RNDN);
485
486 if (c || mpfr_cmp (z2, z))
487 {
488 printf ("Error in mpfr_div rnd=MPFR_RNDN\n");
489 printf ("got "); mpfr_dump (z2);
490 printf ("instead of "); mpfr_dump (z);
491 printf ("inex flag = %d, expected 0\n", c);
492 exit (1);
493 }
494 }
495
496 /* x has still precision 1000, z precision 10, and tmp prec 850 */
497 mpfr_set_prec (z2, 9);
498 for (k = 1; k < KMAX; k++)
499 {
500 mpfr_urandomb (z, RANDS);
501 do
502 {
503 mpfr_urandomb (tmp, RANDS);
504 }
505 while (mpfr_cmp_ui (tmp, 0) == 0);
506 mpfr_mul (x, z, tmp, MPFR_RNDN); /* exact */
507 c = test_div (z2, x, tmp, MPFR_RNDN);
508 /* since z2 has one less bit that z, either the division is exact
509 if z is representable on 9 bits, or we have an even round case */
510
511 c2 = get_inexact (z2, x, tmp);
512 if ((mpfr_cmp (z2, z) == 0 && c) || inex_cmp (c, c2))
513 {
514 printf ("Error in mpfr_div rnd=MPFR_RNDN\n");
515 printf ("got "); mpfr_dump (z2);
516 printf ("instead of "); mpfr_dump (z);
517 printf ("inex flag = %d, expected %d\n", c, c2);
518 exit (1);
519 }
520 else if (c == 2)
521 {
522 mpfr_nexttoinf (z);
523 if (mpfr_cmp(z2, z))
524 {
525 printf ("Error in mpfr_div [even rnd?] rnd=MPFR_RNDN\n");
526 printf ("Dividing ");
527 printf ("got "); mpfr_dump (z2);
528 printf ("instead of "); mpfr_dump (z);
529 printf ("inex flag = %d\n", 1);
530 exit (1);
531 }
532 }
533 else if (c == -2)
534 {
535 mpfr_nexttozero (z);
536 if (mpfr_cmp(z2, z))
537 {
538 printf ("Error in mpfr_div [even rnd?] rnd=MPFR_RNDN\n");
539 printf ("Dividing ");
540 printf ("got "); mpfr_dump (z2);
541 printf ("instead of "); mpfr_dump (z);
542 printf ("inex flag = %d\n", 1);
543 exit (1);
544 }
545 }
546 }
547
548 mpfr_set_prec(x, 1000);
549 mpfr_set_prec(y, 100);
550 mpfr_set_prec(tmp, 850);
551 mpfr_set_prec(z, 10);
552 mpfr_set_prec(z2, 10);
553
554 /* almost exact divisions */
555 for (k = 1; k < KMAX; k++)
556 {
557 do
558 {
559 mpfr_urandomb (z, RANDS);
560 }
561 while (mpfr_cmp_ui (z, 0) == 0);
562 do
563 {
564 mpfr_urandomb (tmp, RANDS);
565 }
566 while (mpfr_cmp_ui (tmp, 0) == 0);
567 mpfr_mul(x, z, tmp, MPFR_RNDN);
568 mpfr_set(y, tmp, MPFR_RNDD);
569 mpfr_nexttoinf (x);
570
571 c = test_div(z2, x, y, MPFR_RNDD);
572 test_div(z3, x, y, MPFR_RNDD);
573 mpfr_set(z, z3, MPFR_RNDD);
574
575 if (c != -1 || mpfr_cmp(z2, z))
576 {
577 printf ("Error in mpfr_div rnd=MPFR_RNDD\n");
578 printf ("got "); mpfr_dump (z2);
579 printf ("instead of "); mpfr_dump (z);
580 printf ("inex flag = %d\n", c);
581 exit (1);
582 }
583
584 mpfr_set (y, tmp, MPFR_RNDU);
585 test_div (z3, x, y, MPFR_RNDU);
586 mpfr_set (z, z3, MPFR_RNDU);
587 c = test_div (z2, x, y, MPFR_RNDU);
588 if (c != 1 || mpfr_cmp (z2, z))
589 {
590 printf ("Error in mpfr_div rnd=MPFR_RNDU\n");
591 printf ("u="); mpfr_dump (x);
592 printf ("v="); mpfr_dump (y);
593 printf ("got "); mpfr_dump (z2);
594 printf ("instead of "); mpfr_dump (z);
595 printf ("inex flag = %d\n", c);
596 exit (1);
597 }
598 }
599
600 mpfr_clear (x);
601 mpfr_clear (y);
602 mpfr_clear (z);
603 mpfr_clear (z2);
604 mpfr_clear (z3);
605 mpfr_clear (tmp);
606 }
607
608 #define MAX_PREC 128
609
610 static void
check_inexact(void)611 check_inexact (void)
612 {
613 mpfr_t x, y, z, u;
614 mpfr_prec_t px, py, pu;
615 int inexact, cmp;
616 mpfr_rnd_t rnd;
617
618 mpfr_init (x);
619 mpfr_init (y);
620 mpfr_init (z);
621 mpfr_init (u);
622
623 mpfr_set_prec (x, 28);
624 mpfr_set_prec (y, 28);
625 mpfr_set_prec (z, 1023);
626 mpfr_set_str_binary (x, "0.1000001001101101111100010011E0");
627 mpfr_set_str (z, "48284762641021308813686974720835219181653367326353400027913400579340343320519877153813133510034402932651132854764198688352364361009429039801248971901380781746767119334993621199563870113045276395603170432175354501451429471578325545278975153148347684600400321033502982713296919861760382863826626093689036010394", 10, MPFR_RNDN);
628 mpfr_div (x, x, z, MPFR_RNDN);
629 mpfr_set_str_binary (y, "0.1111001011001101001001111100E-1023");
630 if (mpfr_cmp (x, y))
631 {
632 printf ("Error in mpfr_div for prec=28, RNDN\n");
633 printf ("Expected "); mpfr_dump (y);
634 printf ("Got "); mpfr_dump (x);
635 exit (1);
636 }
637
638 mpfr_set_prec (x, 53);
639 mpfr_set_str_binary (x, "0.11101100110010100011011000000100001111011111110010101E0");
640 mpfr_set_prec (u, 127);
641 mpfr_set_str_binary (u, "0.1000001100110110110101110110101101111000110000001111111110000000011111001010110100110010111111111101000001011011101011101101000E-2");
642 mpfr_set_prec (y, 95);
643 inexact = test_div (y, x, u, MPFR_RNDN);
644 if (inexact != (cmp = get_inexact (y, x, u)))
645 {
646 printf ("Wrong inexact flag (0): expected %d, got %d\n", cmp, inexact);
647 printf ("x="); mpfr_out_str (stdout, 10, 99, x, MPFR_RNDN); printf ("\n");
648 printf ("u="); mpfr_out_str (stdout, 10, 99, u, MPFR_RNDN); printf ("\n");
649 printf ("y="); mpfr_out_str (stdout, 10, 99, y, MPFR_RNDN); printf ("\n");
650 exit (1);
651 }
652
653 mpfr_set_prec (x, 33);
654 mpfr_set_str_binary (x, "0.101111100011011101010011101100001E0");
655 mpfr_set_prec (u, 2);
656 mpfr_set_str_binary (u, "0.1E0");
657 mpfr_set_prec (y, 28);
658 inexact = test_div (y, x, u, MPFR_RNDN);
659 if (inexact >= 0)
660 {
661 printf ("Wrong inexact flag (1): expected -1, got %d\n",
662 inexact);
663 exit (1);
664 }
665
666 mpfr_set_prec (x, 129);
667 mpfr_set_str_binary (x, "0.111110101111001100000101011100101100110011011101010001000110110101100101000010000001110110100001101010001010100010001111001101010E-2");
668 mpfr_set_prec (u, 15);
669 mpfr_set_str_binary (u, "0.101101000001100E-1");
670 mpfr_set_prec (y, 92);
671 inexact = test_div (y, x, u, MPFR_RNDN);
672 if (inexact <= 0)
673 {
674 printf ("Wrong inexact flag for rnd=MPFR_RNDN(1): expected 1, got %d\n",
675 inexact);
676 mpfr_dump (x);
677 mpfr_dump (u);
678 mpfr_dump (y);
679 exit (1);
680 }
681
682 for (px=2; px<MAX_PREC; px++)
683 {
684 mpfr_set_prec (x, px);
685 mpfr_urandomb (x, RANDS);
686 for (pu=2; pu<=MAX_PREC; pu++)
687 {
688 mpfr_set_prec (u, pu);
689 do { mpfr_urandomb (u, RANDS); } while (mpfr_cmp_ui (u, 0) == 0);
690 {
691 py = MPFR_PREC_MIN + (randlimb () % (MAX_PREC - MPFR_PREC_MIN));
692 mpfr_set_prec (y, py);
693 mpfr_set_prec (z, py + pu);
694 {
695 /* inexact is undefined for RNDF */
696 rnd = RND_RAND_NO_RNDF ();
697 inexact = test_div (y, x, u, rnd);
698 if (mpfr_mul (z, y, u, rnd))
699 {
700 printf ("z <- y * u should be exact\n");
701 exit (1);
702 }
703 cmp = mpfr_cmp (z, x);
704 if (((inexact == 0) && (cmp != 0)) ||
705 ((inexact > 0) && (cmp <= 0)) ||
706 ((inexact < 0) && (cmp >= 0)))
707 {
708 printf ("Wrong inexact flag for rnd=%s\n",
709 mpfr_print_rnd_mode(rnd));
710 printf ("expected %d, got %d\n", cmp, inexact);
711 printf ("x="); mpfr_dump (x);
712 printf ("u="); mpfr_dump (u);
713 printf ("y="); mpfr_dump (y);
714 printf ("y*u="); mpfr_dump (z);
715 exit (1);
716 }
717 }
718 }
719 }
720 }
721
722 mpfr_clear (x);
723 mpfr_clear (y);
724 mpfr_clear (z);
725 mpfr_clear (u);
726 }
727
728 static void
check_special(void)729 check_special (void)
730 {
731 mpfr_t a, d, q;
732 mpfr_exp_t emax, emin;
733 int i;
734
735 mpfr_init2 (a, 100L);
736 mpfr_init2 (d, 100L);
737 mpfr_init2 (q, 100L);
738
739 /* 1/nan == nan */
740 mpfr_set_ui (a, 1L, MPFR_RNDN);
741 MPFR_SET_NAN (d);
742 mpfr_clear_flags ();
743 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
744 MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_NAN);
745
746 /* nan/1 == nan */
747 MPFR_SET_NAN (a);
748 mpfr_set_ui (d, 1L, MPFR_RNDN);
749 mpfr_clear_flags ();
750 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
751 MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_NAN);
752
753 /* +inf/1 == +inf */
754 MPFR_SET_INF (a);
755 MPFR_SET_POS (a);
756 mpfr_set_ui (d, 1L, MPFR_RNDN);
757 mpfr_clear_flags ();
758 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
759 MPFR_ASSERTN (mpfr_inf_p (q));
760 MPFR_ASSERTN (mpfr_sgn (q) > 0);
761 MPFR_ASSERTN (__gmpfr_flags == 0);
762
763 /* +inf/-1 == -inf */
764 MPFR_SET_INF (a);
765 MPFR_SET_POS (a);
766 mpfr_set_si (d, -1, MPFR_RNDN);
767 mpfr_clear_flags ();
768 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
769 MPFR_ASSERTN (mpfr_inf_p (q));
770 MPFR_ASSERTN (mpfr_sgn (q) < 0);
771 MPFR_ASSERTN (__gmpfr_flags == 0);
772
773 /* -inf/1 == -inf */
774 MPFR_SET_INF (a);
775 MPFR_SET_NEG (a);
776 mpfr_set_ui (d, 1L, MPFR_RNDN);
777 mpfr_clear_flags ();
778 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
779 MPFR_ASSERTN (mpfr_inf_p (q));
780 MPFR_ASSERTN (mpfr_sgn (q) < 0);
781 MPFR_ASSERTN (__gmpfr_flags == 0);
782
783 /* -inf/-1 == +inf */
784 MPFR_SET_INF (a);
785 MPFR_SET_NEG (a);
786 mpfr_set_si (d, -1, MPFR_RNDN);
787 mpfr_clear_flags ();
788 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
789 MPFR_ASSERTN (mpfr_inf_p (q));
790 MPFR_ASSERTN (mpfr_sgn (q) > 0);
791 MPFR_ASSERTN (__gmpfr_flags == 0);
792
793 /* 1/+inf == +0 */
794 mpfr_set_ui (a, 1L, MPFR_RNDN);
795 MPFR_SET_INF (d);
796 MPFR_SET_POS (d);
797 mpfr_clear_flags ();
798 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
799 MPFR_ASSERTN (MPFR_IS_ZERO (q));
800 MPFR_ASSERTN (MPFR_IS_POS (q));
801 MPFR_ASSERTN (__gmpfr_flags == 0);
802
803 /* 1/-inf == -0 */
804 mpfr_set_ui (a, 1L, MPFR_RNDN);
805 MPFR_SET_INF (d);
806 MPFR_SET_NEG (d);
807 mpfr_clear_flags ();
808 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
809 MPFR_ASSERTN (MPFR_IS_ZERO (q));
810 MPFR_ASSERTN (MPFR_IS_NEG (q));
811 MPFR_ASSERTN (__gmpfr_flags == 0);
812
813 /* -1/+inf == -0 */
814 mpfr_set_si (a, -1, MPFR_RNDN);
815 MPFR_SET_INF (d);
816 MPFR_SET_POS (d);
817 mpfr_clear_flags ();
818 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
819 MPFR_ASSERTN (MPFR_IS_ZERO (q));
820 MPFR_ASSERTN (MPFR_IS_NEG (q));
821 MPFR_ASSERTN (__gmpfr_flags == 0);
822
823 /* -1/-inf == +0 */
824 mpfr_set_si (a, -1, MPFR_RNDN);
825 MPFR_SET_INF (d);
826 MPFR_SET_NEG (d);
827 mpfr_clear_flags ();
828 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
829 MPFR_ASSERTN (MPFR_IS_ZERO (q));
830 MPFR_ASSERTN (MPFR_IS_POS (q));
831 MPFR_ASSERTN (__gmpfr_flags == 0);
832
833 /* 0/0 == nan */
834 mpfr_set_ui (a, 0L, MPFR_RNDN);
835 mpfr_set_ui (d, 0L, MPFR_RNDN);
836 mpfr_clear_flags ();
837 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
838 MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_NAN);
839
840 /* +inf/+inf == nan */
841 MPFR_SET_INF (a);
842 MPFR_SET_POS (a);
843 MPFR_SET_INF (d);
844 MPFR_SET_POS (d);
845 mpfr_clear_flags ();
846 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
847 MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_NAN);
848
849 /* 1/+0 = +inf */
850 mpfr_set_ui (a, 1, MPFR_RNDZ);
851 mpfr_set_ui (d, 0, MPFR_RNDZ);
852 mpfr_clear_flags ();
853 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
854 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) > 0);
855 MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_DIVBY0);
856
857 /* 1/-0 = -inf */
858 mpfr_set_ui (a, 1, MPFR_RNDZ);
859 mpfr_set_ui (d, 0, MPFR_RNDZ);
860 mpfr_neg (d, d, MPFR_RNDZ);
861 mpfr_clear_flags ();
862 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
863 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) < 0);
864 MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_DIVBY0);
865
866 /* -1/+0 = -inf */
867 mpfr_set_si (a, -1, MPFR_RNDZ);
868 mpfr_set_ui (d, 0, MPFR_RNDZ);
869 mpfr_clear_flags ();
870 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
871 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) < 0);
872 MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_DIVBY0);
873
874 /* -1/-0 = +inf */
875 mpfr_set_si (a, -1, MPFR_RNDZ);
876 mpfr_set_ui (d, 0, MPFR_RNDZ);
877 mpfr_neg (d, d, MPFR_RNDZ);
878 mpfr_clear_flags ();
879 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
880 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) > 0);
881 MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_DIVBY0);
882
883 /* +inf/+0 = +inf */
884 MPFR_SET_INF (a);
885 MPFR_SET_POS (a);
886 mpfr_set_ui (d, 0, MPFR_RNDZ);
887 mpfr_clear_flags ();
888 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
889 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) > 0);
890 MPFR_ASSERTN (__gmpfr_flags == 0);
891
892 /* +inf/-0 = -inf */
893 MPFR_SET_INF (a);
894 MPFR_SET_POS (a);
895 mpfr_set_ui (d, 0, MPFR_RNDZ);
896 mpfr_neg (d, d, MPFR_RNDZ);
897 mpfr_clear_flags ();
898 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
899 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) < 0);
900 MPFR_ASSERTN (__gmpfr_flags == 0);
901
902 /* -inf/+0 = -inf */
903 MPFR_SET_INF (a);
904 MPFR_SET_NEG (a);
905 mpfr_set_ui (d, 0, MPFR_RNDZ);
906 mpfr_clear_flags ();
907 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
908 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) < 0);
909 MPFR_ASSERTN (__gmpfr_flags == 0);
910
911 /* -inf/-0 = +inf */
912 MPFR_SET_INF (a);
913 MPFR_SET_NEG (a);
914 mpfr_set_ui (d, 0, MPFR_RNDZ);
915 mpfr_neg (d, d, MPFR_RNDZ);
916 mpfr_clear_flags ();
917 MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
918 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) > 0);
919 MPFR_ASSERTN (__gmpfr_flags == 0);
920
921 /* check overflow */
922 emax = mpfr_get_emax ();
923 set_emax (1);
924 mpfr_set_ui (a, 1, MPFR_RNDZ);
925 mpfr_set_ui (d, 1, MPFR_RNDZ);
926 mpfr_div_2ui (d, d, 1, MPFR_RNDZ);
927 mpfr_clear_flags ();
928 test_div (q, a, d, MPFR_RNDU); /* 1 / 0.5 = 2 -> overflow */
929 MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) > 0);
930 MPFR_ASSERTN (__gmpfr_flags == (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT));
931 set_emax (emax);
932
933 /* check underflow */
934 emin = mpfr_get_emin ();
935 set_emin (-1);
936 mpfr_set_ui (a, 1, MPFR_RNDZ);
937 mpfr_div_2ui (a, a, 2, MPFR_RNDZ);
938 mpfr_set_prec (d, mpfr_get_prec (q) + 8);
939 for (i = -1; i <= 1; i++)
940 {
941 int sign;
942
943 /* Test 2^(-2) / (+/- (2 + eps)), with eps < 0, eps = 0, eps > 0.
944 -> underflow.
945 With div.c r5513, this test fails for eps > 0 in MPFR_RNDN. */
946 mpfr_set_ui (d, 2, MPFR_RNDZ);
947 if (i < 0)
948 mpfr_nextbelow (d);
949 if (i > 0)
950 mpfr_nextabove (d);
951 for (sign = 0; sign <= 1; sign++)
952 {
953 mpfr_clear_flags ();
954 test_div (q, a, d, MPFR_RNDZ); /* result = 0 */
955 MPFR_ASSERTN (__gmpfr_flags ==
956 (MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT));
957 MPFR_ASSERTN (sign ? MPFR_IS_NEG (q) : MPFR_IS_POS (q));
958 MPFR_ASSERTN (MPFR_IS_ZERO (q));
959 mpfr_clear_flags ();
960 test_div (q, a, d, MPFR_RNDN); /* result = 0 iff eps >= 0 */
961 MPFR_ASSERTN (__gmpfr_flags ==
962 (MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT));
963 MPFR_ASSERTN (sign ? MPFR_IS_NEG (q) : MPFR_IS_POS (q));
964 if (i < 0)
965 mpfr_nexttozero (q);
966 MPFR_ASSERTN (MPFR_IS_ZERO (q));
967 mpfr_neg (d, d, MPFR_RNDN);
968 }
969 }
970 set_emin (emin);
971
972 mpfr_clear (a);
973 mpfr_clear (d);
974 mpfr_clear (q);
975 }
976
977 static void
consistency(void)978 consistency (void)
979 {
980 mpfr_t x, y, z1, z2;
981 int i;
982
983 mpfr_inits (x, y, z1, z2, (mpfr_ptr) 0);
984
985 for (i = 0; i < 10000; i++)
986 {
987 mpfr_rnd_t rnd;
988 mpfr_prec_t px, py, pz, p;
989 int inex1, inex2;
990
991 /* inex is undefined for RNDF */
992 rnd = RND_RAND_NO_RNDF ();
993 px = (randlimb () % 256) + 2;
994 py = (randlimb () % 128) + 2;
995 pz = (randlimb () % 256) + 2;
996 mpfr_set_prec (x, px);
997 mpfr_set_prec (y, py);
998 mpfr_set_prec (z1, pz);
999 mpfr_set_prec (z2, pz);
1000 mpfr_urandomb (x, RANDS);
1001 do
1002 mpfr_urandomb (y, RANDS);
1003 while (mpfr_zero_p (y));
1004 inex1 = mpfr_div (z1, x, y, rnd);
1005 MPFR_ASSERTN (!MPFR_IS_NAN (z1));
1006 p = MAX (MAX (px, py), pz);
1007 if (mpfr_prec_round (x, p, MPFR_RNDN) != 0 ||
1008 mpfr_prec_round (y, p, MPFR_RNDN) != 0)
1009 {
1010 printf ("mpfr_prec_round error for i = %d\n", i);
1011 exit (1);
1012 }
1013 inex2 = mpfr_div (z2, x, y, rnd);
1014 MPFR_ASSERTN (!MPFR_IS_NAN (z2));
1015 if (inex1 != inex2 || mpfr_cmp (z1, z2) != 0)
1016 {
1017 printf ("Consistency error for i = %d, rnd = %s\n", i,
1018 mpfr_print_rnd_mode (rnd));
1019 printf ("inex1=%d inex2=%d\n", inex1, inex2);
1020 printf ("z1="); mpfr_dump (z1);
1021 printf ("z2="); mpfr_dump (z2);
1022 exit (1);
1023 }
1024 }
1025
1026 mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0);
1027 }
1028
1029 /* Reported by Carl Witty on 2007-06-03 */
1030 static void
test_20070603(void)1031 test_20070603 (void)
1032 {
1033 mpfr_t n, d, q, c;
1034
1035 mpfr_init2 (n, 128);
1036 mpfr_init2 (d, 128);
1037 mpfr_init2 (q, 31);
1038 mpfr_init2 (c, 31);
1039
1040 mpfr_set_str (n, "10384593717069655257060992206846485", 10, MPFR_RNDN);
1041 mpfr_set_str (d, "10384593717069655257060992206847132", 10, MPFR_RNDN);
1042 mpfr_div (q, n, d, MPFR_RNDU);
1043
1044 mpfr_set_ui (c, 1, MPFR_RNDN);
1045 if (mpfr_cmp (q, c) != 0)
1046 {
1047 printf ("Error in test_20070603\nGot ");
1048 mpfr_dump (q);
1049 printf ("instead of ");
1050 mpfr_dump (c);
1051 exit (1);
1052 }
1053
1054 /* same for 64-bit machines */
1055 mpfr_set_prec (n, 256);
1056 mpfr_set_prec (d, 256);
1057 mpfr_set_prec (q, 63);
1058 mpfr_set_str (n, "822752278660603021077484591278675252491367930877209729029898240", 10, MPFR_RNDN);
1059 mpfr_set_str (d, "822752278660603021077484591278675252491367930877212507873738752", 10, MPFR_RNDN);
1060 mpfr_div (q, n, d, MPFR_RNDU);
1061 if (mpfr_cmp (q, c) != 0)
1062 {
1063 printf ("Error in test_20070603\nGot ");
1064 mpfr_dump (q);
1065 printf ("instead of ");
1066 mpfr_dump (c);
1067 exit (1);
1068 }
1069
1070 mpfr_clear (n);
1071 mpfr_clear (d);
1072 mpfr_clear (q);
1073 mpfr_clear (c);
1074 }
1075
1076 /* Bug found while adding tests for mpfr_cot */
1077 static void
test_20070628(void)1078 test_20070628 (void)
1079 {
1080 mpfr_exp_t old_emax;
1081 mpfr_t x, y;
1082 int inex, err = 0;
1083
1084 old_emax = mpfr_get_emax ();
1085 set_emax (256);
1086
1087 mpfr_inits2 (53, x, y, (mpfr_ptr) 0);
1088 mpfr_set_si (x, -1, MPFR_RNDN);
1089 mpfr_set_si_2exp (y, 1, -256, MPFR_RNDN);
1090 mpfr_clear_flags ();
1091 inex = mpfr_div (x, x, y, MPFR_RNDD);
1092 if (MPFR_IS_POS (x) || ! mpfr_inf_p (x))
1093 {
1094 printf ("Error in test_20070628: expected -Inf, got\n");
1095 mpfr_dump (x);
1096 err++;
1097 }
1098 if (inex >= 0)
1099 {
1100 printf ("Error in test_20070628: expected inex < 0, got %d\n", inex);
1101 err++;
1102 }
1103 if (! mpfr_overflow_p ())
1104 {
1105 printf ("Error in test_20070628: overflow flag is not set\n");
1106 err++;
1107 }
1108 mpfr_clears (x, y, (mpfr_ptr) 0);
1109 set_emax (old_emax);
1110 }
1111
1112 /* Bug in mpfr_divhigh_n_basecase when all limbs of q (except the most
1113 significant one) are B-1 where B=2^GMP_NUMB_BITS. Since we truncate
1114 the divisor at each step, it might happen at some point that
1115 (np[n-1],np[n-2]) > (d1,d0), and not only the equality.
1116 Reported by Ricky Farr
1117 <https://sympa.inria.fr/sympa/arc/mpfr/2015-10/msg00023.html>
1118 To get a failure, a MPFR_DIVHIGH_TAB entry below the MPFR_DIV_THRESHOLD
1119 limit must have a value 0. With most mparam.h files, this cannot occur. To
1120 make the bug appear, one can configure MPFR with -DMPFR_TUNE_COVERAGE. */
1121 static void
test_20151023(void)1122 test_20151023 (void)
1123 {
1124 mpfr_prec_t p;
1125 mpfr_t n, d, q, q0;
1126 int inex, i;
1127
1128 for (p = GMP_NUMB_BITS; p <= 2000; p++)
1129 {
1130 mpfr_init2 (n, 2*p);
1131 mpfr_init2 (d, p);
1132 mpfr_init2 (q, p);
1133 mpfr_init2 (q0, GMP_NUMB_BITS);
1134
1135 /* generate a random divisor of p bits */
1136 do
1137 mpfr_urandomb (d, RANDS);
1138 while (mpfr_zero_p (d));
1139 /* generate a random non-zero quotient of GMP_NUMB_BITS bits */
1140 do
1141 mpfr_urandomb (q0, RANDS);
1142 while (mpfr_zero_p (q0));
1143 /* zero-pad the quotient to p bits */
1144 inex = mpfr_prec_round (q0, p, MPFR_RNDN);
1145 MPFR_ASSERTN(inex == 0);
1146
1147 for (i = 0; i < 3; i++)
1148 {
1149 /* i=0: try with the original quotient xxx000...000
1150 i=1: try with the original quotient minus one ulp
1151 i=2: try with the original quotient plus one ulp */
1152 if (i == 1)
1153 mpfr_nextbelow (q0);
1154 else if (i == 2)
1155 {
1156 mpfr_nextabove (q0);
1157 mpfr_nextabove (q0);
1158 }
1159
1160 inex = mpfr_mul (n, d, q0, MPFR_RNDN);
1161 MPFR_ASSERTN(inex == 0);
1162 mpfr_nextabove (n);
1163 mpfr_div (q, n, d, MPFR_RNDN);
1164 if (! mpfr_equal_p (q, q0))
1165 {
1166 printf ("Error in test_20151023 for p=%ld, rnd=RNDN, i=%d\n",
1167 (long) p, i);
1168 printf ("n="); mpfr_dump (n);
1169 printf ("d="); mpfr_dump (d);
1170 printf ("expected q0="); mpfr_dump (q0);
1171 printf ("got q="); mpfr_dump (q);
1172 exit (1);
1173 }
1174
1175 inex = mpfr_mul (n, d, q0, MPFR_RNDN);
1176 MPFR_ASSERTN(inex == 0);
1177 mpfr_nextbelow (n);
1178 mpfr_div (q, n, d, MPFR_RNDN);
1179 MPFR_ASSERTN(mpfr_cmp (q, q0) == 0);
1180 }
1181
1182 mpfr_clear (n);
1183 mpfr_clear (d);
1184 mpfr_clear (q);
1185 mpfr_clear (q0);
1186 }
1187 }
1188
1189 /* test a random division of p+extra bits divided by p+extra bits,
1190 with quotient of p bits only, where the p+extra bit approximation
1191 of the quotient is very near a rounding frontier. */
1192 static void
test_bad_aux(mpfr_prec_t p,mpfr_prec_t extra)1193 test_bad_aux (mpfr_prec_t p, mpfr_prec_t extra)
1194 {
1195 mpfr_t u, v, w, q0, q;
1196
1197 mpfr_init2 (u, p + extra);
1198 mpfr_init2 (v, p + extra);
1199 mpfr_init2 (w, p + extra);
1200 mpfr_init2 (q0, p);
1201 mpfr_init2 (q, p);
1202 do mpfr_urandomb (q0, RANDS); while (mpfr_zero_p (q0));
1203 do mpfr_urandomb (v, RANDS); while (mpfr_zero_p (v));
1204
1205 mpfr_set (w, q0, MPFR_RNDN); /* exact */
1206 mpfr_nextabove (w); /* now w > q0 */
1207 mpfr_mul (u, v, w, MPFR_RNDU); /* thus u > v*q0 */
1208 mpfr_div (q, u, v, MPFR_RNDU); /* should have q > q0 */
1209 MPFR_ASSERTN (mpfr_cmp (q, q0) > 0);
1210 mpfr_div (q, u, v, MPFR_RNDZ); /* should have q = q0 */
1211 MPFR_ASSERTN (mpfr_cmp (q, q0) == 0);
1212
1213 mpfr_set (w, q0, MPFR_RNDN); /* exact */
1214 mpfr_nextbelow (w); /* now w < q0 */
1215 mpfr_mul (u, v, w, MPFR_RNDZ); /* thus u < v*q0 */
1216 mpfr_div (q, u, v, MPFR_RNDZ); /* should have q < q0 */
1217 MPFR_ASSERTN (mpfr_cmp (q, q0) < 0);
1218 mpfr_div (q, u, v, MPFR_RNDU); /* should have q = q0 */
1219 MPFR_ASSERTN (mpfr_cmp (q, q0) == 0);
1220
1221 mpfr_clear (u);
1222 mpfr_clear (v);
1223 mpfr_clear (w);
1224 mpfr_clear (q0);
1225 mpfr_clear (q);
1226 }
1227
1228 static void
test_bad(void)1229 test_bad (void)
1230 {
1231 mpfr_prec_t p, extra;
1232
1233 for (p = MPFR_PREC_MIN; p <= 1024; p += 17)
1234 for (extra = 2; extra <= 64; extra++)
1235 test_bad_aux (p, extra);
1236 }
1237
1238 #define TEST_FUNCTION test_div
1239 #define TWO_ARGS
1240 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), randlimb () % 100, RANDS)
1241 #include "tgeneric.c"
1242
1243 static void
test_extreme(void)1244 test_extreme (void)
1245 {
1246 mpfr_t x, y, z;
1247 mpfr_exp_t emin, emax;
1248 mpfr_prec_t p[4] = { 8, 32, 64, 256 };
1249 int xi, yi, zi, j, r;
1250 unsigned int flags, ex_flags;
1251
1252 emin = mpfr_get_emin ();
1253 emax = mpfr_get_emax ();
1254
1255 set_emin (MPFR_EMIN_MIN);
1256 set_emax (MPFR_EMAX_MAX);
1257
1258 for (xi = 0; xi < 4; xi++)
1259 {
1260 mpfr_init2 (x, p[xi]);
1261 mpfr_setmax (x, MPFR_EMAX_MAX);
1262 MPFR_ASSERTN (mpfr_check (x));
1263 for (yi = 0; yi < 4; yi++)
1264 {
1265 mpfr_init2 (y, p[yi]);
1266 mpfr_setmin (y, MPFR_EMIN_MIN);
1267 for (j = 0; j < 2; j++)
1268 {
1269 MPFR_ASSERTN (mpfr_check (y));
1270 for (zi = 0; zi < 4; zi++)
1271 {
1272 mpfr_init2 (z, p[zi]);
1273 RND_LOOP (r)
1274 {
1275 mpfr_clear_flags ();
1276 mpfr_div (z, x, y, (mpfr_rnd_t) r);
1277 flags = __gmpfr_flags;
1278 MPFR_ASSERTN (mpfr_check (z));
1279 ex_flags = MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT;
1280 if (flags != ex_flags)
1281 {
1282 printf ("Bad flags in test_extreme on z = a/b"
1283 " with %s and\n",
1284 mpfr_print_rnd_mode ((mpfr_rnd_t) r));
1285 printf ("a = ");
1286 mpfr_dump (x);
1287 printf ("b = ");
1288 mpfr_dump (y);
1289 printf ("Expected flags:");
1290 flags_out (ex_flags);
1291 printf ("Got flags: ");
1292 flags_out (flags);
1293 printf ("z = ");
1294 mpfr_dump (z);
1295 exit (1);
1296 }
1297 mpfr_clear_flags ();
1298 mpfr_div (z, y, x, (mpfr_rnd_t) r);
1299 flags = __gmpfr_flags;
1300 MPFR_ASSERTN (mpfr_check (z));
1301 ex_flags = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT;
1302 if (flags != ex_flags)
1303 {
1304 printf ("Bad flags in test_extreme on z = a/b"
1305 " with %s and\n",
1306 mpfr_print_rnd_mode ((mpfr_rnd_t) r));
1307 printf ("a = ");
1308 mpfr_dump (y);
1309 printf ("b = ");
1310 mpfr_dump (x);
1311 printf ("Expected flags:");
1312 flags_out (ex_flags);
1313 printf ("Got flags: ");
1314 flags_out (flags);
1315 printf ("z = ");
1316 mpfr_dump (z);
1317 exit (1);
1318 }
1319 }
1320 mpfr_clear (z);
1321 } /* zi */
1322 mpfr_nextabove (y);
1323 } /* j */
1324 mpfr_clear (y);
1325 } /* yi */
1326 mpfr_clear (x);
1327 } /* xi */
1328
1329 set_emin (emin);
1330 set_emax (emax);
1331 }
1332
1333 static void
testall_rndf(mpfr_prec_t pmax)1334 testall_rndf (mpfr_prec_t pmax)
1335 {
1336 mpfr_t a, b, c, d;
1337 mpfr_prec_t pa, pb, pc;
1338
1339 for (pa = MPFR_PREC_MIN; pa <= pmax; pa++)
1340 {
1341 mpfr_init2 (a, pa);
1342 mpfr_init2 (d, pa);
1343 for (pb = MPFR_PREC_MIN; pb <= pmax; pb++)
1344 {
1345 mpfr_init2 (b, pb);
1346 mpfr_set_ui (b, 1, MPFR_RNDN);
1347 while (mpfr_cmp_ui (b, 2) < 0)
1348 {
1349 for (pc = MPFR_PREC_MIN; pc <= pmax; pc++)
1350 {
1351 mpfr_init2 (c, pc);
1352 mpfr_set_ui (c, 1, MPFR_RNDN);
1353 while (mpfr_cmp_ui (c, 2) < 0)
1354 {
1355 mpfr_div (a, b, c, MPFR_RNDF);
1356 mpfr_div (d, b, c, MPFR_RNDD);
1357 if (!mpfr_equal_p (a, d))
1358 {
1359 mpfr_div (d, b, c, MPFR_RNDU);
1360 if (!mpfr_equal_p (a, d))
1361 {
1362 printf ("Error: mpfr_div(a,b,c,RNDF) does not "
1363 "match RNDD/RNDU\n");
1364 printf ("b="); mpfr_dump (b);
1365 printf ("c="); mpfr_dump (c);
1366 printf ("a="); mpfr_dump (a);
1367 exit (1);
1368 }
1369 }
1370 mpfr_nextabove (c);
1371 }
1372 mpfr_clear (c);
1373 }
1374 mpfr_nextabove (b);
1375 }
1376 mpfr_clear (b);
1377 }
1378 mpfr_clear (a);
1379 mpfr_clear (d);
1380 }
1381 }
1382
1383 static void
test_mpfr_divsp2(void)1384 test_mpfr_divsp2 (void)
1385 {
1386 mpfr_t u, v, q;
1387
1388 /* test to exercise r2 = v1 in mpfr_divsp2 */
1389 mpfr_init2 (u, 128);
1390 mpfr_init2 (v, 128);
1391 mpfr_init2 (q, 83);
1392
1393 mpfr_set_str (u, "286677858044426991425771529092412636160", 10, MPFR_RNDN);
1394 mpfr_set_str (v, "241810647971575979588130185988987264768", 10, MPFR_RNDN);
1395 mpfr_div (q, u, v, MPFR_RNDN);
1396 mpfr_set_str (u, "5732952910203749289426944", 10, MPFR_RNDN);
1397 mpfr_div_2ui (u, u, 82, MPFR_RNDN);
1398 MPFR_ASSERTN(mpfr_equal_p (q, u));
1399
1400 mpfr_clear (u);
1401 mpfr_clear (v);
1402 mpfr_clear (q);
1403 }
1404
1405 /* Assertion failure in r10769 with --enable-assert --enable-gmp-internals
1406 (same failure in tatan on a similar example). */
1407 static void
test_20160831(void)1408 test_20160831 (void)
1409 {
1410 mpfr_t u, v, q;
1411
1412 mpfr_inits2 (124, u, v, q, (mpfr_ptr) 0);
1413
1414 mpfr_set_ui (u, 1, MPFR_RNDN);
1415 mpfr_set_str (v, "0x40000000000000005", 16, MPFR_RNDN);
1416 mpfr_div (q, u, v, MPFR_RNDN);
1417 mpfr_set_str (u, "0xfffffffffffffffecp-134", 16, MPFR_RNDN);
1418 MPFR_ASSERTN (mpfr_equal_p (q, u));
1419
1420 mpfr_set_prec (u, 128);
1421 mpfr_set_prec (v, 128);
1422 mpfr_set_str (u, "186127091671619245460026015088243485690", 10, MPFR_RNDN);
1423 mpfr_set_str (v, "205987256581218236405412302590110119580", 10, MPFR_RNDN);
1424 mpfr_div (q, u, v, MPFR_RNDN);
1425 mpfr_set_str (u, "19217137613667309953639458782352244736", 10, MPFR_RNDN);
1426 mpfr_div_2ui (u, u, 124, MPFR_RNDN);
1427 MPFR_ASSERTN (mpfr_equal_p (q, u));
1428
1429 mpfr_clears (u, v, q, (mpfr_ptr) 0);
1430 }
1431
1432 /* With r11138, on a 64-bit machine:
1433 div.c:128: MPFR assertion failed: qx >= __gmpfr_emin
1434 */
1435 static void
test_20170104(void)1436 test_20170104 (void)
1437 {
1438 mpfr_t u, v, q;
1439 mpfr_exp_t emin;
1440
1441 emin = mpfr_get_emin ();
1442 set_emin (-42);
1443
1444 mpfr_init2 (u, 12);
1445 mpfr_init2 (v, 12);
1446 mpfr_init2 (q, 11);
1447 mpfr_set_str_binary (u, "0.111111111110E-29");
1448 mpfr_set_str_binary (v, "0.111111111111E14");
1449 mpfr_div (q, u, v, MPFR_RNDN);
1450 mpfr_clears (u, v, q, (mpfr_ptr) 0);
1451
1452 set_emin (emin);
1453 }
1454
1455 /* With r11140, on a 64-bit machine with GMP_CHECK_RANDOMIZE=1484406128:
1456 Consistency error for i = 2577
1457 */
1458 static void
test_20170105(void)1459 test_20170105 (void)
1460 {
1461 mpfr_t x, y, z, t;
1462
1463 mpfr_init2 (x, 138);
1464 mpfr_init2 (y, 6);
1465 mpfr_init2 (z, 128);
1466 mpfr_init2 (t, 128);
1467 mpfr_set_str_binary (x, "0.100110111001001000101111010010011101111110111111110001110100000001110111010100111010100011101010110000010100000011100100110101101011000000E-6");
1468 mpfr_set_str_binary (y, "0.100100E-2");
1469 /* up to exponents, x/y is exactly 367625553447399614694201910705139062483,
1470 which has 129 bits, thus we are in the round-to-nearest-even case, and since
1471 the penultimate bit of x/y is 1, we should round upwards */
1472 mpfr_set_str_binary (t, "0.10001010010010010000110110010110111111111100011011101010000000000110101000010001011110011011010000111010000000001100101101101010E-3");
1473 mpfr_div (z, x, y, MPFR_RNDN);
1474 MPFR_ASSERTN(mpfr_equal_p (z, t));
1475 mpfr_clears (x, y, z, t, (mpfr_ptr) 0);
1476 }
1477
1478 /* The real cause of the mpfr_ttanh failure from the non-regression test
1479 added in tests/ttanh.c@11993 was actually due to a bug in mpfr_div, as
1480 this can be seen by comparing the logs of the 3.1 branch and the trunk
1481 r11992 with MPFR_LOG_ALL=1 MPFR_LOG_PREC=50 on this particular test
1482 (this was noticed because adding this test to the 3.1 branch did not
1483 yield a failure like in the trunk, though the mpfr_ttanh code did not
1484 change until r11993). This was the bug actually fixed in r12002.
1485 */
1486 static void
test_20171219(void)1487 test_20171219 (void)
1488 {
1489 mpfr_t x, y, z, t;
1490
1491 mpfr_inits2 (126, x, y, z, t, (mpfr_ptr) 0);
1492 mpfr_set_str_binary (x, "0.111000000000000111100000000000011110000000000001111000000000000111100000000000011110000000000001111000000000000111100000000000E1");
1493 /* x = 36347266450315671364380109803814927 / 2^114 */
1494 mpfr_add_ui (y, x, 2, MPFR_RNDN);
1495 /* y = 77885641318594292392624080437575695 / 2^114 */
1496 mpfr_div (z, x, y, MPFR_RNDN);
1497 mpfr_set_ui_2exp (t, 3823, -13, MPFR_RNDN);
1498 MPFR_ASSERTN (mpfr_equal_p (z, t));
1499 mpfr_clears (x, y, z, t, (mpfr_ptr) 0);
1500 }
1501
1502 #if !defined(MPFR_GENERIC_ABI) && GMP_NUMB_BITS == 64
1503 /* exercise mpfr_div2_approx */
1504 static void
test_mpfr_div2_approx(unsigned long n)1505 test_mpfr_div2_approx (unsigned long n)
1506 {
1507 mpfr_t x, y, z;
1508
1509 mpfr_init2 (x, 113);
1510 mpfr_init2 (y, 113);
1511 mpfr_init2 (z, 113);
1512 while (n--)
1513 {
1514 mpfr_urandomb (x, RANDS);
1515 mpfr_urandomb (y, RANDS);
1516 mpfr_div (z, x, y, MPFR_RNDN);
1517 }
1518 mpfr_clear (x);
1519 mpfr_clear (y);
1520 mpfr_clear (z);
1521 }
1522 #endif
1523
1524 /* bug found in ttan with GMP_CHECK_RANDOMIZE=1514257254 */
1525 static void
bug20171218(void)1526 bug20171218 (void)
1527 {
1528 mpfr_t s, c;
1529 mpfr_init2 (s, 124);
1530 mpfr_init2 (c, 124);
1531 mpfr_set_str_binary (s, "-0.1110000111100001111000011110000111100001111000011110000111100001111000011110000111100001111000011110000111100001111000011110E0");
1532 mpfr_set_str_binary (c, "0.1111000011110000111100001111000011110000111100001111000011110000111100001111000011110000111100001111000011110000111100001111E-1");
1533 mpfr_div (c, s, c, MPFR_RNDN);
1534 mpfr_set_str_binary (s, "-1.111000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000");
1535 MPFR_ASSERTN(mpfr_equal_p (c, s));
1536 mpfr_clear (s);
1537 mpfr_clear (c);
1538 }
1539
1540 /* Extended test based on a bug found with flint-arb test suite with a
1541 32-bit ABI: https://bugs.debian.org/cgi-bin/bugreport.cgi?bug=888459
1542 Division of the form: (1 - 2^(-pa)) / (1 - 2^(-pb)).
1543 The result is compared to the one obtained by increasing the precision of
1544 the divisor (without changing its value, so that both results should be
1545 equal). For all of these tests, a failure may occur in r12126 only with
1546 pb=GMP_NUMB_BITS and MPFR_RNDN (and some particular values of pa and pc).
1547 This bug was introduced by r9086, where mpfr_div uses mpfr_div_ui when
1548 the divisor has only one limb.
1549 */
1550 static void
bug20180126(void)1551 bug20180126 (void)
1552 {
1553 mpfr_t a, b1, b2, c1, c2;
1554 int pa, i, j, pc, sa, sb, r, inex1, inex2;
1555
1556 for (pa = 100; pa < 800; pa += 11)
1557 for (i = 1; i <= 4; i++)
1558 for (j = -2; j <= 2; j++)
1559 {
1560 int pb = GMP_NUMB_BITS * i + j;
1561
1562 if (pb > pa)
1563 continue;
1564
1565 mpfr_inits2 (pa, a, b1, (mpfr_ptr) 0);
1566 mpfr_inits2 (pb, b2, (mpfr_ptr) 0);
1567
1568 mpfr_set_ui (a, 1, MPFR_RNDN);
1569 mpfr_nextbelow (a); /* 1 - 2^(-pa) */
1570 mpfr_set_ui (b2, 1, MPFR_RNDN);
1571 mpfr_nextbelow (b2); /* 1 - 2^(-pb) */
1572 inex1 = mpfr_set (b1, b2, MPFR_RNDN);
1573 MPFR_ASSERTN (inex1 == 0);
1574
1575 for (pc = 32; pc <= 320; pc += 32)
1576 {
1577 mpfr_inits2 (pc, c1, c2, (mpfr_ptr) 0);
1578
1579 for (sa = 0; sa < 2; sa++)
1580 {
1581 for (sb = 0; sb < 2; sb++)
1582 {
1583 RND_LOOP_NO_RNDF (r)
1584 {
1585 MPFR_ASSERTN (mpfr_equal_p (b1, b2));
1586 inex1 = mpfr_div (c1, a, b1, (mpfr_rnd_t) r);
1587 inex2 = mpfr_div (c2, a, b2, (mpfr_rnd_t) r);
1588
1589 if (! mpfr_equal_p (c1, c2) ||
1590 ! SAME_SIGN (inex1, inex2))
1591 {
1592 printf ("Error in bug20180126 for "
1593 "pa=%d pb=%d pc=%d sa=%d sb=%d %s\n",
1594 pa, pb, pc, sa, sb,
1595 mpfr_print_rnd_mode ((mpfr_rnd_t) r));
1596 printf ("inex1 = %d, c1 = ", inex1);
1597 mpfr_dump (c1);
1598 printf ("inex2 = %d, c2 = ", inex2);
1599 mpfr_dump (c2);
1600 exit (1);
1601 }
1602 }
1603
1604 mpfr_neg (b1, b1, MPFR_RNDN);
1605 mpfr_neg (b2, b2, MPFR_RNDN);
1606 } /* sb */
1607
1608 mpfr_neg (a, a, MPFR_RNDN);
1609 } /* sa */
1610
1611 mpfr_clears (c1, c2, (mpfr_ptr) 0);
1612 } /* pc */
1613
1614 mpfr_clears (a, b1, b2, (mpfr_ptr) 0);
1615 } /* j */
1616 }
1617
1618 static void
coverage(mpfr_prec_t pmax)1619 coverage (mpfr_prec_t pmax)
1620 {
1621 mpfr_prec_t p;
1622
1623 for (p = MPFR_PREC_MIN; p <= pmax; p++)
1624 {
1625 int inex;
1626 mpfr_t q, u, v;
1627
1628 mpfr_init2 (q, p);
1629 mpfr_init2 (u, p);
1630 mpfr_init2 (v, p);
1631
1632 /* exercise case qx < emin */
1633 mpfr_set_ui_2exp (u, 1, mpfr_get_emin (), MPFR_RNDN);
1634 mpfr_set_ui (v, 4, MPFR_RNDN);
1635
1636 mpfr_clear_flags ();
1637 /* u/v = 2^(emin-2), should be rounded to +0 for RNDN */
1638 inex = mpfr_div (q, u, v, MPFR_RNDN);
1639 MPFR_ASSERTN(inex < 0);
1640 MPFR_ASSERTN(mpfr_zero_p (q) && mpfr_signbit (q) == 0);
1641 MPFR_ASSERTN(mpfr_underflow_p ());
1642
1643 mpfr_clear_flags ();
1644 /* u/v = 2^(emin-2), should be rounded to 2^(emin-1) for RNDU */
1645 inex = mpfr_div (q, u, v, MPFR_RNDU);
1646 MPFR_ASSERTN(inex > 0);
1647 MPFR_ASSERTN(mpfr_cmp_ui_2exp (q, 1, mpfr_get_emin () - 1) == 0);
1648 MPFR_ASSERTN(mpfr_underflow_p ());
1649
1650 mpfr_clear_flags ();
1651 /* u/v = 2^(emin-2), should be rounded to +0 for RNDZ */
1652 inex = mpfr_div (q, u, v, MPFR_RNDZ);
1653 MPFR_ASSERTN(inex < 0);
1654 MPFR_ASSERTN(mpfr_zero_p (q) && mpfr_signbit (q) == 0);
1655 MPFR_ASSERTN(mpfr_underflow_p ());
1656
1657 if (p == 1)
1658 goto end_of_loop;
1659
1660 mpfr_set_ui_2exp (u, 1, mpfr_get_emin (), MPFR_RNDN);
1661 mpfr_nextbelow (u); /* u = (1-2^(-p))*2^emin */
1662 mpfr_set_ui (v, 2, MPFR_RNDN);
1663
1664 mpfr_clear_flags ();
1665 /* u/v = (1-2^(-p))*2^(emin-1), will round to 2^(emin-1) for RNDN */
1666 inex = mpfr_div (q, u, v, MPFR_RNDN);
1667 MPFR_ASSERTN(inex > 0);
1668 MPFR_ASSERTN(mpfr_cmp_ui_2exp (q, 1, mpfr_get_emin () - 1) == 0);
1669 MPFR_ASSERTN(mpfr_underflow_p ());
1670
1671 mpfr_clear_flags ();
1672 /* u/v should round to 2^(emin-1) for RNDU */
1673 inex = mpfr_div (q, u, v, MPFR_RNDU);
1674 MPFR_ASSERTN(inex > 0);
1675 MPFR_ASSERTN(mpfr_cmp_ui_2exp (q, 1, mpfr_get_emin () - 1) == 0);
1676 MPFR_ASSERTN(mpfr_underflow_p ());
1677
1678 mpfr_clear_flags ();
1679 /* u/v should round to +0 for RNDZ */
1680 inex = mpfr_div (q, u, v, MPFR_RNDZ);
1681 MPFR_ASSERTN(inex < 0);
1682 MPFR_ASSERTN(mpfr_zero_p (q) && mpfr_signbit (q) == 0);
1683 MPFR_ASSERTN(mpfr_underflow_p ());
1684
1685 end_of_loop:
1686 mpfr_clear (q);
1687 mpfr_clear (u);
1688 mpfr_clear (v);
1689 }
1690 }
1691
1692 /* coverage for case usize >= n + n in Mulders' algorithm */
1693 static void
coverage2(void)1694 coverage2 (void)
1695 {
1696 mpfr_prec_t p;
1697 mpfr_t q, u, v, t, w;
1698 int inex, inex2;
1699
1700 p = MPFR_DIV_THRESHOLD * GMP_NUMB_BITS;
1701 mpfr_init2 (q, p);
1702 mpfr_init2 (u, 2 * p + 3 * GMP_NUMB_BITS);
1703 mpfr_init2 (v, p);
1704 do mpfr_urandomb (u, RANDS); while (mpfr_zero_p (u));
1705 do mpfr_urandomb (v, RANDS); while (mpfr_zero_p (v));
1706 inex = mpfr_div (q, u, v, MPFR_RNDN);
1707 mpfr_init2 (t, mpfr_get_prec (u));
1708 mpfr_init2 (w, mpfr_get_prec (u));
1709 inex2 = mpfr_mul (t, q, v, MPFR_RNDN);
1710 MPFR_ASSERTN(inex2 == 0);
1711 if (inex == 0) /* check q*v = u */
1712 MPFR_ASSERTN(mpfr_equal_p (u, t));
1713 else
1714 {
1715 if (inex > 0)
1716 mpfr_nextbelow (q);
1717 else
1718 mpfr_nextabove (q);
1719 inex2 = mpfr_mul (w, q, v, MPFR_RNDN);
1720 MPFR_ASSERTN(inex2 == 0);
1721 inex2 = mpfr_sub (t, t, u, MPFR_RNDN);
1722 MPFR_ASSERTN(inex2 == 0);
1723 inex2 = mpfr_sub (w, w, u, MPFR_RNDN);
1724 MPFR_ASSERTN(inex2 == 0);
1725 MPFR_ASSERTN(mpfr_cmpabs (t, w) <= 0);
1726 if (mpfr_cmpabs (t, w) == 0) /* even rule: significand of q should now
1727 be odd */
1728 MPFR_ASSERTN(mpfr_min_prec (q) == mpfr_get_prec (q));
1729 }
1730
1731 mpfr_clear (q);
1732 mpfr_clear (u);
1733 mpfr_clear (v);
1734 mpfr_clear (t);
1735 mpfr_clear (w);
1736 }
1737
1738 int
main(int argc,char * argv[])1739 main (int argc, char *argv[])
1740 {
1741 tests_start_mpfr ();
1742
1743 coverage (1024);
1744 coverage2 ();
1745 bug20180126 ();
1746 bug20171218 ();
1747 testall_rndf (9);
1748 test_20170105 ();
1749 check_inexact ();
1750 check_hard ();
1751 check_special ();
1752 check_lowr ();
1753 check_float (); /* checks single precision */
1754 check_double ();
1755 check_convergence ();
1756 check_64 ();
1757
1758 check4("4.0","4.503599627370496e15", MPFR_RNDZ, 62,
1759 "0.10000000000000000000000000000000000000000000000000000000000000E-49");
1760 check4("1.0","2.10263340267725788209e+187", MPFR_RNDU, 65,
1761 "0.11010011111001101011111001100111110100000001101001111100111000000E-622");
1762 check4("2.44394909079968374564e-150", "2.10263340267725788209e+187",MPFR_RNDU,
1763 65,
1764 "0.11010011111001101011111001100111110100000001101001111100111000000E-1119");
1765
1766 consistency ();
1767 test_20070603 ();
1768 test_20070628 ();
1769 test_20151023 ();
1770 test_20160831 ();
1771 test_20170104 ();
1772 test_20171219 ();
1773 test_generic (MPFR_PREC_MIN, 800, 50);
1774 test_bad ();
1775 test_extreme ();
1776 test_mpfr_divsp2 ();
1777 #if !defined(MPFR_GENERIC_ABI) && GMP_NUMB_BITS == 64
1778 test_mpfr_div2_approx (1000000);
1779 #endif
1780
1781 tests_end_mpfr ();
1782 return 0;
1783 }
1784