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 } 321