計算階乘的另一些有趣的演算法

  一個正整數n的階乘就是前n個正整數的乘積，我們通常需要n-1次乘法操作來算出精確的值。不像等差數列求和、a的n次冪之類的東西，目前求階乘還沒有什麼巨牛無比的高效演算法，我們所能做的僅僅是做一些小的最佳化。

更少的乘法運算次數？
    在高精度運算中，乘法計算的速度遠遠慢於加減法，因此我們有必要減少乘法運算的次數。下面我將做一個非常簡單的變換，使得計算階乘只需要n/2次乘法。繼續看下去之前，你能自己想到這個演算法來嗎？

    我們可以把一個數的階乘轉換為若干個平方差的積。例如，假如我想求9!，我可以把前9個正整數的乘積寫成這個樣子：
   1 * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9
= (5-4) * (5-3) * (5-2) * (5-1) * 5 * (5+1) * (5+2) * (5+3) * (5+4)
= (5-1) * (5+1) * (5-2) * (5+2) * (5-3) * (5+3) * (5-4) * (5+4) * 5
= (5^2 - 1^2) * (5^2 - 2^2) * (5^2 - 3^2) * (5^2 - 4^2) * 5
    注意到一個有趣的事實：上面的四個平方差算出來分別是24, 21, 16, 9，它們之間的差正好是連續的奇數（因為n^2等於前n個正奇數的和）。因此，我們可以用初始數(n/2)^2不斷減去一個個的正奇數，求出所有n/2個平方差，再用n/2次乘法把它們乘起來。這種演算法實現起來非常簡單，並且（當n不大時）同樣只需要單精確度乘高精度，但需要的乘法次數大大減少了。假設我們已經有了一個高精度類，求n!只需要下面幾句話：
long h=n/2, q=h*h;
long r = (n&1)==1 ? 2*q*n : 2*q;
f = LargeInteger.create(r);
for(int d=1; d<n-2; d+=2)
   f = f.multiply(q-=d);

更少的總運算次數？
    盡量提取階乘中的因子2，我們可以得到另一種階乘運算的最佳化方法。這很可能是不需要分解質因數的階乘演算法中最快的一種。
    假如我們需要計算20!，我們可以把20拆成若干組正奇數的乘積：

  1 * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19 * 20
= 1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19 * 2 * 4 * 6 * 8 * 10 * 12 * 14 * 16 * 18 * 20
= 1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19 * 1 * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 2^10
= 1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19 * 1 * 3 * 5 * 7 * 9 * 2 * 4 * 6 * 8 * 10 * 2^10
= 1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19 * 1 * 3 * 5 * 7 * 9 * 1 * 2 * 3 * 4 * 5 * 2^15
= 1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19 * 1 * 3 * 5 * 7 * 9 * 1 * 3 * 5 * 2 * 4 * 2^15
= 1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19 * 1 * 3 * 5 * 7 * 9 * 1 * 3 * 5 * 1 * 2 * 2^17
= 1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19 * 1 * 3 * 5 * 7 * 9 * 1 * 3 * 5 * 1 * 2^18

    只需要一次累乘就可以求到每一組奇數的乘積，最後再花費log(n)次乘法把它們全部乘起來。最後的那個2^18也可以二分計算出來。真正的代碼還有很多細節上的最佳化，另外還借用了遞迴使得操作變得更加簡便。你可以在本文最後附的那個連結裡去找Split-Recursive演算法。

還能再快一點嗎？
    繼續擴充上面的演算法，我們可以想到，如果把每個數的質因數都分解出來，並且統計每種質因子有多少個，我們就可以多次使用二分求冪，再把它們的結果乘起來。注意這裡並不是真的要老老實實地去分解每個數的質因子。對於每個質數x，我們可以很快算出前n個正整數一共包含有多少個質因子x（記得如何求n!末尾有多少個0麼）。這種演算法的效率相當高，已經能夠滿足大多數人的需要了。

