1 /* tsum -- test file for the list summation function 2 3 Copyright 2004-2016 Free Software Foundation, Inc. 4 Contributed by the AriC and Caramba 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 mpfr_prec_t prec = MPFR_PREC_MIN; 130 131 /* Init stuff */ 132 tab = (mpfr_t *) tests_allocate (n * sizeof (mpfr_t)); 133 for (i = 0; i < n; i++) 134 mpfr_init2 (tab[i], f); 135 tabtmp = (mpfr_ptr *) tests_allocate (n * sizeof(mpfr_ptr)); 136 perm = (mpfr_srcptr *) tests_allocate (n * sizeof(mpfr_srcptr)); 137 138 for (i = 0; i < n; i++) 139 { 140 mpfr_urandomb (tab[i], RANDS); 141 tabtmp[i] = tab[i]; 142 } 143 144 mpfr_sum_sort ((mpfr_srcptr *)tabtmp, n, perm, &prec); 145 146 if (check_is_sorted (n, perm) == 0) 147 { 148 printf ("mpfr_sum_sort incorrect.\n"); 149 for (i = 0; i < n; i++) 150 mpfr_dump (perm[i]); 151 exit (1); 152 } 153 154 /* Clear stuff */ 155 for (i = 0; i < n; i++) 156 mpfr_clear (tab[i]); 157 tests_free (tab, n * sizeof (mpfr_t)); 158 tests_free (tabtmp, n * sizeof(mpfr_ptr)); 159 tests_free (perm, n * sizeof(mpfr_srcptr)); 160 } 161 162 static void 163 test_sum (mpfr_prec_t f, unsigned long n) 164 { 165 mpfr_t sum, real_sum, real_non_rounded; 166 mpfr_t *tab; 167 unsigned long i; 168 int rnd_mode; 169 170 /* Init */ 171 tab = (mpfr_t *) tests_allocate (n * sizeof(mpfr_t)); 172 for (i = 0; i < n; i++) 173 mpfr_init2 (tab[i], f); 174 mpfr_inits2 (f, sum, real_sum, real_non_rounded, (mpfr_ptr) 0); 175 176 /* First Uniform */ 177 for (i = 0; i < n; i++) 178 mpfr_urandomb (tab[i], RANDS); 179 algo_exact (real_non_rounded, tab, n, f); 180 for (rnd_mode = 0; rnd_mode < MPFR_RND_MAX; rnd_mode++) 181 { 182 sum_tab (sum, tab, n, (mpfr_rnd_t) rnd_mode); 183 mpfr_set (real_sum, real_non_rounded, (mpfr_rnd_t) rnd_mode); 184 if (mpfr_cmp (real_sum, sum) != 0) 185 { 186 printf ("mpfr_sum incorrect.\n"); 187 mpfr_dump (real_sum); 188 mpfr_dump (sum); 189 exit (1); 190 } 191 } 192 193 /* Then non uniform */ 194 for (i = 0; i < n; i++) 195 { 196 mpfr_urandomb (tab[i], RANDS); 197 if (! mpfr_zero_p (tab[i])) 198 mpfr_set_exp (tab[i], randlimb () % 1000); 199 } 200 algo_exact (real_non_rounded, tab, n, f); 201 for (rnd_mode = 0; rnd_mode < MPFR_RND_MAX; rnd_mode++) 202 { 203 sum_tab (sum, tab, n, (mpfr_rnd_t) rnd_mode); 204 mpfr_set (real_sum, real_non_rounded, (mpfr_rnd_t) rnd_mode); 205 if (mpfr_cmp (real_sum, sum) != 0) 206 { 207 printf ("mpfr_sum incorrect.\n"); 208 mpfr_dump (real_sum); 209 mpfr_dump (sum); 210 exit (1); 211 } 212 } 213 214 /* Clear stuff */ 215 for (i = 0; i < n; i++) 216 mpfr_clear (tab[i]); 217 mpfr_clears (sum, real_sum, real_non_rounded, (mpfr_ptr) 0); 218 tests_free (tab, n * sizeof(mpfr_t)); 219 } 220 221 static 222 void check_special (void) 223 { 224 mpfr_t tab[3], r; 225 mpfr_ptr tabp[3]; 226 int i; 227 228 mpfr_inits (tab[0], tab[1], tab[2], r, (mpfr_ptr) 0); 229 tabp[0] = tab[0]; 230 tabp[1] = tab[1]; 231 tabp[2] = tab[2]; 232 233 i = mpfr_sum (r, tabp, 0, MPFR_RNDN); 234 if (!MPFR_IS_ZERO (r) || !MPFR_IS_POS (r) || i != 0) 235 { 236 printf ("Special case n==0 failed!\n"); 237 exit (1); 238 } 239 240 mpfr_set_ui (tab[0], 42, MPFR_RNDN); 241 i = mpfr_sum (r, tabp, 1, MPFR_RNDN); 242 if (mpfr_cmp_ui (r, 42) || i != 0) 243 { 244 printf ("Special case n==1 failed!\n"); 245 exit (1); 246 } 247 248 mpfr_set_ui (tab[1], 17, MPFR_RNDN); 249 MPFR_SET_NAN (tab[2]); 250 i = mpfr_sum (r, tabp, 3, MPFR_RNDN); 251 if (!MPFR_IS_NAN (r) || i != 0) 252 { 253 printf ("Special case NAN failed!\n"); 254 exit (1); 255 } 256 257 MPFR_SET_INF (tab[2]); 258 MPFR_SET_POS (tab[2]); 259 i = mpfr_sum (r, tabp, 3, MPFR_RNDN); 260 if (!MPFR_IS_INF (r) || !MPFR_IS_POS (r) || i != 0) 261 { 262 printf ("Special case +INF failed!\n"); 263 exit (1); 264 } 265 266 MPFR_SET_INF (tab[2]); 267 MPFR_SET_NEG (tab[2]); 268 i = mpfr_sum (r, tabp, 3, MPFR_RNDN); 269 if (!MPFR_IS_INF (r) || !MPFR_IS_NEG (r) || i != 0) 270 { 271 printf ("Special case -INF failed!\n"); 272 exit (1); 273 } 274 275 MPFR_SET_ZERO (tab[1]); 276 i = mpfr_sum (r, tabp, 2, MPFR_RNDN); 277 if (mpfr_cmp_ui (r, 42) || i != 0) 278 { 279 printf ("Special case 42+0 failed!\n"); 280 exit (1); 281 } 282 283 MPFR_SET_NAN (tab[0]); 284 i = mpfr_sum (r, tabp, 3, MPFR_RNDN); 285 if (!MPFR_IS_NAN (r) || i != 0) 286 { 287 printf ("Special case NAN+0+-INF failed!\n"); 288 exit (1); 289 } 290 291 mpfr_set_inf (tab[0], 1); 292 mpfr_set_ui (tab[1], 59, MPFR_RNDN); 293 mpfr_set_inf (tab[2], -1); 294 i = mpfr_sum (r, tabp, 3, MPFR_RNDN); 295 if (!MPFR_IS_NAN (r) || i != 0) 296 { 297 printf ("Special case +INF + 59 +-INF failed!\n"); 298 exit (1); 299 } 300 301 mpfr_clears (tab[0], tab[1], tab[2], r, (mpfr_ptr) 0); 302 } 303 304 /* bug reported by Joseph S. Myers on 2013-10-27 305 https://sympa.inria.fr/sympa/arc/mpfr/2013-10/msg00015.html */ 306 static void 307 bug20131027 (void) 308 { 309 mpfr_t r, t[4]; 310 mpfr_ptr p[4]; 311 char *s[4] = { 312 "0x1p1000", 313 "-0x0.fffffffffffff80000000000000001p1000", 314 "-0x1p947", 315 "0x1p880" 316 }; 317 int i; 318 319 mpfr_init2 (r, 53); 320 for (i = 0; i < 4; i++) 321 { 322 mpfr_init2 (t[i], i == 0 ? 53 : 1000); 323 mpfr_set_str (t[i], s[i], 0, MPFR_RNDN); 324 p[i] = t[i]; 325 } 326 mpfr_sum (r, p, 4, MPFR_RNDN); 327 328 if (MPFR_NOTZERO (r)) 329 { 330 printf ("mpfr_sum incorrect in bug20131027: expected 0, got\n"); 331 mpfr_dump (r); 332 exit (1); 333 } 334 335 for (i = 0; i < 4; i++) 336 mpfr_clear (t[i]); 337 mpfr_clear (r); 338 } 339 340 int 341 main (void) 342 { 343 mpfr_prec_t p; 344 unsigned long n; 345 346 tests_start_mpfr (); 347 348 check_special (); 349 bug20131027 (); 350 test_sort (1764, 1026); 351 for (p = 2 ; p < 444 ; p += 17) 352 for (n = 2 ; n < 1026 ; n += 42 + p) 353 test_sum (p, n); 354 355 tests_end_mpfr (); 356 return 0; 357 } 358