Home | History | Annotate | Line # | Download | only in tests
tgmpop.c revision 1.1
      1  1.1  mrg /* Test file for mpfr_add_[q,z], mpfr_sub_[q,z], mpfr_div_[q,z], mpfr_mul_[q,z]
      2  1.1  mrg    and mpfr_cmp_[q,z]
      3  1.1  mrg 
      4  1.1  mrg Copyright 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
      5  1.1  mrg Contributed by the Arenaire and Cacao projects, INRIA.
      6  1.1  mrg 
      7  1.1  mrg This file is part of the GNU MPFR Library.
      8  1.1  mrg 
      9  1.1  mrg The GNU MPFR Library is free software; you can redistribute it and/or modify
     10  1.1  mrg it under the terms of the GNU Lesser General Public License as published by
     11  1.1  mrg the Free Software Foundation; either version 3 of the License, or (at your
     12  1.1  mrg option) any later version.
     13  1.1  mrg 
     14  1.1  mrg The GNU MPFR Library is distributed in the hope that it will be useful, but
     15  1.1  mrg WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
     16  1.1  mrg or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
     17  1.1  mrg License for more details.
     18  1.1  mrg 
     19  1.1  mrg You should have received a copy of the GNU Lesser General Public License
     20  1.1  mrg along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
     21  1.1  mrg http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
     22  1.1  mrg 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
     23  1.1  mrg 
     24  1.1  mrg #include <stdio.h>
     25  1.1  mrg #include <stdlib.h>
     26  1.1  mrg #include "mpfr-test.h"
     27  1.1  mrg 
     28  1.1  mrg #define CHECK_FOR(str, cond)                                      \
     29  1.1  mrg if ((cond) == 0) {                                                \
     30  1.1  mrg   printf ("Special case error " str ". Inexact value=%d\n", res); \
     31  1.1  mrg   printf ("Get "); mpfr_dump (y);                                 \
     32  1.1  mrg   printf ("X=  "); mpfr_dump (x);                                 \
     33  1.1  mrg   printf ("Z=  "); mpz_dump (mpq_numref(z));                      \
     34  1.1  mrg   printf ("   /"); mpz_dump (mpq_denref(z));                      \
     35  1.1  mrg   exit (1);                                                       \
     36  1.1  mrg }
     37  1.1  mrg 
     38  1.1  mrg static void
     39  1.1  mrg special (void)
     40  1.1  mrg {
     41  1.1  mrg   mpfr_t x, y;
     42  1.1  mrg   mpq_t z;
     43  1.1  mrg   int res = 0;
     44  1.1  mrg 
     45  1.1  mrg   mpfr_init (x);
     46  1.1  mrg   mpfr_init (y);
     47  1.1  mrg   mpq_init (z);
     48  1.1  mrg 
     49  1.1  mrg   /* cancellation in mpfr_add_q */
     50  1.1  mrg   mpfr_set_prec (x, 60);
     51  1.1  mrg   mpfr_set_prec (y, 20);
     52  1.1  mrg   mpz_set_str (mpq_numref (z), "-187207494", 10);
     53  1.1  mrg   mpz_set_str (mpq_denref (z), "5721", 10);
     54  1.1  mrg   mpfr_set_str_binary (x, "11111111101001011011100101100011011110010011100010000100001E-44");
     55  1.1  mrg   mpfr_add_q (y, x, z, MPFR_RNDN);
     56  1.1  mrg   CHECK_FOR ("cancelation in add_q", mpfr_cmp_ui_2exp (y, 256783, -64) == 0);
     57  1.1  mrg 
     58  1.1  mrg   mpfr_set_prec (x, 19);
     59  1.1  mrg   mpfr_set_str_binary (x, "0.1011110101110011100E0");
     60  1.1  mrg   mpz_set_str (mpq_numref (z), "187207494", 10);
     61  1.1  mrg   mpz_set_str (mpq_denref (z), "5721", 10);
     62  1.1  mrg   mpfr_set_prec (y, 29);
     63  1.1  mrg   mpfr_add_q (y, x, z, MPFR_RNDD);
     64  1.1  mrg   mpfr_set_prec (x, 29);
     65  1.1  mrg   mpfr_set_str_binary (x, "11111111101001110011010001001E-14");
     66  1.1  mrg   CHECK_FOR ("cancelation in add_q", mpfr_cmp (x,y) == 0);
     67  1.1  mrg 
     68  1.1  mrg   /* Inf */
     69  1.1  mrg   mpfr_set_inf (x, 1);
     70  1.1  mrg   mpz_set_str (mpq_numref (z), "395877315", 10);
     71  1.1  mrg   mpz_set_str (mpq_denref (z), "3508975966", 10);
     72  1.1  mrg   mpfr_set_prec (y, 118);
     73  1.1  mrg   mpfr_add_q (y, x, z, MPFR_RNDU);
     74  1.1  mrg   CHECK_FOR ("inf", mpfr_inf_p (y) && mpfr_sgn (y) > 0);
     75  1.1  mrg   mpfr_sub_q (y, x, z, MPFR_RNDU);
     76  1.1  mrg   CHECK_FOR ("inf", mpfr_inf_p (y) && mpfr_sgn (y) > 0);
     77  1.1  mrg 
     78  1.1  mrg   /* Nan */
     79  1.1  mrg   MPFR_SET_NAN (x);
     80  1.1  mrg   mpfr_add_q (y, x, z, MPFR_RNDU);
     81  1.1  mrg   CHECK_FOR ("nan", mpfr_nan_p (y));
     82  1.1  mrg   mpfr_sub_q (y, x, z, MPFR_RNDU);
     83  1.1  mrg   CHECK_FOR ("nan", mpfr_nan_p (y));
     84  1.1  mrg 
     85  1.1  mrg   /* Exact value */
     86  1.1  mrg   mpfr_set_prec (x, 60);
     87  1.1  mrg   mpfr_set_prec (y, 60);
     88  1.1  mrg   mpfr_set_str1 (x, "0.5");
     89  1.1  mrg   mpz_set_str (mpq_numref (z), "3", 10);
     90  1.1  mrg   mpz_set_str (mpq_denref (z), "2", 10);
     91  1.1  mrg   res = mpfr_add_q (y, x, z, MPFR_RNDU);
     92  1.1  mrg   CHECK_FOR ("0.5+3/2", mpfr_cmp_ui(y, 2)==0 && res==0);
     93  1.1  mrg   res = mpfr_sub_q (y, x, z, MPFR_RNDU);
     94  1.1  mrg   CHECK_FOR ("0.5-3/2", mpfr_cmp_si(y, -1)==0 && res==0);
     95  1.1  mrg 
     96  1.1  mrg   /* Inf Rationnal */
     97  1.1  mrg   mpq_set_ui (z, 1, 0);
     98  1.1  mrg   mpfr_set_str1 (x, "0.5");
     99  1.1  mrg   res = mpfr_add_q (y, x, z, MPFR_RNDN);
    100  1.1  mrg   CHECK_FOR ("0.5+1/0", mpfr_inf_p (y) && MPFR_SIGN (y) > 0 && res == 0);
    101  1.1  mrg   res = mpfr_sub_q (y, x, z, MPFR_RNDN);
    102  1.1  mrg   CHECK_FOR ("0.5-1/0", mpfr_inf_p (y) && MPFR_SIGN (y) < 0 && res == 0);
    103  1.1  mrg   mpq_set_si (z, -1, 0);
    104  1.1  mrg   res = mpfr_add_q (y, x, z, MPFR_RNDN);
    105  1.1  mrg   CHECK_FOR ("0.5+ -1/0", mpfr_inf_p (y) && MPFR_SIGN (y) < 0 && res == 0);
    106  1.1  mrg   res = mpfr_sub_q (y, x, z, MPFR_RNDN);
    107  1.1  mrg   CHECK_FOR ("0.5- -1/0", mpfr_inf_p (y) && MPFR_SIGN (y) > 0 && res == 0);
    108  1.1  mrg   res = mpfr_div_q (y, x, z, MPFR_RNDN);
    109  1.1  mrg   CHECK_FOR ("0.5 / (-1/0)", mpfr_zero_p (y) && MPFR_SIGN (y) < 0 && res == 0);
    110  1.1  mrg 
    111  1.1  mrg   /* 0 */
    112  1.1  mrg   mpq_set_ui (z, 0, 1);
    113  1.1  mrg   mpfr_set_ui (x, 42, MPFR_RNDN);
    114  1.1  mrg   res = mpfr_add_q (y, x, z, MPFR_RNDN);
    115  1.1  mrg   CHECK_FOR ("42+0/1", mpfr_cmp_ui (y, 42) == 0 && res == 0);
    116  1.1  mrg   res = mpfr_sub_q (y, x, z, MPFR_RNDN);
    117  1.1  mrg   CHECK_FOR ("42-0/1", mpfr_cmp_ui (y, 42) == 0 && res == 0);
    118  1.1  mrg   res = mpfr_mul_q (y, x, z, MPFR_RNDN);
    119  1.1  mrg   CHECK_FOR ("42*0/1", mpfr_zero_p (y) && MPFR_SIGN (y) > 0 && res == 0);
    120  1.1  mrg   res = mpfr_div_q (y, x, z, MPFR_RNDN);
    121  1.1  mrg   CHECK_FOR ("42/(0/1)", mpfr_inf_p (y) && MPFR_SIGN (y) > 0 && res == 0);
    122  1.1  mrg 
    123  1.1  mrg 
    124  1.1  mrg   mpq_clear (z);
    125  1.1  mrg   mpfr_clear (x);
    126  1.1  mrg   mpfr_clear (y);
    127  1.1  mrg }
    128  1.1  mrg 
    129  1.1  mrg static void
    130  1.1  mrg check_for_zero (void)
    131  1.1  mrg {
    132  1.1  mrg   /* Check that 0 is unsigned! */
    133  1.1  mrg   mpq_t q;
    134  1.1  mrg   mpz_t z;
    135  1.1  mrg   mpfr_t x;
    136  1.1  mrg   int r;
    137  1.1  mrg   mpfr_sign_t i;
    138  1.1  mrg 
    139  1.1  mrg   mpfr_init (x);
    140  1.1  mrg   mpz_init (z);
    141  1.1  mrg   mpq_init (q);
    142  1.1  mrg 
    143  1.1  mrg   mpz_set_ui (z, 0);
    144  1.1  mrg   mpq_set_ui (q, 0, 1);
    145  1.1  mrg 
    146  1.1  mrg   MPFR_SET_ZERO (x);
    147  1.1  mrg   RND_LOOP (r)
    148  1.1  mrg     {
    149  1.1  mrg       for (i = MPFR_SIGN_NEG ; i <= MPFR_SIGN_POS ;
    150  1.1  mrg            i+=MPFR_SIGN_POS-MPFR_SIGN_NEG)
    151  1.1  mrg         {
    152  1.1  mrg           MPFR_SET_SIGN(x, i);
    153  1.1  mrg           mpfr_add_z (x, x, z, (mpfr_rnd_t) r);
    154  1.1  mrg           if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i)
    155  1.1  mrg             {
    156  1.1  mrg               printf("GMP Zero errors for add_z & rnd=%s & s=%d\n",
    157  1.1  mrg                      mpfr_print_rnd_mode ((mpfr_rnd_t) r), i);
    158  1.1  mrg               mpfr_dump (x);
    159  1.1  mrg               exit (1);
    160  1.1  mrg             }
    161  1.1  mrg           mpfr_sub_z (x, x, z, (mpfr_rnd_t) r);
    162  1.1  mrg           if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i)
    163  1.1  mrg             {
    164  1.1  mrg               printf("GMP Zero errors for sub_z & rnd=%s & s=%d\n",
    165  1.1  mrg                      mpfr_print_rnd_mode ((mpfr_rnd_t) r), i);
    166  1.1  mrg               mpfr_dump (x);
    167  1.1  mrg               exit (1);
    168  1.1  mrg             }
    169  1.1  mrg           mpfr_mul_z (x, x, z, (mpfr_rnd_t) r);
    170  1.1  mrg           if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i)
    171  1.1  mrg             {
    172  1.1  mrg               printf("GMP Zero errors for mul_z & rnd=%s & s=%d\n",
    173  1.1  mrg                      mpfr_print_rnd_mode ((mpfr_rnd_t) r), i);
    174  1.1  mrg               mpfr_dump (x);
    175  1.1  mrg               exit (1);
    176  1.1  mrg             }
    177  1.1  mrg           mpfr_add_q (x, x, q, (mpfr_rnd_t) r);
    178  1.1  mrg           if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i)
    179  1.1  mrg             {
    180  1.1  mrg               printf("GMP Zero errors for add_q & rnd=%s & s=%d\n",
    181  1.1  mrg                      mpfr_print_rnd_mode ((mpfr_rnd_t) r), i);
    182  1.1  mrg               mpfr_dump (x);
    183  1.1  mrg               exit (1);
    184  1.1  mrg             }
    185  1.1  mrg           mpfr_sub_q (x, x, q, (mpfr_rnd_t) r);
    186  1.1  mrg           if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i)
    187  1.1  mrg             {
    188  1.1  mrg               printf("GMP Zero errors for sub_q & rnd=%s & s=%d\n",
    189  1.1  mrg                      mpfr_print_rnd_mode ((mpfr_rnd_t) r), i);
    190  1.1  mrg               mpfr_dump (x);
    191  1.1  mrg               exit (1);
    192  1.1  mrg              }
    193  1.1  mrg         }
    194  1.1  mrg     }
    195  1.1  mrg 
    196  1.1  mrg   mpq_clear (q);
    197  1.1  mrg   mpz_clear (z);
    198  1.1  mrg   mpfr_clear (x);
    199  1.1  mrg }
    200  1.1  mrg 
    201  1.1  mrg static void
    202  1.1  mrg test_cmp_z (mpfr_prec_t pmin, mpfr_prec_t pmax, int nmax)
    203  1.1  mrg {
    204  1.1  mrg   mpfr_t x, z;
    205  1.1  mrg   mpz_t  y;
    206  1.1  mrg   mpfr_prec_t p;
    207  1.1  mrg   int res1, res2;
    208  1.1  mrg   int n;
    209  1.1  mrg 
    210  1.1  mrg   mpfr_init (x);
    211  1.1  mrg   mpfr_init2 (z, MPFR_PREC_MIN);
    212  1.1  mrg   mpz_init (y);
    213  1.1  mrg   for(p=pmin ; p < pmax ; p++)
    214  1.1  mrg     {
    215  1.1  mrg       mpfr_set_prec (x, p);
    216  1.1  mrg       for ( n = 0; n < nmax ; n++)
    217  1.1  mrg         {
    218  1.1  mrg           mpfr_urandomb (x, RANDS);
    219  1.1  mrg           mpz_urandomb  (y, RANDS, 1024);
    220  1.1  mrg           if (!MPFR_IS_SINGULAR (x))
    221  1.1  mrg             {
    222  1.1  mrg               mpfr_sub_z (z, x, y, MPFR_RNDN);
    223  1.1  mrg               res1 = mpfr_sgn (z);
    224  1.1  mrg               res2 = mpfr_cmp_z (x, y);
    225  1.1  mrg               if (res1 != res2)
    226  1.1  mrg                 {
    227  1.1  mrg                   printf("Error for mpfr_cmp_z: res=%d sub_z gives %d\n",
    228  1.1  mrg                          res2, res1);
    229  1.1  mrg                   exit (1);
    230  1.1  mrg                 }
    231  1.1  mrg             }
    232  1.1  mrg         }
    233  1.1  mrg     }
    234  1.1  mrg   mpz_clear (y);
    235  1.1  mrg   mpfr_clear (x);
    236  1.1  mrg   mpfr_clear (z);
    237  1.1  mrg }
    238  1.1  mrg 
    239  1.1  mrg static void
    240  1.1  mrg test_cmp_q (mpfr_prec_t pmin, mpfr_prec_t pmax, int nmax)
    241  1.1  mrg {
    242  1.1  mrg   mpfr_t x, z;
    243  1.1  mrg   mpq_t  y;
    244  1.1  mrg   mpfr_prec_t p;
    245  1.1  mrg   int res1, res2;
    246  1.1  mrg   int n;
    247  1.1  mrg 
    248  1.1  mrg   mpfr_init (x);
    249  1.1  mrg   mpfr_init2 (z, MPFR_PREC_MIN);
    250  1.1  mrg   mpq_init (y);
    251  1.1  mrg   for(p=pmin ; p < pmax ; p++)
    252  1.1  mrg     {
    253  1.1  mrg       mpfr_set_prec (x, p);
    254  1.1  mrg       for (n = 0 ; n < nmax ; n++)
    255  1.1  mrg         {
    256  1.1  mrg           mpfr_urandomb (x, RANDS);
    257  1.1  mrg           mpq_set_ui (y, randlimb (), randlimb() );
    258  1.1  mrg           if (!MPFR_IS_SINGULAR (x))
    259  1.1  mrg             {
    260  1.1  mrg               mpfr_sub_q (z, x, y, MPFR_RNDN);
    261  1.1  mrg               res1 = mpfr_sgn (z);
    262  1.1  mrg               res2 = mpfr_cmp_q (x, y);
    263  1.1  mrg               if (res1 != res2)
    264  1.1  mrg                 {
    265  1.1  mrg                   printf("Error for mpfr_cmp_q: res=%d sub_z gives %d\n",
    266  1.1  mrg                          res2, res1);
    267  1.1  mrg                   exit (1);
    268  1.1  mrg                 }
    269  1.1  mrg             }
    270  1.1  mrg         }
    271  1.1  mrg     }
    272  1.1  mrg   mpq_clear (y);
    273  1.1  mrg   mpfr_clear (x);
    274  1.1  mrg   mpfr_clear (z);
    275  1.1  mrg }
    276  1.1  mrg 
    277  1.1  mrg static void
    278  1.1  mrg test_cmp_f (mpfr_prec_t pmin, mpfr_prec_t pmax, int nmax)
    279  1.1  mrg {
    280  1.1  mrg   mpfr_t x, z;
    281  1.1  mrg   mpf_t  y;
    282  1.1  mrg   mpfr_prec_t p;
    283  1.1  mrg   int res1, res2;
    284  1.1  mrg   int n;
    285  1.1  mrg 
    286  1.1  mrg   mpfr_init (x);
    287  1.1  mrg   mpfr_init2 (z, pmax+GMP_NUMB_BITS);
    288  1.1  mrg   mpf_init2 (y, MPFR_PREC_MIN);
    289  1.1  mrg   for(p=pmin ; p < pmax ; p+=3)
    290  1.1  mrg     {
    291  1.1  mrg       mpfr_set_prec (x, p);
    292  1.1  mrg       mpf_set_prec (y, p);
    293  1.1  mrg       for ( n = 0; n < nmax ; n++)
    294  1.1  mrg         {
    295  1.1  mrg           mpfr_urandomb (x, RANDS);
    296  1.1  mrg           mpf_urandomb  (y, RANDS, p);
    297  1.1  mrg           if (!MPFR_IS_SINGULAR (x))
    298  1.1  mrg             {
    299  1.1  mrg               mpfr_set_f (z, y, MPFR_RNDN);
    300  1.1  mrg               mpfr_sub   (z, x, z, MPFR_RNDN);
    301  1.1  mrg               res1 = mpfr_sgn (z);
    302  1.1  mrg               res2 = mpfr_cmp_f (x, y);
    303  1.1  mrg               if (res1 != res2)
    304  1.1  mrg                 {
    305  1.1  mrg                   printf("Error for mpfr_cmp_f: res=%d sub gives %d\n",
    306  1.1  mrg                          res2, res1);
    307  1.1  mrg                   exit (1);
    308  1.1  mrg                 }
    309  1.1  mrg             }
    310  1.1  mrg         }
    311  1.1  mrg     }
    312  1.1  mrg   mpf_clear (y);
    313  1.1  mrg   mpfr_clear (x);
    314  1.1  mrg   mpfr_clear (z);
    315  1.1  mrg }
    316  1.1  mrg 
    317  1.1  mrg static void
    318  1.1  mrg test_specialz (int (*mpfr_func)(mpfr_ptr, mpfr_srcptr, mpz_srcptr, mpfr_rnd_t),
    319  1.1  mrg                void (*mpz_func)(mpz_ptr, mpz_srcptr, mpz_srcptr),
    320  1.1  mrg                const char *op)
    321  1.1  mrg {
    322  1.1  mrg   mpfr_t x1, x2;
    323  1.1  mrg   mpz_t  z1, z2;
    324  1.1  mrg   int res;
    325  1.1  mrg 
    326  1.1  mrg   mpfr_inits2 (128, x1, x2, (mpfr_ptr) 0);
    327  1.1  mrg   mpz_init (z1); mpz_init(z2);
    328  1.1  mrg   mpz_fac_ui (z1, 19); /* 19!+1 fits perfectly in a 128 bits mantissa */
    329  1.1  mrg   mpz_add_ui (z1, z1, 1);
    330  1.1  mrg   mpz_fac_ui (z2, 20); /* 20!+1 fits perfectly in a 128 bits mantissa */
    331  1.1  mrg   mpz_add_ui (z2, z2, 1);
    332  1.1  mrg 
    333  1.1  mrg   res = mpfr_set_z(x1, z1, MPFR_RNDN);
    334  1.1  mrg   if (res)
    335  1.1  mrg     {
    336  1.1  mrg       printf("Specialz %s: set_z1 error\n", op);
    337  1.1  mrg       exit(1);
    338  1.1  mrg     }
    339  1.1  mrg   mpfr_set_z (x2, z2, MPFR_RNDN);
    340  1.1  mrg   if (res)
    341  1.1  mrg     {
    342  1.1  mrg       printf("Specialz %s: set_z2 error\n", op);
    343  1.1  mrg       exit(1);
    344  1.1  mrg     }
    345  1.1  mrg 
    346  1.1  mrg   /* (19!+1) * (20!+1) fits in a 128 bits number */
    347  1.1  mrg   res = mpfr_func(x1, x1, z2, MPFR_RNDN);
    348  1.1  mrg   if (res)
    349  1.1  mrg     {
    350  1.1  mrg       printf("Specialz %s: wrong inexact flag.\n", op);
    351  1.1  mrg       exit(1);
    352  1.1  mrg     }
    353  1.1  mrg   mpz_func(z1, z1, z2);
    354  1.1  mrg   res = mpfr_set_z (x2, z1, MPFR_RNDN);
    355  1.1  mrg   if (res)
    356  1.1  mrg     {
    357  1.1  mrg       printf("Specialz %s: set_z2 error\n", op);
    358  1.1  mrg       exit(1);
    359  1.1  mrg     }
    360  1.1  mrg   if (mpfr_cmp(x1, x2))
    361  1.1  mrg     {
    362  1.1  mrg       printf("Specialz %s: results differ.\nx1=", op);
    363  1.1  mrg       mpfr_print_binary(x1);
    364  1.1  mrg       printf("\nx2=");
    365  1.1  mrg       mpfr_print_binary(x2);
    366  1.1  mrg       printf ("\nZ2=");
    367  1.1  mrg       mpz_out_str (stdout, 2, z1);
    368  1.1  mrg       putchar('\n');
    369  1.1  mrg       exit(1);
    370  1.1  mrg     }
    371  1.1  mrg 
    372  1.1  mrg   mpz_set_ui (z1, 1);
    373  1.1  mrg   mpz_set_ui (z2, 0);
    374  1.1  mrg   mpfr_set_ui (x1, 1, MPFR_RNDN);
    375  1.1  mrg   mpz_func (z1, z1, z2);
    376  1.1  mrg   res = mpfr_func(x1, x1, z2, MPFR_RNDN);
    377  1.1  mrg   mpfr_set_z (x2, z1, MPFR_RNDN);
    378  1.1  mrg   if (mpfr_cmp(x1, x2))
    379  1.1  mrg     {
    380  1.1  mrg       printf("Specialz %s: results differ(2).\nx1=", op);
    381  1.1  mrg       mpfr_print_binary(x1);
    382  1.1  mrg       printf("\nx2=");
    383  1.1  mrg       mpfr_print_binary(x2);
    384  1.1  mrg       putchar('\n');
    385  1.1  mrg       exit(1);
    386  1.1  mrg     }
    387  1.1  mrg 
    388  1.1  mrg   mpz_clear (z1); mpz_clear(z2);
    389  1.1  mrg   mpfr_clears (x1, x2, (mpfr_ptr) 0);
    390  1.1  mrg }
    391  1.1  mrg 
    392  1.1  mrg static void
    393  1.1  mrg test_genericz (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int N,
    394  1.1  mrg                int (*func)(mpfr_ptr, mpfr_srcptr, mpz_srcptr, mpfr_rnd_t),
    395  1.1  mrg                const char *op)
    396  1.1  mrg {
    397  1.1  mrg   mpfr_prec_t prec;
    398  1.1  mrg   mpfr_t arg1, dst_big, dst_small, tmp;
    399  1.1  mrg   mpz_t  arg2;
    400  1.1  mrg   mpfr_rnd_t rnd;
    401  1.1  mrg   int inexact, compare, compare2;
    402  1.1  mrg   unsigned int n;
    403  1.1  mrg 
    404  1.1  mrg   mpfr_inits (arg1, dst_big, dst_small, tmp, (mpfr_ptr) 0);
    405  1.1  mrg   mpz_init (arg2);
    406  1.1  mrg 
    407  1.1  mrg   for (prec = p0; prec <= p1; prec++)
    408  1.1  mrg     {
    409  1.1  mrg       mpfr_set_prec (arg1, prec);
    410  1.1  mrg       mpfr_set_prec (tmp, prec);
    411  1.1  mrg       mpfr_set_prec (dst_small, prec);
    412  1.1  mrg 
    413  1.1  mrg       for (n=0; n<N; n++)
    414  1.1  mrg         {
    415  1.1  mrg           mpfr_urandomb (arg1, RANDS);
    416  1.1  mrg           mpz_urandomb (arg2, RANDS, 1024);
    417  1.1  mrg           rnd = RND_RAND ();
    418  1.1  mrg           mpfr_set_prec (dst_big, 2*prec);
    419  1.1  mrg           compare = func(dst_big, arg1, arg2, rnd);
    420  1.1  mrg           if (mpfr_can_round (dst_big, 2*prec, rnd, rnd, prec))
    421  1.1  mrg             {
    422  1.1  mrg               mpfr_set (tmp, dst_big, rnd);
    423  1.1  mrg               inexact = func(dst_small, arg1, arg2, rnd);
    424  1.1  mrg               if (mpfr_cmp (tmp, dst_small))
    425  1.1  mrg                 {
    426  1.1  mrg                   printf ("Results differ for prec=%u rnd_mode=%s and %s_z:\n"
    427  1.1  mrg                           "arg1=",
    428  1.1  mrg                           (unsigned) prec, mpfr_print_rnd_mode (rnd), op);
    429  1.1  mrg                   mpfr_print_binary (arg1);
    430  1.1  mrg                   printf("\narg2=");
    431  1.1  mrg                   mpz_out_str (stdout, 2, arg2);
    432  1.1  mrg                   printf ("\ngot      ");
    433  1.1  mrg                   mpfr_dump (dst_small);
    434  1.1  mrg                   printf ("expected ");
    435  1.1  mrg                   mpfr_dump (tmp);
    436  1.1  mrg                   printf ("approx   ");
    437  1.1  mrg                   mpfr_dump (dst_big);
    438  1.1  mrg                   exit (1);
    439  1.1  mrg                 }
    440  1.1  mrg               compare2 = mpfr_cmp (tmp, dst_big);
    441  1.1  mrg               /* if rounding to nearest, cannot know the sign of t - f(x)
    442  1.1  mrg                  because of composed rounding: y = o(f(x)) and t = o(y) */
    443  1.1  mrg               if (compare * compare2 >= 0)
    444  1.1  mrg                 compare = compare + compare2;
    445  1.1  mrg               else
    446  1.1  mrg                 compare = inexact; /* cannot determine sign(t-f(x)) */
    447  1.1  mrg               if (((inexact == 0) && (compare != 0)) ||
    448  1.1  mrg                   ((inexact > 0) && (compare <= 0)) ||
    449  1.1  mrg                   ((inexact < 0) && (compare >= 0)))
    450  1.1  mrg                 {
    451  1.1  mrg                   printf ("Wrong inexact flag for rnd=%s and %s_z:\n"
    452  1.1  mrg                           "expected %d, got %d\n",
    453  1.1  mrg                           mpfr_print_rnd_mode (rnd), op, compare, inexact);
    454  1.1  mrg                   printf ("\narg1="); mpfr_print_binary (arg1);
    455  1.1  mrg                   printf ("\narg2="); mpz_out_str(stdout, 2, arg2);
    456  1.1  mrg                   printf ("\ndstl="); mpfr_print_binary (dst_big);
    457  1.1  mrg                   printf ("\ndsts="); mpfr_print_binary (dst_small);
    458  1.1  mrg                   printf ("\ntmp ="); mpfr_dump (tmp);
    459  1.1  mrg                   exit (1);
    460  1.1  mrg                 }
    461  1.1  mrg             }
    462  1.1  mrg         }
    463  1.1  mrg     }
    464  1.1  mrg 
    465  1.1  mrg   mpz_clear (arg2);
    466  1.1  mrg   mpfr_clears (arg1, dst_big, dst_small, tmp, (mpfr_ptr) 0);
    467  1.1  mrg }
    468  1.1  mrg 
    469  1.1  mrg static void
    470  1.1  mrg test_genericq (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int N,
    471  1.1  mrg                int (*func)(mpfr_ptr, mpfr_srcptr, mpq_srcptr, mpfr_rnd_t),
    472  1.1  mrg                const char *op)
    473  1.1  mrg {
    474  1.1  mrg   mpfr_prec_t prec;
    475  1.1  mrg   mpfr_t arg1, dst_big, dst_small, tmp;
    476  1.1  mrg   mpq_t  arg2;
    477  1.1  mrg   mpfr_rnd_t rnd;
    478  1.1  mrg   int inexact, compare, compare2;
    479  1.1  mrg   unsigned int n;
    480  1.1  mrg 
    481  1.1  mrg   mpfr_inits (arg1, dst_big, dst_small, tmp, (mpfr_ptr) 0);
    482  1.1  mrg   mpq_init (arg2);
    483  1.1  mrg 
    484  1.1  mrg   for (prec = p0; prec <= p1; prec++)
    485  1.1  mrg     {
    486  1.1  mrg       mpfr_set_prec (arg1, prec);
    487  1.1  mrg       mpfr_set_prec (tmp, prec);
    488  1.1  mrg       mpfr_set_prec (dst_small, prec);
    489  1.1  mrg 
    490  1.1  mrg       for (n=0; n<N; n++)
    491  1.1  mrg         {
    492  1.1  mrg           mpfr_urandomb (arg1, RANDS);
    493  1.1  mrg           mpq_set_ui (arg2, randlimb (), randlimb() );
    494  1.1  mrg           mpq_canonicalize (arg2);
    495  1.1  mrg           rnd = RND_RAND ();
    496  1.1  mrg           mpfr_set_prec (dst_big, prec+10);
    497  1.1  mrg           compare = func(dst_big, arg1, arg2, rnd);
    498  1.1  mrg           if (mpfr_can_round (dst_big, prec+10, rnd, rnd, prec))
    499  1.1  mrg             {
    500  1.1  mrg               mpfr_set (tmp, dst_big, rnd);
    501  1.1  mrg               inexact = func(dst_small, arg1, arg2, rnd);
    502  1.1  mrg               if (mpfr_cmp (tmp, dst_small))
    503  1.1  mrg                 {
    504  1.1  mrg                   printf ("Results differ for prec=%u rnd_mode=%s and %s_q:\n"
    505  1.1  mrg                           "arg1=",
    506  1.1  mrg                           (unsigned) prec, mpfr_print_rnd_mode (rnd), op);
    507  1.1  mrg                   mpfr_print_binary (arg1);
    508  1.1  mrg                   printf("\narg2=");
    509  1.1  mrg                   mpq_out_str(stdout, 2, arg2);
    510  1.1  mrg                   printf ("\ngot      ");
    511  1.1  mrg                   mpfr_print_binary (dst_small);
    512  1.1  mrg                   printf ("\nexpected ");
    513  1.1  mrg                   mpfr_print_binary (tmp);
    514  1.1  mrg                   printf ("\napprox  ");
    515  1.1  mrg                   mpfr_print_binary (dst_big);
    516  1.1  mrg                   putchar('\n');
    517  1.1  mrg                   exit (1);
    518  1.1  mrg                 }
    519  1.1  mrg               compare2 = mpfr_cmp (tmp, dst_big);
    520  1.1  mrg               /* if rounding to nearest, cannot know the sign of t - f(x)
    521  1.1  mrg                  because of composed rounding: y = o(f(x)) and t = o(y) */
    522  1.1  mrg               if (compare * compare2 >= 0)
    523  1.1  mrg                 compare = compare + compare2;
    524  1.1  mrg               else
    525  1.1  mrg                 compare = inexact; /* cannot determine sign(t-f(x)) */
    526  1.1  mrg               if (((inexact == 0) && (compare != 0)) ||
    527  1.1  mrg                   ((inexact > 0) && (compare <= 0)) ||
    528  1.1  mrg                   ((inexact < 0) && (compare >= 0)))
    529  1.1  mrg                 {
    530  1.1  mrg                   printf ("Wrong inexact flag for rnd=%s and %s_q:\n"
    531  1.1  mrg                           "expected %d, got %d",
    532  1.1  mrg                           mpfr_print_rnd_mode (rnd), op, compare, inexact);
    533  1.1  mrg                   printf ("\narg1="); mpfr_print_binary (arg1);
    534  1.1  mrg                   printf ("\narg2="); mpq_out_str(stdout, 2, arg2);
    535  1.1  mrg                   printf ("\ndstl="); mpfr_print_binary (dst_big);
    536  1.1  mrg                   printf ("\ndsts="); mpfr_print_binary (dst_small);
    537  1.1  mrg                   printf ("\ntmp ="); mpfr_print_binary (tmp);
    538  1.1  mrg                   putchar('\n');
    539  1.1  mrg                   exit (1);
    540  1.1  mrg                 }
    541  1.1  mrg             }
    542  1.1  mrg         }
    543  1.1  mrg     }
    544  1.1  mrg 
    545  1.1  mrg   mpq_clear (arg2);
    546  1.1  mrg   mpfr_clears (arg1, dst_big, dst_small, tmp, (mpfr_ptr) 0);
    547  1.1  mrg }
    548  1.1  mrg 
    549  1.1  mrg static void
    550  1.1  mrg test_specialq (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int N,
    551  1.1  mrg                int (*mpfr_func)(mpfr_ptr, mpfr_srcptr, mpq_srcptr, mpfr_rnd_t),
    552  1.1  mrg                void (*mpq_func)(mpq_ptr, mpq_srcptr, mpq_srcptr),
    553  1.1  mrg                const char *op)
    554  1.1  mrg {
    555  1.1  mrg   mpfr_t fra, frb, frq;
    556  1.1  mrg   mpq_t  q1, q2, qr;
    557  1.1  mrg   unsigned int n;
    558  1.1  mrg   mpfr_prec_t prec;
    559  1.1  mrg 
    560  1.1  mrg   for (prec = p0 ; prec < p1 ; prec++)
    561  1.1  mrg     {
    562  1.1  mrg       mpfr_inits2 (prec, fra, frb, frq, (mpfr_ptr) 0);
    563  1.1  mrg       mpq_init (q1); mpq_init(q2); mpq_init (qr);
    564  1.1  mrg 
    565  1.1  mrg       for( n = 0 ; n < N ; n++)
    566  1.1  mrg         {
    567  1.1  mrg           mpq_set_ui(q1, randlimb(), randlimb() );
    568  1.1  mrg           mpq_set_ui(q2, randlimb(), randlimb() );
    569  1.1  mrg           mpq_canonicalize (q1);
    570  1.1  mrg           mpq_canonicalize (q2);
    571  1.1  mrg           mpq_func (qr, q1, q2);
    572  1.1  mrg           mpfr_set_q (fra, q1, MPFR_RNDD);
    573  1.1  mrg           mpfr_func (fra, fra, q2, MPFR_RNDD);
    574  1.1  mrg           mpfr_set_q (frb, q1, MPFR_RNDU);
    575  1.1  mrg           mpfr_func (frb, frb, q2, MPFR_RNDU);
    576  1.1  mrg           mpfr_set_q (frq, qr, MPFR_RNDN);
    577  1.1  mrg           /* We should have fra <= qr <= frb */
    578  1.1  mrg           if ( (mpfr_cmp(fra, frq) > 0) || (mpfr_cmp (frq, frb) > 0))
    579  1.1  mrg             {
    580  1.1  mrg               printf("Range error for prec=%lu and %s", prec, op);
    581  1.1  mrg               printf ("\nq1="); mpq_out_str(stdout, 2, q1);
    582  1.1  mrg               printf ("\nq2="); mpq_out_str(stdout, 2, q2);
    583  1.1  mrg               printf ("\nfr_dn="); mpfr_print_binary (fra);
    584  1.1  mrg               printf ("\nfr_q ="); mpfr_print_binary (frq);
    585  1.1  mrg               printf ("\nfr_up="); mpfr_print_binary (frb);
    586  1.1  mrg               putchar('\n');
    587  1.1  mrg               exit (1);
    588  1.1  mrg             }
    589  1.1  mrg         }
    590  1.1  mrg 
    591  1.1  mrg       mpq_clear (q1); mpq_clear (q2); mpq_clear (qr);
    592  1.1  mrg       mpfr_clears (fra, frb, frq, (mpfr_ptr) 0);
    593  1.1  mrg     }
    594  1.1  mrg }
    595  1.1  mrg 
    596  1.1  mrg int
    597  1.1  mrg main (int argc, char *argv[])
    598  1.1  mrg {
    599  1.1  mrg   tests_start_mpfr ();
    600  1.1  mrg 
    601  1.1  mrg   special ();
    602  1.1  mrg 
    603  1.1  mrg   test_specialz (mpfr_add_z, mpz_add, "add");
    604  1.1  mrg   test_specialz (mpfr_sub_z, mpz_sub, "sub");
    605  1.1  mrg   test_specialz (mpfr_mul_z, mpz_mul, "mul");
    606  1.1  mrg   test_genericz (2, 100, 100, mpfr_add_z, "add");
    607  1.1  mrg   test_genericz (2, 100, 100, mpfr_sub_z, "sub");
    608  1.1  mrg   test_genericz (2, 100, 100, mpfr_mul_z, "mul");
    609  1.1  mrg   test_genericz (2, 100, 100, mpfr_div_z, "div");
    610  1.1  mrg 
    611  1.1  mrg   test_genericq (2, 100, 100, mpfr_add_q, "add");
    612  1.1  mrg   test_genericq (2, 100, 100, mpfr_sub_q, "sub");
    613  1.1  mrg   test_genericq (2, 100, 100, mpfr_mul_q, "mul");
    614  1.1  mrg   test_genericq (2, 100, 100, mpfr_div_q, "div");
    615  1.1  mrg   test_specialq (2, 100, 100, mpfr_mul_q, mpq_mul, "mul");
    616  1.1  mrg   test_specialq (2, 100, 100, mpfr_div_q, mpq_div, "div");
    617  1.1  mrg   test_specialq (2, 100, 100, mpfr_add_q, mpq_add, "add");
    618  1.1  mrg   test_specialq (2, 100, 100, mpfr_sub_q, mpq_sub, "sub");
    619  1.1  mrg 
    620  1.1  mrg   test_cmp_z (2, 100, 100);
    621  1.1  mrg   test_cmp_q (2, 100, 100);
    622  1.1  mrg   test_cmp_f (2, 100, 100);
    623  1.1  mrg 
    624  1.1  mrg   check_for_zero ();
    625  1.1  mrg 
    626  1.1  mrg   tests_end_mpfr ();
    627  1.1  mrg   return 0;
    628  1.1  mrg }
    629  1.1  mrg 
    630