From d751d9d8edc2d1ea099674efecb04a5033428511 Mon Sep 17 00:00:00 2001
From: Mattias Andrée <maandree@kth.se>
Date: Fri, 13 May 2016 10:32:10 +0200
Subject: How you would calculate factorials efficiently
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Signed-off-by: Mattias Andrée <maandree@kth.se>
---
 doc/not-implemented.tex | 40 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 40 insertions(+)

(limited to 'doc/not-implemented.tex')

diff --git a/doc/not-implemented.tex b/doc/not-implemented.tex
index 3934802..68c5d98 100644
--- a/doc/not-implemented.tex
+++ b/doc/not-implemented.tex
@@ -174,6 +174,46 @@ important function for many combinatorial
 applications, therefore it may be implemented
 in the future if the demand is high enough.
 
+An efficient, yet not optimal, implementation
+of factorials that about halves the number of
+required multiplications compared to the naïve
+method can be derived from the observation
+
+\vspace{1em}
+\( \displaystyle{
+    n! = n!! ~ \lfloor n / 2 \rfloor! ~ 2^{\lfloor n / 2 \rfloor}
+}\), $n$ odd.
+\vspace{1em}
+
+\noindent
+The resulting algorithm can be expressed
+
+\begin{alltt}
+   void
+   fact(z_t r, uint64_t n)
+   \{
+       z_t p, f, two;
+       uint64_t *ns, s = 1, i = 1;
+       zinit(p), zinit(f), zinit(two);
+       zseti(r, 1), zseti(p, 1), zseti(f, n), zseti(two, 2);
+       ns = alloca(zbits(f) * sizeof(*ns));
+       while (n > 1) \{
+           if (n & 1) \{
+               ns[i++] = n;
+               s += n >>= 1;
+           \} else \{
+               zmul(r, r, (zsetu(f, n), f));
+               n -= 1;
+           \}
+       \}
+       for (zseti(f, 1); i-- > 0; zmul(r, r, p);)
+           for (n = ns[i]; zcmpu(f, n); zadd(f, f, two))
+               zmul(p, p, f);
+       zlsh(r, r, s);
+       zfree(two), zfree(f), zfree(p);
+   \}
+\end{alltt}
+
 
 \subsection{Subfactorial}
 \label{sec:Subfactorial}
-- 
cgit v1.2.3-70-g09d2

