mirror of
https://github.com/PDP-10/its.git
synced 2026-02-01 06:22:14 +00:00
Resolves #990: Added SHARE; EIGEN to the build.
Updates the way translated macysma source files is bulid to be more efficient and clean.
This commit is contained in:
Eric Swenson
@@ -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"
|
||||
|
||||
333
src/share/eigen.8
Executable file
333
src/share/eigen.8
Executable file
@@ -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)$
|
||||
Reference in New Issue
Block a user