[Octave-cvsupdate] SF.net SVN: octave:[9375] trunk/octave-forge/extra/quaternion_oo/inst/@ quaternion/mpower.m

[<par...@us...>]

Revision: 9375
          http://octave.svn.sourceforge.net/octave/?rev=9375&view=rev
Author:   paramaniac
Date:     2011-12-13 15:34:25 +0000 (Tue, 13 Dec 2011)
Log Message:
-----------
quaternion: beef up mpower

Modified Paths:
--------------
    trunk/octave-forge/extra/quaternion_oo/inst/@quaternion/mpower.m

Modified: trunk/octave-forge/extra/quaternion_oo/inst/@quaternion/mpower.m
===================================================================
--- trunk/octave-forge/extra/quaternion_oo/inst/@quaternion/mpower.m	2011-12-13 10:33:03 UTC (rev 9374)
+++ trunk/octave-forge/extra/quaternion_oo/inst/@quaternion/mpower.m	2011-12-13 15:34:25 UTC (rev 9375)
@@ -1,4 +1,4 @@
-## Copyright (C) 2010   Lukas F. Reichlin
+## Copyright (C) 2010, 2011   Lukas F. Reichlin
 ##
 ## This program is free software: you can redistribute it and/or modify
 ## it under the terms of the GNU General Public License as published by
@@ -18,36 +18,38 @@
 
 ## Author: Lukas Reichlin <luk...@gm...>
 ## Created: May 2010
-## Version: 0.1
+## Version: 0.2
 
-function a = mpower (a, b)
+function q = mpower (a, b)
 
-  if isa (a, "quaternion")
+  [r, c] = size (a);
+  
+  if (r != c)
+    error ("quaternion: mpower: quaternion matrix must be square");
+  endif
+  
+  if (r == 1 && c == 1)                 # a scalar, b?
+    q = a .^ b;                         # b could be a quaternion
+  elseif (is_real_array (b) && isscalar (b) && fix (b) == b)
+    e = fix (abs (b));
+    switch (sign (b))
+      case -1                           # q^-e
+        a = inv (a);
+        q = a;
+      case 0                            # q^0
+        q = eye (r);                    # alternative: q = quaternion (eye (r))
+        return;
+      case 1;                           # q^e
+        q = a;
+    endswitch  
+    for k = 2 : e
+      q *= a;                           # improvement?: q^8 = ((q^2)^2)^2, q^9 = (((q^2)^2)^2)*q
+    endfor
+  else
+    error ("quaternion: mpower: case not implemented yet");
+    q = expm (logm (a) * b);            # don't know whether this formula is correct
+  endif
 
-    if isscalar (a.w)                   # power takes care of checking type of b
-      a = a .^ b;
-    else
-
-      if fix (b) == b                    # only integers are poorly implemented
-        [n m] = size (a.w);
-        w = eye (n,m);
-        x = zeros (n,m);
-        q = quaternion(w,x,x,x);
-        while b--
-          q *= a;
-        end
-        a = q;
-      else
-        error ("quaternion:devel", ...
-                  "quaternion: power: implemented for integer exponents only");
-        a = expm(b * logm (a));
-      endif
-
-    end
-
-  else
-    error ("quaternion:invalidArgument", "base must be a quaternion.");
-  end
   ## TODO: - q1 ^ q2
   ##       - arrays
 

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.