diff options
| author | Shashank | 2017-05-29 12:40:26 +0530 |
|---|---|---|
| committer | Shashank | 2017-05-29 12:40:26 +0530 |
| commit | 0345245e860375a32c9a437c4a9d9cae807134e9 (patch) | |
| tree | ad51ecbfa7bcd3cc5f09834f1bb8c08feaa526a4 /modules/polynomials/macros | |
| download | scilab_for_xcos_on_cloud-0345245e860375a32c9a437c4a9d9cae807134e9.tar.gz scilab_for_xcos_on_cloud-0345245e860375a32c9a437c4a9d9cae807134e9.tar.bz2 scilab_for_xcos_on_cloud-0345245e860375a32c9a437c4a9d9cae807134e9.zip | |
CMSCOPE changed
Diffstat (limited to 'modules/polynomials/macros')
61 files changed, 1393 insertions, 0 deletions
diff --git a/modules/polynomials/macros/buildmacros.bat b/modules/polynomials/macros/buildmacros.bat new file mode 100755 index 000000000..64bdd6d01 --- /dev/null +++ b/modules/polynomials/macros/buildmacros.bat @@ -0,0 +1,11 @@ +rem Scilab ( http://mwww.scilab.org/ ) - This file is part of Scilab +rem Copyright (C) ????-2008 - INRIA - Allan CORNET +rem +rem This file must be used under the terms of the CeCILL. +rem This source file is licensed as described in the file COPYING, which +rem you should have received as part of this distribution. The terms +rem are also available at +rem http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +@..\..\..\bin\scilex -nwni -ns -e exec('buildmacros.sce');quit;
\ No newline at end of file diff --git a/modules/polynomials/macros/buildmacros.sce b/modules/polynomials/macros/buildmacros.sce new file mode 100755 index 000000000..3b3d3c72e --- /dev/null +++ b/modules/polynomials/macros/buildmacros.sce @@ -0,0 +1,17 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) 2005-2008 - INRIA - Allan CORNET +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +//------------------------------------ +if (isdef("genlib") == %f) then + exec(SCI+"/modules/functions/scripts/buildmacros/loadgenlib.sce"); +end +//------------------------------------ +genlib("polynomialslib","SCI/modules/polynomials/macros",%f,%t); +//------------------------------------ diff --git a/modules/polynomials/macros/chepol.bin b/modules/polynomials/macros/chepol.bin Binary files differnew file mode 100755 index 000000000..7c573d9f2 --- /dev/null +++ b/modules/polynomials/macros/chepol.bin diff --git a/modules/polynomials/macros/chepol.sci b/modules/polynomials/macros/chepol.sci new file mode 100755 index 000000000..522170482 --- /dev/null +++ b/modules/polynomials/macros/chepol.sci @@ -0,0 +1,30 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) INRIA - F.D +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + +function Tn=chepol(n,var) + //Chebychev polynomial + // n :Polynomial order + // var :Polynomial variable (character string) + // Tn :Polynomial in var + // + //! + if n==0 then + Tn=poly(1,var,"coeff"), + elseif n==1 then + Tn=poly(0,var); + else + T0=poly(1,var,"coeff"); + T1=poly(0,var) + for k=2:n + Tn=2*poly(0,var)*T1-T0 + [T1,T0]=(Tn,T1); + end + end + +endfunction diff --git a/modules/polynomials/macros/cleanmacros.bat b/modules/polynomials/macros/cleanmacros.bat new file mode 100755 index 000000000..6294097b5 --- /dev/null +++ b/modules/polynomials/macros/cleanmacros.bat @@ -0,0 +1,13 @@ +rem Scilab ( http://mwww.scilab.org/ ) - This file is part of Scilab +rem Copyright (C) ????-2008 - INRIA - Allan CORNET +rem +rem This file must be used under the terms of the CeCILL. +rem This source file is licensed as described in the file COPYING, which +rem you should have received as part of this distribution. The terms +rem are also available at +rem http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +@del *.bin 2>NUL +@del lib 2>NUL +@del names 2>NUL
\ No newline at end of file diff --git a/modules/polynomials/macros/cmndred.bin b/modules/polynomials/macros/cmndred.bin Binary files differnew file mode 100755 index 000000000..84f8da2f8 --- /dev/null +++ b/modules/polynomials/macros/cmndred.bin diff --git a/modules/polynomials/macros/cmndred.sci b/modules/polynomials/macros/cmndred.sci new file mode 100755 index 000000000..ba80d9f82 --- /dev/null +++ b/modules/polynomials/macros/cmndred.sci @@ -0,0 +1,25 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function [n,d]=cmndred(num,den) + //Syntax: [num,den]=cmndred(num,den) + // + //Given the transfert matrix defined by num./den, cmndred computes + //polynomial matrix n and a common denominator d such that + //n/d=num./den + //! + [m1,n1]=size(num) + d=1;for dk=matrix(den,1,m1*n1),d=lcm([d,dk]),end + for l=1:m1 + for k=1:n1 + n(l,k)=num(l,k)*pdiv(d,den(l,k)); + end; + end; +endfunction diff --git a/modules/polynomials/macros/coffg.bin b/modules/polynomials/macros/coffg.bin Binary files differnew file mode 100755 index 000000000..c33255677 --- /dev/null +++ b/modules/polynomials/macros/coffg.bin diff --git a/modules/polynomials/macros/coffg.sci b/modules/polynomials/macros/coffg.sci new file mode 100755 index 000000000..edba924a7 --- /dev/null +++ b/modules/polynomials/macros/coffg.sci @@ -0,0 +1,49 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA - Francois DELBECQUE +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function [Ns,d]=coffg(Fs) + // [Ns,d]=coffg(Fs) computes Fs^-1 where Fs is a polynomial + // or rational matrix by co-factors method. + // d = common denominator; Ns = numerator (matrix polynomial) + // Fs inverse = Ns/d. + // (Be patient...results are generally reliable) + // See also det, detr, invr, penlaur, glever, leverrier + //! + // + 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 diff --git a/modules/polynomials/macros/colcompr.bin b/modules/polynomials/macros/colcompr.bin Binary files differnew file mode 100755 index 000000000..4af35776e --- /dev/null +++ b/modules/polynomials/macros/colcompr.bin diff --git a/modules/polynomials/macros/colcompr.sci b/modules/polynomials/macros/colcompr.sci new file mode 100755 index 000000000..4a67d7adf --- /dev/null +++ b/modules/polynomials/macros/colcompr.sci @@ -0,0 +1,22 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function [Y,rk,Ac]=colcompr(A); + //[Y,rk,Ac]=colcompr(A); + //column compression of polynomial matrix A + //(left compression) + //Y = right unimodular base + //rk = normal rank of A + //Ac = A*Y + //see rowcompr + //! + [m,n]=size(A); + [Ac,Y,rk]=htrianr(A); +endfunction diff --git a/modules/polynomials/macros/denom.bin b/modules/polynomials/macros/denom.bin Binary files differnew file mode 100755 index 000000000..277b04611 --- /dev/null +++ b/modules/polynomials/macros/denom.bin diff --git a/modules/polynomials/macros/denom.sci b/modules/polynomials/macros/denom.sci new file mode 100755 index 000000000..856b5be14 --- /dev/null +++ b/modules/polynomials/macros/denom.sci @@ -0,0 +1,29 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// ... +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function den=denom(r) + //returns the denominator of a rational matrix + //%Syntax: den=denom(r) + //with + //r: rational function matrix (may be polynomial or scalar matrix) + //den: polynomial matrix + //! + 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 diff --git a/modules/polynomials/macros/derivat.bin b/modules/polynomials/macros/derivat.bin Binary files differnew file mode 100755 index 000000000..caa4e5f8c --- /dev/null +++ b/modules/polynomials/macros/derivat.bin diff --git a/modules/polynomials/macros/derivat.sci b/modules/polynomials/macros/derivat.sci new file mode 100755 index 000000000..56ca58a66 --- /dev/null +++ b/modules/polynomials/macros/derivat.sci @@ -0,0 +1,63 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + +function [p] = derivat(p) + //pd=derivat(p) computes the derivative of the polynomial or rational + //function marix relative to the dummy variable + // + + select typeof(p) + + case "constant" then + p = 0*p; // includes case p == [] + + case "polynomial" then + var = varn(p); // dummy variable + [m,n]=size(p); + deg = degree(p); + degmax = max(deg); + + if (degmax < n*m) + x = poly(0, var); + pd = zeros(p) * x; + for i = 1:degmax + pd = pd + coeff(p,i) * i * x^(i-1); + end + p = pd; + else + for i=1:m + for j=1:n + pij=p(i,j); + nij=deg(i,j); + 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 + end + + case "rational" then + num = p.num; + den = p.den; + + num = derivat(num) .* den - num .* derivat(den); + den = den.^2; + + 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 diff --git a/modules/polynomials/macros/determ.bin b/modules/polynomials/macros/determ.bin Binary files differnew file mode 100755 index 000000000..e2e663d4b --- /dev/null +++ b/modules/polynomials/macros/determ.bin diff --git a/modules/polynomials/macros/determ.sci b/modules/polynomials/macros/determ.sci new file mode 100755 index 000000000..efa7a90ca --- /dev/null +++ b/modules/polynomials/macros/determ.sci @@ -0,0 +1,106 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA - Francois DELBECQUE +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +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 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*max(degree(W))+1; + + if argn(2)==1 then + k=1; + while k < maj, + k=2*k; + end + 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 diff --git a/modules/polynomials/macros/detr.bin b/modules/polynomials/macros/detr.bin Binary files differnew file mode 100755 index 000000000..94005dbc0 --- /dev/null +++ b/modules/polynomials/macros/detr.bin diff --git a/modules/polynomials/macros/detr.sci b/modules/polynomials/macros/detr.sci new file mode 100755 index 000000000..da7058c11 --- /dev/null +++ b/modules/polynomials/macros/detr.sci @@ -0,0 +1,37 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function [d]=detr(h) + //[d]=detr(h) computes the determinant of a polynomial or + //rational function matrix h using Leverrier's method + //! + rhs = argn(2); + if rhs <> 1 then + error(msprintf(gettext("%s: Wrong number of input argument: %d expected.\n"), "detr", 1)); + end + + tof=typeof(h) + if or(tof==["polynomial","constant", "rational"]) 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 + 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 diff --git a/modules/polynomials/macros/diophant.bin b/modules/polynomials/macros/diophant.bin Binary files differnew file mode 100755 index 000000000..a4c1fbd96 --- /dev/null +++ b/modules/polynomials/macros/diophant.bin diff --git a/modules/polynomials/macros/diophant.sci b/modules/polynomials/macros/diophant.sci new file mode 100755 index 000000000..d052415bc --- /dev/null +++ b/modules/polynomials/macros/diophant.sci @@ -0,0 +1,30 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function [x,err]=diophant(p1p2,b) + //solves diophantine equation p1*x1+p2*x2=b + //with p1p2 a polynomial vector [p1 p2] + //b polynomial + //x polynomial vector [x1;x2] + //if the equation is uncompatible err=||p1*x1+p2*x2-b||/||b|| + //else err=0 + //! + p1=p1p2(1);p2=p1p2(2) + [x,u]=bezout(p1,p2) + p1=u(2,2);p2=u(1,2)// + if degree(x)==0 then + x=b*u(:,1) + err=0 + else + [r,q]=pdiv(b,x) + err=norm(coeff(b-x*q),2)/norm(coeff(b),2) + x=q*u(:,1) + end +endfunction diff --git a/modules/polynomials/macros/factors.bin b/modules/polynomials/macros/factors.bin Binary files differnew file mode 100755 index 000000000..8974f7871 --- /dev/null +++ b/modules/polynomials/macros/factors.bin diff --git a/modules/polynomials/macros/factors.sci b/modules/polynomials/macros/factors.sci new file mode 100755 index 000000000..76447ba45 --- /dev/null +++ b/modules/polynomials/macros/factors.sci @@ -0,0 +1,55 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function [lnum,lden,g]=factors(P,flag) + //Given a polynomial or rational P, returns in list lnum polynomials of + //degree 1 or two which are the factors of numerators of P. + // 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 + 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 + 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 diff --git a/modules/polynomials/macros/gcd.bin b/modules/polynomials/macros/gcd.bin Binary files differnew file mode 100755 index 000000000..66145fe9b --- /dev/null +++ b/modules/polynomials/macros/gcd.bin diff --git a/modules/polynomials/macros/gcd.sci b/modules/polynomials/macros/gcd.sci new file mode 100755 index 000000000..ee9f47a80 --- /dev/null +++ b/modules/polynomials/macros/gcd.sci @@ -0,0 +1,58 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function [x, uu] = gcd(p) + //Given a polynomial vector p, [pgcd,u]=gcd(p) computes the gcd + //of components and a unimodular matrix (with polynomial inverse) u, + //with minimal degree such that [p1 p2]*u=[0 ... 0 pgcd] + //! + + if type(p)<>1 & type(p)<>2 & type(p)<>8 then + error(msprintf(_("%s: Wrong type for argument #%d: Integer array or Polynomial expected.\n"), "gcd", 1)); + end + + if type(p)==1 then + if floor(p)<>p then + error(msprintf(_("%s: Wrong type for argument #%d: Integer array or Polynomial expected.\n"), "gcd", 1)); + else + p = iconvert(p,4); + end + end + + [lhs,rhs]=argn(0) + if type(p)==8 then + if lhs==2 then [x,uu]=%i_gcd(p), else x=%i_gcd(p), end + return + end + + [m, n] = size(p) + mn = m*n + p = matrix(p, 1, mn) + x = p(1); + uu = 1 + for l = 2:mn, + [x, u] = bezout(x, p(l)), + if lhs==2 then + uu = [uu(:, 1:l-2) uu(:, l-1)*u(1, [2 1])]; uu(l, l-1:l) = u(2, [2 1]); + end + end, + if lhs==1 then return end + for l = mn:-1:2 + pivot = uu(l, l-1); + for k = l:mn + [r, q] = pdiv(uu(l, k), pivot) + if coeff(q)<>0 then + uu(1:l-1, k) = uu(1:l-1, k)-q*uu(1:l-1, l-1) + uu(l, k) = r; + end + end + end + +endfunction diff --git a/modules/polynomials/macros/hermit.bin b/modules/polynomials/macros/hermit.bin Binary files differnew file mode 100755 index 000000000..23a29226c --- /dev/null +++ b/modules/polynomials/macros/hermit.bin diff --git a/modules/polynomials/macros/hermit.sci b/modules/polynomials/macros/hermit.sci new file mode 100755 index 000000000..150f23abe --- /dev/null +++ b/modules/polynomials/macros/hermit.sci @@ -0,0 +1,39 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function [a,u]=hermit(a) + //[A,U]=hermit(A) + //Hermite form: U is an unimodular matrix such that A*U is + //triangular. The output value of A is A*U. + //Warning: Experimental version + //! + 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 +endfunction diff --git a/modules/polynomials/macros/horner.bin b/modules/polynomials/macros/horner.bin Binary files differnew file mode 100755 index 000000000..a1c927b64 --- /dev/null +++ b/modules/polynomials/macros/horner.bin diff --git a/modules/polynomials/macros/horner.sci b/modules/polynomials/macros/horner.sci new file mode 100755 index 000000000..506ada334 --- /dev/null +++ b/modules/polynomials/macros/horner.sci @@ -0,0 +1,120 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function [r] = horner(p,x) + // horner(P,x) evaluates the polynomial or rational matrix P = P(s) + // when the variable s of the polynomial is replaced by x + // x can be a scalar or polynomial or rational matrix. + // Example: bilinear transform; Assume P = P(s) is a rational matrix + // then the rational matrix P((1+s)/(1-s)) is obtained by + // horner(P,(1+s)/(1-s)); + // To evaluate a rational matrix at given frequencies use + // preferably the freq primitive ; + // See also: freq, repfreq. + // Improvements: + // Special cases aded to improve efficiency: + // - p = row vector, x = column vector + // - p = column vector, x = row vector + // - x = scalar + //! + // + if (argn(2) <> 2) then + error(msprintf(gettext("%s: Wrong number of input argument(s): %d expected.\n"),"horner",2)) + end + + if (size(x, "*") == 0 | size(p, "*") == 0) then + r = [] + return + end + + tp = type(p) + + if (tp <= 2) then + // tp <= 2 <=> matrix of reals, complexes or polynomials + [m,n] = size(p) + + if (m == -1) then + indef=%t, m=1, n=1, p=p+0 + else + indef=%f + end + + [mx,nx] = size(x) + + if (m*n == 1) then + // special case: p = 1x1 polynomial, x = matrix + cp = coeff(p) + r = cp($) * ones(x) + for (k = degree(p) : -1 : 1) + r = r .* x + cp(k) + end + + elseif (n*mx == 1) + // p = one column, x = one row + nd = max(degree(p)); + r = zeros(p) * x; + for (k = nd : -1: 0) + c = coeff(p, k); + r = r .* (ones(p) * x) + c * ones(x); + end + + elseif (m*nx == 1) + // p = one row, x = one column + nd = max(degree(p)); + r = x * zeros(p); + for (k = nd : -1: 0) + c = coeff(p, k); + r = r .* (x * ones(p))+ ones(x) * c; + end + + elseif (mx*nx == 1) + // p = matrix, x = scalar + nd = max(degree(p)); + r = zeros(p); + for (k = nd : -1: 0) + c = coeff(p, k); + r = r * x + c; + end + + else + // other cases + r = [] + for (l = 1 : m) + rk = [] + for (k = 1 : n) + plk = p(l,k) + d = degree(plk) + rlk = coeff(plk,d) * ones(x); // for the case horner(1,x) + for (kk = 1 : d) + rlk = rlk .* x + coeff(plk,d-kk) + end + rk = [rk, rlk] + end + r = [r; rk] + end + end + + if (indef) then + r = r * eye() + end + + elseif (typeof(p) == "rational") then + r = horner(p(2),x) ./ horner(p(3),x) + + elseif (tp == 129) then + // implicit polynomial for indexing + 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 diff --git a/modules/polynomials/macros/hrmt.bin b/modules/polynomials/macros/hrmt.bin Binary files differnew file mode 100755 index 000000000..439eab760 --- /dev/null +++ b/modules/polynomials/macros/hrmt.bin diff --git a/modules/polynomials/macros/hrmt.sci b/modules/polynomials/macros/hrmt.sci new file mode 100755 index 000000000..9cff1de0c --- /dev/null +++ b/modules/polynomials/macros/hrmt.sci @@ -0,0 +1,26 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function [pg,U]=hrmt(v) + // Syntax: [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(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 diff --git a/modules/polynomials/macros/htrianr.bin b/modules/polynomials/macros/htrianr.bin Binary files differnew file mode 100755 index 000000000..ba6b58256 --- /dev/null +++ b/modules/polynomials/macros/htrianr.bin diff --git a/modules/polynomials/macros/htrianr.sci b/modules/polynomials/macros/htrianr.sci new file mode 100755 index 000000000..422a1d454 --- /dev/null +++ b/modules/polynomials/macros/htrianr.sci @@ -0,0 +1,55 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function [A,U,rk]=htrianr(A) + //[A,U,rk]=htrianr(a) + //triangularization of polynomial matrix A. A is [m,n], m<=n. + //U=right unimodular basis + //the output value of A equals A*U + //rk=normal rank of A + //Warning: there is an elimination of neglectable terms + //! + 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:max((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 + 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 max(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; +endfunction diff --git a/modules/polynomials/macros/inv_coeff.bin b/modules/polynomials/macros/inv_coeff.bin Binary files differnew file mode 100755 index 000000000..2455a7546 --- /dev/null +++ b/modules/polynomials/macros/inv_coeff.bin diff --git a/modules/polynomials/macros/inv_coeff.sci b/modules/polynomials/macros/inv_coeff.sci new file mode 100755 index 000000000..f79a1db93 --- /dev/null +++ b/modules/polynomials/macros/inv_coeff.sci @@ -0,0 +1,29 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function [P]=inv_coeff(c,d,name) + // inverse function of coeff + rhs=argn(2); + if rhs <= 2 ; name = "x";end + + [n,m]=size(c); + if rhs <= 1 ; d = (m/n-1) ; end + if d==m-1 then + P=[]; + for l=1:n, P=[P;poly(c(l,:),name,"coeff")];end + return, + end + if modulo(m,d+1) <> 0 then + error(msprintf(_("%s: incompatible input arguments %d and %d\n"),"inv_coeff",1,2)) + end + p=poly(0,name); + P=p.^(0:d)'; + P=c*(P.*.eye(m/(d+1),m/(d+1))) +endfunction diff --git a/modules/polynomials/macros/invr.bin b/modules/polynomials/macros/invr.bin Binary files differnew file mode 100755 index 000000000..b26737dc1 --- /dev/null +++ b/modules/polynomials/macros/invr.bin diff --git a/modules/polynomials/macros/invr.sci b/modules/polynomials/macros/invr.sci new file mode 100755 index 000000000..7df115487 --- /dev/null +++ b/modules/polynomials/macros/invr.sci @@ -0,0 +1,94 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function [f,d]=invr(h,flag) + //if h is a scalar, polynomial or rational fonction matrix, invr + //computes h^(-1). + //! + 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=max(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 + 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; + 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; +endfunction diff --git a/modules/polynomials/macros/lcm.bin b/modules/polynomials/macros/lcm.bin Binary files differnew file mode 100755 index 000000000..c6d7405f9 --- /dev/null +++ b/modules/polynomials/macros/lcm.bin diff --git a/modules/polynomials/macros/lcm.sci b/modules/polynomials/macros/lcm.sci new file mode 100755 index 000000000..f07900011 --- /dev/null +++ b/modules/polynomials/macros/lcm.sci @@ -0,0 +1,45 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function [p, fact] = lcm(p) + //p=lcm(p) computes the lcm of polynomial vector p + //[pp,fact]=lcm(p) computes besides the vector fact of factors + //such that p.*fact=pp*ones(p) + //! + + if type(p)<>1 & type(p)<>2 & type(p)<>8 then + error(msprintf(_("%s: Wrong type for argument #%d: Integer array or Polynomial expected.\n"), "lcm", 1)); + end + + if type(p)==1 then + if floor(p)<>p then + error(msprintf(_("%s: Wrong type for argument #%d: Integer array or Polynomial expected.\n"), "lcm", 1)); + else + p = iconvert(p,4); + end + end + + if type(p)==8 then + if argn(1)==2 then [p, fact] = %i_lcm(p), else p = %i_lcm(p), end + return + end + + [m, n] = size(p), + p = matrix(p, m*n, 1), + p0 = p(1); fact = 1; + for l = 2:m*n, + [u, v] = simp(p0, p(l)), + p0 = p0*v, + fact = [v*fact, u], + end, + fact = matrix(fact, m, n), + p = p0; + +endfunction diff --git a/modules/polynomials/macros/lcmdiag.bin b/modules/polynomials/macros/lcmdiag.bin Binary files differnew file mode 100755 index 000000000..0a9a822d9 --- /dev/null +++ b/modules/polynomials/macros/lcmdiag.bin diff --git a/modules/polynomials/macros/lcmdiag.sci b/modules/polynomials/macros/lcmdiag.sci new file mode 100755 index 000000000..28d82ed5a --- /dev/null +++ b/modules/polynomials/macros/lcmdiag.sci @@ -0,0 +1,49 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function [N,D]=lcmdiag(H,flag) + //returns N and diagonal D + //such that: + // flag='row' => H=D^(-1)*N and D(k,k)=lcm of kth row of H('den') + // flag='col' => H=N*D^(-1) and D(k,k)=lcm of kth col of H('den') + // default flag = 'col' + + if (type(H)<>2 & type(H)<>16) then + error(msprintf(_("%s: Wrong type for argument #%d: Real Matrix of Polynomials expected.\n"), "lcmdiag", 1)); + end + + Num = H("num"); + Den = H("den"); + [m, n] = size(H); + D = []; N = []; + [lhs, rhs] = argn(0); + if rhs==1 then + flag = "col"; + elseif type(flag)<>10 then + error(msprintf(_("%s: Wrong type for argument #%d: String expected.\n"), "lcmdiag", 2)); + end + select flag + case "row" + for k = 1:m + [pgcd, fact] = lcm(Den(k, :)); //H(k,:)=(Num(k,:).*fact)/pgcd + D = diag([diag(D); pgcd]); + N = [N; Num(k, :).*fact]; + end + case "col" + for k = 1:n + [pgcd, fact] = lcm(Den(:, k)); + D = diag([diag(D); pgcd]); + N = [N, Num(:, k).*fact]; + end + else + error(msprintf(_("%s: Wrong value for argument #%d: %s or %s expected.\n"), "lcmdiag", 2, "row", "col")); + end + +endfunction diff --git a/modules/polynomials/macros/lib b/modules/polynomials/macros/lib Binary files differnew file mode 100755 index 000000000..5d6472147 --- /dev/null +++ b/modules/polynomials/macros/lib diff --git a/modules/polynomials/macros/names b/modules/polynomials/macros/names new file mode 100755 index 000000000..a6315f673 --- /dev/null +++ b/modules/polynomials/macros/names @@ -0,0 +1,28 @@ +chepol +cmndred +coffg +colcompr +denom +derivat +determ +detr +diophant +factors +gcd +hermit +horner +hrmt +htrianr +inv_coeff +invr +lcm +lcmdiag +numer +pdiv +pfactors +pol2des +pol2str +polfact +rowcompr +sylm +systmat diff --git a/modules/polynomials/macros/numer.bin b/modules/polynomials/macros/numer.bin Binary files differnew file mode 100755 index 000000000..af0ca3fa6 --- /dev/null +++ b/modules/polynomials/macros/numer.bin diff --git a/modules/polynomials/macros/numer.sci b/modules/polynomials/macros/numer.sci new file mode 100755 index 000000000..299c1cbe7 --- /dev/null +++ b/modules/polynomials/macros/numer.sci @@ -0,0 +1,26 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function num=numer(r) + //returns the numerator num of a rational function matrix r (r may be + //also a scalar or polynomial matrix + //! + 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 + +endfunction diff --git a/modules/polynomials/macros/pdiv.bin b/modules/polynomials/macros/pdiv.bin Binary files differnew file mode 100755 index 000000000..fce62aff6 --- /dev/null +++ b/modules/polynomials/macros/pdiv.bin diff --git a/modules/polynomials/macros/pdiv.sci b/modules/polynomials/macros/pdiv.sci new file mode 100755 index 000000000..2c95716ec --- /dev/null +++ b/modules/polynomials/macros/pdiv.sci @@ -0,0 +1,42 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function [R,Q]=pdiv(P1,P2) + // Element wise euclidan division of a polynomial matrix + // by a polynomial + // This is just a loop for the primitive pppdiv + //! + [lhs,rhs]=argn(0); + [n,m]=size(P1); + [n1,m1]=size(P2); + // Special case for constant matrices + if type(P1)==1&type(P2)==1 then + Q=P1./P2;R=0*P1; + if lhs==1 then R=Q; end + return; + end + R=[],Q=[] + if n1==1 & m1==1 then + for l=1:n, + for k=1:m, + [rlk,qlk]=pppdiv(P1(l,k),P2),R(l,k)=rlk;Q(l,k)=qlk; + end; + end + if lhs==1 then R=Q;end + + return; + end + for l=1:n, + for k=1:m, + [rlk,qlk]=pppdiv(P1(l,k),P2(l,k)),R(l,k)=rlk;Q(l,k)=qlk; + end; + end + if lhs==1 then R=Q; end +endfunction diff --git a/modules/polynomials/macros/pfactors.bin b/modules/polynomials/macros/pfactors.bin Binary files differnew file mode 100755 index 000000000..0240ec103 --- /dev/null +++ b/modules/polynomials/macros/pfactors.bin diff --git a/modules/polynomials/macros/pfactors.sci b/modules/polynomials/macros/pfactors.sci new file mode 100755 index 000000000..a807bd940 --- /dev/null +++ b/modules/polynomials/macros/pfactors.sci @@ -0,0 +1,84 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function [resn,g]=pfactors(pol,flag) + // Given polynomial pol returns in list resn polynomials of + // degree 1 or two which are the factors of pol. + // One has pol= g times product of entries of resn + // if flag=='c' unstable roots are reflected vs the imaginary axis + // if flag=='d' unstable roots are reflected vs unit circle + w=roots(pol); + n=size(w,"*"); + if n==0 then resn=list();g=coeff(pol);return;end + co=coeff(pol);g=co(n+1); + resn=list(); + [LHS,RHS]=argn(0); + if RHS==1 then flag=[];end + if flag==[] then RHS=1;end + if RHS==1 then + kk=1;k=1; + while %T + if abs(imag(w(kk)))<=%eps then + resn(k)=poly(w(kk),varn(pol)); + kk=kk+1;k=k+1; + if kk>n then return;end + end + if abs(imag(w(kk)))>%eps then + resn(k)=real(poly([w(kk),w(kk+1)],varn(pol))); + kk=kk+2;k=k+1; + if kk>n then return;end + end + end + end //RHS=1 + if RHS==2 then + kk=1;k=1; + if flag=="c" then + while %T + if abs(imag(w(kk)))<=%eps then + resn(k)=poly(-abs(w(kk)),varn(pol)); + kk=kk+1;k=k+1; + if kk>n then return;end + end + if abs(imag(w(kk)))>%eps then + if real(w(kk))<0 then + resn(k)=real(poly([w(kk),w(kk+1)],varn(pol))); + else ; + resn(k)=real(poly([-w(kk),-w(kk+1)],varn(pol))); + end + kk=kk+2;k=k+1; + if kk>n then return;end + end + end + end //'c' + if flag=="d" then + while %T + wkk=w(kk); + if abs(imag(wkk))<=%eps then + [themin,which]=min([abs(wkk),1/(abs(wkk))]); + if which==2 then g=-g*real(wkk);end + resn(k)=poly(sign(real(wkk))*themin,varn(pol)); + kk=kk+1;k=k+1; + if kk>n then return;end + end + if abs(imag(wkk))>%eps then + if abs(wkk)<1 then + resn(k)=real(poly([wkk,w(kk+1)],varn(pol))); + else ; + // g=g*wkk*w(kk+1); w(kk+1)= conj(wkk) + g=g*abs(wkk)^2; + zp=[wkk,w(kk+1)];resn(k)=real(poly(ones(zp)./zp,varn(pol))); + end + kk=kk+2;k=k+1; + if kk>n then return;end + end + end + end //'d' + end //RHS=2 +endfunction diff --git a/modules/polynomials/macros/pol2des.bin b/modules/polynomials/macros/pol2des.bin Binary files differnew file mode 100755 index 000000000..7a46e002f --- /dev/null +++ b/modules/polynomials/macros/pol2des.bin diff --git a/modules/polynomials/macros/pol2des.sci b/modules/polynomials/macros/pol2des.sci new file mode 100755 index 000000000..699f1af1e --- /dev/null +++ b/modules/polynomials/macros/pol2des.sci @@ -0,0 +1,21 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function [N,B,C]=pol2des(Ds) + // Given the polynomial matrix Ds= D_0 +D_1 s +D_2 s^2 + ... +D_k s^k, + // pol2des returns three matrices N,B,C (with N nilpotent) such that + // Ds = C (sN-Eye)^-1 B + //! + if type(Ds)==1 then Ds=Ds+0*poly(0,"s")*Ds;end + dg=max(degree(Ds))+1; + [nout,nin]=size(Ds); + [Sl]=markp2ss(coeff(Ds),dg,nout,nin); + N=Sl(2);B=-Sl(3);C=Sl(4) +endfunction diff --git a/modules/polynomials/macros/pol2str.bin b/modules/polynomials/macros/pol2str.bin Binary files differnew file mode 100755 index 000000000..eeeaadc4d --- /dev/null +++ b/modules/polynomials/macros/pol2str.bin diff --git a/modules/polynomials/macros/pol2str.sci b/modules/polynomials/macros/pol2str.sci new file mode 100755 index 000000000..05b47fcd4 --- /dev/null +++ b/modules/polynomials/macros/pol2str.sci @@ -0,0 +1,56 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) 2013 - Samuel GOUGEON +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + +function S = pol2str(p) + if type(p) <> 2 & (typeof(p)~="hypermat" | type(p(:))~=2) then + msg = _("%s: Wrong type for input argument #%d: Polynomial expected.\n") + error(msprintf(msg, "pol2str", 1)); + end + sp = size(p) + p = p(:) + d = max(degree(p)) + c = coeff(p) + cr = real(c) + ci = imag(c) + S = stripblanks(string(c)) + k = find(c==0) + S(k) = "" + k = find(cr~=0 & ci~=0) + k2 = find(k>length(p)) + S(k(k2)) = "+("+S(k(k2))+")" + k = find((cr>0 & ci==0) | (cr==0 & ci>0)) + S(k) = "+"+S(k) + vn = varn(p) + clear ci cr k2 + for j = 2:(d+1) + if j==2 then + e = vn + else + e = vn+msprintf("^%d",j-1) + end + km1 = find(c(:,j)==-1) + if km1~=[] + S(km1,j) = "-"+e + end + kp1 = find(c(:,j)==1) + if kp1~=[] + S(kp1,j) = "+"+e + end + ko = find(c(:,j)~=-1 & c(:,j)~=0 & c(:,j)~=1) + if ko~=[] + S(ko,j) = S(ko,j)+"*"+e + end + end + S = strcat(S,"","c") + k = find(part(S,1)=="+") + S(k) = part(S(k),2:$) + k = find(S=="") + S(k) = "0" + S = matrix(S,sp) +endfunction diff --git a/modules/polynomials/macros/polfact.bin b/modules/polynomials/macros/polfact.bin Binary files differnew file mode 100755 index 000000000..82284998b --- /dev/null +++ b/modules/polynomials/macros/polfact.bin diff --git a/modules/polynomials/macros/polfact.sci b/modules/polynomials/macros/polfact.sci new file mode 100755 index 000000000..225c7e852 --- /dev/null +++ b/modules/polynomials/macros/polfact.sci @@ -0,0 +1,46 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function [f]=polfact(p) + // Minmal factors of p + // f=polfact(p) + // + // p : polynomial + // f : vector [f0 f1 ... fn] such that p=prod(f) + // - f0 constant + // - fi polynomial + //! + // + 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]=gsort(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 +endfunction diff --git a/modules/polynomials/macros/rowcompr.bin b/modules/polynomials/macros/rowcompr.bin Binary files differnew file mode 100755 index 000000000..986a5e861 --- /dev/null +++ b/modules/polynomials/macros/rowcompr.bin diff --git a/modules/polynomials/macros/rowcompr.sci b/modules/polynomials/macros/rowcompr.sci new file mode 100755 index 000000000..2acdfd125 --- /dev/null +++ b/modules/polynomials/macros/rowcompr.sci @@ -0,0 +1,25 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function [X,rk,Ac]=rowcompr(A) + //[X,rk,Ac]=rowcompr(A) + //row compression of polynomial matrix A (toward the bottom) + //X left polynomial unimodular base + //rk=normal rank of A + //Ac = X*A + //Warning: elimination of neglected terms + //! + [n,m]=size(A); + [Ac,U,rk]=htrianr(A'); + Ac=Ac'; + X=U'; + X=X(n:-1:1,:) + Ac=Ac(n:-1:1,:) +endfunction diff --git a/modules/polynomials/macros/sylm.bin b/modules/polynomials/macros/sylm.bin Binary files differnew file mode 100755 index 000000000..d40fa8643 --- /dev/null +++ b/modules/polynomials/macros/sylm.bin diff --git a/modules/polynomials/macros/sylm.sci b/modules/polynomials/macros/sylm.sci new file mode 100755 index 000000000..69ecc64b0 --- /dev/null +++ b/modules/polynomials/macros/sylm.sci @@ -0,0 +1,29 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function [s]=sylm(a,b) + //[s]=sylm(a,b) gives the Sylvester matrix associated to polynomials + //a and b, i.e. the matrix s such that: + // coeff( a*x + b*y )' = s * [coeff(x)';coeff(y)'] + //dimension of s is equal to degree(a)+degree(b) + //If a and b are coprime polynomials + //(rank(sylm(a,b))=degree(a)+degree(b)) the instructions + // u = sylm(a,b) \ eye(na+nb,1) + // x = poly(u(1:nb),'z','coeff') + // y = poly(u(nb+1:na+nb),'z','coeff') + //compute Bezout factors x et y of minimal degree de degre minimal + //such that a*x+b*y=1 + //! + na=degree(a);a=coeff(a)'; + nb=degree(b);b=coeff(b)'; + s(na+nb,na+nb)=0; + for i=1:nb,s(i:na+i,i)=a,end + for i=1:na,s(i:nb+i,nb+i)=b,end +endfunction diff --git a/modules/polynomials/macros/systmat.bin b/modules/polynomials/macros/systmat.bin Binary files differnew file mode 100755 index 000000000..3e1f16562 --- /dev/null +++ b/modules/polynomials/macros/systmat.bin diff --git a/modules/polynomials/macros/systmat.sci b/modules/polynomials/macros/systmat.sci new file mode 100755 index 000000000..123a16b03 --- /dev/null +++ b/modules/polynomials/macros/systmat.sci @@ -0,0 +1,34 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +function [Sm]=systmat(Sl); + // System matrix of the linear system Sl (syslin list) + // in state-space form. + // Sm = [-sE + A B; + // [ C D] + // To get the zeros use det or detr (for square systems) + //! + ty=Sl(1); + ty=ty(1); + if ty=="lss" then + if Sl(7)=="d" then + s=poly(0,"z"); + else + s=poly(0,"s"); + end + Sm=[-s*eye(Sl(2))+Sl(2),Sl(3);Sl(4),Sl(5)]; + return + end + if part(ty,1)=="d" then + s=poly(0,"s"); + Sm=[-s*Sl(6)+Sl(2),Sl(3);Sl(4),Sl(5)]; + return + end +endfunction |