1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) INRIA
3 //
4 // This file must be used under the terms of the CeCILL.
5 // This source file is licensed as described in the file COPYING, which
6 // you should have received as part of this distribution.  The terms
7 // are also available at
8 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10 function f=%p_j_s(p,s)
11 // %p_j_s(p,s)  computes p.^s for p polynomial matrix in special cases
12 //!
14 if s==[] then f=[],return,end
15 if  or(imag(s)<>0)|or(int(s)<>s) then error('%p_j_s: integer power only'),end
16 [m,n]=size(p)
17 [ms,ns]=size(s)
18 if ms==1&ns==1 then
19   if s<0 then
20     if or(abs(coeff(p(:)))*ones(max(degree(p))+1,1)==0) then
21       error(27)
22     end
23     f=rlist(ones(p),p.^(-s),[])
24   else // this case is in fact hard coded
25     f=p.^s
26   end
27 elseif m==1&n==1 then // Element wise exponentiation p.^s with p "scalar"
28   kp=find(s>=0)
29   kn=find(s<0)
30   num=ones(s)
31   den=ones(s)
32   num(kp)=p.^s(kp)
33   p=1/p
34   num(kn)=p(2).^(-s(kn))
35   den(kn)=p(3).^(-s(kn))
36   f=rlist(num,den,[])
37 elseif ms==m&ns==n then  // Element wise exponentiation
38   p=p(:);s=s(:);
39   kp=find(s>=0)
40   kn=find(s<0)
41   num=p
42   den=ones(s)
43   num(kp)=num(kp).^s(kp)
44   if or(abs(coeff(p(kn)))*ones(max(degree(p(kn)))+1,1)==0) then
45     error(27)
46   end
47   num(kn)=ones(p(kn))
48   den(kn)=p(kn).^(-s(kn))
49   f = rlist(matrix(num,n,m),matrix(den,n,m),[])
50 else
51   error(30)
52 end
54 endfunction