1/* Test file for mpfr_div.
2
3Copyright 1999, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
4Contributed by the Arenaire and Cacao projects, INRIA.
5
6This file is part of the GNU MPFR Library.
7
8The GNU MPFR Library is free software; you can redistribute it and/or modify
9it under the terms of the GNU Lesser General Public License as published by
10the Free Software Foundation; either version 3 of the License, or (at your
11option) any later version.
12
13The GNU MPFR Library is distributed in the hope that it will be useful, but
14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
16License for more details.
17
18You should have received a copy of the GNU Lesser General Public License
19along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23#include <stdio.h>
24#include <stdlib.h>
25
26#include "mpfr-test.h"
27
28#ifdef CHECK_EXTERNAL
29static int
30test_div (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
31{
32  int res;
33  int ok = rnd_mode == MPFR_RNDN && mpfr_number_p (b) && mpfr_number_p (c);
34  if (ok)
35    {
36      mpfr_print_raw (b);
37      printf (" ");
38      mpfr_print_raw (c);
39    }
40  res = mpfr_div (a, b, c, rnd_mode);
41  if (ok)
42    {
43      printf (" ");
44      mpfr_print_raw (a);
45      printf ("\n");
46    }
47  return res;
48}
49#else
50#define test_div mpfr_div
51#endif
52
53#define check53(n, d, rnd, res) check4(n, d, rnd, 53, res)
54
55/* return 0 iff a and b are of the same sign */
56static int
57inex_cmp (int a, int b)
58{
59  if (a > 0)
60    return (b > 0) ? 0 : 1;
61  else if (a == 0)
62    return (b == 0) ? 0 : 1;
63  else
64    return (b < 0) ? 0 : 1;
65}
66
67static void
68check4 (const char *Ns, const char *Ds, mpfr_rnd_t rnd_mode, int p,
69        const char *Qs)
70{
71  mpfr_t q, n, d;
72
73  mpfr_inits2 (p, q, n, d, (mpfr_ptr) 0);
74  mpfr_set_str1 (n, Ns);
75  mpfr_set_str1 (d, Ds);
76  test_div(q, n, d, rnd_mode);
77  if (mpfr_cmp_str (q, Qs, ((p==53) ? 10 : 2), MPFR_RNDN) )
78    {
79      printf ("mpfr_div failed for n=%s, d=%s, p=%d, rnd_mode=%s\n",
80              Ns, Ds, p, mpfr_print_rnd_mode (rnd_mode));
81      printf ("got      ");mpfr_print_binary(q);
82      mpfr_set_str (q, Qs, ((p==53) ? 10 : 2), MPFR_RNDN);
83      printf("\nexpected "); mpfr_print_binary(q);
84      putchar('\n');
85      exit (1);
86    }
87  mpfr_clears (q, n, d, (mpfr_ptr) 0);
88}
89
90static void
91check24 (const char *Ns, const char *Ds, mpfr_rnd_t rnd_mode, const char *Qs)
92{
93  mpfr_t q, n, d;
94
95  mpfr_inits2 (24, q, n, d, (mpfr_ptr) 0);
96
97  mpfr_set_str1 (n, Ns);
98  mpfr_set_str1 (d, Ds);
99  test_div(q, n, d, rnd_mode);
100  if (mpfr_cmp_str1 (q, Qs) )
101    {
102      printf ("mpfr_div failed for n=%s, d=%s, prec=24, rnd_mode=%s\n",
103             Ns, Ds, mpfr_print_rnd_mode(rnd_mode));
104      printf ("expected quotient is %s, got ", Qs);
105      mpfr_out_str(stdout,10,0,q, MPFR_RNDN); putchar('\n');
106      exit (1);
107    }
108  mpfr_clears (q, n, d, (mpfr_ptr) 0);
109}
110
111/* the following examples come from the paper "Number-theoretic Test
112   Generation for Directed Rounding" from Michael Parks, Table 2 */
113static void
114check_float(void)
115{
116  check24("70368760954880.0", "8388609.0", MPFR_RNDN, "8.388609e6");
117  check24("140737479966720.0", "16777213.0", MPFR_RNDN, "8.388609e6");
118  check24("70368777732096.0", "8388611.0", MPFR_RNDN, "8.388609e6");
119  check24("105553133043712.0", "12582911.0", MPFR_RNDN, "8.38861e6");
120  /* the exponent for the following example was forgotten in
121     the Arith'14 version of Parks' paper */
122  check24 ("12582913.0", "12582910.0", MPFR_RNDN, "1.000000238");
123  check24 ("105553124655104.0", "12582910.0", MPFR_RNDN, "8388610.0");
124  check24("140737479966720.0", "8388609.0", MPFR_RNDN, "1.6777213e7");
125  check24("70368777732096.0", "8388609.0", MPFR_RNDN, "8.388611e6");
126  check24("105553133043712.0", "8388610.0", MPFR_RNDN, "1.2582911e7");
127  check24("105553124655104.0", "8388610.0", MPFR_RNDN, "1.258291e7");
128
129  check24("70368760954880.0", "8388609.0", MPFR_RNDZ, "8.388608e6");
130  check24("140737479966720.0", "16777213.0", MPFR_RNDZ, "8.388609e6");
131  check24("70368777732096.0", "8388611.0", MPFR_RNDZ, "8.388608e6");
132  check24("105553133043712.0", "12582911.0", MPFR_RNDZ, "8.38861e6");
133  check24("12582913.0", "12582910.0", MPFR_RNDZ, "1.000000238");
134  check24 ("105553124655104.0", "12582910.0", MPFR_RNDZ, "8388610.0");
135  check24("140737479966720.0", "8388609.0", MPFR_RNDZ, "1.6777213e7");
136  check24("70368777732096.0", "8388609.0", MPFR_RNDZ, "8.38861e6");
137  check24("105553133043712.0", "8388610.0", MPFR_RNDZ, "1.2582911e7");
138  check24("105553124655104.0", "8388610.0", MPFR_RNDZ, "1.258291e7");
139
140  check24("70368760954880.0", "8388609.0", MPFR_RNDU, "8.388609e6");
141  check24("140737479966720.0", "16777213.0", MPFR_RNDU, "8.38861e6");
142  check24("70368777732096.0", "8388611.0", MPFR_RNDU, "8.388609e6");
143  check24("105553133043712.0", "12582911.0", MPFR_RNDU, "8.388611e6");
144  check24("12582913.0", "12582910.0", MPFR_RNDU, "1.000000357");
145  check24 ("105553124655104.0", "12582910.0", MPFR_RNDU, "8388611.0");
146  check24("140737479966720.0", "8388609.0", MPFR_RNDU, "1.6777214e7");
147  check24("70368777732096.0", "8388609.0", MPFR_RNDU, "8.388611e6");
148  check24("105553133043712.0", "8388610.0", MPFR_RNDU, "1.2582912e7");
149  check24("105553124655104.0", "8388610.0", MPFR_RNDU, "1.2582911e7");
150
151  check24("70368760954880.0", "8388609.0", MPFR_RNDD, "8.388608e6");
152  check24("140737479966720.0", "16777213.0", MPFR_RNDD, "8.388609e6");
153  check24("70368777732096.0", "8388611.0", MPFR_RNDD, "8.388608e6");
154  check24("105553133043712.0", "12582911.0", MPFR_RNDD, "8.38861e6");
155  check24("12582913.0", "12582910.0", MPFR_RNDD, "1.000000238");
156  check24 ("105553124655104.0", "12582910.0", MPFR_RNDD, "8388610.0");
157  check24("140737479966720.0", "8388609.0", MPFR_RNDD, "1.6777213e7");
158  check24("70368777732096.0", "8388609.0", MPFR_RNDD, "8.38861e6");
159  check24("105553133043712.0", "8388610.0", MPFR_RNDD, "1.2582911e7");
160  check24("105553124655104.0", "8388610.0", MPFR_RNDD, "1.258291e7");
161
162  check24("70368760954880.0", "8388609.0", MPFR_RNDA, "8.388609e6");
163}
164
165static void
166check_double(void)
167{
168  check53("0.0", "1.0", MPFR_RNDZ, "0.0");
169  check53("-7.4988969224688591e63", "4.8816866450288732e306", MPFR_RNDD,
170          "-1.5361282826510687291e-243");
171  check53("-1.33225773037748601769e+199", "3.63449540676937123913e+79",
172          MPFR_RNDZ, "-3.6655920045905428978e119");
173  check53("9.89438396044940256501e-134", "5.93472984109987421717e-67",MPFR_RNDU,
174          "1.6672003992376663654e-67");
175  check53("9.89438396044940256501e-134", "5.93472984109987421717e-67",MPFR_RNDA,
176          "1.6672003992376663654e-67");
177  check53("9.89438396044940256501e-134", "-5.93472984109987421717e-67",
178          MPFR_RNDU, "-1.6672003992376663654e-67");
179  check53("-4.53063926135729747564e-308", "7.02293374921793516813e-84",
180          MPFR_RNDD, "-6.4512060388748850857e-225");
181  check53("6.25089225176473806123e-01","-2.35527154824420243364e-230",
182          MPFR_RNDD, "-2.6540006635008291192e229");
183  check53("6.25089225176473806123e-01","-2.35527154824420243364e-230",
184          MPFR_RNDA, "-2.6540006635008291192e229");
185  check53("6.52308934689126e15", "-1.62063546601505417497e273", MPFR_RNDN,
186          "-4.0250194961676020848e-258");
187  check53("1.04636807108079349236e-189", "3.72295730823253012954e-292",
188          MPFR_RNDZ, "2.810583051186143125e102");
189  /* problems found by Kevin under HP-PA */
190  check53 ("2.861044553323177e-136", "-1.1120354257068143e+45", MPFR_RNDZ,
191           "-2.5727998292003016e-181");
192  check53 ("-4.0559157245809205e-127", "-1.1237723844524865e+77", MPFR_RNDN,
193           "3.6091968273068081e-204");
194  check53 ("-1.8177943561493235e-93", "-8.51233984260364e-104", MPFR_RNDU,
195           "2.1354814184595821e+10");
196}
197
198static void
199check_64(void)
200{
201  mpfr_t x,y,z;
202
203  mpfr_inits2 (64, x, y, z, (mpfr_ptr) 0);
204
205  mpfr_set_str_binary(x, "1.00100100110110101001010010101111000001011100100101010000000000E54");
206  mpfr_set_str_binary(y, "1.00000000000000000000000000000000000000000000000000000000000000E584");
207  test_div(z, x, y, MPFR_RNDU);
208  if (mpfr_cmp_str (z, "0.1001001001101101010010100101011110000010111001001010100000000000E-529", 2, MPFR_RNDN))
209    {
210      printf("Error for tdiv for MPFR_RNDU and p=64\nx=");
211      mpfr_print_binary(x);
212      printf("\ny=");
213      mpfr_print_binary(y);
214      printf("\ngot      ");
215      mpfr_print_binary(z);
216      printf("\nexpected 0.1001001001101101010010100101011110000010111001001010100000000000E-529\n");
217      exit(1);
218    }
219
220  mpfr_clears (x, y, z, (mpfr_ptr) 0);
221}
222
223static void
224check_convergence (void)
225{
226  mpfr_t x, y; int i, j;
227
228  mpfr_init2(x, 130);
229  mpfr_set_str_binary(x, "0.1011111101011010101000001010011111101000011100011101010011111011000011001010000000111100100111110011001010110100100001001000111001E6944");
230  mpfr_init2(y, 130);
231  mpfr_set_ui(y, 5, MPFR_RNDN);
232  test_div(x, x, y, MPFR_RNDD); /* exact division */
233
234  mpfr_set_prec(x, 64);
235  mpfr_set_prec(y, 64);
236  mpfr_set_str_binary(x, "0.10010010011011010100101001010111100000101110010010101E55");
237  mpfr_set_str_binary(y, "0.1E585");
238  test_div(x, x, y, MPFR_RNDN);
239  mpfr_set_str_binary(y, "0.10010010011011010100101001010111100000101110010010101E-529");
240  if (mpfr_cmp (x, y))
241    {
242      printf ("Error in mpfr_div for prec=64, rnd=MPFR_RNDN\n");
243      printf ("got        "); mpfr_print_binary(x); puts ("");
244      printf ("instead of "); mpfr_print_binary(y); puts ("");
245      exit(1);
246    }
247
248  for (i=32; i<=64; i+=32)
249    {
250      mpfr_set_prec(x, i);
251      mpfr_set_prec(y, i);
252      mpfr_set_ui(x, 1, MPFR_RNDN);
253      RND_LOOP(j)
254        {
255          mpfr_set_ui (y, 1, MPFR_RNDN);
256          test_div (y, x, y, (mpfr_rnd_t) j);
257          if (mpfr_cmp_ui (y, 1))
258            {
259              printf ("mpfr_div failed for x=1.0, y=1.0, prec=%d rnd=%s\n",
260                      i, mpfr_print_rnd_mode ((mpfr_rnd_t) j));
261              printf ("got "); mpfr_print_binary(y); puts ("");
262              exit (1);
263            }
264        }
265    }
266
267  mpfr_clear (x);
268  mpfr_clear (y);
269}
270
271#define KMAX 10000
272
273/* given y = o(x/u), x, u, find the inexact flag by
274   multiplying y by u */
275static int
276get_inexact (mpfr_t y, mpfr_t x, mpfr_t u)
277{
278  mpfr_t xx;
279  int inex;
280  mpfr_init2 (xx, mpfr_get_prec (y) + mpfr_get_prec (u));
281  mpfr_mul (xx, y, u, MPFR_RNDN); /* exact */
282  inex = mpfr_cmp (xx, x);
283  mpfr_clear (xx);
284  return inex;
285}
286
287static void
288check_hard (void)
289{
290  mpfr_t u, v, q, q2;
291  mpfr_prec_t precu, precv, precq;
292  int rnd;
293  int inex, inex2, i, j;
294
295  mpfr_init (q);
296  mpfr_init (q2);
297  mpfr_init (u);
298  mpfr_init (v);
299
300  for (precq = MPFR_PREC_MIN; precq <= 64; precq ++)
301    {
302      mpfr_set_prec (q, precq);
303      mpfr_set_prec (q2, precq + 1);
304      for (j = 0; j < 2; j++)
305        {
306          if (j == 0)
307            {
308              do
309                {
310                  mpfr_urandomb (q2, RANDS);
311                }
312              while (mpfr_cmp_ui (q2, 0) == 0);
313            }
314          else /* use q2=1 */
315            mpfr_set_ui (q2, 1, MPFR_RNDN);
316      for (precv = precq; precv <= 10 * precq; precv += precq)
317        {
318          mpfr_set_prec (v, precv);
319          do
320            {
321              mpfr_urandomb (v, RANDS);
322            }
323          while (mpfr_cmp_ui (v, 0) == 0);
324          for (precu = precq; precu <= 10 * precq; precu += precq)
325            {
326              mpfr_set_prec (u, precu);
327              mpfr_mul (u, v, q2, MPFR_RNDN);
328              mpfr_nextbelow (u);
329              for (i = 0; i <= 2; i++)
330                {
331                  RND_LOOP(rnd)
332                    {
333                      inex = test_div (q, u, v, (mpfr_rnd_t) rnd);
334                      inex2 = get_inexact (q, u, v);
335                      if (inex_cmp (inex, inex2))
336                        {
337                          printf ("Wrong inexact flag for rnd=%s: expected %d, got %d\n",
338                                  mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), inex2, inex);
339                          printf ("u=  "); mpfr_dump (u);
340                          printf ("v=  "); mpfr_dump (v);
341                          printf ("q=  "); mpfr_dump (q);
342                          mpfr_set_prec (q2, precq + precv);
343                          mpfr_mul (q2, q, v, MPFR_RNDN);
344                          printf ("q*v="); mpfr_dump (q2);
345                          exit (1);
346                        }
347                    }
348                  mpfr_nextabove (u);
349                }
350            }
351        }
352        }
353    }
354
355  mpfr_clear (q);
356  mpfr_clear (q2);
357  mpfr_clear (u);
358  mpfr_clear (v);
359}
360
361static void
362check_lowr (void)
363{
364  mpfr_t x, y, z, z2, z3, tmp;
365  int k, c, c2;
366
367
368  mpfr_init2 (x, 1000);
369  mpfr_init2 (y, 100);
370  mpfr_init2 (tmp, 850);
371  mpfr_init2 (z, 10);
372  mpfr_init2 (z2, 10);
373  mpfr_init2 (z3, 50);
374
375  for (k = 1; k < KMAX; k++)
376    {
377      do
378        {
379          mpfr_urandomb (z, RANDS);
380        }
381      while (mpfr_cmp_ui (z, 0) == 0);
382      do
383        {
384          mpfr_urandomb (tmp, RANDS);
385        }
386      while (mpfr_cmp_ui (tmp, 0) == 0);
387      mpfr_mul (x, z, tmp, MPFR_RNDN); /* exact */
388      c = test_div (z2, x, tmp, MPFR_RNDN);
389
390      if (c || mpfr_cmp (z2, z))
391        {
392          printf ("Error in mpfr_div rnd=MPFR_RNDN\n");
393          printf ("got        "); mpfr_print_binary(z2); puts ("");
394          printf ("instead of "); mpfr_print_binary(z); puts ("");
395          printf ("inex flag = %d, expected 0\n", c);
396          exit (1);
397        }
398    }
399
400  /* x has still precision 1000, z precision 10, and tmp prec 850 */
401  mpfr_set_prec (z2, 9);
402  for (k = 1; k < KMAX; k++)
403    {
404      mpfr_urandomb (z, RANDS);
405      do
406        {
407          mpfr_urandomb (tmp, RANDS);
408        }
409      while (mpfr_cmp_ui (tmp, 0) == 0);
410      mpfr_mul (x, z, tmp, MPFR_RNDN); /* exact */
411      c = test_div (z2, x, tmp, MPFR_RNDN);
412      /* since z2 has one less bit that z, either the division is exact
413         if z is representable on 9 bits, or we have an even round case */
414
415      c2 = get_inexact (z2, x, tmp);
416      if ((mpfr_cmp (z2, z) == 0 && c) || inex_cmp (c, c2))
417        {
418          printf ("Error in mpfr_div rnd=MPFR_RNDN\n");
419          printf ("got        "); mpfr_print_binary(z2); puts ("");
420          printf ("instead of "); mpfr_print_binary(z); puts ("");
421          printf ("inex flag = %d, expected %d\n", c, c2);
422          exit (1);
423        }
424      else if (c == 2)
425        {
426          mpfr_nexttoinf (z);
427          if (mpfr_cmp(z2, z))
428            {
429              printf ("Error in mpfr_div [even rnd?] rnd=MPFR_RNDN\n");
430              printf ("Dividing ");
431              printf ("got        "); mpfr_print_binary(z2); puts ("");
432              printf ("instead of "); mpfr_print_binary(z); puts ("");
433              printf ("inex flag = %d\n", 1);
434              exit (1);
435            }
436        }
437      else if (c == -2)
438        {
439          mpfr_nexttozero (z);
440          if (mpfr_cmp(z2, z))
441            {
442              printf ("Error in mpfr_div [even rnd?] rnd=MPFR_RNDN\n");
443              printf ("Dividing ");
444              printf ("got        "); mpfr_print_binary(z2); puts ("");
445              printf ("instead of "); mpfr_print_binary(z); puts ("");
446              printf ("inex flag = %d\n", 1);
447              exit (1);
448            }
449        }
450    }
451
452  mpfr_set_prec(x, 1000);
453  mpfr_set_prec(y, 100);
454  mpfr_set_prec(tmp, 850);
455  mpfr_set_prec(z, 10);
456  mpfr_set_prec(z2, 10);
457
458  /* almost exact divisions */
459  for (k = 1; k < KMAX; k++)
460    {
461      do
462        {
463          mpfr_urandomb (z, RANDS);
464        }
465      while (mpfr_cmp_ui (z, 0) == 0);
466      do
467        {
468          mpfr_urandomb (tmp, RANDS);
469        }
470      while (mpfr_cmp_ui (tmp, 0) == 0);
471      mpfr_mul(x, z, tmp, MPFR_RNDN);
472      mpfr_set(y, tmp, MPFR_RNDD);
473      mpfr_nexttoinf (x);
474
475      c = test_div(z2, x, y, MPFR_RNDD);
476      test_div(z3, x, y, MPFR_RNDD);
477      mpfr_set(z, z3, MPFR_RNDD);
478
479      if (c != -1 || mpfr_cmp(z2, z))
480        {
481          printf ("Error in mpfr_div rnd=MPFR_RNDD\n");
482          printf ("got        "); mpfr_print_binary(z2); puts ("");
483          printf ("instead of "); mpfr_print_binary(z); puts ("");
484          printf ("inex flag = %d\n", c);
485          exit (1);
486        }
487
488      mpfr_set (y, tmp, MPFR_RNDU);
489      test_div (z3, x, y, MPFR_RNDU);
490      mpfr_set (z, z3, MPFR_RNDU);
491      c = test_div (z2, x, y, MPFR_RNDU);
492      if (c != 1 || mpfr_cmp (z2, z))
493        {
494          printf ("Error in mpfr_div rnd=MPFR_RNDU\n");
495          printf ("u="); mpfr_dump (x);
496          printf ("v="); mpfr_dump (y);
497          printf ("got        "); mpfr_print_binary (z2); puts ("");
498          printf ("instead of "); mpfr_print_binary (z); puts ("");
499          printf ("inex flag = %d\n", c);
500          exit (1);
501        }
502    }
503
504  mpfr_clear (x);
505  mpfr_clear (y);
506  mpfr_clear (z);
507  mpfr_clear (z2);
508  mpfr_clear (z3);
509  mpfr_clear (tmp);
510}
511
512#define MAX_PREC 128
513
514static void
515check_inexact (void)
516{
517  mpfr_t x, y, z, u;
518  mpfr_prec_t px, py, pu;
519  int inexact, cmp;
520  mpfr_rnd_t rnd;
521
522  mpfr_init (x);
523  mpfr_init (y);
524  mpfr_init (z);
525  mpfr_init (u);
526
527  mpfr_set_prec (x, 28);
528  mpfr_set_prec (y, 28);
529  mpfr_set_prec (z, 1023);
530  mpfr_set_str_binary (x, "0.1000001001101101111100010011E0");
531  mpfr_set_str (z, "48284762641021308813686974720835219181653367326353400027913400579340343320519877153813133510034402932651132854764198688352364361009429039801248971901380781746767119334993621199563870113045276395603170432175354501451429471578325545278975153148347684600400321033502982713296919861760382863826626093689036010394", 10, MPFR_RNDN);
532  mpfr_div (x, x, z, MPFR_RNDN);
533  mpfr_set_str_binary (y, "0.1111001011001101001001111100E-1023");
534  if (mpfr_cmp (x, y))
535    {
536      printf ("Error in mpfr_div for prec=28, RNDN\n");
537      printf ("Expected "); mpfr_dump (y);
538      printf ("Got      "); mpfr_dump (x);
539      exit (1);
540    }
541
542  mpfr_set_prec (x, 53);
543  mpfr_set_str_binary (x, "0.11101100110010100011011000000100001111011111110010101E0");
544  mpfr_set_prec (u, 127);
545  mpfr_set_str_binary (u, "0.1000001100110110110101110110101101111000110000001111111110000000011111001010110100110010111111111101000001011011101011101101000E-2");
546  mpfr_set_prec (y, 95);
547  inexact = test_div (y, x, u, MPFR_RNDN);
548  if (inexact != (cmp = get_inexact (y, x, u)))
549    {
550      printf ("Wrong inexact flag (0): expected %d, got %d\n", cmp, inexact);
551      printf ("x="); mpfr_out_str (stdout, 10, 99, x, MPFR_RNDN); printf ("\n");
552      printf ("u="); mpfr_out_str (stdout, 10, 99, u, MPFR_RNDN); printf ("\n");
553      printf ("y="); mpfr_out_str (stdout, 10, 99, y, MPFR_RNDN); printf ("\n");
554      exit (1);
555    }
556
557  mpfr_set_prec (x, 33);
558  mpfr_set_str_binary (x, "0.101111100011011101010011101100001E0");
559  mpfr_set_prec (u, 2);
560  mpfr_set_str_binary (u, "0.1E0");
561  mpfr_set_prec (y, 28);
562  if ((inexact = test_div (y, x, u, MPFR_RNDN) >= 0))
563    {
564      printf ("Wrong inexact flag (1): expected -1, got %d\n",
565              inexact);
566      exit (1);
567    }
568
569  mpfr_set_prec (x, 129);
570  mpfr_set_str_binary (x, "0.111110101111001100000101011100101100110011011101010001000110110101100101000010000001110110100001101010001010100010001111001101010E-2");
571  mpfr_set_prec (u, 15);
572  mpfr_set_str_binary (u, "0.101101000001100E-1");
573  mpfr_set_prec (y, 92);
574  if ((inexact = test_div (y, x, u, MPFR_RNDN)) <= 0)
575    {
576      printf ("Wrong inexact flag for rnd=MPFR_RNDN(1): expected 1, got %d\n",
577              inexact);
578      mpfr_dump (x);
579      mpfr_dump (u);
580      mpfr_dump (y);
581      exit (1);
582    }
583
584  for (px=2; px<MAX_PREC; px++)
585    {
586      mpfr_set_prec (x, px);
587      mpfr_urandomb (x, RANDS);
588      for (pu=2; pu<=MAX_PREC; pu++)
589        {
590          mpfr_set_prec (u, pu);
591          do { mpfr_urandomb (u, RANDS); } while (mpfr_cmp_ui (u, 0) == 0);
592            {
593              py = MPFR_PREC_MIN + (randlimb () % (MAX_PREC - MPFR_PREC_MIN));
594              mpfr_set_prec (y, py);
595              mpfr_set_prec (z, py + pu);
596                {
597                  rnd = RND_RAND ();
598                  inexact = test_div (y, x, u, rnd);
599                  if (mpfr_mul (z, y, u, rnd))
600                    {
601                      printf ("z <- y * u should be exact\n");
602                      exit (1);
603                    }
604                  cmp = mpfr_cmp (z, x);
605                  if (((inexact == 0) && (cmp != 0)) ||
606                      ((inexact > 0) && (cmp <= 0)) ||
607                      ((inexact < 0) && (cmp >= 0)))
608                    {
609                      printf ("Wrong inexact flag for rnd=%s\n",
610                              mpfr_print_rnd_mode(rnd));
611                      printf ("expected %d, got %d\n", cmp, inexact);
612                      printf ("x="); mpfr_print_binary (x); puts ("");
613                      printf ("u="); mpfr_print_binary (u); puts ("");
614                      printf ("y="); mpfr_print_binary (y); puts ("");
615                      printf ("y*u="); mpfr_print_binary (z); puts ("");
616                      exit (1);
617                    }
618                }
619            }
620        }
621    }
622
623  mpfr_clear (x);
624  mpfr_clear (y);
625  mpfr_clear (z);
626  mpfr_clear (u);
627}
628
629static void
630check_nan (void)
631{
632  mpfr_t  a, d, q;
633  mpfr_exp_t emax, emin;
634  int i;
635
636  mpfr_init2 (a, 100L);
637  mpfr_init2 (d, 100L);
638  mpfr_init2 (q, 100L);
639
640  /* 1/nan == nan */
641  mpfr_set_ui (a, 1L, MPFR_RNDN);
642  MPFR_SET_NAN (d);
643  MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
644  MPFR_ASSERTN (mpfr_nan_p (q));
645
646  /* nan/1 == nan */
647  MPFR_SET_NAN (a);
648  mpfr_set_ui (d, 1L, MPFR_RNDN);
649  MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
650  MPFR_ASSERTN (mpfr_nan_p (q));
651
652  /* +inf/1 == +inf */
653  MPFR_SET_INF (a);
654  MPFR_SET_POS (a);
655  mpfr_set_ui (d, 1L, MPFR_RNDN);
656  MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
657  MPFR_ASSERTN (mpfr_inf_p (q));
658  MPFR_ASSERTN (mpfr_sgn (q) > 0);
659
660  /* 1/+inf == 0 */
661  mpfr_set_ui (a, 1L, MPFR_RNDN);
662  MPFR_SET_INF (d);
663  MPFR_SET_POS (d);
664  MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
665  MPFR_ASSERTN (mpfr_number_p (q));
666  MPFR_ASSERTN (mpfr_sgn (q) == 0);
667
668  /* 0/0 == nan */
669  mpfr_set_ui (a, 0L, MPFR_RNDN);
670  mpfr_set_ui (d, 0L, MPFR_RNDN);
671  MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
672  MPFR_ASSERTN (mpfr_nan_p (q));
673
674  /* +inf/+inf == nan */
675  MPFR_SET_INF (a);
676  MPFR_SET_POS (a);
677  MPFR_SET_INF (d);
678  MPFR_SET_POS (d);
679  MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
680  MPFR_ASSERTN (mpfr_nan_p (q));
681
682  /* 1/+0 = +Inf */
683  mpfr_set_ui (a, 1, MPFR_RNDZ);
684  mpfr_set_ui (d, 0, MPFR_RNDZ);
685  MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
686  MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) > 0);
687
688  /* 1/-0 = -Inf */
689  mpfr_set_ui (a, 1, MPFR_RNDZ);
690  mpfr_set_ui (d, 0, MPFR_RNDZ);
691  mpfr_neg (d, d, MPFR_RNDZ);
692  MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
693  MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) < 0);
694
695  /* -1/+0 = -Inf */
696  mpfr_set_si (a, -1, MPFR_RNDZ);
697  mpfr_set_ui (d, 0, MPFR_RNDZ);
698  MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
699  MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) < 0);
700
701  /* -1/-0 = +Inf */
702  mpfr_set_si (a, -1, MPFR_RNDZ);
703  mpfr_set_ui (d, 0, MPFR_RNDZ);
704  mpfr_neg (d, d, MPFR_RNDZ);
705  MPFR_ASSERTN (test_div (q, a, d, MPFR_RNDZ) == 0); /* exact */
706  MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) > 0);
707
708  /* check overflow */
709  emax = mpfr_get_emax ();
710  set_emax (1);
711  mpfr_set_ui (a, 1, MPFR_RNDZ);
712  mpfr_set_ui (d, 1, MPFR_RNDZ);
713  mpfr_div_2exp (d, d, 1, MPFR_RNDZ);
714  test_div (q, a, d, MPFR_RNDU); /* 1 / 0.5 = 2 -> overflow */
715  MPFR_ASSERTN (mpfr_inf_p (q) && mpfr_sgn (q) > 0);
716  set_emax (emax);
717
718  /* check underflow */
719  emin = mpfr_get_emin ();
720  set_emin (-1);
721  mpfr_set_ui (a, 1, MPFR_RNDZ);
722  mpfr_div_2exp (a, a, 2, MPFR_RNDZ);
723  mpfr_set_prec (d, mpfr_get_prec (q) + 8);
724  for (i = -1; i <= 1; i++)
725    {
726      int sign;
727
728      /* Test 2^(-2) / (+/- (2 + eps)), with eps < 0, eps = 0, eps > 0.
729         -> underflow.
730         With div.c r5513, this test fails for eps > 0 in MPFR_RNDN. */
731      mpfr_set_ui (d, 2, MPFR_RNDZ);
732      if (i < 0)
733        mpfr_nextbelow (d);
734      if (i > 0)
735        mpfr_nextabove (d);
736      for (sign = 0; sign <= 1; sign++)
737        {
738          mpfr_clear_flags ();
739          test_div (q, a, d, MPFR_RNDZ); /* result = 0 */
740          MPFR_ASSERTN (mpfr_underflow_p ());
741          MPFR_ASSERTN (sign ? MPFR_IS_NEG (q) : MPFR_IS_POS (q));
742          MPFR_ASSERTN (MPFR_IS_ZERO (q));
743          mpfr_clear_flags ();
744          test_div (q, a, d, MPFR_RNDN); /* result = 0 iff eps >= 0 */
745          MPFR_ASSERTN (mpfr_underflow_p ());
746          MPFR_ASSERTN (sign ? MPFR_IS_NEG (q) : MPFR_IS_POS (q));
747          if (i < 0)
748            mpfr_nexttozero (q);
749          MPFR_ASSERTN (MPFR_IS_ZERO (q));
750          mpfr_neg (d, d, MPFR_RNDN);
751        }
752    }
753  set_emin (emin);
754
755  mpfr_clear (a);
756  mpfr_clear (d);
757  mpfr_clear (q);
758}
759
760static void
761consistency (void)
762{
763  mpfr_t x, y, z1, z2;
764  int i;
765
766  mpfr_inits (x, y, z1, z2, (mpfr_ptr) 0);
767
768  for (i = 0; i < 10000; i++)
769    {
770      mpfr_rnd_t rnd;
771      mpfr_prec_t px, py, pz, p;
772      int inex1, inex2;
773
774      rnd = RND_RAND ();
775      px = (randlimb () % 256) + 2;
776      py = (randlimb () % 128) + 2;
777      pz = (randlimb () % 256) + 2;
778      mpfr_set_prec (x, px);
779      mpfr_set_prec (y, py);
780      mpfr_set_prec (z1, pz);
781      mpfr_set_prec (z2, pz);
782      mpfr_urandomb (x, RANDS);
783      do
784        mpfr_urandomb (y, RANDS);
785      while (mpfr_zero_p (y));
786      inex1 = mpfr_div (z1, x, y, rnd);
787      MPFR_ASSERTN (!MPFR_IS_NAN (z1));
788      p = MAX (MAX (px, py), pz);
789      if (mpfr_prec_round (x, p, MPFR_RNDN) != 0 ||
790          mpfr_prec_round (y, p, MPFR_RNDN) != 0)
791        {
792          printf ("mpfr_prec_round error for i = %d\n", i);
793          exit (1);
794        }
795      inex2 = mpfr_div (z2, x, y, rnd);
796      MPFR_ASSERTN (!MPFR_IS_NAN (z2));
797      if (inex1 != inex2 || mpfr_cmp (z1, z2) != 0)
798        {
799          printf ("Consistency error for i = %d\n", i);
800          exit (1);
801        }
802    }
803
804  mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0);
805}
806
807/* Reported by Carl Witty on 2007-06-03 */
808static void
809test_20070603 (void)
810{
811  mpfr_t n, d, q, c;
812
813  mpfr_init2 (n, 128);
814  mpfr_init2 (d, 128);
815  mpfr_init2 (q, 31);
816  mpfr_init2 (c, 31);
817
818  mpfr_set_str (n, "10384593717069655257060992206846485", 10, MPFR_RNDN);
819  mpfr_set_str (d, "10384593717069655257060992206847132", 10, MPFR_RNDN);
820  mpfr_div (q, n, d, MPFR_RNDU);
821
822  mpfr_set_ui (c, 1, MPFR_RNDN);
823  if (mpfr_cmp (q, c) != 0)
824    {
825      printf ("Error in test_20070603\nGot        ");
826      mpfr_dump (q);
827      printf ("instead of ");
828      mpfr_dump (c);
829      exit (1);
830    }
831
832  /* same for 64-bit machines */
833  mpfr_set_prec (n, 256);
834  mpfr_set_prec (d, 256);
835  mpfr_set_prec (q, 63);
836  mpfr_set_str (n, "822752278660603021077484591278675252491367930877209729029898240", 10, MPFR_RNDN);
837  mpfr_set_str (d, "822752278660603021077484591278675252491367930877212507873738752", 10, MPFR_RNDN);
838  mpfr_div (q, n, d, MPFR_RNDU);
839  if (mpfr_cmp (q, c) != 0)
840    {
841      printf ("Error in test_20070603\nGot        ");
842      mpfr_dump (q);
843      printf ("instead of ");
844      mpfr_dump (c);
845      exit (1);
846    }
847
848  mpfr_clear (n);
849  mpfr_clear (d);
850  mpfr_clear (q);
851  mpfr_clear (c);
852}
853
854/* Bug found while adding tests for mpfr_cot */
855static void
856test_20070628 (void)
857{
858  mpfr_exp_t old_emax;
859  mpfr_t x, y;
860  int inex, err = 0;
861
862  old_emax = mpfr_get_emax ();
863
864  if (mpfr_set_emax (256))
865    {
866      printf ("Can't change exponent range\n");
867      exit (1);
868    }
869
870  mpfr_inits2 (53, x, y, (mpfr_ptr) 0);
871  mpfr_set_si (x, -1, MPFR_RNDN);
872  mpfr_set_si_2exp (y, 1, -256, MPFR_RNDN);
873  mpfr_clear_flags ();
874  inex = mpfr_div (x, x, y, MPFR_RNDD);
875  if (MPFR_SIGN (x) >= 0 || ! mpfr_inf_p (x))
876    {
877      printf ("Error in test_20070628: expected -Inf, got\n");
878      mpfr_dump (x);
879      err++;
880    }
881  if (inex >= 0)
882    {
883      printf ("Error in test_20070628: expected inex < 0, got %d\n", inex);
884      err++;
885    }
886  if (! mpfr_overflow_p ())
887    {
888      printf ("Error in test_20070628: overflow flag is not set\n");
889      err++;
890    }
891  mpfr_clears (x, y, (mpfr_ptr) 0);
892  mpfr_set_emax (old_emax);
893}
894
895#define TEST_FUNCTION test_div
896#define TWO_ARGS
897#define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), randlimb () % 100, RANDS)
898#include "tgeneric.c"
899
900int
901main (int argc, char *argv[])
902{
903  tests_start_mpfr ();
904
905  check_inexact ();
906  check_hard ();
907  check_nan ();
908  check_lowr ();
909  check_float (); /* checks single precision */
910  check_double ();
911  check_convergence ();
912  check_64 ();
913
914  check4("4.0","4.503599627370496e15", MPFR_RNDZ, 62,
915   "0.10000000000000000000000000000000000000000000000000000000000000E-49");
916  check4("1.0","2.10263340267725788209e+187", MPFR_RNDU, 65,
917   "0.11010011111001101011111001100111110100000001101001111100111000000E-622");
918  check4("2.44394909079968374564e-150", "2.10263340267725788209e+187",MPFR_RNDU,
919         65,
920  "0.11010011111001101011111001100111110100000001101001111100111000000E-1119");
921
922  consistency ();
923  test_20070603 ();
924  test_20070628 ();
925  test_generic (2, 800, 50);
926
927  tests_end_mpfr ();
928  return 0;
929}
930