Artifact 982e724f720c6105407baca5038749b4c41393e4d3a5f22815b15d9757fc63ff:

• Executable file r37/packages/fide/hurwp.red — part of check-in [f2fda60abd] at 2011-09-02 18:13:33 on branch master — Some historical releases purely for archival purposes

  git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 9924) [annotate] [blame] [check-ins using] [more...]

• Executable file r38/packages/fide/hurwp.red — part of check-in [f2fda60abd] at 2011-09-02 18:13:33 on branch master — Some historical releases purely for archival purposes

  git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 9924) [annotate] [blame] [check-ins using]

module hurwp;

% Author: R. Liska.

% Version REDUCE 3.6     05/1991

global '(ofl!* mlist!*)$
fluid '(!*exp !*gcd)$

flag('(tcon),'matflg)$
put('tcon,'msimpfn,'tcon)$
put('tcon,'rtypefn,'getrtypecar)$

procedure tcon u$
% Calculates complex conjugate and transpose matrix
begin
  scalar v,b$
  v:=matsm list('tp,u)$
  for each a in v do
    <<b:=a$
      while b do
        <<rplaca(b,quotsq(subf(numr car b, '((i minus i))),
                          !*f2q denr car b))$
          b:=cdr b >> >>$
  return v
end$

algebraic$
korder lam$

symbolic$

global '(positive!* userpos!* userneg!* !*pfactor)$
!*pfactor:=nil$

procedure positivep u$
% U is prefix form. Procedure tests if U>0, eventually writes this
% condition and puts U into POSITIVE!*. If U<=0 then returns NIL,
% if U>0 then T, in other cases 'COND.
% If it does not know if U>0 and program is running in interactive
% mode it asks user if U>0 and return value is based on user reply.
if numberp u then
    if u>0 then t
      else nil
  else if eqcar(u,'!*sq) and fixp caadr u and fixp cdadr u then
           if caadr u*cdadr u>0 then t
             else nil
  else
    begin
      scalar x,exp$
      exp:=!*exp$
      if !*pfactor and
             member('factor,mlist!*) then
           <<!*exp:=nil$
             u:=aeval list('ppfactor,u) >>$
      u:=prepsq!* simp u$
      !*exp:=exp$
      x:=if terminalp() and null ofl!* then
             begin
               scalar y,z$
               prin2!* "Is it true, that  "$
               maprin u$
               prin2!* "  >  0    ?"$
             a:prin2!* "  Reply (Y/N/?)"$
               terpri!* t$
               y:=read()$
               if y eq 'y then <<z:=t$ userpos!*:=u . userpos!* >>
                 else if y eq 'n
                  then <<z:=nil$ userneg!*:=u . userneg!*>>
                 else if y eq '? then z:='cond
                 else go to a$
               return z
             end
           else 'cond$
      if x eq 'cond then
          if null positive!* then positive!*:= list (1 . u)
            else positive!* := ((caar positive!* + 1) . u) . positive!*$
      return x
    end$

global'(hconds!*)$

algebraic$
array cof(20),fcof(20)$
share hconds!*$

procedure ppfactor x$
begin
  scalar d,n,de$
  d:= old_factorize(num x)$
  n:=for each a in d product a$
  if den x=1 then return n$
  d:= old_factorize(den x)$
  de:=for each a in d product a$
  return (n/de)
end$

procedure hurwitzp u$
% U is a polynomial in LAM. Procedure tests if it is Hurwitz polynomial
% i.e. for all its rools LAMI holds RE(LAMI)<0.
% Returned values: YES  - definitely yes
%                  NO   - definitely no
%                  COND - if conditions holds (all members of POSITIVE!*
%                         are >0)
if im u=0 then rehurwp u
  else cohurwp u$

symbolic$

procedure coef1(u,v,a)$
begin
  scalar lco,l$
  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
end$

procedure rehurwp u$
begin
  scalar deg,hurp,gcd$
  gcd:=!*gcd$
  !*gcd:=t$
  deg:=coef1(car u,'lam,'cof)$
  if deg=0 then return typerr(u,"Polynomial in LAM")$
  positive!* := userpos!* := userneg!* := nil$
  if deg <= 2 then
      <<for i:=0:deg do setel(list('cof,i),
                              aeval list('quotient,
                                         getel list('cof,i),
                                         getel list('cof,deg)))$
        if deg=1 then hurp:=positivep getel list('cof,0)
          else if deg=2 then hurp:=
                   if positivep getel list('cof,0) and
                      positivep getel list('cof,1) then
                        if positive!* then 'cond
                          else t
                      else if positive!* then 'cond    % added 08/08/91
                      else nil$
        printcond(nil) >>
    else hurp:=rehurwp1 deg$
  !*gcd:=gcd$
  return rethurp hurp
end$

procedure rethurp hurp$
<<hconds!*:= 'list . if positive!* then
                         for each a in positive!* collect cdr a
                       else nil$
  !*k2q(if null hurp then 'no
          else if null positive!* then 'yes
          else 'cond) >>$

put('rehurwp,'simpfn,'rehurwp)$

procedure cohurwp u$
begin
  scalar deg$
  u:=reval list('sub,'(equal lam (times i lam)),car u)$
  deg:=coef1(u,'lam,'cof)$
  if deg=0 then return typerr(u,"Polynomial in LAM")$
  positive!* := userpos!* := userneg!* :=nil$
  if aeval list('im,getel list('cof,deg))=0 then
        for j:= 0:deg do setel(list('cof,j),
                              aeval list('times,'i,getel list('cof,j)))$
  return rethurp cohurwp1 (deg)
end$

put('cohurwp,'simpfn,'cohurwp)$

procedure rehurwp1 deg$
begin
  scalar i,bai,bdi,x,lich,sud,bsud,matr,hmat,csud,clich,dsud,dlich$
a:i:=deg$
  csud:=clich:=nil$
  bsud:=t$
b:x:=positivep getel list('cof,i)$
  if null x then go to c
    else if x eq t then bai:=t
    else if x eq 'cond then
      if i=deg or i=0 then <<csud:=caar positive!* . csud$
                             clich:=caar positive!* . clich >>
        else if bsud then csud:=caar positive!* . csud
        else clich:=caar positive!* . clich$
  i:=i-1$
  bsud:=not bsud$
  if i>=0 then go to b$
  go to d$
  % Change of sign AI = - AI
c:if bai or bdi then go to n
    else bai:=t$
  for i:=0:deg do setel(list('cof,i),
                        aeval list('minus,getel list('cof,i)))$
  go to a$
  % Checking DI > 0 - Hurwitz determinants
  % Splitting to odd and even coeffs. AI, A0 is coeff. by LAM**DEG
d:bsud:=t$
  for i:=deg step -1 until 0 do
    <<x:=simp getel list('cof,i)$
      if bsud then sud:=x . sud
        else lich:=x . lich$
      bsud:=not bsud >>$
  sud:=reverse sud$
  lich:=reverse lich$
  % Filling of SUD and LICH on the length DEG by zeroes from right
  sud:=filzero(sud,deg)$
  lich:=filzero(lich,deg)$
  dsud:=dlich:=nil$
  matr:=nil$
  i:=1$
  bsud:=nil$
d1:matr:=nconc(matr,list lich)$
  lich:=(nil . 1) . lich$
d2:hmat:=cutmat(matr,i)$
  x:=mk!*sq detq hmat$
  x:=positivep x$  % Necessary to add storing of odd and even DIs
  if null x then
      if bsud then go to n
        else go to c
    else if x eq t and not bsud then bdi:=t
    else if x eq 'cond then
      if bsud then dsud:=caar positive!* . dsud
        else dlich:=caar positive!* . dlich$
  i:=i+1$
  bsud:=not bsud$
  if i>deg then go to k$
  if not bsud then go to d1$
  matr:=nconc(matr,list sud)$
  sud:=(nil . 1) . sud$
  go to d2$
n:return nil$
k:if null positive!* or ((null csud or null clich) and
                         (null dsud or null dlich))
      then return <<printuser()$ t>>$
  prin2t "If we denote:"$
  printcond(t)$
  printdef('c1,clich:=reverse clich)$
  printdef('c2,csud:=reverse csud)$
  printdef('d1,dlich:=reverse dlich)$
  printdef('d2,dsud:=reverse dsud)$
  prin2t "Necessary and sufficient conditions are:"$
  prin2t if null csud or null clich then         "  (D1)  OR  (D2)"
           else if null dsud or null dlich then  "  (C1)  OR  (C2)"
           else "  (  (C1)  OR  (C2)  )  AND  (  (D1)  OR  (D2)  )"$
  printuser()$
  return 'cond
end$

procedure printcond(x)$
<<if not x then
      prin2t "Necessary and sufficient conditions are:"$
  positive!*:=reverse positive!*$
  for each a in positive!* do
    <<if x then <<prin2!* " ("$
                  prin2!* car a$
                  prin2!* ")  " >>$
      maprin cdr a$
      prin2!* "  >  0"$
      terpri!* t >>$
  if not x then printuser() >>$

procedure printuser()$
if userpos!* or userneg!* then
    <<prin2t"You have explicitly stated:"$
      for each a in userpos!* do <<maprin a$
                                   prin2!* "  >  0"$
                                   terpri!* t >>$
      for each a in userneg!* do <<maprin a$
                                   prin2!* "  <=  0"$
                                   terpri!* t >> >>$

procedure printdef(x,y)$
if y then
    <<prin2!* " ("$
      prin2!* x$
      prin2!* ")   ("$
      prin2!* car y$
      prin2!* ")"$
      if cdr y then for each a in cdr y do
          <<prin2!* " AND ("$
            prin2!* a$
            prin2!* ")" >>$
      terpri!* t >>$

procedure filzero(x,n)$
% Adds zeros (in S.Q. form) to the list X from right on the length N
begin
  scalar y,i$
  y:=x$
  i:=1$
  if null x then return typerr(x,"Empty list")$
  while cdr y do
    <<y:=cdr y$
      i:=i+1>>$
  while i<n do
    <<rplacd(y,list(nil . 1))$
      y:=cdr y$
      i:=i+1 >>$
  return x
end$

procedure cutmat(x,n)$
% From each member of list X, i.e. row of a matrix, remains
% the first N elements
for each a in x collect cutrow(a,n)$

procedure cutrow(y,n)$
begin
  scalar i,z,zz$
  i:=1$
  z:=list car y$
  zz:=z$
  y:=cdr y$
  while i<n do
    <<rplacd(zz,list car y)$
      y:=cdr y$
      zz:=cdr zz$
      i:=i+1 >>$
  return z
end$

procedure cohurwp1 (deg)$
begin
  scalar k,x,y,ak,bk,akk,bkk,matr,hmat$
  % Splitting on RE and IM part
  for j:=0:deg do
    <<x:=getel list('cof,j)$
      y:=simp list('re,x)$
      x:=resimp simp list('quotient,list('difference,x,mk!*sq y),'i)$
      ak:=x . ak$
      bk:=y . bk >>$
  % Construction of coeffs. AI, BI
  positive!*:=userpos!*:=userneg!*:=nil$
  akk:=filzero(ak,2*deg)$
  bkk:=filzero(bk,2*deg)$
  k:=2$
d1:matr:=nconc(matr,list akk)$
  matr:=nconc(matr,list bkk)$
  akk:=(nil . 1) . akk$
  bkk:=(nil . 1) . bkk$
  hmat:=cutmat(matr,k)$
  x:=mk!*sq detq hmat$
  x:=positivep x$
  if null x then go to n$
  if k=2*deg then go to ko$
  k:=k+2$
  go to d1$
n:return nil$
ko:printcond(nil)$
  return t
end$

endmodule;

end;