1
0
mirror of https://github.com/PDP-10/its.git synced 2026-02-23 15:52:10 +00:00

Resolves #1012: Add SHARE;EIGEN DEMO and SHARE;EIGEN USAGE.

This commit is contained in:
Eric Swenson
2018-06-22 12:48:55 -07:00
parent 15f4c90806
commit 1d48b73016
2 changed files with 650 additions and 0 deletions

404
src/share/eigen.demo Executable file
View File

@@ -0,0 +1,404 @@
/* THIS IS THE FILE EIGEN DEMO DSK:SHARE;.
(YOU CAN BATCH OR DEMO THIS FILE, I.E. BATCH(EIGEN,DEMO,DSK,SHARE);, OR
DEMO(EIGEN,DEMO,DSK,SHARE);. NOTE THAT IN THE DEMO MODE YOU HAVE TO HIT THE
SPACE KEY AFTER EACH STEP...)
THE FUNCTIONS IN THE NEW EIGEN PACKAGE ARE DEMONSTRATED HERE. THE DESCRIPTION
OF THE FUNCTIONS CAN BE FOUND IN THE FILE EIGEN USAGE DSK:SHARE;, THE
SOURCE CODE IS ON THE FILE EIGEN > DSK:SHARE; AND THE FASTLOAD FILE IS
EIGEN FASL DSK:SHARE;. ( YOU CAN LOAD THIS ONE USING MACSYMA'S LOADFILE
COMMAND, I.E. LOADFILE(EIGEN,FASL,DSK,SHARE);.)
WE START WITH LOADING THE EIGEN PACKAGE : */
IF NOT STATUS(FEATURE,EIGEN) THEN LOADFILE(EIGEN,FASL,DSK,SHARE);
/* LET US START WITH THE FIRST FUNCTION. (SEE THE DESCRIPTIONS...)
FIRST LET'S DEFINE A COMPLEX VARIABLE... */
Z:A+%I*B;
/* THE CONJUGATE FUNCTION SIMPLY RETURNS THE COMPLEX CONJUGATE OF ITS
ARGUMENT... */
CONJUGATE(Z);
/* NOTE THAT Z COULD BE A MATRIX, A LIST, ETC... */
Z:MATRIX([%I,0],[0,1+%I]);
CONJUGATE(Z);
/* THE NEXT FUNCTION CALCULATES THE INNER PRODUCT OF TWO LISTS...*/
LIST1:[A,B,C,D];
LIST2:[F,G,H,I];
INNERPRODUCT(LIST1,LIST2);
/* THE ELEMENTS OF THE LISTS COULD BE COMPLEX ALSO... */
LIST1:[A+%I*B,C+%I*D];
INNERPRODUCT(LIST1,LIST1);
/* THE NEXT FUNCTION TAKES A LIST AS ITS ARGUMENT AND RETURNS A UNIT
LIST... */
LIST1:[1,1,1,1,1];
UNITVECTOR(LIST1);
LIST2:[1,%I,1-%I,1+%I];
UNITVECTOR(LIST2);
/* THE NEXT FUNCTION TAKES A LIST AS ITS ARGUMENT AND RETURNS A COLUMN
VECTOR... */
LIST1:[A,B,C,D];
COLUMNVECTOR(LIST1);
/* THE NEXT FUNCTION TAKES A LIST OF LISTS AS ITS ARGUMENT AND
ORTHOGONALIZES THEM USING THE GRAM-SCHMIDT ALGORITHM...*/
LISTOFLISTS:[[1,2,3,4],[0,5,4,7],[4,5,6,7],[0,0,1,0]];
/* NOTE THAT THE LISTS IN THIS LIST OF LISTS ARE NOT ORTHOGONAL TO EACH
OTHER... */
INNERPRODUCT([1,2,3,4],[0,5,4,7]);
INNERPRODUCT([1,2,3,4],[4,5,6,7]);
/* BUT AFTER APPLYING THE GRAMSCHMIDT FUNCTION... */
ORTHOGONALLISTS:GRAMSCHMIDT(LISTOFLISTS);
INNERPRODUCT(PART(ORTHOGONALLISTS,1),PART(ORTHOGONALLISTS,2));
INNERPRODUCT(PART(ORTHOGONALLISTS,2),PART(ORTHOGONALLISTS,3));
/* NOTE THAT ORHTOGONALLISTS CONTAINS INTEGERS THAT ARE FACTORED.
IF YOU DO NOT LIKE THIS FORM, YOU CAN SIMPLY RATSIMP THE RESULT : */
RATSIMP(ORTHOGONALLISTS);
/* THE NEXT FUNCTION TAKES A MATRIX AS ITS ARGUMENT AND RETURNS THE
EIGENVALUES OF THAT MATRIX... */
MATRIX1:MATRIX([M1,0,0,0,M5],[0,M2,0,0,M5],[0,0,M3,0,M5],[0,0,0,M4,M5],[0,0,0,0,0]);
/* THIS IS THE MATRIX THAT CAUSED A LOT OF TROUBLE FOR THE OLD EIGEN
PACKAGE... IT TOOK ~170 SECONDS TO FIND THE EIGEN VECTORS OF THIS
MATRIX... YOU SHOULD BE ABLE TO DO IT IN YOUR HEAD IN ABOUT 20 SECONDS
... THE NEW EIGEN PACKAGE HANDLES IT IN ABOUT 10 SECONDS... ANYWAY,
LET'S KEEP GOING... */
EIGENVALUES(MATRIX1);
/* THE FIRST SUBLIST IN THE ANSWER IS THE EIGENVALUES, SECOND LIST IS
THEIR MULTIPLICITIES IN THE CORRESPONDING ORDER...
THE NEXT FUNCTION TAKES A MATRIX AS ITS ARGUMENT AND RETURNS THE
EIGEN VALUES AND THE EIGEN VECTORS OF THAT MATRIX... */
EIGENVECTORS(MATRIX1);
/* FIRST SUBLIST IN THE ANSWER IS THE OUTPUT OF THE EIGENVALUES COMMAND
THE OTHERS ARE THE EIGEN VECTORS CORRESPONDING TO THOSE EIGEN VALUES...
NOTICE THAT THIS COMMAND IS MORE POWERFUL THAN THE EIGENVALUES COMMAND
BECAUSE IT DETERMINES BOTH THE EIGEN VALUES AND THE EIGEN VECTORS...
IF YOU ALREADY KNOW THE EIGEN VALUES, YOU CAN SET THE KNOWNEIGVALS FLAG
TO TRUE AND THE GLOBAL VARIABLE LISTEIGVALS TO THE LIST OF EIGEN
VALUES... THIS WILL MAKE THE EXECUTION OF EIGENVECTORS COMMAND FASTER
BECAUSE IT DOESN'T HAVE TO FIND THE EIGEN VALUES ITSELF... */
MATRIX2:MATRIX([1,2,3,4],[0,3,4,5],[0,0,5,6],[0,0,0,9]);
/* THE NEXT FUNCTION TAKES A MATRIX AS ITS ARGUMENT AND RETURNS THE
EIGENVALUES AND THE UNIT EIGEN VECTORS OF THAT MATRIX... */
UNITEIGENVECTORS(MATRIX2);
/* IF YOU ALREADY KNOW THE EIGENVECTORS YOU CAN SET THE FLAG
KNOWNEIGVECTS TO TRUE AND THE GLOBAL VARIABLE LISTEIGVECTS TO THE
LIST OF THE EIGEN VECTORS...
THE NEXT FUNCTION TAKES A MATRIX AS ITS ARGUMENT AND RETURNS THE EIGEN
VALUES AND THE UNIT EIGEN VECTORS OF THAT MATRIX. IN ADDITION IF
THE FLAG NONDIAGONALIZABLE IS FALSE,TWO GLOBAL MATRICES LEFTMATRIX AND
RIGHTMATRIX WILL BE GENERATED. THESE MATRICES HAVE THE PROPERTY THAT
LEFTMATRIX.(MATRIX).RIGHTMATRIX IS A DIAGONAL MATRIX WITH THE EIGEN
VALUES OF THE (MATRIX) ON THE DIAGONAL... */
SIMILARITYTRANSFORM(MATRIX1)$
NONDIAGONALIZABLE;
RATSIMP(LEFTMATRIX.MATRIX1.RIGHTMATRIX);
/* NOW THAT YOU KNOW HOW TO USE THE EIGEN PACKAGE, HERE ARE SOME
EXAMPLES ABOUT HOW NOT TO USE IT.
CONSIDER THE FOLLOWING MATRIX : */
MATRIX3:MATRIX([1,0],[0,1]);
/* AS YOU'VE UNDOUBTEDLY NOTICED, THIS IS THE 2*2 IDENTITY MATRIX.
LET'S FIND THE EIGEN VALUES AND THE EIGEN VECTORS OF THIS MATRIX...
*/
EIGENVECTORS(MATRIX3);
/* "NOTHING SPECIAL HAPPENED", YOU SAY. EVERYONE KNOWS WHAT THE EIGEN
VALUES AND THE EIGEN VECTORS OF THE IDENTITY MATRIX ARE, RIGHT?
RIGHT. NOW CONSIDER THE FOLLOWING MATRIX : */
MATRIX4:MATRIX([1,E],[E,1]);
/* LET E>0, BUT AS SMALL AS YOU CAN IMAGINE. SAY 10^(-100).
LET'S FIND THE EIGEN VALUES AND THE EIGEN VECTORS OF THIS MATRIX :
*/
EIGENVECTORS(MATRIX4);
/* SINCE E~10^(-100), THE EIGEN VALUES OF MATRIX4 ARE EQUAL TO THE
EIGEN VALUES OF MATRIX3 TO A VERY GOOD ACCURACY. BUT, LOOK
AT THE EIGEN VECTORS!!! EIGEN VECTORS OF MATRIX4 ARE NOWHERE
NEAR THE EIGEN VECTORS OF MATRIX3. THERE IS ANGLE OF %PI/4
BETWEEN THE CORRESPONDING EIGEN VECTORS. SO, ONE LEARNS
ANOTHER FACT OF LIFE :
MATRICES WHICH HAVE APPROXIMATELY THE SAME EIGEN VALUES DO NOT
HAVE APPROXIMATELY THE SAME EIGEN VECTORS IN GENERAL.
THIS EXAMPLE MAY SEEM ARTIFICIAL TO YOU, BUT IT IS NOT IF YOU THINK
A LITTLE BIT MORE ABOUT IT. SO, PLEASE BE CAREFUL WHEN YOU
APPROXIMATE THE ENTRIES OF WHATEVER MATRIX YOU HAVE. YOU MAY
GET GOOD APPROXIMATIONS TO ITS EIGEN VALUES, HOWEVER THE EIGEN
VECTORS YOU GET MAY BE ENTIRELY SPURIOUS( OR SOME MAY BE CORRECT,
BUT SOME OTHERS MAY BE TOTALLY WRONG ).
NOW, HERE IS ANOTHER SAD STORY :
LET'S TAKE A LOOK AT THE FOLLOWING MATRIX : */
MATRIX5:MATRIX([5/2,50-25*%I],[50+25*%I,2505/2]);
/* NICE LOOKING MATRIX, ISN'T IT? AS USUAL, WE WILL FIND THE EIGEN
VALUES AND THE EIGEN VECTORS OF IT : */
EIGENVECTORS(MATRIX5);
/* WELL, HERE THEY ARE. SUPPOSE THAT THIS WAS NOT WHAT YOU WANTED.
INSTEAD OF THOSE SQRT(70)'S, YOU WANT THE NUMERICAL VALUES OF
EVERYTHING. ONE WAY OF DOING THIS IS TO SET THE FLAG "NUMER"
TO TRUE AND USE THE EIGENVECTORS COMMAND AGAIN : */
NUMER:TRUE;
EIGENVECTORS(MATRIX5);
/* OOOPS!!! WHAT HAPPENED?? WE GOT THE EIGEN VALUES, BUT THERE ARE
NO EIGENVECTORS. NONSENSE, THERE MUST BE A BUG IN EIGEN, RIGHT?
WRONG. THERE IS NO BUG IN EIGEN. WE HAVE DONE SOMETHING WHICH
WE SHOULD NOT HAVE DONE. LET ME EXPLAIN :
WHEN ONE IS SOLVING FOR THE EIGEN VECTORS, ONE HAS TO FIND THE
SOLUTION TO HOMOGENEOUS EQUATIONS LIKE : */
EQUATION1:A*X+B*Y=0;
EQUATION2:C*X+D*Y=0;
/* IN ORDER FOR THIS SET OF EQUATIONS TO HAVE A SOLUTION OTHER THAN
THE TRIVIAL SOLUTION ( THE ONE IN WHICH X=0 AND Y=0 ), THE
DETERMINANT OF THE COEFFICIENTS ( IN THIS CASE A*D-B*C ) SHOULD
VANISH. EXACTLY. IF THE DETERMINANT DOES NOT VANISH THE ONLY
SOLUTION WILL BE THE TRIVIAL SOLUTION AND WE WILL GET NO EIGEN
VECTORS. DURING THIS DEMO, I DID NOT SET A,B,C,D TO ANY
PARTICULAR VALUES. LET'S SEE WHAT HAPPENS WHEN WE TRY TO SOLVE
THE SET ABOVE : */
ALGSYS([EQUATION1,EQUATION2],[X,Y]);
/* YOU SEE? THE INFAMOUS TRIVIAL SOLUTION. NOW LET ME SET A,B,C,D
TO SOME NUMERICAL VALUES : */
A:4;
B:6;
C:2;
D:3;
A*D-B*C;
EQUATION1:EV(EQUATION1);
EQUATION2:EV(EQUATION2);
ALGSYS([EQUATION1,EQUATION2],[X,Y]);
/* NOW WE HAVE A NONTRIVIAL SOLUTION WITH ONE ARBITRARY CONSTANT.
( %R(SOMETHING) ). WHAT HAPPENED IN THE PREVIOUS CASE IS THAT
THE NUMERICAL ERRORS CAUSED THE DETERMINANT NOT TO VANISH, HENCE
ALGSYS GAVE THE TRIVIAL SOLUTION AND WE GOT NO EIGEN VECTORS.
IF YOU WANT A NUMERICAL ANSWER, FIRST CALCULATE IT EXACTLY,
THEN SET "NUMER" TO TRUE AND EVALUATE THE ANSWER. */
NUMER:FALSE;
NOTNUMERICAL:EIGENVECTORS(MATRIX5);
NUMER:TRUE;
EV(NOTNUMERICAL);
/* YOU SEE, IT WORKS NOW. ACTUALLY, IF YOU HAVE A MATRIX WITH
NUMERICAL ENTRIES AND YOU CAN LIVE WITH REASONABLY ACCURATE
ANSWERS, THERE ARE MUCH BETTER (FASTER) PROGRAMS. ASK SOMEBODY
ABOUT THE IMSL ROUTINES ON THE SHARE DIRECTORY...
THIS IS ALL... IF YOU THINK THAT THE NAMES OF THE FUNCTIONS ARE TOO
LONG, THERE ARE SHORTER NAMES FOR THEM AND THEY ARE GIVEN IN THE FILE
EIGEN USAGE DSK:SHARE;. GOOD LUCK!!!!!!!!!!!!!...... YEKTA */
/* Part II. Matrices, over a finite field, Zp. If you are a physicist, or
engineer, your exposure to matrices has been limited to matrices
with entries in the reals or the complex domain. You probably know
that spiralling sinks are associated with complex eigenvalues with
length (called modulus) less than one. Spiralling sources similarly
are associated with complex eigenvalues with length >1.
Over a finite field, the geometric intuition is not there, but
what you think should happen does. Eigenvalues and eigenvectors
obey the same definitions. The only problem is that not all the
eigenvalues might lie in the field.
Finite fields can be found useful as a quantum model--
there are a finite number of energy states, and the energies
cycle back after an electromagnetic emission. This is only
one hypothetical model. I am sure that more physical models
could be found.
In any event, mathematically these are well-defined
quantities, which behave much better than floating point
realities. Either enjoy them as they are, or take it as
a challenge to find them a physical reality!!!
--ncs */
/* */
reset()$
/* A MATRIX THAT I WORKED WITH AND FOUND SOME NICE THEOREMS ABOUT IS
THE PASCAL MATRIX J.
HERE, I WORK WITH THE 3 X 3 MATRIX MODULO 3, AND THE 5 X 5 MATRIX
MODULO 5. */
m[i,j]:= binomial(i+j-2,j-1)$
pascal[k]:= genmatrix(m,k,k);
PASCAL[3];
jordanform(pascal[3]);
modulus:3$
jordanform(pascal[3]);
/* By introducing the modulus 3, the matrix pascal[3],
has a totally different Jordan form. Notice the off-diagonal
terms. */
/* NOW LET'S TRY PASCAL[5] WITH MODULUS 5. */
modulus:5$
PASCAL[5];
jordanform(pascal[5]);
/* THERE JUST AREN'T ENOUGH ROOTS IN MODULULO 5. IN A SUITABLE LARGER
FIELD ALL THE ROOTS WOULD LIE.
HERE WE COULD TAKE THAT FIELD TO BE THE REALS BY modulus:false$
THEN jordanform(pascal[5]); */
/* Knowing that the roots are 1, w, and w^2, the cube roots of unity, with */
/* multiplicities 1, 2, and 2, I use EIGENVECTOR and macsyma's */
/* ALGEBRAIC INTEGER capability */
tellrat(w^2+w+1)$
knowneigvals:true$
listeigvals:[[1,w,w^2],[1,2,2]]$
eigenvectors(PASCAL[5]);
/* Suppose we are interested in Finite Fermat Transforms (FFT).
In modulus 17, the second Fermat number, a fourth root is 4.
So for convolutions of length 4, the following transformation
matrix is important. */
kill(m)$
m[i,j]:= v^((i-1)*(j-1))$
/* HERE IS THE STRUCTURE OF THE FFT. NOTICE IT IS THE SAME FFT
AS COOLEY-TUKEY'S FFT OVER THE COMPLEX FIELD. A relevant
reference is E. Landau "NUMBER THEORY", especially Schur's
proof of Gauss's sum. */
genmatrix(m,4,4);
modulus:5$
kill(m)$
m[i,j]:=2^((i-1)*(j-1))$
a:matrixmap(fullratsimp,genmatrix(m,4,4));
jordanform(a);
/* MUCH MORE INFORMATION CAN BE OBTAINED IN THIS EXAMPLE BY USE OF
THE DOT OPERATOR. THE FOLLOWING IS a^2 */
matrixmap(fullratsimp,a.a);
/* another cute example. The Vandermonde matrix, evaluated at the
fourth root of unity mod 17, and dimensioned
as a 4 x 4 matrix */
kill(m)$
m[i,j]:= 4^(i-1)$
modulus:17$
genmatrix(m,4,4);
jordanform(%);
/* now isn't that unbelievable!!!! THERE ARE FOUR EIGENVALUES, ALL ZERO,
TWO JORDAN BLOCKS OF SIZE ONE,
AND ONE JORDAN BLOCK OF SIZE TWO. */
/* SO THIS IS ANOTHER INTERESTING EXAMPLE. */
kill(m)$
m[i,j]:= i+j$
modulus:5$
matrixmap(fullratsimp,genmatrix(m,4,4));
jordanform(%);
/* THE OLD RANK WASN'T A VERY VERSATILE COMMAND. IT COULD TELL THE RANK OF
MANY MATRICES, BUT WHOEVER WROTE
IT EXPECTED THAT THE CONSUMER WOULD ONLY INPUT INTEGERS. THUS RANK WOULD
SOMETIMES GIVE THE WRONG ANSWER...
TO REMEDY THIS ILL, I WROTE THE COMMAND PRANK WHICH USES THE SAME ALGORITHM
AS THE OLD RANK, HOWEVER
IT WRAPS THE ZERO-EQUIVALENCE STEP IN A LAMBDA-WRAPPER. THIS ALLOWS THE
CONSUMER FREEDOM OF CHOICE.*/
FOOBAR:MATRIX([SIN(X),1-COS(X)],[COS(X)+1,SIN(X)]);
/* THIS EXAMPLE DATES FROM 1982. RANK GIVES THE WRONG ANSWER AS CAN BE SEEN */
RANK(FOOBAR);
/* NOW LETS TRY pRANK!! */
RANKPRO:LAMBDA([X],FULLRATSIMP(TRIGREDUCE(X)))$
PRANK(FOOBAR);
/* SO BY TWIDDLING WITH rankpro YOU SHOULD HAVE COMPLETE FREEDOM OF CHOICE
AS TO YOUR SIMPLIFIERS.
HOWEVER, IT SHOULD BE MENTIONED THAT D. RICHARDSON PROVED IN 1968 THAT
ZERO-EQUIVALENCE IS
FORMALLY UNSOLVABLE. SEE "SOME UNSOLVABLE PROBLEMS INVOLVING ELEMENTARY
FUNCTIONS OF A REAL VARIABLE"
J. SYMBOLIC LOGIC,33,1968,PP.514-520 */
/* So that's all for now*/