Serge Steer [Thu, 14 Aug 2008 09:11:58 +0000 (09:11 +0000)]
14 files changed:

index 0d3d0ef..7471583 100644 (file)
@@ -17,29 +17,33 @@ function [Ns,d]=coffg(Fs)
// See also det, detr, invr, penlaur, glever, leverrier
//!
//
-if type(Fs)==2 then
-   [n,np]=size(Fs);
-   if n<>np then error('First argument to coffg must be square!');end
-   d=det(Fs) // common denominator
-   n1=n;
-   for kk=1:n1,for l=1:n1,
-     signe=(-1)^(kk+l);
-     col=[1:kk-1,kk+1:n1];row=[1:l-1,l+1:n1];
-     Ns(kk,l)=-signe*det(Fs(row,col))
-    end;end
-    Ns=-Ns;
-elseif type(Fs)==16|type(Fs)==15
-   [n,np]=size(Fs);
-   if n<>np then error('First argument to coffg must be square!');end
-   d=det(Fs) // common denominator
-   n1=n;
-   for kk=1:n1,for l=1:n1,
-     signe=(-1)^(kk+l);
-     col=[1:kk-1,kk+1:n1];row=[1:l-1,l+1:n1];
-     Ns(kk,l)=-signe*det(Fs(row,col))
-   end;end
-   Ns=-Ns;
-else
-error('First argument to coffg is invalid');
-end
+  if or(typeof(Fs)==['polynomial' 'constant']) then
+    [n,np]=size(Fs);
+    if n<>np then
+      error(msprintf(gettext("%s: Wrong size for input argument #%d: A square matrix expected.\n"),"coffg",1))
+    end
+    d=det(Fs) // common denominator
+    n1=n;
+    for kk=1:n1,for l=1:n1,
+       signe=(-1)^(kk+l);
+       col=[1:kk-1,kk+1:n1];row=[1:l-1,l+1:n1];
+       Ns(kk,l)=-signe*det(Fs(row,col))
+      end;end
+      Ns=-Ns;
+  elseif typeof(Fs)=="rational" then
+    [n,np]=size(Fs);
+    if n<>np then
+      error(msprintf(gettext("%s: Wrong size for input argument #%d: A square matrix expected.\n"),"coffg",1))
+    end
+    d=det(Fs) // common denominator
+    n1=n;
+    for kk=1:n1,for l=1:n1,
+       signe=(-1)^(kk+l);
+       col=[1:kk-1,kk+1:n1];row=[1:l-1,l+1:n1];
+       Ns(kk,l)=-signe*det(Fs(row,col))
+      end;end
+      Ns=-Ns;
+  else
+    error(msprintf(gettext("%s: Wrong type for input argument #%d: A floating point number or polynomial or rational fraction array expected.\n"),"detr",1))
+  end
endfunction
index 6fa26a9..efe2ac5 100644 (file)
@@ -17,22 +17,14 @@ function den=denom(r)
//r: rational function matrix (may be polynomial or scalar matrix)
//den: polynomial matrix
//!
-select type(r)
-case 1 then
-  den=ones(r);
-case 2 then
-  den=ones(r);
-
-//-compat next case retained for list/tlist compatibility
-case 15 then
-  r1=r(1);
-  if r1(1)<>'r' then error(92,1),end
-  den=r(3)
-case 16 then
-  r1=r(1);
-  if r1(1)<>'r' then error(92,1),end
-  den=r(3)
-else
-  error(92,1)
-end
+  select typeof(r)
+  case 'constant' then
+    den=ones(r);
+  case 'polynomial' then
+    den=ones(r);
+  case 'rational' then
+    den=r.den
+  else
+    error(msprintf(gettext("%s: Wrong type for input argument #%d: A floating point number or polynomial or rational fraction array expected.\n"),"denom",1))
+  end
endfunction
index d255a02..33f6a08 100644 (file)
@@ -12,39 +12,33 @@ function p=derivat(p)
//pd=derivat(p)  computes the derivative of the polynomial or rational
//function marix relative to the dummy variable
//!
-t=type(p)
-if t==1 then p=0*p,return,end
-if t==2 then
-  [m,n]=size(p);var=varn(p);
-  for i=1:m
-    for j=1:n
-      pij=p(i,j);nij=degree(pij);
-      if nij==0 then
-       p(i,j)=0
-      else
-       pij=coeff(pij).*(0:nij),p(i,j)=poly(pij(2:nij+1),var,'c')
-      end;
-    end;
-  end;
-  return
-end;
-
-if t==16 then
-  p1=p(1);
-  if p1(1)=='r' then
-    num=p(2);den=p(3)
+  select typeof(p)
+  case "constant" then
+    p=0*p
+  case "polynomial" then
+    [m,n]=size(p);var=varn(p);
+    for i=1:m
+      for j=1:n
+       pij=p(i,j);nij=degree(pij);
+       if nij==0 then
+         p(i,j)=0;
+       else
+         pij=coeff(pij).*(0:nij);
+         p(i,j)=poly(pij(2:nij+1),var,'c');
+       end
+      end
+    end
+  case "rational" then
+    num=p.num;den=p.den
[m,n]=size(num)
for i=1:m
for j=1:n
-       num(i,j)=derivat(num(i,j))*den(i,j)...
-           -num(i,j)*derivat(den(i,j))
-       den(i,j)=den(i,j)**2
-      end;
-    end;
-//    p=syslin(p(4),num,den)
-    p(2)=num;p(3)=den;
-    return
-  end;
-end;
-error('incorrect data type')
+       num(i,j)=derivat(num(i,j))*den(i,j)-num(i,j)*derivat(den(i,j));
+       den(i,j)=den(i,j)**2;
+      end
+    end
+    p.num=num;p.den=den;
+  else
+    error(msprintf(gettext("%s: Wrong type for input argument #%d: A floating point number or polynomial or rational fraction array expected.\n"),"derivat",1))
+  end
endfunction
index a04d212..dfa428b 100644 (file)

