Боюсь, я нарушил правило никогда не писать код настолько умно, насколько это возможно. У меня происходит бесконечная рекурсия при попытке вычислить определитель матрицы 8x8 с использованием расширения Лапласа.
Программа пытается выполнить проективное преобразование окружности, проекцию которой определяют 4 исходные точки и 4 точки назначения.
Основная идея взята из этой страницы, на которой показана используемая матрица 8x8. решить для коэффициентов матрицы проекции. Затем я попытался реализовать необходимые функции из определений в википедии.
xform.ps (использует библиотеку G для функций 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}]
Может ли кто-нибудь определить, что я делаю неправильно, или есть более простой способ вычисления миноров и детерминантов, который позволил бы более чисто избежать моей проблемы?
