title: Fast exact bigint factorial
body: <p>I have a fixed-point bignumber library and want to implement fast factorial with no precision loss.</p>
  <p>After some math tricks on paper I got this formula:</p>
  <pre><code>(4N)!=((2N)!).((2N)!).{ (2N+1).(2N+3).(2N+5)...(4N-1) }.(2^N)/(N!)
  </code></pre>
  <p>This is already pretty fast, and with some programming tricks the complexity nears <code>~ O(log(n))</code>.</p>
  <p>To be clear, my current implementation is this:</p>
  <pre><code>//---------------------------------------------------------------------------
  longnum fact(const DWORD &amp;x,longnum &amp;h) // h return (x&gt;&gt;1)! to speed up computation
      {
      if (x==0) { h=1; return  1; }
      if (x==1) { h=1; return  1; }
      if (x==2) { h=1; return  2; }
      if (x==3) { h=1; return  6; }
      if (x==4) { h=2; return 24; }
      int N4,N2,N,i; longnum c,q;
      N=(x&gt;&gt;2);
      N2=N&lt;&lt;1;
      N4=N&lt;&lt;2;
      h=fact(N2,q);                                          // get 2N! and N!
      c=h*h; for (i=(N2+1)|1;i&lt;=N4;i+=2) c*=i; c/=q;         // c= ((2N!)^2)*T1 / N!
      for (i=N4+1;i&lt;=x;i++) c*=i; c.round(); c&lt;&lt;=N  ;        // convert 4N! -&gt; x!, cut off precision losses
      for (i=(N2+1)|1,N2=x&gt;&gt;1;i&lt;=N2;i++) h*=i; h.round();    // convert 2N! -&gt; (x/2)!, cut off precision losses
      return c;
      }
  //---------------------------------------------------------------------------
  longnum fact(const DWORD &amp;x)
      {
      longnum tmp;
      return fact(x,tmp);
      }
  //---------------------------------------------------------------------------
  </code></pre>
  <p>Now my question:</p>
  <ol>
  <li><p>Is there a <strong>fast way to obtain</strong> <code>N!</code> from this <strong>term:</strong> <code>T1 = { (2N+1).(2N+3).(2N+5)...(4N-1) }</code>?</p>
  <p>Already answered.</p>
  </li>
  </ol>
  <p>So to be clear, I need to extract this <strong>unknown term:</strong></p>
  <pre><code>T2 = (4N)! / (((2N)!).((2N)!))
  </code></pre>
  <p>so:</p>
  <pre><code>(4N)! = (((2N)!).((2N)!)).T2
  </code></pre>
  <p>This would help a lot because then it would not be needed to compute <code>.../(N!)</code> for factorial.</p>
  <p>The <code>T1</code> term is always integer-decomposable to this:</p>
  <pre><code>T1 = T2 * N!
  </code></pre>
  <p>Finally, it hit me :) I have done a little program for primes decomposition of factorials and then suddenly all becomes much clearer:</p>
  <pre><code>4! =  2!.2!.(2^1).(3^1) = 24
  8! =  4!.4!.(2^1).(5^1).(7^1) = 40320
  12! =  6!.6!.(2^2).(3^1).(7^1).(11^1) = 479001600
  16! =  8!.8!.(2^1).(3^2).(5^1).(11^1).(13^1) = 20922789888000
  20! =  10!.10!.(2^2).(11^1).(13^1).(17^1).(19^1) = 2432902008176640000
  24! =  12!.12!.(2^2).(7^1).(13^1).(17^1).(19^1).(23^1) = 620448401733239439360000
  28! =  14!.14!.(2^3).(3^3).(5^2).(17^1).(19^1).(23^1) = 304888344611713860501504000000
  32! =  16!.16!.(2^1).(3^2).(5^1).(17^1).(19^1).(23^1).(29^1).(31^1) = 263130836933693530167218012160000000
  36! =  18!.18!.(2^2).(3^1).(5^2).(7^1).(11^1).(19^1).(23^1).(29^1).(31^1) = 371993326789901217467999448150835200000000
  40! =  20!.20!.(2^2).(3^2).(5^1).(7^1).(11^1).(13^1).(23^1).(29^1).(31^1).(37^1) = 815915283247897734345611269596115894272000000000
  </code></pre>
  <p>After analyzing the prime exponents of the <code>T2</code> term (the rest after half factorials ^ 2) I derive the formula for them:</p>
  <pre><code>T2(4N) = multiplication(i=2,3,5,7,11,13,17,...) of ( i ^ sum(j=1,2,3,4,5,...) of (4N/(i^j))-(2N/(i^j)) )
  </code></pre>
  <ul>
  <li>where multiplication is through all <code>primes &lt;= 4N</code></li>
  <li>where sumation is until <code>i^j &lt;= 4N</code></li>
  </ul>
  <p>The problem is that the divisions <code>4N/(i^j)</code> and <code>2N/(i^j)</code>  must be done in <strong>integer math</strong> so they <strong>cannot be simplified easily</strong>.</p>
  <p>So I have another question:</p>
  <ol>
  <li><p><strong>How can I compute this: <code>exponent(i) = sum(j=1,2,3,4,5,...) of (N/(i^j))</code> effectively?</strong></p>
  <p><code>i</code> is any prime where <code>i&lt;=N</code>. It should be easy.</p>
  </li>
  </ol>
  <p>Now I calculate the exponent <code>e</code> for prime <code>i</code> inside the <code>T2(N)</code> term like this (but this is too complex for my taste):</p>
  <pre><code>for (e=0,a=N/i,b=(N&gt;&gt;1)/i;(a)||(b);e+=a-b-b,a/=i,b/=i);
  </code></pre>
  <p>... I will try implement <code>T2</code> into <code>fact(x)</code> and compare speeds ...</p>