reduce-historical/r38/packages/symmetry/symatvec.red

module symatvec;

% Symmetry 

% Author : Karin Gatermann
%         Konrad-Zuse-Zentrum fuer
%         Informationstechnik Berlin
%         Heilbronner Str. 10
%         W-1000 Berlin 31
%         Germany
%         Email: Gatermann@sc.ZIB-Berlin.de


% symatvec.red

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  functions for matrix vector operations
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

symbolic procedure gen!+can!+bas(dimension);
% returns the canonical basis of R^dimension as a vector list
begin
scalar eins,nullsq,i,j,ll;
   eins:=(1 ./ 1);
   nullsq:=(nil ./ 1);
   ll:= for i:=1:dimension collect
           for j:=1:dimension collect
              if i=j then eins else nullsq;
   return ll; 
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  matrix functions
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

symbolic procedure alg!+matrix!+p(mat1);
% returns true if the matrix is a matrix from algebraic level
begin
scalar len,elem;
  if length(mat1)<1 then rederr("should be a matrix");
  if not(car (mat1) = 'mat) then rederr("should be a matrix");
  mat1:=cdr mat1;
  if length(mat1)<1 then rederr("should be a matrix");
  len:=length(car mat1);
  for each elem in cdr mat1 do
    if not(length(elem)=len) then rederr("should be a matrix");
  return t;
end;

symbolic procedure matrix!+p(mat1);
% returns true if the matrix is a matrix in internal structure
begin
scalar dimension,z,res;
  if length(mat1)<1 then return nil;
  dimension:=length(car mat1);
  res:=t;
  for each z in cdr mat1 do
    if not(dimension = length(z)) then res:=nil;
  return res;
end;

symbolic procedure squared!+matrix!+p(mat1);
% returns true if the matrix is a matrix in internal structure
begin
  if (matrix!+p(mat1) and (get!+row!+nr(mat1) = get!+col!+nr(mat1)))
      then return t;
end;

symbolic procedure equal!+matrices!+p(mat1,mat2);
% returns true if the matrices are equal ( internal structure)
begin
scalar s,z,helpp,mathelp,sum,rulesum,rule1,rule2;
   if (same!+dim!+squared!+p(mat1,mat2)) then
       <<
           mathelp:=
             mk!+mat!+plus!+mat(mat1,
                 mk!+scal!+mult!+mat((-1 ./ 1),mat2));
           sum:=(nil ./ 1);
           for each z in mathelp do
                for each s in z do
                  if !*complex then
                     sum:=addsq(sum,multsq(s,mk!+conjugate!+sq s)) else
                     sum:=addsq(sum,multsq(s,s));
      %      print!-sq(sum);
      rulesum:=change!+sq!+to!+algnull(sum);
      if rulesum = 0 then helpp:=t else helpp:=nil;
 %     print!-sq(simp rulesum);
%           if null(numr(simp prepsq(sum))) then helpp:=t 
% else helpp:=nil;
       >> else helpp:=nil;
   return helpp;
end;

symbolic procedure get!+row!+nr(mat1);
% returns the number of rows
begin
   return length(mat1);
end;

symbolic procedure get!+col!+nr(mat1);
% returns the number of columns
begin
   return length(car mat1);
end;

symbolic procedure get!+mat!+entry(mat1,z,s);
% returns the matrix element in row z and column s
begin
   return nth(nth(mat1,z),s);
end;

symbolic procedure same!+dim!+squared!+p(mat1,mat2);
% returns true if the matrices are both squared matrices 
% of the same dimension
% (internal structur)
begin
  if (squared!+matrix!+p(mat1) and squared!+matrix!+p(mat2) and
       (get!+row!+nr(mat1) = get!+row!+nr(mat1))) 
          then return t;
end;

symbolic procedure mk!+transpose!+matrix(mat1);
% returns the transposed matrix (internal structure)
begin
scalar z,s,tpmat1;
  if not(matrix!+p(mat1)) then rederr("no matrix in transpose");
  tpmat1:=for z:=1:get!+col!+nr(mat1) collect
           for s:=1:get!+row!+nr(mat1) collect
                get!+mat!+entry(mat1,s,z);
  return tpmat1
end;

symbolic procedure mk!+conjugate!+matrix(mat1);
% returns the matrix with conjugate elements (internal structure)
begin
scalar z,s,tpmat1;
  if not(matrix!+p(mat1)) then rederr("no matrix in conjugate matrix");
  tpmat1:=for z:=1:get!+row!+nr(mat1) collect
           for s:=1:get!+col!+nr(mat1) collect
              mk!+conjugate!+sq(get!+mat!+entry(mat1,z,s));
  return tpmat1
end;

symbolic procedure mk!+hermitean!+matrix(mat1);
% returns the transposed matrix (internal structure)
begin
   if !*complex then
   return mk!+conjugate!+matrix(mk!+transpose!+matrix(mat1)) else
   return mk!+transpose!+matrix(mat1);
end;

symbolic procedure unitarian!+p(mat1);
% returns true if matrix is orthogonal or unitarian resp.
begin
scalar mathermit,unitmat1;
  mathermit:=mk!+mat!+mult!+mat(mk!+hermitean!+matrix(mat1),mat1);
  unitmat1:=mk!+unit!+mat(get!+row!+nr(mat1));
  if equal!+matrices!+p(mathermit,unitmat1) then return t;
end;

symbolic procedure mk!+mat!+mult!+mat(mat1,mat2);
% returns a matrix= matrix1*matrix2 (internal structure)
begin
scalar dims1,dimz1,dims2,s,z,res,sum,k;
  if not(matrix!+p(mat1)) then rederr("no matrix in mult");
  if not(matrix!+p(mat2)) then rederr("no matrix in mult");
  dims1:=get!+col!+nr(mat1);
  dimz1:=get!+row!+nr(mat1);
  dims2:=get!+col!+nr( mat2);
  if not(dims1 = get!+row!+nr(mat2)) then 
     rederr("matrices can not be multiplied");
  res:=for z:=1:dimz1 collect
         for s:=1:dims2 collect
           <<
              sum:=(nil ./ 1);
              for k:=1:dims1 do 
               sum:=addsq(sum,
                      multsq(
                       get!+mat!+entry(mat1,z,k),
                       get!+mat!+entry(mat2,k,s)
                             )
                          );
              sum:=subs2 sum where !*sub2=t;
              sum
           >>;
   return res;
end;

symbolic procedure mk!+mat!+plus!+mat(mat1,mat2);
% returns a matrix= matrix1 + matrix2 (internal structure)
begin
scalar dims,dimz,s,z,res,sum;
  if not(matrix!+p(mat1)) then rederr("no matrix in add");
  if not(matrix!+p(mat2)) then rederr("no matrix in add");
  dims:=get!+col!+nr(mat1);
  dimz:=get!+row!+nr(mat1);
  if not(dims = get!+col!+nr(mat2)) then
          rederr("wrong dimensions in add");
  if not(dimz = get!+row!+nr(mat2)) then
          rederr("wrong dimensions in add");
  res:=for z:=1:dimz collect
          for s:=1:dims collect
           <<
            sum:=addsq(
                   get!+mat!+entry(mat1,z,s),
                   get!+mat!+entry(mat2,z,s)
                 );
              sum:=subs2 sum where !*sub2=t;
              sum      
           >>;
  return res;
end;

symbolic procedure mk!+mat!*mat!*mat(mat1,mat2,mat3);
% returns a matrix= matrix1*matrix2*matrix3 (internal structure)
begin
scalar res;
  res:= mk!+mat!+mult!+mat(mat1,mat2);
  return mk!+mat!+mult!+mat(res,mat3);
end;

symbolic procedure add!+two!+mats(mat1,mat2);
% returns a matrix=( matrix1, matrix2 )(internal structure)
begin
scalar dimz,z,res;
  if not(matrix!+p(mat1)) then rederr("no matrix in add");
  if not(matrix!+p(mat2)) then rederr("no matrix in add");
  dimz:=get!+row!+nr(mat1);
  if not(dimz = get!+row!+nr(mat2)) then rederr("wrong dim in add");
  res:=for z:=1:dimz collect
      append(nth(mat1,z),nth(mat2,z));
  return res;
end;

symbolic procedure mk!+scal!+mult!+mat(scal1,mat1);
% returns a matrix= scalar*matrix (internal structure)
begin
scalar res,z,s,prod;
  if not(matrix!+p(mat1)) then rederr("no matrix in add");
  res:=for each z in mat1 collect
         for each s in z collect
           <<
              prod:=multsq(scal1,s);
              prod:=subs2 prod where !*sub2=t;
              prod
           >>;
  return res;
end;

symbolic procedure mk!+trace(mat1);
% returns the trace of the matrix (internal structure)
begin
scalar spurx,s;
  if not(squared!+matrix!+p(mat1)) then
          rederr("no square matrix in add");
  spurx :=(nil ./ 1);
  for s:=1:get!+row!+nr(mat1) do
     spurx :=addsq(spurx,get!+mat!+entry(mat1,s,s));
   spurx :=subs2 spurx where !*sub2=t;
  return spurx
end;

symbolic procedure mk!+block!+diagonal!+mat(mats);
% returns a blockdiagonal matrix from 
% a list of matrices (internal structure)
begin
  if length(mats)<1 then rederr("no list in mkdiagonalmats");
  if length(mats)=1 then return car mats else
     return fill!+zeros(car mats,mk!+block!+diagonal!+mat(cdr(mats)));
end;

symbolic procedure fill!+zeros(mat1,mat2);
% returns a blockdiagonal matrix from 2 matrices (internal structure)
begin
scalar nullmat1,nullmat2;
  nullmat1:=mk!+null!+mat(get!+row!+nr(mat2),get!+col!+nr(mat1));
  nullmat2:=mk!+null!+mat(get!+row!+nr(mat1),get!+col!+nr(mat2));
  return append(add!+two!+mats(mat1,nullmat2),
                    add!+two!+mats(nullmat1,mat2));
end;

symbolic procedure mk!+outer!+mat(innermat);
% returns a matrix for algebraic level
begin
 scalar res,s,z;
 if not(matrix!+p(innermat)) then rederr("no matrix in mkoutermat");
 res:= for each z in innermat collect
        for each s in z collect
            prepsq s;
 return append(list('mat),res);
end;

symbolic procedure mk!+inner!+mat(outermat);
% returns a matrix in internal structure
begin
   scalar res,s,z;
   res:= for each z in cdr outermat collect
          for each s in z collect
             simp s;
   if matrix!+p(res) then return res else 
        rederr("incorrect input in mkinnermat");
end;

symbolic procedure mk!+resimp!+mat(innermat);
% returns a matrix in internal structure
begin
   scalar res,s,z;
   res:= for each z in innermat collect
          for each s in z collect
             resimp s;
   return res;
end;

symbolic procedure mk!+null!+mat(dimz,dims);
% returns a matrix of zeros in internal structure
begin
scalar nullsq,s,z,res;
   nullsq:=(nil ./ 1);
   res:=for z:=1:dimz collect
           for s:=1:dims collect  nullsq;
  return res;
end;

symbolic procedure mk!+unit!+mat(dimension);
% returns a squared unit matrix in internal structure
begin
   return gen!+can!+bas(dimension);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  vector functions
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

symbolic procedure vector!+p(vector1);
% returns the length of a vector
% vector -- list of sqs
begin
  if length(vector1)>0 then return t;
end;

symbolic procedure get!+vec!+dim(vector1);
% returns the length of a vector
% vector -- list of sqs
begin
  return length(vector1);
end;

symbolic procedure get!+vec!+entry(vector1,elem);
% returns the length of a vector
% vector -- list of sqs
begin
  return nth(vector1,elem);
end;

symbolic procedure mk!+mat!+mult!+vec(mat1,vector1);
% returns a vector= matrix*vector (internal structure)
begin
scalar z;
  return for each z in mat1 collect
           mk!+real!+inner!+product(z,vector1);
end;

symbolic procedure mk!+scal!+mult!+vec(scal1,vector1);
% returns a vector= scalar*vector (internal structure)
begin
scalar entry,res,h;
  res:=for each entry in vector1 collect
     <<
        h:=multsq(scal1,entry);
        h:=subs2 h where !*sub2=t;
        h
     >>;
  return res;
end;

symbolic procedure mk!+vec!+add!+vec(vector1,vector2);
% returns a vector= vector1+vector2 (internal structure)
begin
scalar ent,res,h;
  res:=for ent:=1:get!+vec!+dim(vector1) collect
       <<
         h:= addsq(get!+vec!+entry(vector1,ent),
                get!+vec!+entry(vector2,ent));
         h:=subs2 h where !*sub2=t;
         h
       >>;
  return res;
end;

symbolic procedure mk!+squared!+norm(vector1);
% returns a scalar= sum vector_i^2 (internal structure)
begin
   return mk!+inner!+product(vector1,vector1);
end;

symbolic procedure my!+nullsq!+p(scal);
% returns true, if ths sq is zero
begin
   if null(numr( scal)) then return t;
end;

symbolic procedure mk!+null!+vec(dimen);
% returns a vector of zeros
begin
scalar nullsq,i,res;
   nullsq:=(nil ./ 1);
   res:=for i:=1:dimen collect nullsq;
   return res;
end;

symbolic procedure mk!+conjugate!+vec(vector1);
% returns a vector of zeros
begin
scalar z,res;
   res:=for each z in vector1  collect mk!+conjugate!+sq(z);
   return res;
end;

symbolic procedure null!+vec!+p(vector1);
% returns a true, if vector is the zero vector
begin
    if my!+nullsq!+p(mk!+squared!+norm(vector1)) then
       return t;
end;

symbolic procedure mk!+normalize!+vector(vector1);
% returns a normalized vector (internal structure)
begin
scalar scalo,vecres;
  scalo:=simp!* {'sqrt, mk!*sq(mk!+squared!+norm(vector1))};
  if my!+nullsq!+p(scalo) then 
     vecres:= mk!+null!+vec(get!+vec!+dim(vector1)) else
      <<
         scalo:=simp prepsq scalo;
         scalo:=quotsq((1 ./ 1),scalo);
         vecres:= mk!+scal!+mult!+vec(scalo,vector1);
      >>;
  return vecres;
end;

symbolic procedure mk!+inner!+product(vector1,vector2);
% returns the inner product of vector1 and vector2 (internal structure)
begin
scalar z,sum,vec2;
  if not(get!+vec!+dim(vector1) = get!+vec!+dim(vector2)) then
        rederr("wrong dimensions in innerproduct");
  sum:=(nil ./ 1);
  if !*complex then vec2:=mk!+conjugate!+vec(vector2) else
    vec2:=vector2;
  for z:=1:get!+vec!+dim(vector1) do 
      sum:=addsq(sum,multsq(
            get!+vec!+entry(vector1,z),
            get!+vec!+entry(vec2,z)
                           )
                );
  sum:=subs2 sum where !*sub2=t;
  return sum;
end;

symbolic procedure mk!+real!+inner!+product(vector1,vector2);
% returns the inner product of vector1 and vector2 (internal structure)
begin
scalar z,sum;
  if not(get!+vec!+dim(vector1) = get!+vec!+dim(vector2)) then
        rederr("wrong dimensions in innerproduct");
  sum:=(nil ./ 1);
  for z:=1:get!+vec!+dim(vector1) do 
      sum:=addsq(sum,multsq(
            get!+vec!+entry(vector1,z),
            get!+vec!+entry(vector2,z)
                           )
                );
  sum:=subs2 sum where !*sub2=t;
  return sum;
end;

symbolic procedure mk!+Gram!+Schmid(vectorlist,vector1);
% returns a vectorlist of orthonormal vectors
% assumptions: vectorlist is orthonormal basis, internal structure
begin
scalar i,orthovec,scalo,vectors1;
  orthovec:=vector1;
  for i:=1:(length(vectorlist)) do
     <<
       scalo:= negsq(mk!+inner!+product(orthovec,nth(vectorlist,i)));
       orthovec:=mk!+vec!+add!+vec(orthovec,
          mk!+scal!+mult!+vec(scalo,nth(vectorlist,i)));
     >>;
  orthovec:=mk!+normalize!+vector(orthovec);
  if null!+vec!+p(orthovec) then 
     vectors1:=vectorlist else
     vectors1:=add!+vector!+to!+list(orthovec,vectorlist);
  return vectors1
end;

symbolic procedure Gram!+Schmid(vectorlist);
% returns a vectorlist of orthonormal vectors
begin
scalar ortholist,i;
  if length(vectorlist)<1 then rederr("error in Gram Schmid");
  if vector!+p(car vectorlist) then    
      ortholist:=nil
        else rederr("strange in Gram-Schmid");
  for i:=1:length(vectorlist) do
        ortholist:=mk!+Gram!+Schmid(ortholist,nth(vectorlist,i));
  return ortholist;
end;

symbolic procedure add!+vector!+to!+list(vector1,vectorlist);
% returns a list of vectors consisting of vectorlist 
% and the vector1 at the end
% internal structure
begin
    return append(vectorlist,list(vector1));
end;

symbolic procedure mk!+internal!+mat(vectorlist);
% returns a matrix consisting of columns 
% equal to the vectors in vectorlist
% internal structure
begin
  return mk!+transpose!+matrix(vectorlist);
end;

symbolic procedure mat!+veclist(mat1);
% returns a vectorlist consisting of the columns of the matrix
% internal structure
begin
  return mk!+transpose!+matrix(mat1);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% some useful functions
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

symbolic procedure change!+sq!+to!+int(scal1);
% scal1 -- sq which is an integer
% result is a nonnegative integer
begin 
  scalar nr;
  nr:=simp!* prepsq scal1;
  if (denr(nr) = 1) then return numr(nr) else
    rederr("no integer in change!+sq!+to!+int");
end;

symbolic procedure change!+int!+to!+sq(scal1);
% scal1 --  integer for example 1 oder 2 oder 3
% result is a sq
begin 
  return (scal1 ./ 1);
end;

symbolic procedure change!+sq!+to!+algnull(scal1);
begin
scalar rulesum,storecomp;
           if !*complex then
              <<
                 storecomp:=t;
                 off complex;
              >> else
              <<
                 storecomp:=nil;
              >>;
          rulesum:=evalwhereexp ({'(list (list 
 (REPLACEBY 
   (COS (!~ X))
   (TIMES
      (QUOTIENT 1 2)
 (PLUS (EXPT E (TIMES I (!~ X))) (EXPT E (MINUS (TIMES I (!~ X))))) ))
 (REPLACEBY 
   (SIN (!~ X))
   (TIMES
      (QUOTIENT 1 (times 2 i))
 (difference (EXPT E (TIMES I (!~ X)))
      (EXPT E (MINUS (TIMES I (!~ X))))) ))

))
, prepsq(scal1)});
  rulesum:=reval rulesum;
  if storecomp then on complex;
 % print!-sq(simp (rulesum));
  return rulesum;
end;

symbolic procedure mk!+conjugate!+sq(mysq);
begin
    return conjsq(mysq);
 %   return subsq(mysq,'(( i . (minus i))));
end;

symbolic procedure mk!+equation(arg1,arg2);
begin
  return list('equal,arg1,arg2);
end;

symbolic procedure outer!+equation!+p(outerlist);
begin
    if eqcar(outerlist, 'equal) then return t
end;

symbolic procedure mk!+outer!+list(innerlist);
begin
  return append (list('list),innerlist)
end;

symbolic procedure mk!+inner!+list(outerlist);
begin
   if outer!+list!+p(outerlist) then return cdr outerlist;
end;

symbolic procedure outer!+list!+p(outerlist);
begin
  if eqcar(outerlist, 'list) then return t
end;

symbolic procedure equal!+lists!+p(ll1,ll2);
begin
  return (list!+in!+list!+p(ll1,ll2) and list!+in!+list!+p(ll2,ll1));
end;

symbolic procedure list!+in!+list!+p(ll1,ll2);
begin
  if length(ll1)=0 then return t else
       return (memq(car ll1,ll2) and list!+in!+list!+p(cdr ll1,ll2));
end;

symbolic procedure print!-matrix(mat1);
begin
  writepri (mkquote mk!+outer!+mat(mat1),'only);
end;

symbolic procedure print!-sq(mysq);
begin
  writepri (mkquote prepsq(mysq),'only);
end;

endmodule;

end;