1 /* Test file for mpfr_factorial. 2 3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4 Contributed by the Arenaire and Cacao 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 <stdio.h> 24 #include <stdlib.h> 25 26 #include "mpfr-test.h" 27 28 #define TEST_FUNCTION mpfr_fac_ui 29 30 static void 31 special (void) 32 { 33 mpfr_t x, y; 34 int inex; 35 36 mpfr_init (x); 37 mpfr_init (y); 38 39 mpfr_set_prec (x, 21); 40 mpfr_set_prec (y, 21); 41 mpfr_fac_ui (x, 119, MPFR_RNDZ); 42 mpfr_set_str_binary (y, "0.101111101110100110110E654"); 43 if (mpfr_cmp (x, y)) 44 { 45 printf ("Error in mpfr_fac_ui (119)\n"); 46 exit (1); 47 } 48 49 mpfr_set_prec (y, 206); 50 inex = mpfr_fac_ui (y, 767, MPFR_RNDN); 51 mpfr_set_prec (x, 206); 52 mpfr_set_str_binary (x, "0.110111100001000001101010010001000111000100000100111000010011100011011111001100011110101000111101101100110001001100110100001001111110000101010000100100011100010011101110000001000010001100010000101001111E6250"); 53 if (mpfr_cmp (x, y)) 54 { 55 printf ("Error in mpfr_fac_ui (767)\n"); 56 exit (1); 57 } 58 if (inex <= 0) 59 { 60 printf ("Wrong flag for mpfr_fac_ui (767)\n"); 61 exit (1); 62 } 63 64 mpfr_set_prec (y, 202); 65 mpfr_fac_ui (y, 69, MPFR_RNDU); 66 67 mpfr_clear (x); 68 mpfr_clear (y); 69 } 70 71 static void 72 test_int (void) 73 { 74 unsigned long n0 = 1, n1 = 80, n; 75 mpz_t f; 76 mpfr_t x, y; 77 mpfr_prec_t prec_f, p; 78 int r; 79 int inex1, inex2; 80 81 mpz_init (f); 82 mpfr_init (x); 83 mpfr_init (y); 84 85 mpz_fac_ui (f, n0 - 1); 86 for (n = n0; n <= n1; n++) 87 { 88 mpz_mul_ui (f, f, n); /* f = n! */ 89 prec_f = mpz_sizeinbase (f, 2) - mpz_scan1 (f, 0); 90 for (p = MPFR_PREC_MIN; p <= prec_f; p++) 91 { 92 mpfr_set_prec (x, p); 93 mpfr_set_prec (y, p); 94 for (r = 0; r < MPFR_RND_MAX; r++) 95 { 96 inex1 = mpfr_fac_ui (x, n, (mpfr_rnd_t) r); 97 inex2 = mpfr_set_z (y, f, (mpfr_rnd_t) r); 98 if (mpfr_cmp (x, y)) 99 { 100 printf ("Error for n=%lu prec=%lu rnd=%s\n", 101 n, (unsigned long) p, mpfr_print_rnd_mode ((mpfr_rnd_t) r)); 102 exit (1); 103 } 104 if ((inex1 < 0 && inex2 >= 0) || (inex1 == 0 && inex2 != 0) 105 || (inex1 > 0 && inex2 <= 0)) 106 { 107 printf ("Wrong inexact flag for n=%lu prec=%lu rnd=%s\n", 108 n, (unsigned long) p, mpfr_print_rnd_mode ((mpfr_rnd_t) r)); 109 exit (1); 110 } 111 } 112 } 113 } 114 115 mpz_clear (f); 116 mpfr_clear (x); 117 mpfr_clear (y); 118 } 119 120 static void 121 overflowed_fac0 (void) 122 { 123 mpfr_t x, y; 124 int inex, rnd, err = 0; 125 mpfr_exp_t old_emax; 126 127 old_emax = mpfr_get_emax (); 128 129 mpfr_init2 (x, 8); 130 mpfr_init2 (y, 8); 131 132 mpfr_set_ui (y, 1, MPFR_RNDN); 133 mpfr_nextbelow (y); 134 set_emax (0); /* 1 is not representable. */ 135 RND_LOOP (rnd) 136 { 137 mpfr_clear_flags (); 138 inex = mpfr_fac_ui (x, 0, (mpfr_rnd_t) rnd); 139 if (! mpfr_overflow_p ()) 140 { 141 printf ("Error in overflowed_fac0 (rnd = %s):\n" 142 " The overflow flag is not set.\n", 143 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 144 err = 1; 145 } 146 if (rnd == MPFR_RNDZ || rnd == MPFR_RNDD) 147 { 148 if (inex >= 0) 149 { 150 printf ("Error in overflowed_fac0 (rnd = %s):\n" 151 " The inexact value must be negative.\n", 152 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 153 err = 1; 154 } 155 if (! mpfr_equal_p (x, y)) 156 { 157 printf ("Error in overflowed_fac0 (rnd = %s):\n" 158 " Got ", mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 159 mpfr_print_binary (x); 160 printf (" instead of 0.11111111E0.\n"); 161 err = 1; 162 } 163 } 164 else 165 { 166 if (inex <= 0) 167 { 168 printf ("Error in overflowed_fac0 (rnd = %s):\n" 169 " The inexact value must be positive.\n", 170 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 171 err = 1; 172 } 173 if (! (mpfr_inf_p (x) && MPFR_SIGN (x) > 0)) 174 { 175 printf ("Error in overflowed_fac0 (rnd = %s):\n" 176 " Got ", mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 177 mpfr_print_binary (x); 178 printf (" instead of +Inf.\n"); 179 err = 1; 180 } 181 } 182 } 183 set_emax (old_emax); 184 185 if (err) 186 exit (1); 187 mpfr_clear (x); 188 mpfr_clear (y); 189 } 190 191 int 192 main (int argc, char *argv[]) 193 { 194 unsigned int prec, err, yprec, n, k, zeros; 195 int rnd; 196 mpfr_t x, y, z, t; 197 int inexact; 198 199 tests_start_mpfr (); 200 201 special (); 202 203 test_int (); 204 205 mpfr_init (x); 206 mpfr_init (y); 207 mpfr_init (z); 208 mpfr_init (t); 209 210 mpfr_fac_ui (y, 0, MPFR_RNDN); 211 212 if (mpfr_cmp_ui (y, 1)) 213 { 214 printf ("mpfr_fac_ui(0) does not give 1\n"); 215 exit (1); 216 } 217 218 for (prec = 2; prec <= 100; prec++) 219 { 220 mpfr_set_prec (x, prec); 221 mpfr_set_prec (z, prec); 222 mpfr_set_prec (t, prec); 223 yprec = prec + 10; 224 mpfr_set_prec (y, yprec); 225 226 for (n = 0; n < 50; n++) 227 for (rnd = 0; rnd < MPFR_RND_MAX; rnd++) 228 { 229 inexact = mpfr_fac_ui (y, n, (mpfr_rnd_t) rnd); 230 err = (rnd == MPFR_RNDN) ? yprec + 1 : yprec; 231 if (mpfr_can_round (y, err, (mpfr_rnd_t) rnd, (mpfr_rnd_t) rnd, prec)) 232 { 233 mpfr_set (t, y, (mpfr_rnd_t) rnd); 234 inexact = mpfr_fac_ui (z, n, (mpfr_rnd_t) rnd); 235 /* fact(n) ends with floor(n/2)+floor(n/4)+... zeros */ 236 for (k=n/2, zeros=0; k; k >>= 1) 237 zeros += k; 238 if (MPFR_EXP(y) <= (mpfr_exp_t) (prec + zeros)) 239 /* result should be exact */ 240 { 241 if (inexact) 242 { 243 printf ("Wrong inexact flag: expected exact\n"); 244 exit (1); 245 } 246 } 247 else /* result is inexact */ 248 { 249 if (!inexact) 250 { 251 printf ("Wrong inexact flag: expected inexact\n"); 252 printf ("n=%u prec=%u\n", n, prec); 253 mpfr_print_binary(z); puts (""); 254 exit (1); 255 } 256 } 257 if (mpfr_cmp (t, z)) 258 { 259 printf ("results differ for x="); 260 mpfr_out_str (stdout, 2, prec, x, MPFR_RNDN); 261 printf (" prec=%u rnd_mode=%s\n", prec, 262 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 263 printf (" got "); 264 mpfr_out_str (stdout, 2, prec, z, MPFR_RNDN); 265 puts (""); 266 printf (" expected "); 267 mpfr_out_str (stdout, 2, prec, t, MPFR_RNDN); 268 puts (""); 269 printf (" approximation was "); 270 mpfr_print_binary (y); 271 puts (""); 272 exit (1); 273 } 274 } 275 } 276 } 277 278 mpfr_clear (x); 279 mpfr_clear (y); 280 mpfr_clear (z); 281 mpfr_clear (t); 282 283 overflowed_fac0 (); 284 285 tests_end_mpfr (); 286 return 0; 287 } 288