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