octave_packages/quaternion-2.0.0/@quaternion/mpower.m

   1 ## Copyright (C) 2010, 2011   Lukas F. Reichlin
   2 ##
   3 ## This program is free software: you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation, either version 3 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program.  If not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## Matrix power operator of quaternions.  Used by Octave for "q^x".
  18
  19 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
  20 ## Created: May 2010
  21 ## Version: 0.2
  22
  23 function q = mpower (a, b)
  24
  25   [r, c] = size (a);
  26
  27   if (r != c)
  28     error ("quaternion: mpower: quaternion matrix must be square");
  29   endif
  30
  31   if (r == 1 && c == 1)                 # a scalar, b?
  32     q = a .^ b;                         # b could be a quaternion
  33   elseif (is_real_array (b) && isscalar (b) && fix (b) == b)
  34     e = fix (abs (b));
  35     switch (sign (b))
  36       case -1                           # q^-e
  37         a = inv (a);
  38         q = a;
  39       case 0                            # q^0
  40         q = eye (r);                    # alternative: q = quaternion (eye (r))
  41         return;
  42       case 1;                           # q^e
  43         q = a;
  44     endswitch
  45     for k = 2 : e
  46       q *= a;                           # improvement?: q^8 = ((q^2)^2)^2, q^9 = (((q^2)^2)^2)*q
  47     endfor
  48   else
  49     error ("quaternion: mpower: case not implemented yet");
  50     q = expm (logm (a) * b);            # don't know whether this formula is correct
  51   endif
  52
  53   ## TODO: - q1 ^ q2
  54   ##       - arrays
  55
  56 endfunction