[~/]kenrou% telnet localhost 2238 Trying 127.0.0.1... Connected to localhost. Escape character is '^]'. Connection closed by foreign host. [~/]kenrou% telnet localhost 2238 Trying 127.0.0.1... Connected to localhost. Escape character is '^]'. Connected to the VAX780 simulator DZ device, line 1 Berkeley 4.1 VAX/UNIX (Amnesia-Vax) login: kenrou Password: Last login: Sat Aug 6 19:52:45 on tty00 Welcome to Berkeley Vax/UNIX (4.1bsd revised 1 Sept. 1981) You have mail. 4.1BSD kenrou% cd work/Programs/vaxima 4.1BSD kenrou% ls kihon.v machin.v primer.txt tensor.v 4.1BSD kenrou% vaxima Vaxima 2.04 Sun Sep 12 15:31:19 1982 (c1) (x+1)**3; 3 (d1) (x + 1) (c2) expand(d1); 3 2 (d2) x + 3 x + 3 x + 1 (c3) load("eigen.v"); Batching the file /usr/mac/share/eigen.v (c4) /* 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] (c5) SSTATUS(FEATURE,EIGEN)$ (c6) HERMITIANMATRIX:FALSE$ (c7) NONDIAGONALIZABLE:FALSE$ (c8) KNOWNEIGVALS:FALSE$ (c9) KNOWNEIGVECTS:FALSE$ (c10) LISTEIGVECTS:[]$ (c11) LISTEIGVALS:[]$ (c12) CONJUGATE(X):=SUBLIS('[%I=-%I],X)$ (c13) INNERPRODUCT(X,Y):=BLOCK([LISTARITH],LISTARITH:TRUE,RATSIMP(CONJUGATE(X).Y))$ (c14) UNITVECTOR(X):=BLOCK([LISTARITH,INTRN],LISTARITH:TRUE,INTRN:INNERPRODUCT(X,X), INTRN:SQRT(INTRN),X/INTRN)$ (c15) COLUMNVECTOR(X):=TRANSPOSE(MATRIX(X))$ (c16) 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)))$ (c17) 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)))$ (c18) 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)))$ (c19) /* 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))$ (c20) 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)))$ (c21) 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)))$ (c22) CONJ(X):=CONJUGATE(X)$ (c23) INPROD(X,Y):=INNERPRODUCT(X,Y)$ (c24) UVECT(X):=UNITVECTOR(X)$ (c25) COVECT(X):=COLUMNVECTOR(X)$ (c26) GSCHMIT(X):=GRAMSCHMIDT(X)$ (c27) EIVALS(MAT):=EIGENVALUES(MAT)$ (c28) EIVECTS(MAT):=EIGENVECTORS(MAT)$ (c29) UEIVECTS(MAT):=UNITEIGENVECTORS(MAT)$ (c30) SIMTRAN(MAT):=SIMILARITYTRANSFORM(MAT)$ Batching done.(d31) /usr/mac/share/eigen.v (c32) load("kihon.v"); Batching the file kihon.v (c33) /* kihon.v */ colvector(A) := transpose(matrix(A))$ (c34) rowvector(A) := matrix(A)$ (c35) 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 )$ [*symbol:342{90%}; list:931{47%}; fixnum:51{2%}; ut:0%] (c36) cp_col_ad(A, k, i, j) := block([B], mode_declare(B, any), B : copymatrix(A), B[i] : B[i] + k * B[j], B )$ (c37) cp_col_sp(A, k, i) := block([B], mode_declare(B, any), B : copymatrix(A), B[i] : k * B[i], B )$ (c38) cp_row_ex(A, i, j) := block([B], mode_declare(B, any), B : cp_col_ex(transpose(A), i, j), transpose(B) )$ (c39) cp_row_ad(A, k, i, j) := block([B], mode_declare(B, any), B : cp_col_ad(transpose(A), k, i, j), transpose(B) )$ (c40) cp_row_sp(A, k, i) := block([B], mode_declare(B, matrix), B : cp_col_sp(transpose(A), k, i), transpose(B) )$ (c41) col_sp(A, k, i) := block( A[i] : k * A[i], A )$ (c42) col_ad(A, k, i, j) := block( A[i] : A[i] + k * A[j], A )$ (c43) col_ex(A, i, j) := block([dummy], dummy : A[i], A[i] : A[j], A[j] : dummy, A )$ (c44) 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 )$ (c45) 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 )$ (c46) 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 )$ (c47) 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) ) )$ (c48) 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 )$ (c49) cp_kaidan(A) := block([CPA], mode_declare(CPA, any), CPA : copymatrix(A), kaidan(CPA) )$ (c50) 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) )$ (c51) 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 )$ (c52) cp_kai2hyo(A) := block([CPA], mode_declare(CPA, any), CPA : copymatrix(A), kai2hyo(CPA) )$ (c53) hyoujyun(A) := kai2hyo(kaidan(A))$ (c54) cp_hyoujyun(A) := kai2hyo(cp_kaidan(A))$ (c55) 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 )$ (c56) 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 )$ (c57) 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 )$ (c58) 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) ) )$ (c59) 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 )$ (c60) 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) ) )$ (c61) 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)[*list:931{54%}; fixnum:51{2%}; ut:78%] ), finish, A )$ (c62) 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) )$ (c63) 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 )$ (c64) 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] )$ (c65) 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) ) )$ (c66) 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 )$ (c67) cp_kamisankaku(A) := block([CPA], mode_declare(CPA, any), CPA : copymatrix(A), kamisankaku(CPA) )$ (c68) 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) ) )$ (c69) 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 )$ (c70) cp_luresolution(A) := block([B], mode_declare(B, any), B : copymatrix(A), luresolution(B) )$ (c71) lurslv(A) := luresolution(A)$ (c72) cp_lurslv(A) := cp_luresolution(A)$ (c73) lfctr(A) := first(cp_luresolution(A))$ (c74) ufctr(A) := first(rest(cp_luresolution(A)))$ (c75) 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 )$ (c76) mydet(A) := gyouretsushiki(A)$ (c77) 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 [] )$ (c78) 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 [] )$ (c79) ltom(a) := listtomatrix(a)$ (c80) mtol(a) := matrixtolist(a)$ (c81) 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) )$ [*list:931{56%}; fixnum:51{2%}; ut:80%] (c82) diagonalize(A) := block([P], mode_declare(P, any), P : diagonalizer(A), ratsimp(ratsimp(P) . A . ratsimp(transpose(P))) )$ (c83) unitarize(A) := block([AL], mode_declare(AL, list), AL : maplist(lambda([x], unitvector(x)), gramschmidt(mtol(transpose(A)))), ratsimp(transpose(ltom(AL))) )$ (c84) gson(AL) := maplist(lambda([x], unitvector(x)), gramschmidt(AL))$ (c85) 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) ) )$ (c86) 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.(d87) kihon.v (c88) load("tensor.v"); Batching the file tensor.v (c89) /* tensor.v */ rectyp(typ, tmp) := block( tmp : first(tmp), if listp(tmp) then ( typ : cons(length(tmp), typ), rectyp(typ, tmp) ) else typ )$ (c90) pseudotype(f) := block([ttyp, ttmp], mode_declare(ttyp, list, ttmp, list), ttmp : copylist(f), ttyp : cons(length(ttmp), []), reverse(rectyp(ttyp, ttmp)) )$ (c91) pseudoproduct(a, b) := block( if not listp(first(a)) then maplist(lambda([x], x * b), a) else maplist(lambda([x], pseudoproduct(x, b)), a) )$ (c92) zerolist(n) := if n = 0 then [] else cons(0, zerolist(n - 1))$ (c93) zerotensor(ttyp) := block( if (not listp(ttyp)) or (ttyp = []) then ( print("the argument is not a type of tensor"), [] ) else if rest(ttyp) = [] then zerolist(first(ttyp)) else pseudoproduct(zerolist(first(ttyp)), zerotensor(rest(ttyp))) )$ (c94) tensortype(f) := block([typ, za, ta], mode_declare(typ, list, za, list, ta, list), typ : pseudotype(f), za : zerotensor(typ), ta : 0 * f, if za = ta then typ else [] )$ (c95) tensorp(a) := block([za, ta], mode_declare(za, list, ta, list), za : zerotensor(pseudotype(a)), ta : 0 * a, if za = ta then true else false )$ (c96) sametensortypep(a, b) := block([za : 0 * a, zb : 0 * b], mode_declare(za, list, zb, list), if za = zb then true else false )$ (c97) tensorrank(a) := length(tensortype(a))$ (c98) tensorproduct(a, b) := block( mode_declare(b, list), if tensorp(a) and tensorp(b) then pseudoproduct(a, b) else [] )$ (c99) infix("oxo")$ (c100) a oxo b := tensorproduct(a, b)$ (c101) exchinext(a) := block([b, c, ni, nip], mode_declare(b, list, c, list, ni, fixnum, nip, fixnum), ni : length(a), nip : length(a[1]), b : [], for j : 1 thru nip do ( c : [], for k : 1 thru ni do c : endcons(a[k][j], c), b : endcons(c, b) ), b )$ (c102) matrixlikep(a) := block( if tensorrank(a) = 2 then true else false )$ (c103) exchnext(a, i) := block([n], mode_declare(n, fixnum), n : tensorrank(a), if (i > n - 1) or (i < 1) then ( print("Index is out of range"), [] ) else if i = 1 then exchinext(a) else maplist(lambda([x], exchnext(x, i - 1)), a) )$ (c104) tensortranspose(a, i, j) := block([b], mode_declare(b, list), if not tensorp(a) then ( print("Object is not a tensor"), [] ) else if i = j then ( print("Invalid indices"), [] ) else if i > j then tensortranspose(a, j, i) else if j = i + 1 then exchnext(a, i) else ( b : copylist(a), for k : i thru j - 1 do b : exchnext(b, k), for k : j - 2 thru i step -1 do b : exchnext(b, k), b ) )$ (c105) exchangetensortype(a, i, j) := tensortranspose(a, i, j)$ (c106) exchttype(a, i, j) := tensortranspose(a, i, j)$ (c107) squarematrixlikep(a) := block( if tensorrank(a) # 2 then false else if length(a) # length(a[1]) then false else true )$ (c108) pseudotrace(a) := block([n, tra], mode_declare(n, fixnum, tra, any), if not squarematrixlikep(a) then ( print("Indices are not square-matrix-like"), [] ) else ( n : length(a), tra : sum(a[i][i], i, 1, n), tra ) )$ (c109) 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 [] )$ (c110) 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 [] )$ (c111) ltom(a) := listtomatrix(a)$ (c112) mtol(a) := matrixtolist(a)$ (c113) tomatrixlast2(a) := block( if squarematrixlikep(a) then listtomatrix(a) else maplist(lambda([x], tomatrixlast2(x)), a) )$ (c114) matrixtoscalar(a) := block([n], mode_declare(n, fixnum), if matrixp(a) then ( n : length(a), sum(a[i][i], i, 1, n) ) else maplist(lambda([x], matrixtoscalar(x)), a) )$ (c115) contractlast2(a) := matrixtoscalar(tomatrixlast2(a))$ (c116) tensorcontract(a, i, j) := block([n, b], mode_declare(n, fixnum, b, list), n : tensorrank(a), if n = 0 then ( print("Object is not a tensor"), [] ) else if (i = j) or (i < 1) or [*list:931{60%}; fixnum:51{2%}; ut:80%] (i > n) or (j < 1) or (j > n) then ( print("Invalid indices"), [] ) else if i > j then tensorcontract(a, j, i) else if i = n - 1 then contractlast2(a) else ( b : copylist(a), if j # n then for p : j thru n - 1 do b : exchnext(b, p), if i # n - 1 then for p : i thru n - 2 do b : exchnext(b, p), b : contractlast2(b), b ) )$ (c117) cntrct(a, i, j) := tensorcontract(a, i, j)$ (c118) simplistp(f) := block( if not listp(first(f)) then true else false )$ (c119) tensorvectorcontract(a, b) := block([c : copylist(a)], mode_declare(c, list), if (tensorrank(b) # 1) or (not tensorp(a)) then ( print("Objects have invalid types"), [] ) else if first(reverse(tensortype(c))) = length(b) then contractlast2(tensorproduct(a, b)) else [] )$ (c120) infix("tvc")$ (c121) a tvc b := tensorvectorcontract(a, b)$ (c122) tensortensorcontract(a, b, i, j) := block([n], mode_declare(n, fixnum), n : tensorrank(a), if (n = 0) or (not tensorp(b)) then ( print("Objects may not be tensor"), [] ) else tensorcontract(tensorproduct(a, b), i, n + j) )$ (c123) ttc(a, b, i, j) := tensortensorcontract(a, b, i, j)$ (c124) rec(xx, acc) := block( if xx = [] then acc else if not listp(xx) then cons(xx, acc) else rec(first(xx), rec(rest(xx), acc)) )$ (c125) lflatten(x) := rec(x, [])$ (c126) flatten1(a) := maplist(lambda([x], lflatten(x)), a)$ (c127) demote(x) := block([tx, tt], mode_declare(tx, list, tt, list), tt : [], tx : tensortranspose(copylist(x), 2, 3), loop, if first(tx) # [] then ( tt : append(tt, flatten1(first(tx))), tx : rest(tx) ), if tx # [] then go(loop), tt )$ (c128) pauli(g1, g2) := demote(tensorproduct(g1, g2))$ (c129) pauli4(g1, g2, g3, g4) := pauli(pauli(g1, g2), pauli(g3, g4))$ (c130) pauli8(g1, g2, g3, g4, g5, g6, g7, g8) := pauli(pauli4(g1, g2, g3, g4), pauli4(g5, g6, g7, g8))$ Batching done.(d131) tensor.v (c132) linenum:0$ (c1) a1 : matrix([1/2],[1],[3/2]); [ 1 ] [ - ] [ 2 ] [ ] (d1) [ 1 ] [ ] [ 3 ] [ - ] [ 2 ] (c2) a2 : matrix([1],[5/4],[3/2]); [ 1 ] [ ] [ 5 ] [ - ] (d2) [ 4 ] [ ] [ 3 ] [ - ] [ 2 ] (c3) a3 : matrix([7/3],[8/3],[3]); [ 7 ] [ - ] [ 3 ] [ ] (d3) [ 8 ] [ - ] [ 3 ] [ ] [ 3 ] (c4) A : addcol(a1, a2, a3); [ 1 7 ] [ - 1 - ] [ 2 3 ] [ ] [ 5 8 ] (d4) [ 1 - - ] [ 4 3 ] [ ] [ 3 3 ] [ - - 3 ] [ 2 2 ] (c5) cp_kaidan(A); [ 14 ] [ 1 2 -- ] [ 3 ] [ ] [ 5 8 ] [ 1 - - ] [ 4 3 ] [ ] [ 3 3 ] [ - - 3 ] [ 2 2 ] [ 14 ] [ 1 2 -- ] [ 3 ] [ ] [ 3 ] [ 0 - - - 2 ] [ 4 ] [ ] [ 3 3 ] [ - - 3 ] [ 2 2 ] [ 14 ] [ 1 2 -- ] [ 3 ] [ ] [ 3 ] [ 0 - - - 2 ] [ 4 ] [ ] [ 3 ] [ 0 - - - 4 ] [ 2 ] [ 14 ] [ 1 2 -- ] [ 3 ] [ ] [ 8 ] [ 0 1 - ] [ 3 ] [ ] [ 3 ] [ 0 - - - 4 ] [ 2 ] [ 2 ] [ 1 0 - - ] [ 3 ] [ ] [ 8 ] [ 0 1 - ] [ 3 ] [ ] [ 3 ] [ 0 - - - 4 ] [ 2 ] [ 2 ] [ 1 0 - - ] [ 3 ] [ ] [ 8 ] [ 0 1 - ] [ 3 ] [ ] [ 0 0 0 ] [ 2 ] [ 1 0 - - ] [ 3 ] [ ] (d5) [ 8 ] [ 0 1 - ] [ 3 ] [ ] [ 0 0 0 ] (c6) rank(A); (d6) 2 (c7) gyouretsushiki(A); [ 1 7 ] [ - 1 - ] [ 2 3 ] [ ] [ 3 ] [ 0 - - - 2 ] [ 4 ] [ ] [ 3 3 ] [ - - 3 ] [ 2 2 ] [ 1 7 ] [ - 1 - ] [ 2 3 ] [ ] [ 3 ] [ 0 - - - 2 ] [ 4 ] [ ] [ 3 ] [ 0 - - - 4 ] [ 2 ] [*list:941{61%}; fixnum:51{2%}; ut:37%] [ 1 7 ] [ - 1 - ] [ 2 3 ] [ ] [ 3 ] [ 0 - - - 2 ] [ 4 ] [ ] [ 0 0 0 ] (d7) 0 (c8) kill(a1, a2, a3, A); (d8) done (c9) a1 : matrix([1],[2],[0]); [ 1 ] [ ] (d9) [ 2 ] [ ] [ 0 ] (c10) a2 : matrix([4],[3],[0]); [ 4 ] [ ] (d10) [ 3 ] [ ] [ 0 ] (c11) a3 : matrix([0],[2],[1]); [ 0 ] [ ] (d11) [ 2 ] [ ] [ 1 ] (c12) A : addcol(A); First argument to ADDCOL must be a matrix (c13) A : addcol(a1, a2, a3); [ 1 4 0 ] [ ] (d13) [ 2 3 2 ] [ ] [ 0 0 1 ] (c14) AL : matrixtolist(transpose(A)); (d14) [[1, 2, 0], [4, 3, 0], [0, 2, 1]] (c15) AAL : gson(AL); /usr/mac/maxsrc/sublis.o being loaded. [fasl /usr/mac/maxsrc/sublis.o] 1 2 2 1 (d15) [[-------, -------, 0], [-------, - -------, 0], [0, 0, 1]] sqrt(5) sqrt(5) sqrt(5) sqrt(5) (c16) AA : transpose(listtomatrix((AAL)); Missing ")" aa : transpose ( listtomatrix ( ( aal ) ) **$** Please rephrase or edit (c16) AA : transpose(listtomatrix(AAL)); [ 1 2 ] [ ------- ------- 0 ] [ sqrt(5) sqrt(5) ] [ ] (d16) [ 2 1 ] [ ------- - ------- 0 ] [ sqrt(5) sqrt(5) ] [ ] [ 0 0 1 ] (c17) kill(a1, a2, a3, A, AL, AAL, AA); (d17) done (c18) A : matrix([2,1,1],[1,2,1],[1,1,2]); [ 2 1 1 ] [ ] (d18) [ 1 2 1 ] [ ] [ 1 1 2 ] (c19) eeq : [gyouretsushiki(A - x * ident(3)) = 0]; [ 2 - x 1 1 ] [ ] [ 1 1 ] [ 0 - x - ----- + 2 1 - ----- ] [ 2 - x 2 - x ] [ ] [ 1 1 2 - x ] [ 2 - x 1 1 ] [ ] [ 1 1 ] [ 0 - x - ----- + 2 1 - ----- ] [ 2 - x 2 - x ] [ ] [ 1 1 ] [ 0 1 - ----- - x - ----- + 2 ] [ 2 - x 2 - x ] [ 2 - x 1 1 ] [ ] [ 1 1 ] [ 0 - x - ----- + 2 1 - ----- ] [ 2 - x 2 - x ] [ ] [ 1 2 ] [ (1 - -----) ] [ 2 - x 1 ] [ 0 0 - x - --------------- - ----- + 2 ] [ 1 2 - x ] [ - x - ----- + 2 ] [ 2 - x ] 1 2 (1 - -----) 1 2 - x 1 (d19) [(2 - x) (- x - ----- + 2) (- x - --------------- - ----- + 2) = 0] 2 - x 1 2 - x - x - ----- + 2 2 - x (c20) ratsimp(eeq); [*list:951{61%}; fixnum:51{2%}; ut:68%] 3 2 (d20) [- x + 6 x - 9 x + 4 = 0] (c21) solve(eeq, [x]); (d21) [x = 1, x = 4] (c22) eigenvectors(A); [autoload rat/result] [fasl /usr/mac/rat/result.o] (d22) [[[4, 1], [1, 2]], [1, 1, 1], [1, 0, - 1], [0, 1, - 1]] (c23) P : diagonalizer(A); [*list:961{60%}; fixnum:51{2%}; ut:38%] [*list:971{60%}; fixnum:51{2%}; ut:66%] [ 1 1 1 ] [ ------- ------- ------- ] [ sqrt(3) sqrt(3) sqrt(3) ] [ ] [ 1 1 ] (d23) [ ------- 0 - ------- ] [ sqrt(2) sqrt(2) ] [ ] [ sqrt(2) sqrt(2) sqrt(2) ] [ - --------- ------- - --------- ] [ 2 sqrt(3) sqrt(3) 2 sqrt(3) ] (c24) Q : transpose(P); [ 1 1 sqrt(2) ] [ ------- ------- - --------- ] [ sqrt(3) sqrt(2) 2 sqrt(3) ] [ ] [ 1 sqrt(2) ] (d24) [ ------- 0 ------- ] [ sqrt(3) sqrt(3) ] [ ] [ 1 1 sqrt(2) ] [ ------- - ------- - --------- ] [ sqrt(3) sqrt(2) 2 sqrt(3) ] (c25) P . A . Q; [ 4 0 0 ] [ ] (d25) [ 0 1 0 ] [ ] [ 0 0 1 ] (c26) DA : diagonalize(A); [*list:971{60%}; fixnum:51{2%}; ut:54%] [ 4 0 0 ] [ ] (d26) [ 0 1 0 ] [ ] [ 0 0 1 ] (c27) DN : matrix([4^n,0,0],[0,1,0],[0,0,1]); [ n ] [ 4 0 0 ] (d27) [ ] [ 0 1 0 ] [ ] [ 0 0 1 ] (c28) AN : ratsimp(Q . DN . P); [ n n n ] [ 4 + 2 4 - 1 4 - 1 ] [ ------ ------ ------ ] [ 3 3 3 ] [ ] [ n n n ] (d28) [ 4 - 1 4 + 2 4 - 1 ] [ ------ ------ ------ ] [ 3 3 3 ] [ ] [ n n n ] [ 4 - 1 4 - 1 4 + 2 ] [ ------ ------ ------ ] [ 3 3 3 ] (c29) kill(A, P, Q, D, DN, AN); (d29) done (c30) f(x) := x^2/(%e)^x; 2 x (d30) f(x) := --- x %e (c31) diff(f(x), x); - x 2 - x (d31) 2 x %e - x %e (c32) ratsimp(%d31); (d32) %d31 (c33) ratsimp(d31); 2 - x (d33) - (x - 2 x) %e (c34) diff(f(x), x, 2); [*list:971{60%}; fixnum:51{2%}; ut:65%] 2 - x - x - x (d34) x %e - 4 x %e + 2 %e (c35) ratsimp(d34); 2 - x (d35) (x - 4 x + 2) %e (c36) kill(f); (d36) done (c37) f(x) := log((%e)^(x^2) + 1); 2 x (d37) f(x) := log(%e + 1) (c38) diff(f(x), x); 2 x 2 x %e (d38) -------- 2 x %e + 1 (c39) diff(f(x), x, 2); 2 2 2 2 2 x 2 x x 4 x %e 4 x %e 2 %e (d39) - ----------- + --------- + -------- 2 2 2 x 2 x x (%e + 1) %e + 1 %e + 1 (c40) ratsimp(d39); 2 2 2 x 2 x 2 %e + (4 x + 2) %e (d40) -------------------------- 2 2 2 x x %e + 2 %e + 1 (c41) kill(f); (d41) done (c42) f(x) := asinh(x); (d42) f(x) := asinh(x) (c43) diff(f(x), x); 1 (d43) ------------ 2 sqrt(x + 1) (c44) f1(x) := 1/sqrt(x^2 + 1); 1 (d44) f1(x) := ------------ 2 sqrt(x + 1) (c45) diff(f(x), x, 2); x (d45) - ----------- 2 3/2 (x + 1) (c46) f2(x) := -x/(x^2+1)^(3/2); x (d46) f2(x) := - ----------- 2 3/2 (x + 1) (c47) diff(f(x), x, 3); 2 3 x 1 (d47) ----------- - ----------- 2 5/2 2 3/2 (x + 1) (x + 1) (c48) f3(x) := 3*x^2/(x^2+1)^(5/2) - 1/(x^2+1)^(3/2); 2 3 x 1 (d48) f3(x) := ----------- - ----------- 2 5/2 2 3/2 (x + 1) (x + 1) (c49) f(0) + f1(0)*x/1! + f2(0)*x^2/2! + f3(0)*x^3/3!; 3 x (d49) x - -- 6 (c50) taylor(asinh(x), x, 0, 3); 3 x (d50)/T/ x - -- + . . . 6 (c51) kill(f,f1,f2,f3); (d51) done (c52) f(x, y) := sin(x) * cos(y); (d52) f(x, y) := sin(x) cos(y) (c53) diff(f(x,y), x); (d53) cos(x) cos(y) (c54) fx(x,y) := cos(x) * cos(y); (d54) fx(x, y) := cos(x) cos(y) (c55) diff(f(x,y), y); (d55) - sin(x) sin(y) (c56) fy(x,y) := -sin(x) * sin(y); (d56) fy(x, y) := - sin(x) sin(y) (c57) z(x,y) := fx(%pi/6, %pi/3)*(x-%pi/6) + fy(%pi/6, %pi/3)*(y-%pi/3) + f(%pi/6, %pi/3); %pi %pi %pi %pi %pi %pi %pi %pi (d57) z(x, y) := fx(---, ---) (x - ---) + fy(---, ---) (y - ---) + f(---, ---) 6 3 6 6 3 3 6 3 (c58) z(x,y); %pi %pi sqrt(3) (y - ---) sqrt(3) (x - ---) 3 6 1 (d58) - ----------------- + ----------------- + - 4 4 4 (c59) kill(f, fx, fy, z); (d59) done (c60) f(x, y) := x * sin(y); (d60) f(x, y) := x sin(y) (c61) integrate(f(x,y), y); (d61) - x cos(y) (c62) integrate(f(x,y), y, 0, %pi/2); [*flonum:12{8%}; list:971{61%}; fixnum:51{2%}; ut:64%] [autoload macrak/rpart] [fasl /usr/mac/macrak/rpart.o] (d62) x (c63) integrtate(x, x, 0, %pi/4); %pi (d63) integrtate(x, x, 0, ---) 4 (c64) integrate(x, x, 0, %pi/4); 2 %pi (d64) ---- 32 (c65) integrate(f(x,y), x); 2 x sin(y) (d65) --------- 2 (c66) integrate(f(x,y), x, 0, %pi/4); [*list:981{61%}; fixnum:51{2%}; ut:74%] 2 %pi sin(y) (d66) ----------- 32 /* (c67), (d67) deleted */ (c68) integrate(((%pi)^2) * sin(y)/32, y, 0, %pi/2); [*list:991{61%}; fixnum:51{2%}; ut:56%] 2 %pi (d68) ---- 32 (c69) ff(x,y) := -x^2*cos(y)/2; 2 x cos(y) (d69) ff(x, y) := - --------- 2 (c70) ff(%pi/4,%pi/2)-ff(%pi/4,0)-ff(0,%pi/2)+ff(0,0); 2 %pi (d70) ---- 32 (c71) integrate(integrate(f(x,y), x, 0, %pi/4), y, 0, %pi/4); [*list:1001{60%}; fixnum:51{2%}; ut:51%] [*list:1001{60%}; fixnum:51{2%}; ut:69%] sqrt(2) 2 (1 - -------) %pi 2 (d71) ------------------ 32 (c72) integrate(integrate(f(x,y), x, 0, %pi/4), y, 0, %pi/2); 2 %pi (d72) ---- 32 (c73) kill(f, ff); (d73) done (c74) f(x,y) := x^2 + y^2; 2 2 (d74) f(x, y) := x + y /* (c75) deleted */ (c76) integrate(f(x,y), y, 0, 3*sqrt(1-x^2/4)); Is (x - 2) (x + 2) positive, negative, or zero? negative; 2 Is x - 4 negative or zero? negative; [*list:1011{60%}; fixnum:51{2%}; ut:65%] 2 2 sqrt(4 - x ) (3 x + 36) (d76) ------------------------ 8 (c77) integrate(sqrt(4-x^2)*(3*x^2+36)/8, x, -2, 2); [*list:1011{60%}; fixnum:51{2%}; ut:92%] 39 %pi (d77) ------ 4 (c78) integrate(sqrt(4-x^2)*(3*x^2+36)/8, x); 2 3/2 2 x 3 x (4 - x ) 39 x sqrt(4 - x ) 78 asin(-) - --------------- + ----------------- 2 4 2 (d78) ------------------------------------------------ 8 (c79) subst([y = 3 * v, x = 2 * u], f(x,y)); 2 2 (d79) 9 v + 4 u (c80) J1 : matrix([2, 0], [0, 3]); [ 2 0 ] (d80) [ ] [ 0 3 ] (c81) JA1 : determinant(J1); (d81) 6 (c82) g(u, v) := (4 * u^2 + 9 * v^2) * JA1; 2 2 (d82) g(u, v) := (4 u + 9 v ) ja1 (c83) subst([u = r * cos(theta), v = r * sin(theta)], g(u, v)); 2 2 2 2 (d83) 6 (9 r sin (theta) + 4 r cos (theta)) (c84) J2 : matrix([cos(theta), -r*sin(theta)], [sin(theta), r*cos(theta)]); [ cos(theta) - r sin(theta) ] (d84) [ ] [ sin(theta) r cos(theta) ] (c85) JA2 : determinant(J2); 2 2 (d85) r sin (theta) + r cos (theta) (c86) trigsimp(JA2); 2 2 (d86) trigsimp(r sin (theta) + r cos (theta)) (c87) h(r, theta) := 6 * (9*r^2*sin(theta)^2 + 4*r^2*cos(theta)^2) * JA2; 2 2 2 2 (d87) h(r, theta) := 6 (9 r sin (theta) + 4 r cos (theta)) ja2 (c88) integrate(integrate(h(r, theta), theta, 0, %pi), r, 0, 1); [*list:1011{60%}; fixnum:51{2%}; ut:15%] [*list:1021{60%}; fixnum:51{2%}; ut:69%] Is r positive, negative, or zero? positive; [*list:1021{60%}; fixnum:51{2%}; ut:52%] [*list:1021{60%}; fixnum:51{2%}; ut:77%] [*list:1031{60%}; fixnum:51{2%}; ut:66%] [*list:1031{60%}; fixnum:51{2%}; ut:66%] 39 %pi (d88) ------ 4 (c89) kill(f, J1, JA1, x, y, u, v, r, theta, J2, JA2); (d89) done (c90) eq1 : [x+2*y+3*z=6, 5*x+6*y+7*z=4, 4*x+4*y+4*z=-2]; [*list:1031{59%}; fixnum:51{2%}; ut:66%] (d90) [3 z + 2 y + x = 6, 7 z + 6 y + 5 x = 4, 4 z + 4 y + 4 x = - 2] (c91) xx : [x, y, z]; (d91) [x, y, z] (c92) A : coefmatrix(eq1, xx); [ 1 2 3 ] [ ] (d92) [ 5 6 7 ] [ ] [ 4 4 4 ] (c93) bb : matrix([6], [4], [-2]); [ 6 ] [ ] (d93) [ 4 ] [ ] [ - 2 ] (c94) Ab : addcol(A, bb); [ 1 2 3 6 ] [ ] (d94) [ 5 6 7 4 ] [ ] [ 4 4 4 - 2 ] (c95) cp_kaidan(Ab); [ 1 2 3 6 ] [ ] [ 0 - 4 - 8 - 26 ] [ ] [ 4 4 4 - 2 ] [ 1 2 3 6 ] [ ] [ 0 - 4 - 8 - 26 ] [ ] [ 0 - 4 - 8 - 26 ] [ 1 2 3 6 ] [ ] [ 13 ] [ 0 1 2 -- ] [ 2 ] [ ] [ 0 - 4 - 8 - 26 ] [ 1 0 - 1 - 7 ] [ ] [ 13 ] [ 0 1 2 -- ] [ 2 ] [ ] [ 0 - 4 - 8 - 26 ] [ 1 0 - 1 - 7 ] [ ] [ 13 ] [ 0 1 2 -- ] [ 2 ] [ ] [ 0 0 0 0 ] [ 1 0 - 1 - 7 ] [ ] [ 13 ] (d95) [ 0 1 2 -- ] [ 2 ] [ ] [ 0 0 0 0 ] (c96) rank(A); (d96) 2 (c97) rank(Ab); (d97) 2 (c98) solve(eq1, xx); Dependent equations eliminated: (3) 4 %r1 - 13 (d98) [[x = %r1 - 7, y = - ----------, z = %r1]] 2 (c99) kill(eq1, A, xx, bb, Ab); (d99) done (c100) A : matrix([42, 51, 24], [30, 36, 15], [29, 34, 12]); [ 42 51 24 ] [ ] (d100) [ 30 36 15 ] [ ] [ 29 34 12 ] (c101) myinv(A); [ 42 51 24 1 0 0 ] [ ] [ 30 36 15 0 1 0 ] [ ] [ 29 34 12 0 0 1 ] [ 17 4 1 ] [ 1 -- - -- 0 0 ] [ 14 7 42 ] [ ] [ 30 36 15 0 1 0 ] [ ] [ 29 34 12 0 0 1 ] [ 17 4 1 ] [ 1 -- - -- 0 0 ] [ 14 7 42 ] [ ] [ 3 15 5 ] [ 0 - - - -- - - 1 0 ] [ 7 7 7 ] [ ] [ 29 34 12 0 0 1 ] [ 17 4 1 ] [ 1 -- - -- 0 0 ] [ 14 7 42 ] [ ] [ 3 15 5 ] [ 0 - - - -- - - 1 0 ] [ 7 7 7 ] [ ] [ 17 32 29 ] [ 0 - -- - -- - -- 0 1 ] [ 14 7 42 ] [ 17 4 1 ] [ 1 -- - -- 0 0 ] [ 14 7 42 ] [ ] [ 5 7 ] [ 0 1 5 - - - 0 ] [ 3 3 ] [ ] [ 17 32 29 ] [ 0 - -- - -- - -- 0 1 ] [ 14 7 42 ] [ 11 17 ] [ 1 0 - -- - 2 -- 0 ] [ 2 6 ] [ ] [ 5 7 ] [ 0 1 5 - - - 0 ] [ 3 3 ] [ ] [ 17 32 29 ] [ 0 - -- - -- - -- 0 1 ] [ 14 7 42 ] [ 11 17 ] [ 1 0 - -- - 2 -- 0 ] [ 2 6 ] [ ] [ 5 7 ] [ 0 1 5 - - - 0 ] [ 3 3 ] [ ] [ 3 4 17 ] [ 0 0 - - - -- 1 ] [ 2 3 6 ] [ 11 17 ] [ 1 0 - -- - 2 -- 0 ] [ 2 6 ] [ ] [ 5 7 ] [ 0 1 5 - - - 0 ] [ 3 3 ] [ ] [ 8 17 2 ] [ 0 0 1 - - -- - ] [ 9 9 3 ] [ 26 68 11 ] [ 1 0 0 -- - -- -- ] [ 9 9 3 ] [ ] [ 5 7 ] [ 0 1 5 - - - 0 ] [ 3 3 ] [ ] [ 8 17 2 ] [ 0 0 1 - - -- - ] [ 9 9 3 ] [*list:1031{60%}; fixnum:51{2%}; ut:59%] [ 26 68 11 ] [ 1 0 0 -- - -- -- ] [ 9 9 3 ] [ ] [ 25 64 10 ] [ 0 1 0 - -- -- - -- ] [ 9 9 3 ] [ ] [ 8 17 2 ] [ 0 0 1 - - -- - ] [ 9 9 3 ] [ 26 68 11 ] [ -- - -- -- ] [ 9 9 3 ] [ ] [ 25 64 10 ] (d101) [ - -- -- - -- ] [ 9 9 3 ] [ ] [ 8 17 2 ] [ - - -- - ] [ 9 9 3 ] (c102) B : matrix([1, -1, 0, 0], [-1, 2, -1, 0], [0, -1, 2, -1], [0, 0, -1, 2]); [ 1 - 1 0 0 ] [ ] [ - 1 2 - 1 0 ] (d102) [ ] [ 0 - 1 2 - 1 ] [ ] [ 0 0 - 1 2 ] (c103) myinv(B); [ 1 - 1 0 0 1 0 0 0 ] [ ] [ - 1 2 - 1 0 0 1 0 0 ] [ ] [ 0 - 1 2 - 1 0 0 1 0 ] [ ] [ 0 0 - 1 2 0 0 0 1 ] [ 1 - 1 0 0 1 0 0 0 ] [ ] [ 0 1 - 1 0 1 1 0 0 ] [ ] [ 0 - 1 2 - 1 0 0 1 0 ] [ ] [ 0 0 - 1 2 0 0 0 1 ] [ 1 0 - 1 0 2 1 0 0 ] [ ] [ 0 1 - 1 0 1 1 0 0 ] [ ] [ 0 - 1 2 - 1 0 0 1 0 ] [ ] [ 0 0 - 1 2 0 0 0 1 ] [ 1 0 - 1 0 2 1 0 0 ] [ ] [ 0 1 - 1 0 1 1 0 0 ] [ ] [ 0 0 1 - 1 1 1 1 0 ] [ ] [ 0 0 - 1 2 0 0 0 1 ] [ 1 0 0 - 1 3 2 1 0 ] [ ] [ 0 1 - 1 0 1 1 0 0 ] [ ] [ 0 0 1 - 1 1 1 1 0 ] [ ] [ 0 0 - 1 2 0 0 0 1 ] [ 1 0 0 - 1 3 2 1 0 ] [ ] [ 0 1 0 - 1 2 2 1 0 ] [ ] [ 0 0 1 - 1 1 1 1 0 ] [ ] [ 0 0 - 1 2 0 0 0 1 ] [ 1 0 0 - 1 3 2 1 0 ] [ ] [ 0 1 0 - 1 2 2 1 0 ] [ ] [ 0 0 1 - 1 1 1 1 0 ] [ ] [ 0 0 0 1 1 1 1 1 ] [ 1 0 0 0 4 3 2 1 ] [ ] [ 0 1 0 - 1 2 2 1 0 ] [ ] [ 0 0 1 - 1 1 1 1 0 ] [ ] [ 0 0 0 1 1 1 1 1 ] [ 1 0 0 0 4 3 2 1 ] [ ] [ 0 1 0 0 3 3 2 1 ] [ ] [ 0 0 1 - 1 1 1 1 0 ] [ ] [ 0 0 0 1 1 1 1 1 ] [ 1 0 0 0 4 3 2 1 ] [ ] [ 0 1 0 0 3 3 2 1 ] [ ] [ 0 0 1 0 2 2 2 1 ] [ ] [ 0 0 0 1 1 1 1 1 ] [ 4 3 2 1 ] [ ] [ 3 3 2 1 ] (d103) [ ] [ 2 2 2 1 ] [ ] [ 1 1 1 1 ] (c104) kill(A, B); (d104) done (c105) quit(); 4.1BSD kenrou% logout Berkeley 4.1 VAX/UNIX (Amnesia-Vax) login: