File: determ.mc

package info (click to toggle)
maxima 5.6-17
  • links: PTS
  • area: main
  • in suites: woody
  • size: 30,572 kB
  • ctags: 47,715
  • sloc: ansic: 154,079; lisp: 147,553; asm: 45,843; tcl: 16,744; sh: 11,057; makefile: 7,198; perl: 1,842; sed: 334; fortran: 24; awk: 5
file content (82 lines) | stat: -rw-r--r-- 2,660 bytes parent folder | download | duplicates (2)
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"$