reduce-historical/r38/packages/fide/charpol.red

module charpol;

% Author: R. Liska.

% Version REDUCE 3.6     05/1991.

fluid '(!*exp !*gcd !*prfourmat)$

switch prfourmat$
!*prfourmat:=t$

procedure coefc1 uu$
begin
  scalar lco,l,u,v,a$
  u:=car uu$
  v:=cadr uu$
  a:=caddr uu$
  lco:=aeval list('coeff,u,v)$
  lco:=cdr lco$
  l:=length lco - 1$
  for i:=0:l do
    <<setel(list(a,i),car lco)$
      lco:=cdr lco >>$
  return (l . 1)
end$

deflist('((coefc1 coefc1)),'simpfn)$

global '(cursym!* coords!* icoords!* unvars!*)$

icoords!*:='(i j k l m n i1 j1 k1 l1 m1 n1)$

flag('(tcon unit charmat ampmat denotemat),'matflg)$

put('unit,'rtypefn,'getrtypecar)$
put('charmat,'rtypefn,'getrtypecar)$
put('ampmat,'rtypefn,'getrtypecar)$
put('denotemat,'rtypefn,'getrtypecar)$

procedure unit u$
generateident length matsm u$

procedure charmat u$
matsm list('difference,list('times,'lam,list('unit,u)),u)$

procedure charpol u$
begin
  scalar x,complexx;
  complexx:=!*complex;
  algebraic on complex;
  x:=simp list('det,list('charmat,carx(u,'charpol)))$
  if null complexx then algebraic off complex;
  return x
end;

put('charpol,'simpfn,'charpol)$

algebraic$
korder lam$

procedure re(x)$
sub(i=0,x)$

procedure im(x)$
(x-re(x))/i$

procedure con(x)$
sub(i=-i,x)$

procedure complexpol x$
begin
  scalar y$
  y:=troot1 x$
  return if im y=0 then y
           else y*con y
end$

procedure troot1 x$
begin
  scalar y$
  y:=x$
  while not(sub(lam=0,y)=y) and sub(lam=1,y)=0 do y:=y/(lam-1)$
  x:=sub(lam=1,y)$
  if not(numberp y or (numberp num y and numberp den y)) then
      write " If  ",re x,"  =  0  and  ",im x,
            "  =  0  , a root of the polynomial is equal to 1."$
  return y
end$

procedure hurw(x)$
% X is a polynomial in LAM, all its roots are |LAMI|<1 <=> for all roots
% of the polynomial HURW(X) holds RE(LAMI)<0.
(lam-1)**deg(num x,lam)*sub(lam=(lam+1)/(lam-1),x)$

symbolic$

procedure unfunc u$
<<unvars!*:=u$
  for each a in u do put(a,'simpfn,'simpiden) >>$

put('unfunc,'stat,'rlis)$

global '(denotation!* denotid!*)$
denotation!*:=nil$
denotid!*:='a$

procedure denotid u$
<<denotid!*:=car u$nil>>$

put('denotid,'stat,'rlis)$

procedure cleardenot$
denotation!*:=nil$

put('cleardenot,'stat,'endstat)$
flag('(cleardenot),'eval)$

algebraic$
array cofpol!*(20)$

procedure denotepol u$
begin
  scalar nco,dco$
  dco:=den u$
  u:=num u$
  nco:=coefc1 (u,lam,cofpol!*)$
  for j:=0:nco do cofpol!*(j):=cofpol!*(j)/dco$
  denotear nco$
  u:=for j:=0:nco sum lam**j*cofpol!*(j)$
  return u
end$

symbolic$

put('denotear,'simpfn,'denotear)$

procedure denotear u$
begin
  scalar nco,x$
  nco:=car u$
  for i:=0:nco do
    <<x:=list('cofpol!*,i)$
      setel(x,mk!*sq denote(getel x,0,i)) >>$
  return (nil .1)
end$

procedure denotemat u$
begin
  scalar i,j,x$
  i:=0$
  x:=for each a in  matsm u collect
       <<i:=i+1$
         j:=0$
         for each b in a collect
           <<j:=j+1$
             denote(mk!*sq b,i,j) >> >>$
  return x
end$

procedure denote(u,i,j)$
% U is prefix form, I,J are integers
begin
  scalar reu,imu,ireu,iimu,eij,fgcd$
  if atom u then return simp u$
  fgcd:=!*gcd$
  !*gcd:=t$
  reu:=simp!* list('re,u)$
  imu:=simp!* list('im,u)$
  !*gcd:=fgcd$
  eij:=append(explode i,explode j)$
  ireu:=intern compress append(append(explode denotid!* ,'(r)),eij)$
  iimu:=intern compress append(append(explode denotid!* ,'(i)),eij)$
  if car reu then insdenot(ireu,reu)$
  if car imu then insdenot(iimu,imu)$
  return simp list('plus,
                   if car reu then ireu
                     else 0,
                   list('times,
                        'i,
                        if car imu then iimu
                          else 0))
end$

procedure insdenot(iden,u)$
denotation!*:=(u . iden) . denotation!*$

procedure prdenot$
for each a in reverse denotation!* do
  assgnpri(list('!*sq,car a,t),list cdr a,'only)$

put('prdenot,'stat,'endstat)$
flag('(prdenot),'eval)$

procedure ampmat u$
begin
  scalar x,i,h1,h0,un,rh1,rh0,ru,ph1,ph0,!*exp,!*gcd,complexx$
  complexx:=!*complex;
  !*exp:=t$
  fouriersubs()$
  u:=car matsm u$
  x:=for each a in coords!* collect if a='t then 0
                                      else list('times,
                                                tcar get(a,'index),
                                                get(a,'wave),
                                                get(a,'step))$
  x:=list('expp,'plus . x)$
  x:=simp x$
  u:=for each a in u collect resimp quotsq(a,x)$
  gonsubs()$
  algebraic on complex;
  u:=for each a in u collect resimp a$
  remfourier()$
a:if null u then go to d$
  ru:=caar u$
  un:=unvars!*$
  i:=1$
b:if un then go to c$
  rh1:=reverse rh1$
  rh0:=reverse rh0$
  h1:=rh1 . h1$
  h0:=rh0 . h0$
  rh0:=rh1:=nil$
  u:=cdr u$
  go to a$
c:rh1:=coefck(ru,list('u1!*,i)) . rh1$
  rh0:=negsq coefck(ru,list('u0!*,i)) . rh0$
  un:=cdr un$
  i:=i+1$
  go to b$
d:h1:=reverse h1$
  h0:=reverse h0$
  if !*prfourmat then
      <<ph1:=for each a in h1 collect
               for each b in a collect prepsq!* b$
        setmatpri('h1,nil . ph1)$
        ph1:=nil$
        ph0:=for each a in h0 collect
               for each b in a collect prepsq!* b$
        setmatpri('h0,nil . ph0)$
        ph0:=nil >>$
  !*gcd:=t;
  x:=if length h1=1 then list list quotsq(caar h0,caar h1)
       else lnrsolve(h1,h0)$
  if null complexx then algebraic off complex;
  return x
end$

procedure coefck(x,y)$
% X is standard form, Y is prefix form, returns coefficient of Y
% appearing in X, i.e. X contains COEFCK(X,Y)*Y
begin
  scalar ky,xs$
  ky:=!*a2k y$
  xs:=car subf(x,list(ky . 0))$
  xs:=addf(x,negf xs)$
  if null xs then return(nil . 1)$
  xs:=quotf1(xs,!*k2f ky)$
  return if null xs then <<msgpri
    ("COEFCK failed on ",list('sq!*,x . 1,nil)," with ",y,'hold)$
                           xread nil$
                           !*f2q xs>>
           else !*f2q xs
end$

procedure simpfour u$
begin
  scalar nrunv,x,ex,arg,mv,cor,incr,lcor$
  nrunv:=get(car u,'nrunvar)$
a:u:=cdr u$
  if null u then go to r$
  arg:=simp car u$
  mv:=mvar car arg$
  if not atom mv or not numberp cdr arg then return msgpri
      ("Bad index ",car u,nil,nil,'hold)$
  cor:=tcar get(mv,'coord)$
  if not(cor member coords!*) then return msgpri
      ("Term ",car u," contains non-coordinate ",mv,'hold)$
  if cor member lcor then return msgpri
      ("Term ",car u," means second appearance of coordinate ",cor,
       'hold)$
  if not(cor='t) and cdr arg>get(cor,'maxden)
    then put(cor,'maxden,cdr arg)$
  lcor:=cor . lcor$
  incr:=addsq(arg,negsq !*k2q mv)$
  if not flagp(cor,'uniform) then return lprie
      ("Non-uniform grids not yet supported")$
  if cor='t then go to ti$
  ex:=list('times,car u,get(cor,'step),get(cor,'wave)) . ex$
  go to a$
ti:if null car incr then x:=list('u0!*,nrunv)
    else if incr= 1 . 1 then x:=list('u1!*,nrunv)
    else return lprie "Scheme is not twostep in time"$
  go to a$
r:for each a in setdiff(coords!*,lcor) do
    if a='t then x:=list('u0!*,nrunv)
      else ex:=list('times,tcar get(a,'index),get(a,'step),get(a,'wave))
                  . ex$
  return simp list('times,x,list('expp,'plus . ex))
end$

procedure fouriersubs$
begin
  scalar x,i$
  for each a in '(expp u1!* u0!*) do put(a,'simpfn,'simpiden)$
  x:=unvars!*$
  i:=1$
a:if null x then go to b$
  put(car x,'nrunvar,i)$
  i:=i+1$
  x:=cdr x$
  go to a$
b:flag(unvars!*,'full)$
  for each a in unvars!* do put(a,'simpfn,'simpfour)$
  for each a in coords!* do
    if not(a='t) then
    <<x:=intern compress append(explode 'h,explode a)$
      put(a,'step,x)$
      if not flagp(a,'uniform) then put(x,'simpfn,'simpiden)$
      x:=intern compress append(explode 'k,explode a)$
      put(a,'wave,x)$
      x:=intern compress append(explode 'a,explode a)$
      put(a,'angle,x)$
      put(a,'maxden,1) >>$
  algebraic for all z,y,v let
    expp((z+y)/v)=expp(z/v)*expp(y/v),
    expp(z+y)=expp z*expp y
end$

procedure gonsubs$
begin
  scalar xx$
  algebraic for all z,y,v clear expp((z+y)/v),expp(z+y)$
  for each a in coords!* do
    if not(a='t) then
      <<xx:=list('quotient,
                 list('times,
                      get(a,'maxden),
                      get(a,'angle)),
                 get(a,'step))$
        setk(get(a,'wave),aeval xx)$
        xx:=list('times,
                 get(a,'wave),
                 get(a,'step))$
        mathprint list('setq,
                       get(a,'angle),
                       if get(a,'maxden)=1 then xx
                         else list('quotient,
                                   xx,
                                   get(a,'maxden))) >>$
  algebraic for all x let expp x=cos x+i*sin x$
  algebraic for all x,n such that numberp n and n>1 let
    sin(n*x)=sin x*cos((n-1)*x)+cos x*sin((n-1)*x),
    cos(n*x)=cos x*cos((n-1)*x)-sin x*sin((n-1)*x)$
  for each a in unvars!* do
    <<put(a,'simpfn,'simpiden)$
      remprop(a,'nrunvar) >>
end$

procedure remfourier$
<<algebraic for all x clear expp x$
  algebraic for all x,n such that numberp n and n>1 clear
    sin(n*x),cos(n*x)$
  for each a in coords!* do
    if not(a='t) then
    <<remprop(a,'step)$
      remprop(remprop(a,'wave),'avalue)$
      remprop(a,'maxden) >> >>$

operator numberp$

endmodule;

end;