階乗の計算(python)

階乗の高速計算ができそうなコードを作ってみました。

def dif(a,b,c,d):
    ac=c-a
    bd=d-b
    return a*bd+d*ac
lis=[1,2,3,4]
lis=[5,6,7,8]
lis=list(range(1,17))
lis2=[]
for i in range(2):#ループの回数。4**i
    lens2=[]
    lens=len(lis)//4
    for j in range(lens):
        dif1=dif(*lis[j*4:j*4+4])
        dif2=lis[j*4]*lis[j*4+1]
        lis2+=[dif2*(dif1+dif2)]
        print(dif1,dif2,dif2*(dif1+dif2))
    lis=lis2

素因数分解の問題を解くために、いろいろな方法を試してみましたが、合同式の階乗が高速に求められれば、変数が揃ったところで0になるからそれでいいや…と思い、もうそろそろ終了にしようかなと。

上のコードはループ1回につき4つの数字の積をまとめて計算します。それぞれの数字が等しい差により連続な数字のリストであるけれど、特にこの条件が必要な訳ではない。
記載のコードでは、1から16の階乗を計算するようになっています。

証明)
a,b,c,dを左が一番小さく、等しい差で順に増加する4つの連続する数字であるとする。
c-a=b-dより、
ab-cd=ab-ad+ad-cd=a(b-d)+d(a-c)である
aからdの積abcdについては、
abcd=ab*cd=ab*(cd-ab+ab)=ab*(ab+a(b-d)+d(a-c))

アウトプットが、

10 2 24
26 30 1680
42 90 11880
58 182 43680
518878080 40320 20922789888000
であり、関数の出力のdif1に規則性が見られるため、もう少し改善が望めそうである。ループ1回につき、階乗のnが4倍ということに一応なります。

速いのかはわからない。
階乗の数字がループ毎に4倍になるので、モンゴメリ法みたいに倍々に増やして、途中で0になったらそこを探索していく方法でいいかな。

合同式の計算でかつn=64^3(=262144)で計算してみると、ざっと0.1秒、まあそんなに早いわけではなさそう。素の階乗計算であれば約2秒でした。

扱いやすい二項係数へと変形して、応用を効かせたら以下の式が成り立つかもしれないので、少し調べてみようかなと。
$${{}_{4m} C_{2n}= {}_{2m} C_{2n}+(n!)^2/(2n!)({}_{4m} C_{n}({}_{3m} C_{n}-{}_{m} C_{n})+ {}_{m} C_{n}({}_{4m} C_{n}-{}_{2m} C_{n}))}$$