Better fix after https://codereview.scilab.org/#/c/14828/
[scilab.git] / scilab / modules / overloading / macros / %s_pow.sci
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.1-en.txt
9
10 function x=%s_pow(a,p)
11     //   g_pow - A^p special cases
12     //%CALLING SEQUENCE
13     //   X=g_pow(A,p)
14     //%PARAMETERS
15     //   A   : constant or square hermitian or diagonalizable matrix
16     //   X   : square matrix
17     //   p   : square matrix or scalar
18     //%DESCRIPTION
19     //This function is called by the operation ^ to compute A^p in special cases
20     // - A scalar and p square matrix
21     // - A square matrix p is not an integer
22     //!
23     [m,n]=size(a)
24     [mp,np]=size(p)
25     if m*n==1&mp==np then  //a^P
26         flag=or(p<>p')
27         r=and(imag(p)==0)&imag(a)==0
28         if ~flag then
29             //Hermitian matrix
30             [u,s]=schur(p);
31             w=a.^diag(s);
32             x=u*diag(a.^diag(s))*u';
33             if r then
34                 x=real(x)
35             end
36         else
37             [s,u,bs]=bdiag(p+0*%i);
38             if max(bs)>1 then
39                 error(msprintf(_("%s: Unable to diagonalize.\n"),"%s_pow"));
40             end
41             w=diag(s);
42             x=u*diag(a.^diag(s))*inv(u);
43         end
44         if r then x=real(x), end
45     elseif m==n&mp*np==1 then  //A^p  p non integer
46         flag=or(a<>a')
47         if ~flag then
48             //Hermitian matrix
49             r=and(imag(a)==0)
50             [u,s]=schur(a);
51             x=u*diag(diag(s).^p)*u';
52             if r then
53                 if s>=0&imag(p)==0 then
54                     x=real(x)
55                 end
56             end
57         else
58             //General matrix
59             r=and(imag(a)==0)
60             [s,u,bs]=bdiag(a+0*%i);
61             if max(bs)>1 then
62                 error(msprintf(_("%s: Unable to diagonalize.\n"),"%s_pow"));
63             end
64             x=u*diag(diag(s).^p)*inv(u);
65         end
66         if int(p)==p & real(p)==p & r then
67             x=real(x);
68         end
69     else
70         error(43)
71     end
72 endfunction