tests/tsum.c

   1 /* tsum -- test file for the list summation function
   2
   3 Copyright 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
   4 Contributed by the AriC and Caramel 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 http://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 <stdlib.h>
  24 #include <stdio.h>
  25 #include "mpfr-test.h"
  26
  27 static int
  28 check_is_sorted (unsigned long n, mpfr_srcptr *perm)
  29 {
  30   unsigned long i;
  31
  32   for (i = 0; i < n - 1; i++)
  33     if (MPFR_GET_EXP(perm[i]) < MPFR_GET_EXP(perm[i+1]))
  34       return 0;
  35   return 1;
  36 }
  37
  38 static int
  39 sum_tab (mpfr_ptr ret, mpfr_t *tab, unsigned long n, mpfr_rnd_t rnd)
  40 {
  41   mpfr_ptr *tabtmp;
  42   unsigned long i;
  43   int inexact;
  44   MPFR_TMP_DECL(marker);
  45
  46   MPFR_TMP_MARK(marker);
  47   tabtmp = (mpfr_ptr *) MPFR_TMP_ALLOC(n * sizeof(mpfr_srcptr));
  48   for (i = 0; i < n; i++)
  49     tabtmp[i] = tab[i];
  50
  51   inexact = mpfr_sum (ret, tabtmp, n, rnd);
  52   MPFR_TMP_FREE(marker);
  53   return inexact;
  54 }
  55
  56
  57 static mpfr_prec_t
  58 get_prec_max (mpfr_t *tab, unsigned long n, mpfr_prec_t f)
  59 {
  60   mpfr_prec_t res;
  61   mpfr_exp_t min, max;
  62   unsigned long i;
  63
  64   i = 0;
  65   while (MPFR_IS_ZERO (tab[i]))
  66     {
  67       i++;
  68       if (i == n)
  69         return MPFR_PREC_MIN;  /* all values are 0 */
  70     }
  71
  72   if (! mpfr_check (tab[i]))
  73     {
  74       printf ("tab[%lu] is not valid.\n", i);
  75       exit (1);
  76     }
  77   MPFR_ASSERTN (MPFR_IS_FP (tab[i]));
  78   min = max = MPFR_GET_EXP(tab[i]);
  79   for (i++; i < n; i++)
  80     {
  81       if (! mpfr_check (tab[i]))
  82         {
  83           printf ("tab[%lu] is not valid.\n", i);
  84           exit (1);
  85         }
  86       MPFR_ASSERTN (MPFR_IS_FP (tab[i]));
  87       if (! MPFR_IS_ZERO (tab[i]))
  88         {
  89           if (MPFR_GET_EXP(tab[i]) > max)
  90             max = MPFR_GET_EXP(tab[i]);
  91           if (MPFR_GET_EXP(tab[i]) < min)
  92             min = MPFR_GET_EXP(tab[i]);
  93         }
  94     }
  95   res = max - min;
  96   res += f;
  97   res += __gmpfr_ceil_log2 (n) + 1;
  98   return res;
  99 }
 100
 101
 102 static void
 103 algo_exact (mpfr_t somme, mpfr_t *tab, unsigned long n, mpfr_prec_t f)
 104 {
 105   unsigned long i;
 106   mpfr_prec_t prec_max;
 107
 108   prec_max = get_prec_max(tab, n, f);
 109   mpfr_set_prec (somme, prec_max);
 110   mpfr_set_ui (somme, 0, MPFR_RNDN);
 111   for (i = 0; i < n; i++)
 112     {
 113       if (mpfr_add(somme, somme, tab[i], MPFR_RNDN))
 114         {
 115           printf ("FIXME: algo_exact is buggy.\n");
 116           exit (1);
 117         }
 118     }
 119 }
 120
 121 /* Test the sorting function */
 122 static void
 123 test_sort (mpfr_prec_t f, unsigned long n)
 124 {
 125   mpfr_t *tab;
 126   mpfr_ptr *tabtmp;
 127   mpfr_srcptr *perm;
 128   unsigned long i;
 129
 130   /* Init stuff */
 131   tab = (mpfr_t *) (*__gmp_allocate_func) (n * sizeof (mpfr_t));
 132   for (i = 0; i < n; i++)
 133     mpfr_init2 (tab[i], f);
 134   tabtmp = (mpfr_ptr *) (*__gmp_allocate_func) (n * sizeof(mpfr_ptr));
 135   perm = (mpfr_srcptr *) (*__gmp_allocate_func) (n * sizeof(mpfr_srcptr));
 136
 137   for (i = 0; i < n; i++)
 138     {
 139       mpfr_urandomb (tab[i], RANDS);
 140       tabtmp[i] = tab[i];
 141     }
 142
 143   mpfr_sum_sort ((mpfr_srcptr *)tabtmp, n, perm);
 144
 145   if (check_is_sorted (n, perm) == 0)
 146     {
 147       printf ("mpfr_sum_sort incorrect.\n");
 148       for (i = 0; i < n; i++)
 149         mpfr_dump (perm[i]);
 150       exit (1);
 151     }
 152
 153   /* Clear stuff */
 154   for (i = 0; i < n; i++)
 155     mpfr_clear (tab[i]);
 156   (*__gmp_free_func) (tab, n * sizeof (mpfr_t));
 157   (*__gmp_free_func) (tabtmp, n * sizeof(mpfr_ptr));
 158   (*__gmp_free_func) (perm, n * sizeof(mpfr_srcptr));
 159 }
 160
 161 static void
 162 test_sum (mpfr_prec_t f, unsigned long n)
 163 {
 164   mpfr_t sum, real_sum, real_non_rounded;
 165   mpfr_t *tab;
 166   unsigned long i;
 167   int rnd_mode;
 168
 169   /* Init */
 170   tab = (mpfr_t *) (*__gmp_allocate_func) (n * sizeof(mpfr_t));
 171   for (i = 0; i < n; i++)
 172     mpfr_init2 (tab[i], f);
 173   mpfr_inits2 (f, sum, real_sum, real_non_rounded, (mpfr_ptr) 0);
 174
 175   /* First Uniform */
 176   for (i = 0; i < n; i++)
 177     mpfr_urandomb (tab[i], RANDS);
 178   algo_exact (real_non_rounded, tab, n, f);
 179   for (rnd_mode = 0; rnd_mode < MPFR_RND_MAX; rnd_mode++)
 180     {
 181       sum_tab (sum, tab, n, (mpfr_rnd_t) rnd_mode);
 182       mpfr_set (real_sum, real_non_rounded, (mpfr_rnd_t) rnd_mode);
 183       if (mpfr_cmp (real_sum, sum) != 0)
 184         {
 185           printf ("mpfr_sum incorrect.\n");
 186           mpfr_dump (real_sum);
 187           mpfr_dump (sum);
 188           exit (1);
 189         }
 190     }
 191
 192   /* Then non uniform */
 193   for (i = 0; i < n; i++)
 194     {
 195       mpfr_urandomb (tab[i], RANDS);
 196       if (! mpfr_zero_p (tab[i]))
 197         mpfr_set_exp (tab[i], randlimb () % 1000);
 198     }
 199   algo_exact (real_non_rounded, tab, n, f);
 200   for (rnd_mode = 0; rnd_mode < MPFR_RND_MAX; rnd_mode++)
 201     {
 202       sum_tab (sum, tab, n, (mpfr_rnd_t) rnd_mode);
 203       mpfr_set (real_sum, real_non_rounded, (mpfr_rnd_t) rnd_mode);
 204       if (mpfr_cmp (real_sum, sum) != 0)
 205         {
 206           printf ("mpfr_sum incorrect.\n");
 207           mpfr_dump (real_sum);
 208           mpfr_dump (sum);
 209           exit (1);
 210         }
 211     }
 212
 213   /* Clear stuff */
 214   for (i = 0; i < n; i++)
 215     mpfr_clear (tab[i]);
 216   mpfr_clears (sum, real_sum, real_non_rounded, (mpfr_ptr) 0);
 217   (*__gmp_free_func) (tab, n * sizeof(mpfr_t));
 218 }
 219
 220 static
 221 void check_special (void)
 222 {
 223   mpfr_t tab[3], r;
 224   mpfr_ptr tabp[3];
 225   int i;
 226
 227   mpfr_inits (tab[0], tab[1], tab[2], r, (mpfr_ptr) 0);
 228   tabp[0] = tab[0];
 229   tabp[1] = tab[1];
 230   tabp[2] = tab[2];
 231
 232   i = mpfr_sum (r, tabp, 0, MPFR_RNDN);
 233   if (!MPFR_IS_ZERO (r) || !MPFR_IS_POS (r) || i != 0)
 234     {
 235       printf ("Special case n==0 failed!\n");
 236       exit (1);
 237     }
 238
 239   mpfr_set_ui (tab[0], 42, MPFR_RNDN);
 240   i = mpfr_sum (r, tabp, 1, MPFR_RNDN);
 241   if (mpfr_cmp_ui (r, 42) || i != 0)
 242     {
 243       printf ("Special case n==1 failed!\n");
 244       exit (1);
 245     }
 246
 247   mpfr_set_ui (tab[1], 17, MPFR_RNDN);
 248   MPFR_SET_NAN (tab[2]);
 249   i = mpfr_sum (r, tabp, 3, MPFR_RNDN);
 250   if (!MPFR_IS_NAN (r) || i != 0)
 251     {
 252       printf ("Special case NAN failed!\n");
 253       exit (1);
 254     }
 255
 256   MPFR_SET_INF (tab[2]);
 257   MPFR_SET_POS (tab[2]);
 258   i = mpfr_sum (r, tabp, 3, MPFR_RNDN);
 259   if (!MPFR_IS_INF (r) || !MPFR_IS_POS (r) || i != 0)
 260     {
 261       printf ("Special case +INF failed!\n");
 262       exit (1);
 263     }
 264
 265   MPFR_SET_INF (tab[2]);
 266   MPFR_SET_NEG (tab[2]);
 267   i = mpfr_sum (r, tabp, 3, MPFR_RNDN);
 268   if (!MPFR_IS_INF (r) || !MPFR_IS_NEG (r) || i != 0)
 269     {
 270       printf ("Special case -INF failed!\n");
 271       exit (1);
 272     }
 273
 274   MPFR_SET_ZERO (tab[1]);
 275   i = mpfr_sum (r, tabp, 2, MPFR_RNDN);
 276   if (mpfr_cmp_ui (r, 42) || i != 0)
 277     {
 278       printf ("Special case 42+0 failed!\n");
 279       exit (1);
 280     }
 281
 282   MPFR_SET_NAN (tab[0]);
 283   i = mpfr_sum (r, tabp, 3, MPFR_RNDN);
 284   if (!MPFR_IS_NAN (r) || i != 0)
 285     {
 286       printf ("Special case NAN+0+-INF failed!\n");
 287       exit (1);
 288     }
 289
 290   mpfr_set_inf (tab[0], 1);
 291   mpfr_set_ui  (tab[1], 59, MPFR_RNDN);
 292   mpfr_set_inf (tab[2], -1);
 293   i = mpfr_sum (r, tabp, 3, MPFR_RNDN);
 294   if (!MPFR_IS_NAN (r) || i != 0)
 295     {
 296       printf ("Special case +INF + 59 +-INF failed!\n");
 297       exit (1);
 298     }
 299
 300   mpfr_clears (tab[0], tab[1], tab[2], r, (mpfr_ptr) 0);
 301 }
 302
 303
 304 int
 305 main (void)
 306 {
 307   mpfr_prec_t p;
 308   unsigned long n;
 309
 310   tests_start_mpfr ();
 311
 312   check_special ();
 313   test_sort (1764, 1026);
 314   for (p = 2 ; p < 444 ; p += 17)
 315     for (n = 2 ; n < 1026 ; n += 42 + p)
 316       test_sum (p, n);
 317
 318   tests_end_mpfr ();
 319   return 0;
 320 }