Edison.X. Blog

K algorithm ( "Super-random" number generator)

這是由 Kunth 於 The Art Of Programming 裡提出說明，當時 Kunth 試著以此演算法建立較好之 PRNG，也強調讀者不必特別研究它，故這部份沒太大興趣者可跳過，流程大致如下，其中 ** 代表次方，如 3**5 代表 3 的 5 次方。(以下特別注意 K1, K2)

INIT: 給一 10 位數 X ，進行以下步驟
K1. 選擇迭代次數： Y:= X/(10 ** 9)，執行 K2~K13 ，Y+1次。
K2. 選擇隨機步驟： Z:= Y/(10 ** 8) mod 10，步驟轉向  K(3+Z)
K3. 確保 X >= 5* 10**9： 若 X < 5000000000 ， X := X + 5000000000
K4. 平方取中法： X := (X*X / 10 ** 5) mod 10
K5. 進行乘法： X:= 1001001001X
K6. 假(偽)補數： 若 X < 10 ** 8 ， X:= X+9814055677；否則 X:=10 ** 10 - X
K7. 互換二半： X 高5位與 X 低5位 (十進制) 交換，即 X:= (X mod 10 ** 5) + (X / 10 **5)
K8. 進行乘法： 同 K5
K9. 減小數字：將 X 之十進制表示之每個非 0 數字減 1
K10. 9999修改：若 X<100000，X:=X*X + 9999；否則 X:=X-9999
K11. 正規化：此時 X 不能為零，若 X < 10 ** 9 ，X:=10X，重覆此步驟
K12. 修正之平方取中法： X:= ( X*(X-1) / 10 ** 5  ) mod ( 10 ** 10)，即中間十位數取代 X
K13. 重複 ? ：若 Y>0，Y:=Y-1，回到步驟 K2；若 Y=0，其 X 即為隨機值。

最終，此算法收斂至 605038420，而以另一數開始時 (原文未提及為何數)，於 7401 個數後，以長度為 3178 之週期做為循環。結論是：隨機亂數並非透過隨機選擇之方法便可生成，應 應用某些理論方可完善實現。

由於 Kunth 所提出之 K演算法，必須有 20 位數(十進制) 之整數型態方可快速進行，而 unsigned long long 最大值也才 19 位數，雖可自行調整試試，但效果未必像 Kunth 所做如此。此處不再加以 implement，唯留下二問題供參考。

K9  步驟，該如何實現？
K11步驟，語意敘述上是一個 while loop，是否有更快之解法？

其實 Kunth 在這篇演算法裡面，有提到不少 div, mod 技巧，如互換高底 5 個數 (十進制)、平方取中法之小技巧，唯此處輕輕帶過，推測 Kunth 可能有其他更快之解法，或認為 "這是常識"，只是我沒想到。

以 unsigned long long，「不考慮 overflow 情況下」 (雖然已知這是不可能的)，基礎板的碼參考如下

unsigned long long K9(unsigned long long X)
{
    unsigned long long multi=1;
    while(multi<=X){                    // K9
       if( (X % multi*10 ) / multi !=0) X-=multi;
       multi*=10;
    }
    return X;
}

unsigned long long K11(unsigned long long X)
{
    while(X<1000000000) X*=10;
    return X;
}

如何更快？此處便留下另一伏筆。