1 /* mpfr_fac_ui -- factorial of a non-negative integer
3 Copyright 2001, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
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.
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.
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. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
26 /* The computation of n! is done by
31 /* FIXME: efficient problems with large arguments; see comments in gamma.c. */
34 mpfr_fac_ui (mpfr_ptr y
, unsigned long int x
, mpfr_rnd_t rnd_mode
)
36 mpfr_t t
; /* Variable of Intermediary Calculation*/
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 */
45 MPFR_SAVE_EXPO_DECL (expo
);
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 */
52 MPFR_SAVE_EXPO_MARK (expo
);
54 /* Initialisation of the Precision */
57 /* compute the size of intermediary variable */
58 Nt
= Ny
+ 2 * MPFR_INT_CEIL_LOG2 (x
) + 7;
60 mpfr_init2 (t
, Nt
); /* initialise of intermediary variable */
63 MPFR_ZIV_INIT (loop
, Nt
);
66 /* compute factorial */
67 inexact
= mpfr_set_ui (t
, 1, rnd
);
68 for (i
= 2 ; i
<= x
; i
++)
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? */
77 err
= Nt
- 1 - MPFR_INT_CEIL_LOG2 (Nt
);
79 round
= !inexact
|| mpfr_can_round (t
, err
, rnd
, MPFR_RNDZ
,
80 Ny
+ (rnd_mode
== MPFR_RNDN
));
82 if (MPFR_LIKELY (round
))
84 /* If inexact = 0, then t is exactly x!, so round is the
86 Otherwise, t != x! since we rounded to zero or away. */
87 round
= mpfr_set (y
, t
, rnd_mode
);
93 else if ((inexact
< 0 && round
<= 0)
94 || (inexact
> 0 && round
>= 0))
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
;
101 MPFR_ZIV_NEXT (loop
, Nt
);
102 mpfr_set_prec (t
, Nt
);
104 MPFR_ZIV_FREE (loop
);
107 MPFR_SAVE_EXPO_FREE (expo
);
108 return mpfr_check_range (y
, inexact
, rnd_mode
);