1/* mpfr_fac_ui -- factorial of a non-negative integer 2 3Copyright 2001, 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#define MPFR_NEED_LONGLONG_H 24#include "mpfr-impl.h" 25 26 /* The computation of n! is done by 27 28 n!=prod^{n}_{i=1}i 29 */ 30 31/* FIXME: efficient problems with large arguments; see comments in gamma.c. */ 32 33int 34mpfr_fac_ui (mpfr_ptr y, unsigned long int x, mpfr_rnd_t rnd_mode) 35{ 36 mpfr_t t; /* Variable of Intermediary Calculation*/ 37 unsigned long i; 38 int round, inexact; 39 40 mpfr_prec_t Ny; /* Precision of output variable */ 41 mpfr_prec_t Nt; /* Precision of Intermediary Calculation variable */ 42 mpfr_prec_t err; /* Precision of error */ 43 44 mpfr_rnd_t rnd; 45 MPFR_SAVE_EXPO_DECL (expo); 46 MPFR_ZIV_DECL (loop); 47 48 /***** test x = 0 and x == 1******/ 49 if (MPFR_UNLIKELY (x <= 1)) 50 return mpfr_set_ui (y, 1, rnd_mode); /* 0! = 1 and 1! = 1 */ 51 52 MPFR_SAVE_EXPO_MARK (expo); 53 54 /* Initialisation of the Precision */ 55 Ny = MPFR_PREC (y); 56 57 /* compute the size of intermediary variable */ 58 Nt = Ny + 2 * MPFR_INT_CEIL_LOG2 (x) + 7; 59 60 mpfr_init2 (t, Nt); /* initialise of intermediary variable */ 61 62 rnd = MPFR_RNDZ; 63 MPFR_ZIV_INIT (loop, Nt); 64 for (;;) 65 { 66 /* compute factorial */ 67 inexact = mpfr_set_ui (t, 1, rnd); 68 for (i = 2 ; i <= x ; i++) 69 { 70 round = mpfr_mul_ui (t, t, i, rnd); 71 /* assume the first inexact product gives the sign 72 of difference: is that always correct? */ 73 if (inexact == 0) 74 inexact = round; 75 } 76 77 err = Nt - 1 - MPFR_INT_CEIL_LOG2 (Nt); 78 79 round = !inexact || mpfr_can_round (t, err, rnd, MPFR_RNDZ, 80 Ny + (rnd_mode == MPFR_RNDN)); 81 82 if (MPFR_LIKELY (round)) 83 { 84 /* If inexact = 0, then t is exactly x!, so round is the 85 correct inexact flag. 86 Otherwise, t != x! since we rounded to zero or away. */ 87 round = mpfr_set (y, t, rnd_mode); 88 if (inexact == 0) 89 { 90 inexact = round; 91 break; 92 } 93 else if ((inexact < 0 && round <= 0) 94 || (inexact > 0 && round >= 0)) 95 break; 96 else /* inexact and round have opposite signs: we cannot 97 compute the inexact flag. Restart using the 98 symmetric rounding. */ 99 rnd = (rnd == MPFR_RNDZ) ? MPFR_RNDU : MPFR_RNDZ; 100 } 101 MPFR_ZIV_NEXT (loop, Nt); 102 mpfr_set_prec (t, Nt); 103 } 104 MPFR_ZIV_FREE (loop); 105 106 mpfr_clear (t); 107 MPFR_SAVE_EXPO_FREE (expo); 108 return mpfr_check_range (y, inexact, rnd_mode); 109} 110 111 112 113 114