archive/share/trash/determ.mac

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