VAX780 simulator V3.8-1
Listening on port 2238 (socket 4)
Modem control activated

Boot
: hp(0,0)vmunix
123060+27528+24628 start 0xF5C
Berkeley VAX/UNIX Version 4.9  Wed Feb 17 15:27:46 PST 1982
real mem  = 8322048
avail mem = 7738368
mcr0 at tr1
mcr1 at tr2
uba0 at tr3
dz0 at uba0 csr 160100 vec 300, ipl 15
mba0 at tr8
hp0 at mba0 drive 0
hp1 at mba0 drive 1
hp2 at mba0 drive 2
hp3 at mba0 drive 3
mba1 at tr9
ht0 at mba1 drive 0
tu0 at ht0 slave 0
tu1 at ht0 slave 1
root on hp0
WARNING: should run interleaved swap with >= 2Mb
Automatic reboot in progress...
Mon Oct 11 18:27:25 GMT 1976
/dev/hp0a: 678 files 4278 blocks 3345 free
/dev/rhp0g: 2590 files 47295 blocks 94283 free
Mon Oct 11 18:27:31 GMT 1976
Mounted /usr on /dev/hp0g
preserving editor files
clearing /tmp
starting daemons: update cron accounting network mail printer.
Mon Oct 11 18:27:32 GMT 1976


Berkeley 4.1 VAX/UNIX (Amnesia-Vax)

login:


--------------------------------------------------------------------------------------------

[~/]kenrou% telnet localhost 2238
Trying 127.0.0.1...
Connected to localhost.
Escape character is '^]'.


Connected to the VAX780 simulator DZ device, line 0



Berkeley 4.1 VAX/UNIX (Amnesia-Vax)

login: kenrou
Password:
Last login: Mon Oct 11 16:11:55 on tty00

Welcome to Berkeley Vax/UNIX (4.1bsd revised 1 Sept. 1981)
4.1BSD kenrou% cd Programs/vaxima
4.1BSD kenrou% ls
kihon.v     primer.txt  tensor.v
4.1BSD kenrou% vaxima

     Vaxima 2.04
  Sun Sep 12 15:31:19 1982

(c1) load("eigen.v");
Batching the file /usr/mac/share/eigen.v

(c2) /* THIS IS THE FILE EIGEN > DSK:SHARE;.  April 1982.  */

EVAL_WHEN(TRANSLATE_FILE,
        MODEDECLARE([HERMITIANMATRIX,NONDIAGONALIZABLE,KNOWNEIGVALS,
        KNOWNEIGVECTS],BOOLEAN,
        [INDEX1,INDEX2,INDEX3,INDEX4,DIMNSN,COUNT,%RNUM],FIXNUM),
        DECLARE([HERMITIANMATRIX,NONDIAGONALIZABLE,KNOWNEIGVALS,
        KNOWNEIGVECTS,LISTEIGVECTS,LISTEIGVALS,%RNUM,LISTARITH,PROGRAMMODE,
        ALGEBRAIC,REALONLY,MULTIPLICITIES,SOLVEEXPLICIT],SPECIAL))$

/usr/mac/transl/trmode.o being loaded.
[fasl /usr/mac/transl/trmode.o]

(c3) SSTATUS(FEATURE,EIGEN)$

(c4) HERMITIANMATRIX:FALSE$

(c5) NONDIAGONALIZABLE:FALSE$

(c6) KNOWNEIGVALS:FALSE$

(c7) KNOWNEIGVECTS:FALSE$

(c8) LISTEIGVECTS:[]$

(c9) LISTEIGVALS:[]$