function res=determ(W,k)

-       // determinant of a polynomial or rational matrix by FFT
-       // W=square polynomial matrix
-       // k=``predicted'' degree of the determinant of W i.e. k is
-       // an integer larger or equal to the actual degree of W.
-       // Method: evaluate the determinant of W for the Fourier frequencies
-       // and apply inverse fft to the coefficients of the determinant.
-       // See also detr
-
-       if W==[] then
-               res=1;
-               return;
-       end;
-
-       if size(W,1)<>size(W,2) then
-               error('argument must be a square matrix'),
-       end
-
-       n1=size(W,1)
-
-       // small cases
-
-       if n1==1 then
-               res=W;
-               return;
-       elseif n1==2 then
-               res = W(1,1)*W(2,2) - W(1,2)*W(2,1);
-               return;
-       end
-
-       //upper bound of the determinant degree
-
-       maj = n1*maxi(degree(W))+1;
-
-       if argn(2)==1 then
-               k=1;
-               while k < maj,
-                       k=2*k;
-               end
-       end
-
-       if n1>8 then
-               write(%io(2),'Computing determinant: Be patient...');
-       end
-
-       // Default Values
-       e=0*ones(k,1);
-       e(2)=1;
-
-       // Paramètres de clean
-       epsa=1.d-10;
-       epsr=0;//no relative rounding
-
-       if k==1 then
-               ksi=1;
-       else
-               ksi=fft(e,-1);
-       end
-
-       fi=[];
-
-       if ~isreal(W,0) then
-               // Cas Complexe
-               for kk=1:k,
-                       fi=[fi,det(horner(W,ksi(kk)))];
-               end
-               Temp0 = poly(fft(fi,1),varn(W),'c');
-               Temp1 = clean(real(Temp0),epsa,epsr)+%i*clean(imag(Temp0),epsa,epsr);
-
-       else
-               // Cas Réel
-               for kk=1:k,fi=[fi,det(freq(W,ones(W),ksi(kk)))];end
-               Temp1 = clean(real(poly(fft(fi,1),varn(W),'c')),epsa,epsr);
-       end
-
-       if argn(2)==1 then
-
-               // Cas où k est défini dans les paramètres d'entrée.
-               // On va maintenant annuler tous les coefficients
-               // dont le degré est supérieur à maj
-
-               Temp2 = coeff(Temp1);
-               for i=1:maj,
-                       Temp2(i) = 0;
-               end
-               res = Temp1 - poly(Temp2,varn(W),"coeff");
-               return;
-
-       else
-               // Cas où k n'est pas défini dans les paramètres d'entrée
-               res = Temp1;
-               return;
-       end
+// determinant of a polynomial or rational matrix by FFT
+// W=square polynomial matrix
+// k=``predicted'' degree of the determinant of W i.e. k is
+// an integer larger or equal to the actual degree of W.
+// Method: evaluate the determinant of W for the Fourier frequencies
+// and apply inverse fft to the coefficients of the determinant.
+// See also detr
+
+  if and(typeof(W)<>['rational','polynomial','constant']) then
+    error(msprintf(gettext("%s: Wrong type for input argument #%d: A floating point number or polynomial or rational fraction array expected.\n"),"determ",1))
+  end
+  if size(W,1)<>size(W,2) then
+     error(msprintf(gettext("%s: Wrong size for input argument #%d: A square matrix expected.\n"),"determ",1))
+  end
+
+  if W==[] then
+    res=1;
+    return;
+  end;
+
+
+  n1=size(W,1)
+
+  // small cases
+
+  if n1==1 then
+    res=W;
+    return;
+  elseif n1==2 then
+    res = W(1,1)*W(2,2) - W(1,2)*W(2,1);
+    return;
+  end
+
+  //upper bound of the determinant degree
+
+  maj = n1*maxi(degree(W))+1;
+
+  if argn(2)==1 then
+    k=1;
+    while k < maj,
+      k=2*k;
+    end
+  end
+
+  if n1>8 then
+    mprintf('Computing determinant: Be patient...');
+  end
+
+  // Default Values
+  e=0*ones(k,1);
+  e(2)=1;
+
+  // Paramètres de clean
+  epsa=1.d-10;
+  epsr=0;//no relative rounding
+
+  if k==1 then
+    ksi=1;
+  else
+    ksi=fft(e,-1);
+  end
+
+  fi=[];
+
+  if ~isreal(W,0) then
+    // Cas Complexe
+    for kk=1:k,
+      fi=[fi,det(horner(W,ksi(kk)))];
+    end
+    Temp0 = poly(fft(fi,1),varn(W),'c');
+    Temp1 = clean(real(Temp0),epsa,epsr)+%i*clean(imag(Temp0),epsa,epsr);
+
+  else
+    // Cas Réel
+    for kk=1:k,fi=[fi,det(freq(W,ones(W),ksi(kk)))];end
+    Temp1 = clean(real(poly(fft(fi,1),varn(W),'c')),epsa,epsr);
+  end
+
+  if argn(2)==1 then
+
+    // Cas où k est défini dans les paramètres d'entrée.
+    // On va maintenant annuler tous les coefficients
+    // dont le degré est supérieur à maj
+
+    Temp2 = coeff(Temp1);
+    for i=1:maj,
+      Temp2(i) = 0;
+    end
+    res = Temp1 - poly(Temp2,varn(W),"coeff");
+    return;
+
+  else
+    // Cas où k n'est pas défini dans les paramètres d'entrée
+    res = Temp1;
+    return;
+  end

