1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
|
/*
?putprop(?quote(gcdivide),?filestrip(?cdr([functs,lisp,dsk,share])),
?quote(?autoload)); */
eval_when([translate,batch,demo,load,loadfile],
load("functs"))$
det&& det(m):=block([n,a],/* local(a), */
n:length(m),
if n < 2 then error ("improper argument:",m),
/* array(a,n,n), */
a:?gensym(),
apply('array,[a,n,n]),
for i thru n do for j thru n do arrayapply(a,[i,j])::m[i,j],
arrayapply(a,[0,0])::1,
catch(for k from 2 step 2 thru n do /* iterate for each pair of rows */
(if k#n then block([c0,l1,l2,u], /* omit calculation on last
iteration for n even */
l1:(if a[k-1,k-1]#0 then /* check for possible pivoting */
false
else /* zero element means pivoting is necessary */
for s from k thru n do /* search for nonzero element */
if a[s,k-1]#0 then /* found pivot */
(for j thru n do /* interchange rows s and k-1 */
block([t],
t:a[k-1,j],
a[k-1,j]:a[s,j],
a[s,j]:t),
return(false)) /* indicate pivot was found */
else /* no pivot found means */
if s=n then throw(0)), /* singular matrix */
l2:(if l1 then
true
else
/* search through rows for nonzero multiplier */
for t from k thru n do
(c0:determinant(matrix(
[a[k-1,k-1],a[k-1,k]],
[a[t,k-1],a[t,k]])),
if c0#0 then
(u:t,
return(false)) /* indicate multiplier was found */
else /* no multiplier means singular matrix */
if t=n then throw(0))),
if l2 then
return(true)
else /* calculate coefficient from multiplier */
(c0:gcdivide(c0,a[k-2,k-2]),
if u#k then /* if multiplier was found in a different row */
for j thru n do block([t], /* interchange rows t and k */
t:a[k,j],
a[k,j]:a[t,j],
a[t,j]:t),
/* iterate through remaining rows */
for i from k+1 thru n do block([c1,c2],
c1:gcdivide(-determinant(matrix(
[a[k-1,k-1],a[k-1,k]],
[a[i,k-1],a[i,k]])),a[k-2,k-2]),
c2:gcdivide(determinant(matrix(
[a[k,k-1],a[k,k]],
[a[i,k-1],a[i,k]])),a[k-2,k-2]),
a[i,k-1]:0,
a[i,k]:0,
/* iterate through remaining columns */
for j from k+1 thru n do
a[i,j]:gcdivide(c0*a[i,j]+c1*a[k,j]+c2*a[k-1,j],
a[k-2,k-2])),
c0:0,
for j from k+1 thru n do
(a[k,j]:0,
a[k-1,j]:0))),
/* omit calculation on last iteration for n odd */
a[k,k]:if k=n-1 then 0 else
gcdivide(determinant(genmatrix(a,k,k,k-1)),a[k-2,k-2]),
a[k-2,k-2]:0,
a[k-1,k-1]:0,
a[k,k-1]:0,
a[k-1,k]:0,
/* on last iteration indicate nonsingular matrix */
if k=n or k=n-1 then return(false)),
a[n,n])) /* if a singular matrix has not been thrown,
then catch the determinant */$
"end of determ.mc"$
|