另一種詭異的階乘演算法：
    這個演算法可能是所有有名字的階乘演算法中最慢的一個了（Additive Moessner演算法），它對一個數列進行重複的累加操作，一次次地計算首碼和，總共將花費O(n^3)次加法操作。但是，令人費解的是，這個簡單的程式為什麼可以輸出前n個正整數的階乘呢？
a[0]:=1;
for i:=1 to n do
begin
   a[i]:=0;
   for j:=n downto 1 do
   begin
      for k:=1 to j do
         a[k]:=a[k]+a[k-1]
      write(a[i],' ');
   end;
end;
    我在網上搜尋相關的東西時找到了另一個有趣的東西。對一個初始時全為1的數列反覆進行這兩個操作：累加求首碼和，然後以1,2,3,...的間隔劃掉其中一部分數（即划去所有位置編號為三角形數的數）形成新的序列。類似的數列操作方法最先由Alfred Moessner提出的，我們這裡不妨把它叫做Moessner數列。你會發現，第n輪操作開始前，數列的第一個數恰好是n! 。看看下面的例子吧：

1 1 1 1 1 1 1 1 1  1  1  1  1  1  1 ...
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ...
x 2 x 4 5 x 7 8 9  x 11 12 13 14  x ...

2 4  5  7  8  9 11 12 13 14 ...
2 6 11 18 26 35 46 58 71 85 ...
x 6  x 18 26  x 46 58 71  x ...

6 18 26 46  58  71 ...
6 24 50 96 154 225 ...
x 24  x 96 154   x ...

24  96 154 ...
24 120 274 ...
x 120  x  ...

120 ...
.....

    當然，發現前面O(n^3)的程式和這個Moessner數列的關聯時我很是吃了一驚：在前面的程式裡，如果你輸出每一次i迴圈末所得到的數列，你會發現輸出的這些數正好就是後面這個問題裡被我們劃掉的數，而它們其實就是第一類Stirling數！
    這到底是為什麼呢？是什麼東西把階乘、第一類Stirling數、Moessner數列和那個O(n^3)的程式聯絡在一起的呢？昨天，我想這個問題想了一天，最後終於想通了。如果把Moessner數列排列成這個樣子，一切就恍然大悟了：

  
    仔細觀察，我們會發現：
    1. 按照Moessner數列的定義，每個數都應該等於它左邊的數和左上方的數的和（這個“左邊”可以跳過若干空格）。例如，35 = 9 + 26，46 = 11 + 35。排成一系列三角形後，每個三角形最右邊一列的數就是被划去的數，它永遠不能參與它下面的那些行的運算。
    2. 設a[n,i,j]表示左起第n個三角形陣列中的第i行右起第j列上的數，則a[n,i,j]=a[n-1,i-1,j]*n + a[n-1,i,j]，例如274=50*5+24。如果遞推時遇到空白位置而它左邊隔若干空格的地方還有數的話，則需要用左邊的數來補，例如18=4*4+2。對於每個三角形的最後一列來說，這個性質實際上就是第一類Stirling數的遞推關係，因此Moessner數列中才會出現第一類Stirling數。
    3. 在第一類Stirling數中，s(n,1)=n! ，也即左起第n個三角形最底端的那個數等於n!。從上面的第二個性質來看，這也是顯然的。
    4. O(n^3)的演算法實際上就是在繪製上面這個圖。每一次j迴圈末，我們得到的序列是第i個三角形中每一行左起第j個數組成的序列。例如，計算第5個三角形內的數時，程式首先累加出1, 11, 46, 96, 120, 120，這樣便算出了a[5]=120，數列的前5個數再次累加即得到1, 12, 58, 154, 274，由此算出a[4]=274。
    第二個性質可以利用第一個性質進行數學歸納法證明，證明很簡單，我就不多說了。現在我儘可能少寫一些繁瑣的細節，節約一些時間用來複習古代漢語。

做人要厚道，
轉貼請註明出處。

查看更多：
http://www.luschny.de/math/factorial/FastFactorialFunctions.htm
http://www.luschny.de/math/factorial/index.html <---- 巨牛，20多種階乘演算法的代碼！