1/* tsum -- test file for the list summation function
2
3Copyright 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4Contributed by the AriC and Caramel 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 <stdlib.h>
24#include <stdio.h>
25#include "mpfr-test.h"
26
27static int
28check_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
38static int
39sum_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
57static mpfr_prec_t
58get_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
102static void
103algo_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 */
122static void
123test_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
161static void
162test_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
220static
221void 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
304int
305main (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}
321