Nice Integer Eigenvectors

[Mmmm]

This program gives you nice integer eigenvectors for 2x2 and 3x3 matrices rather than the normalised ones given by the EIGENVV program (or the seemingly random mess that the HP39gii gives out).

EIGEN() will output a list of a matrix of column eigenvectors and its diagonalisation if possible or a vector of eigenvalues if not possible.
EIGVALS() outputs a vector of eigenvalues
EIGVECT() outputs a single eigenvector matrix (vectors as columns).

it only works for matrices with real integer eigenvalues. If it can't solve the problem it falls back to EIGENVV (ie if it is non-integer, complex or not 2x2 or 3x3).

eg: Comparison of output: EIGENVV() Vs EIGEN()
   

It will also solve some problems that EIGENVV fails on.
eg: EIGENVV() and EIGENVAL() Vs EIGEN() and EIGVALS()
   

Code:

EIGVALS(matrix); //declare_EIGVALS_here_to_preserve
                            //function_order_in_user_menu

SIMPLIFY(vector)
BEGIN
LOCAL p,size,elements,factor;

size:=SIZE(vector);
elements:=size(1);

//find common factor of vector elements
factor:=vector(1);
FOR p FROM 2 TO elements DO
  factor:=igcd(vector(p),factor);
END;

//simplify
IF factor≠0 THEN
  vector:=vector/factor;
END;

//make first element non-zero
FOR p FROM 1 TO elements DO
  IF vector(p)≠0 THEN
    IF vector(p)≠ABS(vector(p)) THEN
      vector:=-1*vector;       
    END;
    BREAK;
  END;
END;
RETURN(vector);
END;

EXPORT EIGEN(matrix)
BEGIN
LOCAL m,vals,n,vect,eigmat,row,dim;
LOCAL col,nonparallel,valsmat,size;

//find dimension of matrix
size:=SIZE(matrix);
dim:=size(1);
IF size(2)≠dim THEN
dim:=0;
END;
  
//initialise vars
eigmat:=IDENMAT(dim);
vect:=MAKEMAT(1,dim);
row:=0;

//populate eigenvalue vector
vals:=EIGVALS(matrix);

//check matrix is compatible with EIGEN
//if not run EIGENVV instead
IF dim<2 OR dim>3 THEN //use EIGENVV instead
  //leave row==0
ELSE
  FOR row FROM 1 TO dim DO
    IF IM(vals(row))≠0 THEN
      BREAK;
    END;
    IF FP(vals(row))≠0 THEN
      BREAK;
    END;
  END;
END;
IF row<(dim+1) THEN
  RETURN(EIGENVV(matrix));
END;

//main loop - find each eigenvector
FOR n FROM 1 TO dim DO
  //find characteristic matrix
  m:=(matrix-vals(n)*IDENMAT(dim));

  //Simplify rows of m
  FOR row FROM 1 TO dim DO
    m(row):=SIMPLIFY(m(row));
  END;

  IF RANK(m)==0 THEN //matrix == kI so eigenvectors=I
    eigmat:=IDENMAT(dim);
    BREAK;
  END;

  //find first non-zero element
  FOR row FROM 1 TO dim DO
    FOR col FROM 1 TO dim DO
      IF m(row,col)≠0 THEN
        BREAK;
      END;
    END;
    IF col<(dim+1) THEN
      BREAK;
    END;    
  END;
  
  IF dim==2 THEN //2x2 eigenvector calculation
    vect(1):=-m(row,2);
    vect(2):=m(row,1);
  ELSE //3x3 eigenvector calculations
    IF RANK(m)==1 THEN //only 1 li row - find a vector perp to m
      nonparallel:=m(row);
      nonparallel(1+(col MOD dim)):=1+m(row,(1+(col MOD dim)));
      vect:=CROSS(m(row),nonparallel);
    ELSE //2 li rows - use cross product
      IF m(row)≠m(row+1) AND m(row+1)≠[0,0,0] THEN
        vect:=CROSS(m(row),m(row+1));
      ELSE
        vect:=CROSS(m(row),m(row+2));
      END;
    END;
  END;
  vect:=SIMPLIFY(vect);
  
  // in case of repeated eigenvalue find vector perp
  // to charmatrix and repeated eigenvector,
  // but only if we have 3 non li row vects in m and dim>2   
  IF n>2 AND dim>2 THEN
    IF vect==eigmat(n-2) AND RANK(m)==1 THEN //find a perp vector
      vect:=CROSS(vect,m(row));
      vect:=SIMPLIFY(vect);
    END;
  END; 
  IF n>1 AND dim>2 THEN
    IF vect==eigmat(n-1) AND RANK(m)==1 THEN
      vect:=CROSS(vect,m(row));
      vect:=SIMPLIFY(vect);
    END;
  END;

//populate eigenvector matrix (rows)
  eigmat(n):=vect;
END;

//zero repeated eigenvectors
FOR n FROM 1 TO dim DO
  IF eigmat(n)==eigmat(1+(n MOD dim)) THEN
    FOR col FROM 1 TO dim DO
      eigmat(MAX(n,(1+(n MOD dim))),col):=0;
    END;
  END;
END;

//change to a column matrix
eigmat:=TRN(eigmat);

//display result
IF DET(eigmat)≠0 THEN
valsmat:=eigmat*matrix*eigmat;
  //diagonalise and show as a list
CONCAT(eigmat,valsmat);
ELSE
   CONCAT(eigmat,vals);
END;

END;

EXPORT EIGVECT(matrix)
BEGIN
//brief function to return the vector matrix only
LOCAL result;
result:=EIGEN(matrix);
result(1);
END;

EXPORT EIGVALS(matrix)
BEGIN
LOCAL a,b,c,d,eqn,size,dim;

//determine the dimensions of the matrix
size:=SIZE(matrix);
dim:=size(1);
IF size(2)≠dim THEN
dim:=0;
END;

//use built in eigenvalue command if matrix not compatible
IF dim<2 OR dim>3 THEN
  RETURN(EIGENVAL(matrix));
END;

IF dim==2 THEN
  //Generate coefficients for characteristic equation (quadratic eqn)
  a:=1;
  b:=-(matrix(1,1)+matrix(2,2));
  c:=DET(matrix);

  //form quadratic coefficient vector
  eqn:=[0,0,0];
  eqn(1):=ROUND(a,8);
  eqn(2):=ROUND(b,8);
  eqn(3):=ROUND(c,8);

ELSE //(if dim==3)
  //Generate coefficients for characteristic equation (cubic eqn)
  a:=-1;
  b:=TRACE(matrix);
  c:=-matrix(1,1)*matrix(3,3)-matrix(1,1)*matrix(2,2)-matrix(2,2)*matrix(3,3)+matrix(1,3)*matrix(3,1)+matrix(1,2)*matrix(2,1)+matrix(2​,3)*matrix(3,2);
  d:=matrix(1,1)*matrix(2,2)*matrix(3,3)+matrix(1,2)*matrix(2,3)*matrix(3,1)+​matrix(1,3)*matrix(2,1)*matrix(3,2)-matrix(1,1)*matrix(2,3)*matrix(3,2)-matrix(3,3)*matrix(1,2)*matrix(2,1)-matrix(2,2)*matrix(1,3)*matrix(3,1);

  //form cubic coefficient vector
  eqn:=[0,0,0,0];
  eqn(1):=ROUND(a,8);
  eqn(2):=ROUND(b,8);
  eqn(3):=ROUND(c,8);
  eqn(4):=ROUND(d,8);
END;
//solve polynomial
ROUND(POLYROOT(eqn),8);
END;