行列式の値を計算するプログラムをPythonで作ってみた。(3)

　良く考えてみたら、昨日のプログラムで直接計算させてみたらいいかもということでやってみました。(^_^;
　この展開方法は帰納的に作ったので何展開か分かりません。公式集でも見ればたぶん何かあるでしょう。(^_^;
　それにしても、最初の頃に比べれば随分速くなりました。PythonはCに比べれば遅いですが、少し改良すればすぐに結果に反映するのでやり甲斐があります。(^_^;

● CalcDetN1.py

# coding: UTF-8
# CalcDetN1.py

import itertools

def Sgn(p):
    r,n = 1,len(p)
    v,w = [0]+list(p),list(range(n+1))
    for i in range(1,n+1):
        v[i]+=1; w[v[i]] = i
    for i in range(1,n):
        j = v[i]
        if j!=i:
            v[w[i]] = j; w[j] = w[i]; r*=-1
    return r

def det3(m):
    r = m[0][0]*(m[1][1]*m[2][2]-m[1][2]*m[2][1])
    r+= m[0][1]*(m[1][2]*m[2][0]-m[1][0]*m[2][2])
    r+= m[0][2]*(m[1][0]*m[2][1]-m[1][1]*m[2][0])
    return r

def det4(m):
    r = (m[0][0]*m[1][1]-m[0][1]*m[1][0])*(m[2][2]*m[3][3]-m[2][3]*m[3][2])
    r-= (m[0][0]*m[1][2]-m[0][2]*m[1][0])*(m[2][1]*m[3][3]-m[2][3]*m[3][1])
    r+= (m[0][0]*m[1][3]-m[0][3]*m[1][0])*(m[2][1]*m[3][2]-m[2][2]*m[3][1])
    r+= (m[0][1]*m[1][2]-m[0][2]*m[1][1])*(m[2][0]*m[3][3]-m[2][3]*m[3][0])
    r-= (m[0][1]*m[1][3]-m[0][3]*m[1][1])*(m[2][0]*m[3][2]-m[2][2]*m[3][0])
    r+= (m[0][2]*m[1][3]-m[0][3]*m[1][2])*(m[2][0]*m[3][1]-m[2][1]*m[3][0])
    return r

def detN(mat):
    n = len(mat)
    if n==1: return mat[0][0]
    if n==2: return mat[0][0]*mat[1][1]-mat[1][0]*mat[0][1]
    if n==3: return det3(mat)
    if n==4: return det4(mat)

    a = [[mat[i][j] for j in range(n//2)  ] for i in range(n)]
    b = [[mat[i][j] for j in range(n//2,n)] for i in range(n)]
    res = 0
    P = range(n)
    for p in itertools.combinations(P,n//2):
        q = list(p)+[k for k in P if k not in p ]
        ma = [a[q[i]] for i in range(n//2)]
        mb = [b[q[i]] for i in range(n//2,n)]
        res+= Sgn(q)*detN(ma)*detN(mb)
    return res

def main():
    import numpy as np
    import random
    from time import time
    tm=time()  # Timer Start

    for n in range(1,10+1):
        mat=[[random.randint(0,99) for x in range(n)] for y in range(n)]
        a = np.array(mat)
        print("%2d : %s | %s"%(n,detN(mat),np.linalg.det(a)))

    tm=time()-tm  # Timer Stop
    print("Runtime : %.3f [sec]"%tm)

if __name__ == '__main__':
    main()

●実行結果

 1 : 86 | 86.0
 2 : 3515 | 3515.0
 3 : 389892 | 389892.0
 4 : -3700583 | -3700583.0
 5 : 58778732 | 58778732.0
 6 : -26582694044 | -26582694044.0
 7 : 3236361776857 | 3.23636177686e+12
 8 : 3687627036878 | 3.68762703688e+12
 9 : -87510485617826102 | -8.75104856178e+16
10 : -998051754221172129 | -9.98051754221e+17
Runtime : 0.153 [sec]

※参考URL
●行列式の値を計算するプログラムをPythonで作ってみた。
●行列式の値を計算するプログラムをPythonで作ってみた。(2)
●行列式の展開式を文字列で求めるプログラムをPythonで作ってみた。

線型代数学 (数学選書 (1))