(c10) CONJUGATE(X):=SUBLIS('[%I=-%I],X)$

(c11) INNERPRODUCT(X,Y):=BLOCK([LISTARITH],LISTARITH:TRUE,RATSIMP(CONJUGATE(X).Y))$

(c12) UNITVECTOR(X):=BLOCK([LISTARITH,INTRN],LISTARITH:TRUE,INTRN:INNERPRODUCT(X,X),
INTRN:SQRT(INTRN),X/INTRN)$

(c13) COLUMNVECTOR(X):=TRANSPOSE(MATRIX(X))$

(c14) GRAMSCHMIDT(X):=
                BLOCK([LISTARITH,DIMNSN,LISTALL,INTERN,COUNT,DENOM,UNIT,INDEX1,
                INDEX2],
                LISTARITH:TRUE,DIMNSN:LENGTH(X),LISTALL:[PART(X,1)],
                COUNT:1,IF DIMNSN=1 THEN RETURN(X)
                ELSE (FOR INDEX1:2 THRU DIMNSN DO
                (UNIT:PART(X,INDEX1),FOR INDEX2 THRU COUNT DO
                (INTERN:PART(LISTALL,INDEX2),DENOM:INNERPRODUCT(INTERN,INTERN),
                UNIT:FACTOR(RATSIMP(UNIT-INNERPRODUCT(INTERN,UNIT)*INTERN/DENOM
                ))),
                COUNT:COUNT+1,LISTALL:ENDCONS(UNIT,LISTALL)),
                RETURN(LISTALL)))$

(c15) EIGENVALUES(MAT):=
                BLOCK([DIMNSN,LISTALL,SOLUTION,MULTIPLICITIES,SOLVEEXPLICIT,
                DUMMY:?GENSYM(),INDEX2],
                LISTALL:[],
                SOLVEEXPLICIT:TRUE,
                DIMNSN:LENGTH(MAT),
                SOLUTION:BLOCK([PROGRAMMODE:TRUE],
                SOLVE(CHARPOLY(MAT,DUMMY),DUMMY)),
                IF SOLUTION=[] THEN 
                (PRINT(" "),PRINT("SOLVE is unable to find the roots of"),
                PRINT("the characteristic polynomial."),
                RETURN(LISTALL))
                ELSE (FOR INDEX2 THRU DIMNSN DO
                (DIMNSN:DIMNSN-PART(MULTIPLICITIES,INDEX2)+1,
                LISTALL:ENDCONS(RHS(PART(SOLUTION,INDEX2)),LISTALL)),
                LISTALL:ENDCONS(MULTIPLICITIES,[LISTALL]),
                RETURN(LISTALL)))$

(c16) EIGENVECTORS(MAT):=
                BLOCK([EQUATIONS,UNKNOWNS,SOLUTION,LISTALL,EIGVALS,DIMNSN,
                COUNT,VECTR,INDEX3,INDEX4,INDEX2,INDEX1,NOTKNWN,MATRX,MMATRX,
                UNIT,MULTIPLICITIES,%RNUM,REALONLY,ALGEBRAIC,INTERM,INTERN],
                LOCAL(NOTKNWN),UNKNOWNS:[],DIMNSN:LENGTH(MAT),COUNT:DIMNSN,
                IF KNOWNEIGVALS THEN EIGVALS:LISTEIGVALS
                ELSE EIGVALS:EIGENVALUES(MAT),
                IF EIGVALS=[] THEN (NONDIAGONALIZABLE:TRUE,RETURN(EIGVALS))
                ELSE (MULTIPLICITIES:PART(EIGVALS,2),
                FOR INDEX1 THRU DIMNSN DO
                UNKNOWNS:ENDCONS(NOTKNWN[INDEX1],UNKNOWNS),
                VECTR:COLUMNVECTOR(UNKNOWNS),MATRX:MAT.VECTR,
                NONDIAGONALIZABLE:FALSE,
                LISTALL:[EIGVALS],REALONLY:FALSE,ALGEBRAIC:TRUE,
                FOR INDEX1 THRU COUNT DO 
                (COUNT:COUNT-PART(MULTIPLICITIES,INDEX1)+1,
                MMATRX:MATRX-PART(EIGVALS,1,INDEX1)*VECTR,
                EQUATIONS:[],
                FOR INDEX2 THRU DIMNSN DO
                EQUATIONS:CONS(MMATRX[INDEX2,1],EQUATIONS),%RNUM:0,
                SOLUTION:ALGSYS(EQUATIONS,UNKNOWNS),
                INTERM:MAP(RHS,SOLUTION[1]),
                UNIT:[],IF %RNUM#PART(MULTIPLICITIES,INDEX1)
                THEN NONDIAGONALIZABLE:TRUE,
                FOR INDEX3 THRU %RNUM DO
                (INTERN:SUBSTVECTK(%RNUM_LIST,INDEX3,INTERM),
                UNIT:APPEND(UNIT,[INTERN])),
                IF UNIT=[] THEN
                (PRINT(" "),PRINT("ALGSYS failure: The eigenvector(s) for the",
                INDEX1,"th eigenvalue will be missing.")),
                IF HERMITIANMATRIX AND %RNUM>1 THEN UNIT:GRAMSCHMIDT(UNIT),
                LISTALL:APPEND(LISTALL,UNIT)),
                RETURN(LISTALL)))$

(c17) /* The first arg is of the form [r1,r2,r3].
   We want to construct [r1=0,r2=1,r3=0] for example. */

SUBSTVECTK(L,N,EXP):=
                BLOCK([SUB_LIST:[],J:0],
                FOR VAR IN L DO
                (J:J+1,SUB_LIST:CONS(VAR = IF J=N THEN 1 ELSE 0,SUB_LIST)),
                SUBLIS(SUB_LIST,EXP))$

(c18) UNITEIGENVECTORS(MAT):=
                BLOCK([LISTUEVEC,LISTALL,INDEX1,UNIT],
                IF KNOWNEIGVECTS THEN LISTUEVEC:LISTEIGVECTS
                ELSE LISTUEVEC:EIGENVECTORS(MAT),
                IF LISTUEVEC=[] THEN RETURN(LISTUEVEC)
                ELSE (LISTALL:[PART(LISTUEVEC,1)],
                FOR INDEX1:2 THRU LENGTH(LISTUEVEC) DO 
                (UNIT:PART(LISTUEVEC,INDEX1),
                UNIT:RATSIMP(UNITVECTOR(UNIT)),
                LISTALL:ENDCONS(UNIT,LISTALL)),
                RETURN(LISTALL)))$

(c19) SIMILARITYTRANSFORM(MAT):=
                BLOCK([LISTVEC,LISTUEVEC],
                LISTUEVEC:UNITEIGENVECTORS(MAT),
                IF NONDIAGONALIZABLE THEN RETURN(LISTUEVEC)
                ELSE (LISTVEC:DELETE(PART(LISTUEVEC,1),LISTUEVEC),
                RIGHTMATRIX:TRANSPOSE(APPLY(MATRIX,LISTVEC)),
                IF HERMITIANMATRIX THEN
                LEFTMATRIX:CONJUGATE(TRANSPOSE(RIGHTMATRIX))
                ELSE LEFTMATRIX:RIGHTMATRIX^^-1,
                RETURN(LISTUEVEC)))$

(c20) CONJ(X):=CONJUGATE(X)$

(c21) INPROD(X,Y):=INNERPRODUCT(X,Y)$

(c22) UVECT(X):=UNITVECTOR(X)$

(c23) COVECT(X):=COLUMNVECTOR(X)$

(c24) GSCHMIT(X):=GRAMSCHMIDT(X)$

(c25) EIVALS(MAT):=EIGENVALUES(MAT)$

(c26) EIVECTS(MAT):=EIGENVECTORS(MAT)$

(c27) UEIVECTS(MAT):=UNITEIGENVECTORS(MAT)$

(c28) SIMTRAN(MAT):=SIMILARITYTRANSFORM(MAT)$

Batching done.(d29)                         /usr/mac/share/eigen.v

(c30) load("kihon.v");
Batching the file kihon.v

(c31) /* kihon.v */

colvector(A) := transpose(matrix(A))$

(c32) rowvector(A) := matrix(A)$

(c33) cp_col_ex(A, i, j) := block([B, dummy],
    mode_declare(B, any, dummy, any),
    B : copymatrix(A),
    dummy : B[i],
    B[i] : B[j],
    B[j] : dummy,
    B
)$

(c34) cp_col_ad(A, k, i, j) := block([B],
    mode_declare(B, any),
    B : copymatrix(A),
    B[i] : B[i] + k * B[j],
    B
)$

(c35) cp_col_sp(A, k, i) := block([B],
    mode_declare(B, any),
    B : copymatrix(A),
    B[i] : k * B[i],
    B
)$
[*symbol:342{90%}; list:931{47%}; fixnum:51{2%}; ut:0%]

(c36) cp_row_ex(A, i, j) := block([B],
    mode_declare(B, any),
    B : cp_col_ex(transpose(A), i, j),
    transpose(B)
)$

(c37) cp_row_ad(A, k, i, j) := block([B],
    mode_declare(B, any),
    B : cp_col_ad(transpose(A), k, i, j),
    transpose(B)
)$

(c38) cp_row_sp(A, k, i) := block([B],
    mode_declare(B, matrix),
    B : cp_col_sp(transpose(A), k, i),
    transpose(B)
)$

(c39) col_sp(A, k, i) := block(
    A[i] : k * A[i],
    A
)$

(c40) col_ad(A, k, i, j) := block(
    A[i] : A[i] + k * A[j],
    A
)$

(c41) col_ex(A, i, j) := block([dummy],
    dummy : A[i],
    A[i] : A[j],
    A[j] : dummy,
    A
)$

(c42) row_sp(A, k, i) := block([len],
    mode_declare(len, fixnum),
    len : length(transpose(A)[1]),
    for p : 1 thru len do
        A[p][i] : k * A[p][i],
    A
)$

(c43) row_ad(A, k, i, j) := block([len],
    mode_declare(len, fixnum),
    len : length(transpose(A)[1]),
    for p : 1 thru len do
        A[p][i] : A[p][i] + k * A[p][j],
    A
)$

(c44) row_ex(A, i, j) := block([len, c],
    mode_declare(len, fixnum, c, any),
    len : length(transpose(A)[1]),
    for p : 1 thru len do (
        c : A[p][i],
        A[p][i] : A[p][j],
        A[p][j] : c
    ),
    A
)$

(c45) kaidan_findnonzero(A, pi, pj, mm) := block([pix],
    mode_declare(pix, fixnum),
    pix : pi,
    jump,
    if not A[pix, pj] = 0 then
        pix
    else if pix = mm then
        (-1)
    else
    (
        pix : pix + 1,
        go(jump)
    )
)$

(c46) kaidan(A) := block([m, n, i, j, ir],
    mode_declare(m, fixnum, n, fixnum, i, fixnum, j, fixnum, ir, fixnum),

    if not matrixp(A) then
    (
        print(A, "is not a matrix"),
        []
    ),

    m : length(transpose(A)[1]),
    n : length(A[1]),
    i : 1,
    j : 1,

    loop,

    ir : kaidan_findnonzero(A, i, j, m),
    if ir = -1 then
    (
        if j = n then
            go(finish)
        else
        (
            j : j + 1,
            go(loop)
        )
    )
    else if ir > i then
    (
        col_ex(A, i, ir),
        print(A)
    ),

    if not A[i, j] = 1 then
    (
        col_sp(A, 1/A[i, j], i),
        print(A)
    ),

    if i > 1 then
        for k : 1 thru i - 1 do
            if not A[k, j] = 0 then
            (
                col_ad(A, -A[k, j], k, i),
                print(A)
            ),

    if i < m then
        for k : i + 1 thru m do
            if not A[k, j] = 0 then
            (
                col_ad(A, -A[k, j], k, i),
                print(A)
            ),
 
    if i < m and j < n then
    (
        i : i + 1,
        j : j + 1,
        go(loop)
    ),

    finish,

    A
)$

(c47) cp_kaidan(A) := block([CPA],
    mode_declare(CPA, any),
    CPA : copymatrix(A),
    kaidan(CPA)
)$

(c48) k2h_findfirst(A, ii, nn) := block([ij],
    mode_declare(ij, fixnum),
    for l : ii thru nn do
        if  A[ii, l] = 1 then
            ij : l
        else
            ij : l + 1,
        if ij <= nn then
            return(ij)
        else
            return(-1)
)$

(c49) kai2hyo(A) := block([m, n, minima],
    mode_declare(m, fixnum, n, fixnum, minima, fixnum),
    if  not matrixp(A) then
    (
        print(A, "is not a matrix!"),
        []
    ),
    n : length(A[1]),
    m : length(transpose(A)[1]),
    if m > n then
        minima : n
    else
        minima : m,
    block([j],
        mode_declare(j, fixnum),
        for i : 1 thru minima do
        (
            j : k2h_findfirst(A, i, n),
            if (not j = -1) and (not i = j) then
            (
                row_ex(A, i, j),
                print(A)
            ),
            if i < n then
                for k : i + 1 thru n do
                    if not A[i, k] = 0 then
                    (
                        row_ad(A, -A[i, k], k, i),
                        print(A)
                    )
        )
    ),
    A
)$

(c50) cp_kai2hyo(A) := block([CPA],
    mode_declare(CPA, any),
    CPA : copymatrix(A),
    kai2hyo(CPA)
)$

(c51) hyoujyun(A) := kai2hyo(kaidan(A))$

(c52) cp_hyoujyun(A) := kai2hyo(cp_kaidan(A))$

(c53) matblock(Mat, i1, i2, j1, j2) := block([A, m, n],
    mode_declare(A, any, m, fixnum, n, fixnum),
    m : i2 - i1,
    n : j2 - j1,
    A : zeromatrix(m + 1, n + 1),
    for i : 1 thru m + 1 do
        for j : 1 thru n + 1 do
            A[i][j] : Mat[i1 + i - 1][j1 + j - 1],
    A
)$

(c54) matreplace(Mat1, Mat2, i1, i2, j1, j2) := block([A, m, n],
    mode_declare(A, any, m, fixnum, n, fixnum),
    m : i2 - i1,
    n : j2 - j1,
    A : copymatrix(Mat1),
    for i : 1 thru m + 1 do
        for j : 1 thru n + 1 do
            A[i1 + i - 1][j1 + j - 1] : Mat2[i][j],
    A
)$

(c55) myinv(A) := block([I, CPA, AI, P, m, n],
    mode_declare(I, any, CPA, any, AI, any, P, any, m, fixnum, n, fixnum),
    m : length(transpose(A)[1]),
    n : length(A[1]),
    if not m = n then
    (
        print(A, "is not a square matrix!"),
        []
    ),
    if not m = rank(A) then
    (
        print(A, "in not a non-singular matrix1"),
        []
    ),
    CPA : copymatrix(A),
    I : copymatrix(ident(m)),
    AI : addcol(CPA, I),
    print(AI),
    kaidan(AI),
    P : matblock(AI, 1, m, n + 1, n + m),
    P
)$

(c56) qp_findnonzero(A, pi, pj, mm) := block([pix],
    mode_declare(pix, fixnum),
    pix : pi,
    loop1,
        if not A[pix, pj] = 0 then
            pix
        else if pix = mm then
            -1
        else
        (
            pix : pix + 1,
            go(loop1)
        )
)$

(c57) queryp(A) := block([P, CPA, m, n, i, j, ir, dummy, c],
    mode_declare(P, any, CPA, any, m, fixnum, n, fixnum,
                 i, fixnum, j, fixnum, ir, fixnum, dummy, any, c, any),
    if not matrixp(A) then
    (
        print(A, "is not a matrix!"),
        []
    ),
    m : length(transpose(A)[1]),
    n : length(A[1]),

    P : copymatrix(ident(m)),
    CPA : copymatrix(A),
    i : 1,
    j : 1,
    loop,
        ir : qp_findnonzero(CPA, i, j, m),
        if ir = -1 then
        (
            if j = n then
                go(finish)
            else
            (
                j : j + 1,
                go(loop)
            )
        )
        else if ir > 1 then
        (
            dummy : CPA[i],
            CPA[i] : CPA[ir],
            CPA[ir] : dummy,
            dummy : P[i],
            P[i] : P[ir],
            P[ir] : dummy,
            print(CPA, "    ", P)
        ),
        c : 1 / CPA[i, j],
        CPA[i] : c * CPA[i],
        P[i] : c * P[i],
        print(CPA, "    ", P),
        if i > 1 then
            for k : 1 thru i - 1 do
                if not CPA[k, j] = 0 then
                (
                    c : CPA[k, j],
                    CPA[k] : CPA[k] - c * CPA[i],
                    P[k] : P[k] - c * P[i],
                    print(CPA, "    ", P)
                ),
        if i < m then
            for k : i + 1 thru m do
                if not CPA[k, j] = 0 then
                (
                    c : CPA[k, j],
                    CPA[k] : CPA[k] - c * CPA[i],
                    P[k] : P[k] - c * P[i],
                    print(CPA, "    ", P)
                ),
        if i < m and j < n then
        (
            i : i + 1,
            j : j + 1,
            go(loop)
        ),
    finish,
    P
)$

(c58) qk_findnonzero(A, pi, pj, mm) := block([pix],
    mode_declare(pix, fixnum),
    pix : pi,
    loop1,
    if not A[pix, pj] = 0 then
        pix
    else if pix = mm then
        -1
    else
    (
        pix : pix + 1,
        go(loop1)
    )
)$

(c59) qaidan(A) := block([m, n, i, j, ir],
    mode_declare(m, fixnum, n, fixnum, i, fixnum, j, fixnum, ir, fixnum),

    if not matrixp(A) then
    (
        print(A, " is not a matrix!"),
        []
    ),
    m : length(transpose(A)[1]),
    n : length(A[1]),
    i : 1,
    j : 1,
    loop,
        ir : qk_findnonzero(A, i, j, m),
        if ir = -1 then
        (
            if j = n then
                go(finish)
            else
            (
                j : j + 1,
                go(loop)
            )
        )
        else if ir > 1 then
            col_ex(A, i, ir),
        if not A[i, j] = 1 then
            col_sp(A, 1/A[i, j], i),
        if i > 1 then
            for k : 1 thru i - 1 do
                if not A[k, j] = 0 then
                    col_ad(A, -A[k, j], k, i),
        if i < m then
            for k : i + 1 thru m do
                if not A[k, j] = 0 then
                    col_ad(A, -A[k, j], k, i),
        if i < m and j < n then
        (
            i : i + 1,
            j : j + 1,
            go(loop)
        ),
    finish,
    A
)$
[*list:931{55%}; fixnum:51{2%}; ut:82%]

(c60) qq_findfirst(A, ii, nn) := block([ij],
    mode_declare(ij, fixnum),
    for l : ii thru nn do
        if A[ii, l] = 1 then
            ij : l
        else
            ij : l + 1,
    if ij <= nn then
        return(ij)
    else
        return(-1)
)$

(c61) queryq(A) := block([m, n, minima, Q, CPA],
    mode_declare(m, fixnum, n, fixnum, minima, fixnum, Q, any, CPA, any),
    if not matrixp(A) then
    (
        print(A, " is not a matrix!"),
        []
    ),
    n : length(A[1]),
    m : length(transpose(A)[1]),
    Q : copymatrix(ident(n)),
    CPA : copymatrix(A),
    if m > n then
        minima : n
    else
        minima : m,

    block([j],
        mode_declare(j, fixnum),
        for i : 1 thru minima do
        (
            j : qq_findfirst(CPA, i, n),
            if (not j = -1) and (not i = j) then
            (
                row_ex(CPA, i, j),
                row_ex(Q, i, j),
                print(Q)
            ),
            if i < n then
                for k : i + 1 thru n do
                    if not CPA[i, k] = 0 then
                    (
                        row_ad(Q, -CPA[i, k], k, i),
                        row_ad(CPA, -CPA[i, k], k, i),
                        print(Q)
                    )
        )
    ),
    Q
)$

(c62) getpq(A) := block([P, Q, B, C],
    mode_declare(P, any, Q, any, B, any, C, any),
    B : copymatrix(A),
    C : qaidan(B),
    P : queryp(A),
    Q : queryq(C),
    [P, Q]
)$

(c63) ks_findnonzero(A, pi, pj, mm) := block([pix],
    mode_declare(pix, fixnum),
    pix : pi,
    loop1,
    if not A[pix, pj] = 0 then
        pix
    else if pix = mm then
        -1
    else
    (
        pix : pix + 1,
        go(loop1)
    )
)$

(c64) kamisankaku(A) := block([m, n, i, ir],
    mode_declare(m, fixnum, n, fixnum, i, fixnum, ir, fixnum),
    if not matrixp(A) then
    (
        print(A, " is not a matrix!"),
        []
    ),
    m : length(transpose(A)[1]),
    n : length(A[1]),
    if not m = n then
    (
        print(A, " is not a square matrix!"),
        []
    ),
    i : 1,
    loop,
        if i = m then
            go(finish),
        if not A[i, i] = 0 then
        (
            for k : i + 1 thru m do
                if not A[k, i] = 0 then
                (
                    col_ad(A, -(A[k, i]/A[i, i]), k, i),
                    print(A)
                )
        )
        else
        (
            ir : ks_findnonzero(A, i, i, m),
            if ir = -1 then
            (
                i : i + 1,
                go(loop)
            )
            else
            (
                col_ad(A, 1, i, ir),
                print(A),
                for k : i + 1 thru m do
                    if not A[k, i] = 0 then
                    (
                        col_ad(A, -(A[k, i]/A[i, i]), k, i),
                        print(A)
                    )
            )
        ),
        if i < m then
        (
            i : i + 1,
            go(loop)
        ),
    finish,
    A
)$

(c65) cp_kamisankaku(A) := block([CPA],
    mode_declare(CPA, any),
    CPA : copymatrix(A),
    kamisankaku(CPA)
)$

(c66) lu_findnonzero(A, pi, pj, mm) := block([pix],
    mode_declare(pix, fixnum),
    pix : pi,
    loop1,
        if not A[pix, pj] = 0 then
            pix
        else if pix = mm then
            -1
        else
        (
            pix : pix + 1,
            go(loop1)
        )
)$

(c67) luresolution(A) := block([P, D, L, LU, I, m, n, i, ir],
    mode_declare(P, any, D, any, L, any, LU, any, I, any,
                 m, fixnum, n, fixnum, i, fixnum, ir, fixnum),
    if not matrixp(A) then
    (
        print(A, " is not a matrix!"),
        []
    ),
    m : length(transpose(A)[1]),
    n : length(A[1]),
    if not m = n then
    (
        print(A, " is not a square-matrix!"),
        []
    ),
    P : copymatrix(ident(n)),
    i : 1,
    loop,
        if i = m then
            go(finish),
        if not A[i, i] = 0 then
        (
            for k : i + 1 thru m do
                if not A[k, i] = 0 then
                (
                    col_ad(P, -(A[k, i]/A[i, i]), k, i),
                    col_ad(A, -(A[k, i]/A[i, i]), k, i)
                )
        )
        else
        (
            ir : lu_findnonzero(A, i, i, m),
            if ir = -1 then
            (
                i : i + 1,
                go(loop)
            )
            else
            (
                col_ad(P, 1, i, ir),
                col_ad(A, 1, i, ir),
                for k : i + 1 thru m do
                    if not A[k, i] = 0 then
                    (
                        col_ad(P, -(A[k, i]/A[i, i]), k, i),
                        col_ad(A, -(A[k, i]/A[i, i]), k, i)
                    )
            )
        ),
        if i < m then
        (
            i : i + 1,
            go(loop)
        ),
    finish,
    I : copymatrix(ident(n)),
    D : addcol(P, I),
    qaidan(D),
    L : matblock(D, 1, n, n + 1, n + n),
    LU : [L, A],
    LU
)$

(c68) cp_luresolution(A) := block([B],
    mode_declare(B, any),
    B : copymatrix(A),
    luresolution(B)
)$

(c69) lurslv(A) := luresolution(A)$

(c70) cp_lurslv(A) := cp_luresolution(A)$

(c71) lfctr(A) := first(cp_luresolution(A))$

(c72) ufctr(A) := first(rest(cp_luresolution(A)))$

(c73) gyouretsushiki(A) := block([B, detm, m, n],
    mode_declare(B, any, detm, any, m, fixnum, n, fixnum),
    B : copymatrix(A),
    m : length(transpose(B)[1]),
    n : length(B[1]),
    if not m = n then
    (
        print(A, "is not a square-matrix!"),
        []
    ),
    kamisankaku(B),
    detm : prod(B[k, k], k, 1, m),
    detm
)$

(c74) mydet(A) := gyouretsushiki(A)$

(c75) listtomatrix(a) := block([tv, m, b],
    mode_declare(tv, any, m, fixnum, b, list),
    if listp(a) then
    (
        b : copylist(a),
        m : length(a),
        tv : matrix(first(b)),
        b : rest(b),

        for i : 1 thru m - 1 do
            tv : addrow(tv, b[i]),
        tv
    )
    else
        []
)$

(c76) matrixtolist(a) := block([tl, m],
    mode_declare(tl, list, m, fixnum),
    if matrixp(a) then
    (
        tl : [],
        m : length(a),
        for i : m thru 1 step -1 do
            tl : cons(a[i], tl),
        tl
    )
    else
        []
)$

(c77) ltom(a) := listtomatrix(a)$

(c78) mtol(a) := matrixtolist(a)$

(c79) diagonalizer(A) := block([B, AL],
    mode_declare(B, any, AL, list),
    AL : maplist(lambda([x], unitvector(x)),
                 gramschmidt(rest(eigenvectors(A)))),
    B : ltom(AL),
    ratsimp(B)
)$

(c80) diagonalize(A) := block([P],
    mode_declare(P, any),
    P : diagonalizer(A),
    ratsimp(ratsimp(P) . A . ratsimp(transpose(P)))
)$
[*list:931{56%}; fixnum:51{2%}; ut:77%]

(c81) unitarize(A) := block([AL],
    mode_declare(AL, list),
    AL : maplist(lambda([x], unitvector(x)),
                 gramschmidt(mtol(transpose(A)))),
    ratsimp(transpose(ltom(AL)))
)$

(c82) gson(AL) := maplist(lambda([x], unitvector(x)), gramschmidt(AL))$

(c83) qs_findnonzero(A, pi, pj, mm) := block([pix],
    mode_declare(pix, fixnum),
    pix : pi,
    loop1,
        if not A[pix, pj] = 0 then
            pix
        else if pix = mm then
            -1
        else
        (
            pix : pix + 1,
            go(loop1)
        )
)$

(c84) qamisankaku(A) := block([m, n, i, ir],
    mode_declare(m, fixnum, n, fixnum, i, fixnum, ir, fixnum),
    if not matrixp(A) then
    (
        print(A, " is not a matrix!"),
        []
    ),
    m : length(transpose(A)[1]),
    n : length(A[1]),
    if not m = n then
    (
        print(A, " is not a square-matrix!"),
        []
    ),
    i : 1,
    loop,
        if i = m then
            go(finish),
        if not A[i, i] = 0 then
            for k : i + 1 thru m do
                col_ad(A, -(A[k, i]/A[i, i]), k, i)
        else
        (
            ir : qs_findnonzero(A, i, i, m),
            if ir = -1 then
            (
                i : i + 1,
                go(loop)
            )
            else
            (
                col_ad(A, 1, i, ir),
                for k : i + 1 thru m do
                    if not A[k, i] = 0 then
                        col_ad(A, -(A[k, i]/A[i, i]), k, i)
            )
        ),
        if i < m then
        (
            i : i + 1,
            go(loop)
        ),
    finish,
    A
)$

Batching done.(d85)                                 kihon.v

(c86) A : matrix([1,0,2],[0,1,2],[2,2,-1]);
                                 [ 1  0   2  ]
                                 [           ]
(d86)                            [ 0  1   2  ]
                                 [           ]
                                 [ 2  2  - 1 ]

(c87) P : diagonalizer(A);
[autoload rat/result]
[fasl /usr/mac/rat/result.o]

/usr/mac/maxsrc/sublis.o being loaded.
[fasl /usr/mac/maxsrc/sublis.o]
[*list:931{58%}; fixnum:51{2%}; ut:71%]
                       [    1         1           2    ]
                       [ -------   -------   - ------- ]
                       [ sqrt(6)   sqrt(6)     sqrt(6) ]
                       [                               ]
                       [    1         1          1     ]
(d87)                  [ -------   -------    -------  ]
                       [ sqrt(3)   sqrt(3)    sqrt(3)  ]
                       [                               ]
                       [    1          1               ]
                       [ -------  - -------      0     ]
                       [ sqrt(2)    sqrt(2)            ]

(c88) D : P . A . transpose(P);
                                 [ - 3  0  0 ]
                                 [           ]
(d88)                            [  0   3  0 ]
                                 [           ]
                                 [  0   0  1 ]

(c89) quit();
4.1BSD kenrou% logout


Berkeley 4.1 VAX/UNIX (Amnesia-Vax)

login: