From 80823faac6955cb42422f5d39e0db17c326e2177 Mon Sep 17 00:00:00 2001 From: Sam Hocevar Date: Tue, 31 Oct 2017 14:19:44 +0100 Subject: [PATCH] Make fast real factorial more generic. This new version allows to compute odd and even double factorials without harming performance for the standard factorial. --- src/math/real.cpp | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/math/real.cpp b/src/math/real.cpp index 7b841be5..c822c208 100644 --- a/src/math/real.cpp +++ b/src/math/real.cpp @@ -861,16 +861,22 @@ template<> real pow(real const &x, real const &y) return is_odd ? exp(y * log(-x)) : -exp(y * log(-x)); } -static real fast_fact(unsigned int x) +/* A fast factorial implementation for small numbers. An optional + * step argument allows to compute double factorials (i.e. with + * only the odd or the even terms. */ +static real fast_fact(unsigned int x, unsigned int step = 1) { - real ret = real::R_1(); + unsigned int start = (x + step - 1) % step + 1; + real ret = start ? real(start) : real(step); uint64_t multiplier = 1; - for (unsigned int i = 1, exponent = 0;;) + for (unsigned int i = start, exponent = 0;;) { - if (i++ >= x) + if (i >= x) return ldexp(ret * multiplier, exponent); + i += step; + /* Accumulate the power of two part in the exponent */ unsigned int tmp = i; while ((tmp & 1) == 0)