计算行列式的无限循环

Question

我担心我违反了永远不要尽可能聪明地编写代码的规则。尝试使用 LaPlace 展开计算 8x8 矩阵的行列式时，我发生了无限递归。

程序尝试执行圆的投影变换，该投影由 4 个源点和 4 个目标点定义。

基本思想来自this page，它显示了一个用于求解投影矩阵系数的 8x8 矩阵。然后我尝试从维基百科的定义中实现所需的功能。

xform.ps（将G library 用于函数matmul transpose：

%!
errordict/rangecheck{pstack()= countexecstack array execstack == quit}put
errordict/dictstackoverflow{countexecstack array execstack dup length 21 sub 20 getinterval == quit}put
4(G/G)run <<

/main {
    300 400  100  0 360 arc  flattenpath

    100 100  500 100  500 500  100 500
    200 200  300 200  300 300  200 300  xformpath

    stroke
}

/xformpath {
    10 dict begin
        {mul neg}
        { *- v3 u3  v2 u2  v1 u1  v0 u0
             y3 x3  y2 x2  y1 x1  y0 x0 }{exch def}forall

        [
            [ x0  y0  1   0   0   0   x0 u0 *-  y0 u0 *- ]
            [ x1  y1  1   0   0   0   x1 u1 *-  y1 u1 *- ]
            [ x2  y2  1   0   0   0   x2 u2 *-  y2 u2 *- ]
            [ x3  y3  1   0   0   0   x3 u3 *-  y3 u3 *- ]
            [ 0   0   0   x0  y0  1   x0 v0 *-  y0 v0 *- ]
            [ 0   0   0   x1  y1  1   x1 v1 *-  y1 v1 *- ]
            [ 0   0   0   x2  y2  1   x2 v2 *-  y2 v2 *- ]
            [ 0   0   0   x3  y3  1   x3 v3 *-  y3 v3 *- ]
        ]
        invertmat
        [ u0 u1 u2 u3 v0 v1 v2 v3 ] transpose
        matmul
        massage-vector-into-matrix
        /P exch def

        [
            { project2 /moveto cvx }
            { project2 /lineto cvx }
            {}
            { /closepath cvx } pathforall
        ] cvx
        newpath exec

    end
}

/massage-vector-into-matrix {
    transpose 0 get
    dup 0 3 getinterval exch
    dup 3 3 getinterval exch
    6 2 getinterval aload pop 1  3 array astore
    3 array astore
}

/project2 {
    hom P matmul het
}

/hom { % x y  ->  [ x y 1 ]
    1  3 array astore
}

/het { % [ x y w ]  ->  x/w y/w
    aload pop dup 3 1 roll div 3 1 roll div exch
}

/invertmat {
    dup  det 1 exch div  dup dup 3 array astore exch
    adjugate  { { mul } sop } vop
}

/det {
    << exch
        /A exch
    >> begin
        A length 2 eq {
            aload pop
            aload pop /d exch def /c exch def
            aload pop /b exch def /a exch def
            a d mul b c mul sub
        }{
            /M A length def  % LaPlace expansion on first column
            /N A 0 get length def
            0
            0 1 M _1 { /i exch def
                i 0 A minor
                i 2 mod 1 eq { neg } if
                add
            } for
        } ifelse
    end
}

/adjugate {
    cofactor transpose
}

/cofactor {
    << exch
        /A exch
        /M 1 index       length
        /N 3 index 0 get length
    >> begin
        [
            0 1 M _1 { /i exch def
                [
                    0 1 N _1 { /j exch def
                        i j A minor
                        i j add 2 mod 1 eq { neg } if       
                    } for
                ]
            } for
        ]
    end
}

/minor {
    3 dict begin
        dup length exch
        dup 0 get length exch
        {A N M n m}{exch def}forall
        [
            0 1 M _1 { /i exch def
                [
                    0 1 N _1 { /j exch def
                        m i eq  n j eq  or  not {
                            A i get j get
                        } if
                    } for
                ]
            } for
        ]

        det
    end
}

/_1 { 1 sub }

>> begin main

错误是/dictstackoverflow。转储 execstack 的尾部给出了这个：

[--%for_pos_int_continue-- {i 2 mod 1 eq {neg} if add} {end} {end} 1 1 7 
{/i exch def i 0 A minor i 2 mod 1 eq {neg} if add} --%for_pos_int_continue-- 
{i 2 mod 1 eq {neg} if add} {end} {end} 1 1 7 {/i exch def i 0 A minor i 
2 mod 1 eq {neg} if add} --%for_pos_int_continue-- {i 2 mod 1 eq {neg} if 
add} {end} {A length 2 eq {aload pop aload pop /d exch def /c exch def aload 
pop /b exch def /a exch def a d mul b c mul sub} {/M A length def /N A 0 ge
t length def 0 0 1 M _1 {/i exch def i 0 A minor i 2 mod 1 eq {neg} if add} 
for} ifelse end}]

谁能发现我做错了什么，或者有没有更简单的方法来计算可以更干净地避免我的问题的未成年人和行列式？

Answer 1

我得到了一个决定因素似乎起作用的函数。当然，它使用 LaPlace 扩展。但是通过扩展第一列，可以使用getinterval 提取所有部分行。所以它沿着第一列组装未成年人，没有做双循环，看起来更容易阅读。

%!
%errordict/rangecheck{pstack()= countexecstack array execstack == quit}put
4(G/G)run
/block { << exch {} forall >> begin main } def
{ %block

main {
    [[3 8][4 6]] determinant ==
    [[ 12 2 3 ][4 5 6][7 8 9]] determinant ==
    [[ 1 0 0 ][0 1 0][0 0 1]] determinant ==
    [[6 1 1][4 -2 5][2 8 7]] determinant ==
    [[1 2 1 0][0 3 1 1][-1 0 3 1][3 1 2 0]] determinant ==
    quit
}

determinant {
    dup length 3 le {
        dup length 0 eq {
            pop 0
        }{
        dup length 1 eq {
            0 get
        }{
        dup length 2 eq {
            aload pop aload pop 3 2 roll aload pop % c d a b
            4 -1 roll mul neg % d a -bc
            3 1 roll mul add % ad-bc
        }{
            dup 1 get 1 2 getinterval 1 index 2 get 1 2 getinterval 2 array astore determinant
            1 index 0 get 0 get mul exch
            dup 0 get 1 2 getinterval 1 index 2 get 1 2 getinterval 2 array astore determinant
            1 index 1 get 0 get mul neg exch
            dup 0 get 1 2 getinterval 1 index 1 get 1 2 getinterval 2 array astore determinant
            exch 2 get 0 get mul  add add
        }ifelse }ifelse }ifelse
    }{
        %dup ==()=
        dup { 0 get } map exch
        { 1  1 index length 1 sub  getinterval } map  % [ ai0 ] [ ai1..n ]
        [ exch
          0  1  2 index length 1 sub { % [ ... ai1..n  i
              2 copy  2 copy  0 exch getinterval 3 1 roll  % [ ... a i a0..11..n ai1..n i
              1 add  1 index length 1 index sub  getinterval
              combine 3 1 roll pop
          } for pop ] %   pstack()=
        { determinant } map % dup ==()=
        { mul } vop
        aload length 1 sub { sub } repeat  %  -+-+ == -(-(-(-)))
    } ifelse
}

combine {
    %    exch [ exch {}forall counttomark 2 add -1 roll {}forall ]
    2 copy length exch length add array
    dup 0 4 index putinterval
    dup 4 -1 roll length 4 -1 roll putinterval
}

map {
    [ 3 1 roll forall ]
}

} block