endfunction
index a298133..132b59f 100644 (file)
@@ -9,31 +9,37 @@

function [d]=detr(h)
-//[d]=detr(h)  computes de determinant of a polynomial or
+//[d]=detr(h)  computes the determinant of a polynomial or
//rational function matrix h using Leverrier's method
//!
-h1=h(1);
-if type(h)< 3 then
-   [m,n]=size(h);
-   if m<>n then error(20),end
-   f=eye(n,n);
-   for k=1:n-1,
-       b=h*f,
-       d=-sum(diag(b))/k
-       f=b+eye(n,n)*d,
-   end
-   d=-sum(diag(h*f))/n;
-              else //
-   if h1(1)<> 'r' then error(44),end
-   [m,n]=size(h(2));
-   if m<>n then error(20),end
-   f=eye(n,n);
-   for k=1:n-1,
-       b=h*f,
-       d=0;for l=1:n,d=d+b(l,l),end,d=-d/k;
-       f=b+eye(n,n)*d,
-   end
-   b=h*f;d=0;for l=1:n,d=d+b(l,l),end;d=-d/n;
-end
-if 2*int(n/2)<>n then d=-d;end
+  tof=typeof(h)
+  if or(tof==['polynomial','constant']) then
+    [m,n]=size(h);
+    if m<>n then
+      error(msprintf(gettext("%s: Wrong size for input argument #%d: A square matrix expected.\n"),"detr",1))
+    end
+    f=eye(n,n);
+    for k=1:n-1,
+      b=h*f,
+      d=-sum(diag(b))/k
+      f=b+eye(n,n)*d,
+    end
+    d=-sum(diag(h*f))/n;
+    if 2*int(n/2)<>n then d=-d;end
+  elseif tof=='rational' then //
+    [m,n]=size(h(2));
+    if m<>n then
+      error(msprintf(gettext("%s: Wrong size for input argument #%d: A square matrix expected.\n"),"detr",1))
+    end
+    f=eye(n,n);
+    for k=1:n-1,
+      b=h*f,
+      d=0;for l=1:n,d=d+b(l,l),end,d=-d/k;
+      f=b+eye(n,n)*d,
+    end
+    b=h*f;d=0;for l=1:n,d=d+b(l,l),end;d=-d/n;
+    if 2*int(n/2)<>n then d=-d;end
+  else
+    error(msprintf(gettext("%s: Wrong type for input argument #%d: A floating point number or polynomial or rational fraction array expected.\n"),"detr",1))
+  end
endfunction
index 50dac40..0f9507a 100644 (file)
@@ -14,21 +14,42 @@ function [lnum,lden,g]=factors(P,flag)
// and in lden the factors of denominator of P. g is the gain.
// if flag=='c' unstable roots are reflected vs the imaginary axis
// if flag=='d' unstable roots are reflected vs unit circle
-[LHS,RHS]=argn(0);
-if RHS==1 then flag=[];end
-if type(P)==2 then [lnum,lden]=pfactors(P,flag);
-  if LHS >=3 then error('invalid LHS');end
-  return;end
-  if type(P)==16 then
-    if typeof(P)=='state-space' then P=ss2tf(P);end
-    [lnum,gn]=pfactors(P('num'),flag);
-    [lden,gd]=pfactors(P('den'),flag);g=gn/gd;
+  [LHS,RHS]=argn(0);
+  if RHS==1 then flag=[];end
+  select typeof(P)
+  case 'polynomial' then
+    if size(P,'*')<>1 then
+      error(msprintf(gettext("%s: Wrong size for input argument #%d: A polynomial expected.\n"),"factors",1))
+    end
+    [lnum,lden]=pfactors(P,flag);
+  case 'rational' then
+    if size(P,'*')<>1 then
+      error(msprintf(gettext("%s: Wrong size for input argument #%d: A rational fraction expected.\n"),"factors",1))
+    end
+    [lnum,gn]=pfactors(P.num,flag);
+    [lden,gd]=pfactors(P.den,flag);
+    g=gn/gd;
if LHS==1 then
num=g;
for k=lnum;num=num.*k;end
den=1;
for k=lden;den=den.*k;end
-      lnum=syslin(P('dt'),num,den);return
+      lnum=syslin(P.dt,num,den);return
+    end
+  case 'state-space' then
+    if size(P,'*')<>1 then
+      error(msprintf(gettext("%s: Wrong size for input argument #%d: A single input, single output system expected.\n"),"factors",1))
+    end
+
+    P=ss2tf(P)
+    [lnum,gn]=pfactors(P.num,flag);
+    [lden,gd]=pfactors(P.den,flag);g=gn/gd;
+    if LHS==1 then
+      num=g;for k=lnum;num=num.*k;end
+      den=1;for k=lden;den=den.*k;end
+      lnum=syslin(P.dt,num,den);return
end
+  else
+    error(msprintf(gettext("%s: Wrong type for input argument #%d: Linear dynamical system or polynomial array expected.\n" ),"factors",1))
end
endfunction
index 62c2e29..66cac08 100644 (file)
@@ -14,20 +14,26 @@ function [a,u]=hermit(a)
//triangular. The output value of A is A*U.
//Warning: Experimental version
//!
-[m,n]=size(a);if m<>n then error('square matrix only!'),end
-[a,u]=htrianr(a)
-for l=n-1:-1:1
- dl(l:n)=degree(a(l,l:n));
- for k=l+1:n
-    if dl(k)>=dl(l) then
-       all=a(l,l);
-       if norm(coeff(all),1) > 1.d-10 then
-       [r,q]=pdiv(a(l,k),a(l,l))
-       if l>1 then a(1:l-1,k)=a(1:l-1,k)-a(1:l-1,l)*q;end
-       a(l,k)=r
-       u(:,k)=u(:,k)-u(:,l)*q
-       end
+  if type(a)>2 then
+    error(msprintf(gettext("%s: Wrong type for input argument: Polynomial array expected.\n"),"hermit"))
+  end
+  [m,n]=size(a);
+  if m<>n then
+    error(msprintf(gettext("%s: Wrong size for input argument #%d: A square matrix expected.\n" ),"hermit",1))
+  end
+  [a,u]=htrianr(a)
+  for l=n-1:-1:1
+    dl(l:n)=degree(a(l,l:n));
+    for k=l+1:n
+      if dl(k)>=dl(l) then
+       all=a(l,l);
+       if norm(coeff(all),1) > 1.d-10 then
+         [r,q]=pdiv(a(l,k),a(l,l))
+         if l>1 then a(1:l-1,k)=a(1:l-1,k)-a(1:l-1,l)*q;end
+         a(l,k)=r
+         u(:,k)=u(:,k)-u(:,l)*q
+       end
+      end
end
- end
-end
+  end
endfunction
index 29eff37..1d719f2 100644 (file)
@@ -20,10 +20,12 @@ function r=horner(p,x)
// See also: freq, repfreq.
//!
//
-  if argn(2)<>2 then error('horner requires two arguments'),end
+  if argn(2)<>2 then
+    error(msprintf(gettext("%s: Wrong number of input argument(s): %d expected.\n"),'horner',2))
+  end
if x==[]|p==[] then r=[],return,end
tp=type(p)
-  if tp<=10 then
+  if tp<=2 then
[m,n]=size(p)
if m==-1 then indef=%t,m=1,n=1,p=p+0;else indef=%f,end
if m*n==1 then //special case for a single polynomial
@@ -47,9 +49,11 @@ function r=horner(p,x)
end
end
if indef then r=r*eye(),end
-  elseif tp==16 then
+  elseif typeof(p)=='rational' then
r=horner(p(2),x)./horner(p(3),x),
elseif tp==129 then
r=horner(p(:),x);r=r(1):r(2):r(3)
+  else
+    error(msprintf(gettext("%s: Unexpected type for input argument #%d.\n"),'horner',1))
end
endfunction
index 3d56667..0fa6008 100644 (file)
@@ -13,12 +13,14 @@ function [pg,U]=hrmt(v)
// Finds unimodular U and pg = gcd of a row of polynomials v
// such that v*U = [pg,0]
//!
-[n,m]=size(v)
-if n>1 then error(60),end
-pg=v(1)
-U=eye(m,m)
-for k=2:m
- [pg,uk]=bezout(pg,v(k))
- U(:,k-1:k)=U(:,k-1:k)*uk(:,[2 1])
-end
+  [n,m]=size(v)
+  if n>1 then
+    error(msprintf(gettext("%s: Wrong size for input argument #%d: A row vector expected.\n"),"hrmt",1))
+  end
+  pg=v(1)
+  U=eye(m,m)
+  for k=2:m
+    [pg,uk]=bezout(pg,v(k))
+    U(:,k-1:k)=U(:,k-1:k)*uk(:,[2 1])
+  end
endfunction
index 757f5c1..a1fa337 100644 (file)
@@ -16,36 +16,40 @@ function [A,U,rk]=htrianr(A)
//rk=normal rank of A
//Warning: there is an elimination of neglectable terms
//!
-A=clean(A);
-[m,n]=size(A);U=eye(n,n);
-l1=n+1;
-for l=m:-1:maxi((m-n),1)
- l1=l1-1;
- if l1<>0 then
-  Al=A(l,1:l1);
-  if norm(coeff(Al),1) > 1.d-10 then
-    [pg,Ul]=hrmt(Al);
-    Ul=clean(Ul,1.d-10);
-    A(l,1:l1)=[0*ones(1,l1-1) pg];
-    U(:,1:l1)=U(:,1:l1)*Ul;
-         if l>1 then
+  if type(A)>2 then
+    error(msprintf(gettext("%s: Wrong type for input argument: Polynomial array expected.\n"),"htrianr"))
+  end
+  A=clean(A);
+  [m,n]=size(A);U=eye(n,n);
+  l1=n+1;
+  for l=m:-1:maxi((m-n),1)
+    l1=l1-1;
+    if l1<>0 then
+      Al=A(l,1:l1);
+      if norm(coeff(Al),1) > 1.d-10 then
+       [pg,Ul]=hrmt(Al);
+       Ul=clean(Ul,1.d-10);
+       A(l,1:l1)=[0*ones(1,l1-1) pg];
+       U(:,1:l1)=U(:,1:l1)*Ul;
+       if l>1 then
A(1:l-1,1:l1)=A(1:l-1,1:l1)*Ul;
end
-  else
-  l1=l1+1
+      else
+       l1=l1+1
+      end
+    end
end
- end
-end
-U=clean(U,1.d-10);
-k0=0;k1=0;tol=norm(coeff(A),1);
-v=[];w=[];
-for k=1:n
-   if maxi(abs(coeff(A(:,k)))) <= sqrt(%eps)*tol then
-      k0=k0+1;v=[v,k];                           else
+  U=clean(U,1.d-10);
+  k0=0;k1=0;tol=norm(coeff(A),1);
+  v=[];w=[];
+  for k=1:n
+    if maxi(abs(coeff(A(:,k)))) <= sqrt(%eps)*tol then
+      k0=k0+1;v=[v,k];
+    else
k1=k1+1,w=[w,k];
-   end
-end
-ww=[v,w];
-A=A(:,ww);U=U(:,ww);
-rk=n-k0;
+    end
+  end
+  ww=[v,w];
+  A=A(:,ww);U=U(:,ww);
+  rk=n-k0;
endfunction
index d971d85..4255f90 100644 (file)
@@ -21,7 +21,7 @@ function [P]=inv_coeff(c,d,name)
return,
end
if modulo(m,d+1) <> 0 then
-    error(msprintf(_("%s: incompatible  input arguments %d and %d"),"inv_coeff",1,2))
+    error(msprintf(_("%s: incompatible input arguments %d and %d"),"inv_coeff",1,2))
end
p=poly(0,name);
P=p.^(0:d)';
index fbdd83f..38fc380 100644 (file)
@@ -12,90 +12,83 @@ function [f,d]=invr(h,flag)
//if h is a scalar, polynomial or rational fonction matrix, invr
//computes h^(-1).
//!
-
-if type(h)==1 then f=inv(h);return;end
-[lhs,rhs]=argn(0);
-if rhs==1 then flag='C';end
-if type(h) == 2 then
-  //     POLYNOMIAL MATRIX
-  [m,n]=size(h);
-  if m<>n then error(20),end
-  ndeg=maxi(degree(h));
-  if ndeg==1 then
-    //           MATRIX PENCIL
-    E=coeff(h,1);A=-coeff(h,0);
-    if norm(E-eye(E),1) < 100*%eps then
-      // sI -A
-      [num,den]=coff(A,varn(h));f=num/den;
-      return
-    end
-    [Bfs,Bis,chis]=glever(E,A,varn(h));
-    f=Bfs/chis - Bis;
-    if lhs==2 then
-      d=lcm(f('den'));
-      f=f*d;f=f('num');
+  if argn(2)==1 then
+    flag='C';
+  end
+  lhs=argn(1)
+  select typeof(h)
+  case 'constant' then
+    f=inv(h);
+  case 'polynomial' then //POLYNOMIAL MATRIX
+    [m,n]=size(h);
+    if m<>n then error(20),end
+    ndeg=maxi(degree(h));
+    if ndeg==1 then   //MATRIX PENCIL
+      E=coeff(h,1);A=-coeff(h,0);
+      if norm(E-eye(E),1) < 100*%eps then
+       // sI -A
+       [num,den]=coff(A,varn(h));f=num/den;
+      else
+       [Bfs,Bis,chis]=glever(E,A,varn(h));
+       f=Bfs/chis - Bis;
+       if lhs==2 then
+         d=lcm(f('den'));
+         f=f*d;f=f('num');
+       end
+      end
+    else // GENERAL POLYNOMIAL MATRIX
+      select flag
+      case 'L'
+       f=eye(n,n);
+       for k=1:n-1,
+         b=h*f,
+         d=-sum(diag(b))/k
+         f=b+eye(n,n)*d,
+       end;
+       d=sum(diag(h*f))/n,
+       if degree(d)==0 then d=coeff(d),end,
+       if lhs==1 then f=f/d;end
+      case 'C'
+       [f,d]=coffg(h);
+       if degree(d)==0 then d=coeff(d),end
+       if lhs==1 then f=f/d;end
+      else
+       error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),..
+                      "invr",2,"''C'',''D''"))
+      end;
end
-    return;
-  end;
-  //   GENERAL POLYNOMIAL MATRIX
-  select flag
-  case 'L'
-    f=eye(n,n);
-    for k=1:n-1,
-      b=h*f,
-      d=-sum(diag(b))/k
-      f=b+eye(n,n)*d,
-    end;
-    d=sum(diag(h*f))/n,
-    if degree(d)==0 then d=coeff(d),end,
-    if lhs==1 then f=f/d;end
-    return;
-  case 'C'
-    [f,d]=coffg(h);
-    if degree(d)==0 then d=coeff(d),end
-    if lhs==1 then f=f/d;end
-    return;
-  end;
-end
-
-//-compat type(h)==15 retained for list/tlist compatibility
-if type(h)==15|type(h)==16 then
-  h1=h(1);
-  if h1(1)<> 'r' then error(44),end
-  [m,n]=size(h(2));
-  if m<>n then error(20),end
-  select flag
-  case 'L'
-    //    Leverrier
-    f=eye(n,n);
-    for k=1:n-1,
-      b=h*f,
-      d=0;for l=1:n,d=d+b(l,l),end,d=-d/k;
-      f=b+eye(n,n)*d,
+  case 'rational' then
+    [m,n]=size(h(2));
+    if m<>n then error(20),end
+    select flag
+    case 'L' //    Leverrier
+      f=eye(n,n);
+      for k=1:n-1,
+       b=h*f,
+       d=0;for l=1:n,d=d+b(l,l),end,d=-d/k;
+       f=b+eye(n,n)*d,
+      end;
+      b=h*f;d=0;for l=1:n,d=d+b(l,l),end;d=d/n,
+      if lhs==1 then f=f/d;end
+    case 'A' //   lcm of all denominator entries
+      denh=lcm(h('den'));
+      Num=h*denh;Num=Num('num');
+      [N,d]=coffg(Num);
+      f=N*denh;
+      if lhs==1 then f=f/d;end
+    case 'C'// default method by polynomial inverse
+      [Nh,Dh]=lcmdiag(h); //h=Nh*inv(Dh); Dh diagonal;
+      [N,d]=coffg(Nh);
+      f=Dh*N;
+      if lhs==1 then f=f/d;end
+    case 'Cof'// cofactors method
+      [f,d]=coffg(h);
+      if lhs==1 then f=f/d;end
+    else
+       error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),..
+                      "invr",2,"''L'',''A'',''C'',''Cof''"))
end;
-    b=h*f;d=0;for l=1:n,d=d+b(l,l),end;d=d/n,
-    if lhs==1 then f=f/d;end
-    return;
-  case 'A'
-    //   lcm of all denominator entries
-    denh=lcm(h('den'));
-    Num=h*denh;Num=Num('num');
-    [N,d]=coffg(Num);
-    f=N*denh;
-    if lhs==1 then f=f/d;end
-    return;
-  case 'C'
-    // default method by polynomial inverse
-    [Nh,Dh]=lcmdiag(h); //h=Nh*inv(Dh); Dh diagonal;
-    [N,d]=coffg(Nh);
-    f=Dh*N;
-    if lhs==1 then f=f/d;end
-    return;
-  case 'Cof'
-    // cofactors method
-    [f,d]=coffg(h);
-    if lhs==1 then f=f/d;end
+  else
+    error(msprintf(gettext("%s: Wrong type for input argument #%d: A floating point number or polynomial or rational fraction array expected.\n"),"invr",1))
end;
-end;
-error('invalid input to invr');
endfunction
@@ -12,21 +12,15 @@ function num=numer(r)
//returns the numerator num of a rational function matrix r (r may be
//also a scalar or polynomial matrix
//!
-r1=r(1);
-select type(r)
-case 1 then
-  num=r;
-case 2 then
-  num=r;
+  select typeof(r)
+  case 'constant' then
+    num=r;
+  case 'polynomial' then
+    num=r;
+  case 'rational' then
+    num=r.num
+  else
+    error(msprintf(gettext("%s: Wrong type for input argument #%d: A floating point number or polynomial or rational fraction array expected.\n"),"numer",1))
+  end

