diff --git a/build/lisp.tcl b/build/lisp.tcl index 0b11532b..bc062744 100644 --- a/build/lisp.tcl +++ b/build/lisp.tcl @@ -471,6 +471,7 @@ respond ";BKPT" "(quit)" respond "*" ":midas maxtul;ts mcl_mcldmp midas\r" respond "*" ":link maxtul;.good. complr,sys;ts complr\r" respond "*" ":link liblsp;gcdemn fasl,lisp;\r" +respond "*" ":link maxtul;ts utmcl,maxtul;ts mcl\r" respond "*" "complr\013" respond "_" "mrg;macros\r" @@ -535,21 +536,17 @@ respond "*" ":link sys3;ts macsym,maxdmp;loser >\r" ### build ctensr for macsyma respond "*" "macsym\013" -respond "(C1)" "translate_file(\"sharem\\;packg >\");" -respond "(C2)" "quit();" -respond ":KILL" "complr\013" -respond "_" "sharem;packg fasl_packg trlisp\r" -respond "_" "\032" -type ":kill\r" +respond "(C1)" "compile_lisp_file(translate_file(\"sharem\\;packg >\")\[2\]);" +respond "(C2)" "compile_lisp_file(translate_file(\"tensor\\;ctensr funcs\")\[2\]);" +respond "Type ALL;" "all;" +respond "Type ALL;" "all;" +respond "(C3)" "quit();" + +### build eigen for macsyma respond "*" "macsym\013" -respond "(C1)" "translate_file(\"tensor\\;ctensr funcs\");" -respond "Type ALL;" "all;" +respond "(C1)" "compile_lisp_file(translate_file(\"share\\;eigen >\")\[2\]);" respond "Type ALL;" "all;" respond "(C2)" "quit();" -respond ":KILL" "complr\013" -respond "_" "share;ctensr fasl_tensor;ctensr trlisp\r" -respond "_" "\032" -type ":kill\r" ### more lisplib stuff respond "*" "complr\013" diff --git a/src/share/eigen.8 b/src/share/eigen.8 new file mode 100755 index 00000000..74379ce1 --- /dev/null +++ b/src/share/eigen.8 @@ -0,0 +1,333 @@ + + +EVAL_WHEN(BATCH,TTYOFF:TRUE)$ +/* -*- Mode: MACSYMA -*- */ +/* THIS IS THE SOURCE CODE FOR THE NEW EIGEN PACKAGE AND IT IS + MACSYMA BATCHABLE.*/ + +/* Copyright (c) 1985 Massachusetts Institute of Technology */ +/* */ +/*This material was developed by the Scheme project at the */ +/*Massachusetts Institute of Technology, Department of */ +/*Electrical Engineering and Computer Science. Permission */ +/*to copy this software, to redistribute it, and to use it */ +/*for any purpose is granted, subject to the following restrictions and */ +/*understandings. */ +/* */ +/*1.Any copy made of this software must include this copyright notice */ +/*in full.*/ +/*2. Users of this software agree to make their best effort (a) */ +/*to return to the MIT Scheme project any improvements or extensions */ +/*that they make, so that these may be included in future releases; */ +/*and (b) to inform MIT of noteworthy uses of this software. */ +/* */ +/*3. All materials developed as a consequence of the use of this */ +/*software shall duly acknowledge such use, in accordance with the */ +/*usual standards of acknowledging credit in academic research. */ +/* */ +/*4.MIT has made no warrantee or representation that the operation */ +/*of this software will be error-free, and MIT is under no obligation */ +/*to provide any services, by way of maintenance, update, or otherwise. */ +/* */ +/*5.In conjunction with products arising from the use of this material, */ +/*there shall be no use of the name of the Massachusetts Institute of */ +/*Technology nor of any adaptation thereof in any advertising, promotional, */ +/*or sales literature without prior written consent from MIT in each case. */ +/* */ +/*Functions written in upper case are authored by Yekta Gursel. Functions*/ +/*written in lower case are authored by Nicholas Strauss. Any comments,*/ +/*criticisms, or bugs should be sent to either Yekta@mc or ncs@mc.*/ + +EVAL_WHEN(TRANSLATE, + MODEDECLARE( + [INDEX1,INDEX2,INDEX3,INDEX4,uu,COUNT,DIMNSN,%RNUM],FIXNUM), + DECLARE([%RNUM],SPECIAL) +)$ + +SSTATUS(FEATURE,EIGEN)$ + +/* Variable Definitions */ + +DEFINE_VARIABLE(HERMITIANMATRIX,FALSE,BOOLEAN)$ +DEFINE_VARIABLE(NONDIAGONALIZABLE,FALSE,BOOLEAN)$ +DEFINE_VARIABLE(KNOWNEIGVALS,FALSE,BOOLEAN)$ +DEFINE_VARIABLE(KNOWNEIGVECTS,FALSE,BOOLEAN)$ + +DEFINE_VARIABLE(LISTEIGVECTS,[],ANY)$ +DEFINE_VARIABLE(LISTEIGVALS,[],ANY)$ +DEFINE_VARIABLE(LEFTMATRIX,'LEFTMATRIX,ANY)$ +DEFINE_VARIABLE(RIGHTMATRIX,'RIGHTMATRIX,ANY)$ + + +define_variable(rankpro,'rankpro,any)$ + +CONJUGATE(X):=SUBLIS('[%I=-%I],X)$ +CONJ(X):=CONJUGATE(X)$ + +INNERPRODUCT(X,Y):=BLOCK([LISTARITH:TRUE],RATSIMP(CONJUGATE(X).Y))$ +INPROD(X,Y):=INNERPRODUCT(X,Y)$ + +UVECT(X):=UNITVECTOR(X)$ + +UNITVECTOR(X):=BLOCK( + [LISTARITH:TRUE], + X/SQRT(INNERPRODUCT(X,X)))$ + +COVECT(X):=COLUMNVECTOR(X)$ + +COLUMNVECTOR(X):=TRANSPOSE(MATRIX(X))$ + +GSCHMIT(X):=GRAMSCHMIDT(X)$ + +GRAMSCHMIDT(X):=BLOCK( + [LISTARITH:TRUE,DIMNSN:LENGTH(X),LISTALL:[INPART(X,1)],INTERN,COUNT:1,DENOM,UNIT], + IF DIMNSN=1 + THEN X + ELSE (FOR INDEX1:2 THRU DIMNSN DO + (UNIT:INPART(X,INDEX1), + FOR INDEX2 THRU COUNT DO + (INTERN:INPART(LISTALL,INDEX2), + DENOM:INNERPRODUCT(INTERN,INTERN), + UNIT:FACTOR(RATSIMP(UNIT-INNERPRODUCT(INTERN,UNIT)*INTERN/DENOM))), + COUNT:COUNT+1, + LISTALL:ENDCONS(UNIT,LISTALL)), + LISTALL))$ + +EIVALS(MAT):=EIGENVALUES(MAT)$ + +EIGENVALUES(MAT):=BLOCK( + [LISTALL:[],SOLVEEXPLICIT:TRUE,DIMNSN:LENGTH(MAT),SOLUTION,MULTIPLICITIES, + DUMMY:?GENSYM()],if modulus=false then ( + SOLUTION:BLOCK([PROGRAMMODE:TRUE], + SOLVE(CHARPOLY(MAT,DUMMY),DUMMY)), + IF SOLUTION=[] + THEN (PRINT(" +SOLVE is unable to find the roots of +the characteristic polynomial."), + RETURN(LISTALL)) + ELSE (FOR INDEX2 THRU DIMNSN DO + (DIMNSN:DIMNSN-PART(MULTIPLICITIES,INDEX2)+1, + LISTALL:ENDCONS(RHS(PART(SOLUTION,INDEX2)),LISTALL)), + ENDCONS(MULTIPLICITIES,[LISTALL]))) + else eigenfinitefield(mat))$ + +EIVECTS(MAT):=EIGENVECTORS(MAT)$ + +EIGENVECTORS(MAT):=BLOCK( + [EQUATIONS,UNKNOWNS:[],SOLUTION,LISTALL,EIGVALS,DIMNSN:LENGTH(MAT), + sumofmultiplicities:0,uu,mult, + COUNT,VECTR,NOTKNWN,MATRX,MMATRX, + UNIT,MULTIPLICITIES,%RNUM,REALONLY:REALONLY,ALGEBRAIC:ALGEBRAIC,INTERM,INTERN], + DECLARE(LISTALL,SPECIAL), + LOCAL(NOTKNWN), + COUNT:DIMNSN, + EIGVALS:IF KNOWNEIGVALS + THEN LISTEIGVALS + ELSE EIGENVALUES(MAT), + IF EIGVALS=[] + THEN (NONDIAGONALIZABLE:TRUE, + RETURN(EIGVALS)) + ELSE (MULTIPLICITIES:INPART(EIGVALS,2),mult:multiplicities, + for uu:1 thru length(multiplicities) do( + sumofmultiplicities:first(mult)+sumofmultiplicities,mult:rest(mult)), + if modulus#false then count:sumofmultiplicities, + 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-INPART(MULTIPLICITIES,INDEX1)+1, + MMATRX:MATRX-INPART(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#INPART(MULTIPLICITIES,INDEX1) + THEN NONDIAGONALIZABLE:TRUE, + FOR INDEX3 THRU %RNUM DO + (INTERN:SUBSTVECTK(%RNUM_LIST,INDEX3,INTERM), + UNIT:APPEND(UNIT,[INTERN])), + IF UNIT=[] + THEN (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)), + LISTALL))$ + +/* 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))$ + +UEIVECTS(MAT):=UNITEIGENVECTORS(MAT)$ + +UNITEIGENVECTORS(MAT):=BLOCK( + [LISTUEVEC:IF KNOWNEIGVECTS + THEN LISTEIGVECTS + ELSE EIGENVECTORS(MAT),LISTALL,UNIT], + DECLARE([LISTALL,LISTUEVEC],SPECIAL), + IF LISTUEVEC=[] + THEN LISTUEVEC + ELSE (LISTALL:[INPART(LISTUEVEC,1)], + FOR INDEX1:2 THRU LENGTH(LISTUEVEC) DO + (UNIT:RATSIMP(UNITVECTOR(INPART(LISTUEVEC,INDEX1))), + LISTALL:ENDCONS(UNIT,LISTALL)), + LISTALL))$ + +SIMTRAN(MAT):=SIMILARITYTRANSFORM(MAT)$ + +SIMILARITYTRANSFORM(MAT):=BLOCK( + [LISTVEC,LISTUEVEC:UNITEIGENVECTORS(MAT)], + DECLARE(LISTUEVEC,SPECIAL), + IF NONDIAGONALIZABLE + THEN LISTUEVEC + ELSE (LISTVEC:DELETE(INPART(LISTUEVEC,1),LISTUEVEC), + RIGHTMATRIX:TRANSPOSE(APPLY(MATRIX,LISTVEC)), + LEFTMATRIX:IF HERMITIANMATRIX + THEN CONJUGATE(TRANSPOSE(RIGHTMATRIX)) + ELSE RIGHTMATRIX^^-1, + LISTUEVEC))$ + + +jrank(j,generalsimp):= + if length(j)<=r then nil else + if generalsimp(determinant(j))#0 + then r:max(r,length(j)) else + for i:1 thru length(j)do( + for k:1 thru length(j)do( + jrank(minor(j,i,k),generalsimp)))$ + +prank(j):=block([r:0], + if rankpro=false then jrank(j,lambda([x],fullratsimp(x))) else + jrank(j,rankpro),r)$ + +rankpro:false$ + +eigenfinite(j):= + Block([f,x,listofeigenvalues,dimsn], + listofeigenvalues:[], + dimsn:length(j), + define(f(x),fullratsimp(charpoly(j,x))), + for i : 1 thru modulus step 1 + do ( if fullratsimp(f(i))=0 then + listofeigenvalues: cons(i, listofeigenvalues)), + listofeigenvalues)$ + + +eigenfinitefield(zorro):= + block([partswitch:true,j,nicepoly,dummy:?gensym(),listofeigenvalues:[], + j:matrixmap(totaldisrep,zorro), + multiplicity:[],monopart], + nicepoly:factor(fullratsimp(charpoly(j,dummy))), + y:1, + if nicepoly=dummy then (listofeigenvalues:cons(0,listofeigenvalues),multiplicity:cons(1,multiplicity)) else ( + if inpart(nicepoly,0)=\^ and freeof(\+,nicepoly) and freeof(\-,nicepoly) + then (listofeigenvalues:cons(0,listofeigenvalues),multiplicity:cons(inpart(nicepoly,2),multiplicity)) + else ( + if not inpart(nicepoly,0)=\* then + (if hipow(nicepoly,dummy) = 1 then ( + if freeof(\^,nicepoly) then + (listofeigenvalues:cons(-inpart(nicepoly,1),listofeigenvalues), + multiplicity:cons(1,multiplicity)) /* multiplicity one */ + else (listofeigenvalues:cons(-inpart(nicepoly,1,1),listofeigenvalues) + ,multiplicity:cons(inpart(nicepoly ,2),multiplicity))) else nil) /* high multiplicity */ + else ( + while inpart(nicepoly,y) # end /* pick poly apart */ + do (monopart:inpart(nicepoly,y), /* monopart a piece */ + + if freeof(\+,monopart) and freeof(\-,monopart) and not freeof(dummy,monopart) then /* monopart zero root */ + (if freeof (\^,monopart) then (listofeigenvalues:cons(0,listofeigenvalues), + multiplicity:cons(1,multiplicity)) + else + (listofeigenvalues:cons(0,listofeigenvalues), /* single term case */ + multiplicity:cons(inpart(monopart,2),multiplicity))) + else + (if hipow(monopart,dummy) = 1 then ( + if freeof(\^,monopart) then /* multiplicity one */ + (listofeigenvalues:cons(-inpart(monopart,1),listofeigenvalues), + multiplicity:cons(1,multiplicity)) + else (listofeigenvalues:cons(-inpart(monopart,1,1),listofeigenvalues), /* high multiplicity */ + multiplicity:cons(inpart(monopart,2),multiplicity))) + + else nil), + y:y+1)))), + [listofeigenvalues,multiplicity])$ + + +jordanblock(x,j):= + block([dimsn,r,b,temp,jtemp,listofblock,c], + dimsn:length(j), + temp:ident(dimsn), + jtemp:j - x*temp, + for i:1 thru dimsn step 1 + do( + temp:temp . jtemp, + if not freeof(\%i,temp) then temp:radcan(temp), + r[i]:prank(temp)), + r[dimsn + 1]:r[dimsn], + listofblock:[dimsn - 2*r[1] + r[2]], + for k:2 thru dimsn step 1 + do(listofblock:cons(r[k+1] - 2*r[k] + r[k-1],listofblock)), + reverse(listofblock))$ + + +prettyjordan(jordanlist,mat):= block([dimsn:length(mat),lamb,place,k, + index,m,errorlist:jordanlist], + index:0, + m:zeromatrix(dimsn,dimsn), + while jordanlist#[] do( + lamb:first(first(jordanlist)), + place:last(first(jordanlist)), + jordanlist:rest(jordanlist), + k:1, + while place#[] do( + for i:1 thru first(place) do( + for j:1 thru k-1 do( + setelmx(lamb,index+j,index+j,m), + setelmx(1,index+j,index+j+1,m)),setelmx(lamb,index+k,index+k,m), + index:index+k), + k:k+1, + place:rest(place))), + if index#dimsn then m:[errorlist,"field not algebraically closed"] + else [],m)$ + +jordanform(j):= + block([grandlist,listofmultiplicity,listofeigenvalues,dimsn,listoflistofblocks,c], + grandlist:eigenvalues(j), + listofeigenvalues: first(grandlist), + listofmultiplicity:last(grandlist), + dimsn:length(j), + listoflistofblocks:[], + c:length(listofeigenvalues), + for i:1 thru c step 1 + do(if first(listofmultiplicity)#1 then + (listoflistofblocks:cons( + [first(listofeigenvalues),jordanblock(first(listofeigenvalues),j)], + listoflistofblocks), + listofeigenvalues:rest(listofeigenvalues), + listofmultiplicity:rest(listofmultiplicity)) + else + (listoflistofblocks:cons( + [first(listofeigenvalues),[1]], + listoflistofblocks), + listofeigenvalues:rest(listofeigenvalues), + listofmultiplicity:rest(listofmultiplicity))), + prettyjordan(listoflistofblocks,j))$ + +EVAL_WHEN(BATCH,TTYOFF:FALSE)$