-//-compat next case retained for list/tlist compatibility
-case 15 then
-  if r1(1)<>'r' then error(92,1),end
-  num=r(2)
-case 16 then
-  if r1(1)<>'r' then error(92,1),end
-  num=r(2)
-else
-  error(92,1)
-end
endfunction
index c1f7746..507f425 100644 (file)
@@ -12,26 +12,35 @@ function [f]=polfact(p)
// Minmal factors of p
// f=polfact(p)
//
-// p : polynomila
+// p : polynomial
// f : vector [f0 f1 ... fn] such that p=prod(f)
//     - f0  constant
//     - fi polynomial
//!
//
-if size(p,'*')<>1 then error('polynomial argument required!'),end
-if norm(imag(coeff(p)))<>0 then error('real case only!'),end
-n=degree(p);f=coeff(p,n);
-if n==0 then return,end
-var=varn(p);
-r=roots(p);[s,k]=sort(abs(r));r=r(k)
-k=1;
-while k<=n do,
-  if imag(r(k))<>0 then
-    f=[f,real(poly(r(k:k+1),var))]
-    k=k+2
-  else
-    f=[f,poly(r(k),var)]
-    k=k+1
+  if type(p)>2 then
+    error(msprintf(gettext("%s: Wrong type for input argument: Polynomial array expected.\n"),'polfact'))
+  end
+  if size(p,'*')<>1 then
+    error(msprintf(gettext("%s: Wrong size for input argument #%d: A polynomial expected.\n"),'polfact',1))
+  end
+
+  if norm(imag(coeff(p)))<>0 then
+    error(msprintf(gettext("%s: Input argument #%d must be real.\n"),'polfact',1))
+  end
+  p=real(p)
+  n=degree(p);f=coeff(p,n);
+  if n==0 then return,end
+  var=varn(p);
+  r=roots(p);[s,k]=sort(abs(r));r=r(k)
+  k=1;
+  while k<=n do,
+    if imag(r(k))<>0 then
+      f=[f,real(poly(r(k:k+1),var))]
+      k=k+2
+    else
+      f=[f,poly(r(k),var)]
+      k=k+1
+    end
end
-end
endfunction