--- a/src/xpow.cc
+++ b/src/xpow.cc
@@ -33,6 +33,7 @@
 #include "CDiagMatrix.h"
 #include "CMatrix.h"
 #include "EIG.h"
+#include "fEIG.h"
 #include "dDiagMatrix.h"
 #include "dMatrix.h"
 #include "mx-cm-cdm.h"
@@ -1260,6 +1261,1218 @@
   return result;
 }
 
+static inline int
+xisint (float x)
+{
+  return (D_NINT (x) == x
+	  && ((x >= 0 && x < INT_MAX)
+	      || (x <= 0 && x > INT_MIN)));
+}
+
+// Safer pow functions.
+//
+//       op2 \ op1:   s   m   cs   cm
+//            +--   +---+---+----+----+
+//   scalar   |     | 1 | 5 |  7 | 11 |
+//                  +---+---+----+----+
+//   matrix         | 2 | * |  8 |  * |
+//                  +---+---+----+----+
+//   complex_scalar | 3 | 6 |  9 | 12 |
+//                  +---+---+----+----+
+//   complex_matrix | 4 | * | 10 |  * |
+//                  +---+---+----+----+
+
+// -*- 1 -*-
+octave_value
+xpow (float a, float b)
+{
+  float retval;
+
+  if (a < 0.0 && static_cast<int> (b) != b)
+    {
+      FloatComplex atmp (a);
+
+      return std::pow (atmp, b);
+    }
+  else
+    retval = std::pow (a, b);
+
+  return retval;
+}
+
+// -*- 2 -*-
+octave_value
+xpow (float a, const FloatMatrix& b)
+{
+  octave_value retval;
+
+  octave_idx_type nr = b.rows ();
+  octave_idx_type nc = b.cols ();
+
+  if (nr == 0 || nc == 0 || nr != nc)
+    error ("for x^A, A must be square");
+  else
+    {
+      FloatEIG b_eig (b);
+
+      if (! error_state)
+	{
+	  FloatComplexColumnVector lambda (b_eig.eigenvalues ());
+	  FloatComplexMatrix Q (b_eig.eigenvectors ());
+
+	  for (octave_idx_type i = 0; i < nr; i++)
+	    {
+	      FloatComplex elt = lambda(i);
+	      if (std::imag (elt) == 0.0)
+		lambda(i) = std::pow (a, std::real (elt));
+	      else
+		lambda(i) = std::pow (a, elt);
+	    }
+	  FloatComplexDiagMatrix D (lambda);
+
+	  retval = FloatComplexMatrix (Q * D * Q.inverse ());
+	}
+      else
+	error ("xpow: matrix diagonalization failed");
+    }
+
+  return retval;
+}
+
+// -*- 3 -*-
+octave_value
+xpow (float a, const FloatComplex& b)
+{
+  FloatComplex result;
+  FloatComplex atmp (a);
+  result = std::pow (atmp, b);
+  return result;
+}
+
+// -*- 4 -*-
+octave_value
+xpow (float a, const FloatComplexMatrix& b)
+{
+  octave_value retval;
+
+  octave_idx_type nr = b.rows ();
+  octave_idx_type nc = b.cols ();
+
+  if (nr == 0 || nc == 0 || nr != nc)
+    error ("for x^A, A must be square");
+  else
+    {
+      FloatEIG b_eig (b);
+
+      if (! error_state)
+	{
+	  FloatComplexColumnVector lambda (b_eig.eigenvalues ());
+	  FloatComplexMatrix Q (b_eig.eigenvectors ());
+
+	  for (octave_idx_type i = 0; i < nr; i++)
+	    {
+	      FloatComplex elt = lambda(i);
+	      if (std::imag (elt) == 0.0)
+		lambda(i) = std::pow (a, std::real (elt));
+	      else
+		lambda(i) = std::pow (a, elt);
+	    }
+	  FloatComplexDiagMatrix D (lambda);
+
+	  retval = FloatComplexMatrix (Q * D * Q.inverse ());
+	}
+      else
+	error ("xpow: matrix diagonalization failed");
+    }
+
+  return retval;
+}
+
+// -*- 5 -*-
+octave_value
+xpow (const FloatMatrix& a, float b)
+{
+  octave_value retval;
+
+  octave_idx_type nr = a.rows ();
+  octave_idx_type nc = a.cols ();
+
+  if (nr == 0 || nc == 0 || nr != nc)
+    error ("for A^b, A must be square");
+  else
+    {
+      if (static_cast<int> (b) == b)
+	{
+	  int btmp = static_cast<int> (b);
+	  if (btmp == 0)
+	    {
+	      retval = FloatDiagMatrix (nr, nr, 1.0);
+	    }
+	  else
+	    {
+	      // Too much copying?
+	      // FIXME -- we shouldn't do this if the exponent is
+	      // large...
+
+	      FloatMatrix atmp;
+	      if (btmp < 0)
+		{
+		  btmp = -btmp;
+
+		  octave_idx_type info;
+		  float rcond = 0.0;
+		  MatrixType mattype (a);
+
+		  atmp = a.inverse (mattype, info, rcond, 1);
+
+		  if (info == -1)
+		    warning ("inverse: matrix singular to machine\
+ precision, rcond = %g", rcond);
+		}
+	      else
+		atmp = a;
+
+	      FloatMatrix result (atmp);
+
+	      btmp--;
+
+	      while (btmp > 0)
+		{
+		  if (btmp & 1)
+		    result = result * atmp;
+
+		  btmp >>= 1;
+
+		  if (btmp > 0)
+		    atmp = atmp * atmp;
+		}
+
+	      retval = result;
+	    }
+	}
+      else
+	{
+	  FloatEIG a_eig (a);
+
+	  if (! error_state)
+	    {
+	      FloatComplexColumnVector lambda (a_eig.eigenvalues ());
+	      FloatComplexMatrix Q (a_eig.eigenvectors ());
+
+	      for (octave_idx_type i = 0; i < nr; i++)
+		lambda(i) = std::pow (lambda(i), b);
+
+	      FloatComplexDiagMatrix D (lambda);
+
+	      retval = FloatComplexMatrix (Q * D * Q.inverse ());
+	    }
+	  else
+	    error ("xpow: matrix diagonalization failed");
+	}
+    }
+
+  return retval;
+}
+
+// -*- 6 -*-
+octave_value
+xpow (const FloatMatrix& a, const FloatComplex& b)
+{
+  octave_value retval;
+
+  octave_idx_type nr = a.rows ();
+  octave_idx_type nc = a.cols ();
+
+  if (nr == 0 || nc == 0 || nr != nc)
+    error ("for A^b, A must be square");
+  else
+    {
+      FloatEIG a_eig (a);
+
+      if (! error_state)
+	{
+	  FloatComplexColumnVector lambda (a_eig.eigenvalues ());
+	  FloatComplexMatrix Q (a_eig.eigenvectors ());
+
+	  for (octave_idx_type i = 0; i < nr; i++)
+	    lambda(i) = std::pow (lambda(i), b);
+
+	  FloatComplexDiagMatrix D (lambda);
+
+	  retval = FloatComplexMatrix (Q * D * Q.inverse ());
+	}
+      else
+	error ("xpow: matrix diagonalization failed");
+    }
+
+  return retval;
+}
+
+// -*- 7 -*-
+octave_value
+xpow (const FloatComplex& a, float b)
+{
+  FloatComplex result;
+
+  if (xisint (b))
+    result = std::pow (a, static_cast<int> (b));
+  else
+    result = std::pow (a, b);
+
+  return result;
+}
+
+// -*- 8 -*-
+octave_value
+xpow (const FloatComplex& a, const FloatMatrix& b)
+{
+  octave_value retval;
+
+  octave_idx_type nr = b.rows ();
+  octave_idx_type nc = b.cols ();
+
+  if (nr == 0 || nc == 0 || nr != nc)
+    error ("for x^A, A must be square");
+  else
+    {
+      FloatEIG b_eig (b);
+
+      if (! error_state)
+	{
+	  FloatComplexColumnVector lambda (b_eig.eigenvalues ());
+	  FloatComplexMatrix Q (b_eig.eigenvectors ());
+
+	  for (octave_idx_type i = 0; i < nr; i++)
+	    {
+	      FloatComplex elt = lambda(i);
+	      if (std::imag (elt) == 0.0)
+		lambda(i) = std::pow (a, std::real (elt));
+	      else
+		lambda(i) = std::pow (a, elt);
+	    }
+	  FloatComplexDiagMatrix D (lambda);
+
+	  retval = FloatComplexMatrix (Q * D * Q.inverse ());
+	}
+      else
+	error ("xpow: matrix diagonalization failed");
+    }
+
+  return retval;
+}
+
+// -*- 9 -*-
+octave_value
+xpow (const FloatComplex& a, const FloatComplex& b)
+{
+  FloatComplex result;
+  result = std::pow (a, b);
+  return result;
+}
+
+// -*- 10 -*-
+octave_value
+xpow (const FloatComplex& a, const FloatComplexMatrix& b)
+{
+  octave_value retval;
+
+  octave_idx_type nr = b.rows ();
+  octave_idx_type nc = b.cols ();
+
+  if (nr == 0 || nc == 0 || nr != nc)
+    error ("for x^A, A must be square");
+  else
+    {
+      FloatEIG b_eig (b);
+
+      if (! error_state)
+	{
+	  FloatComplexColumnVector lambda (b_eig.eigenvalues ());
+	  FloatComplexMatrix Q (b_eig.eigenvectors ());
+
+	  for (octave_idx_type i = 0; i < nr; i++)
+	    {
+	      FloatComplex elt = lambda(i);
+	      if (std::imag (elt) == 0.0)
+		lambda(i) = std::pow (a, std::real (elt));
+	      else
+		lambda(i) = std::pow (a, elt);
+	    }
+	  FloatComplexDiagMatrix D (lambda);
+
+	  retval = FloatComplexMatrix (Q * D * Q.inverse ());
+	}
+      else
+	error ("xpow: matrix diagonalization failed");
+    }
+
+  return retval;
+}
+
+// -*- 11 -*-
+octave_value
+xpow (const FloatComplexMatrix& a, float b)
+{
+  octave_value retval;
+
+  octave_idx_type nr = a.rows ();
+  octave_idx_type nc = a.cols ();
+
+  if (nr == 0 || nc == 0 || nr != nc)
+    error ("for A^b, A must be square");
+  else
+    {
+      if (static_cast<int> (b) == b)
+	{
+	  int btmp = static_cast<int> (b);
+	  if (btmp == 0)
+	    {
+	      retval = FloatDiagMatrix (nr, nr, 1.0);
+	    }
+	  else
+	    {
+	      // Too much copying?
+	      // FIXME -- we shouldn't do this if the exponent is
+	      // large...
+
+	      FloatComplexMatrix atmp;
+	      if (btmp < 0)
+		{
+		  btmp = -btmp;
+
+		  octave_idx_type info;
+		  float rcond = 0.0;
+		  MatrixType mattype (a);
+
+		  atmp = a.inverse (mattype, info, rcond, 1);
+
+		  if (info == -1)
+		    warning ("inverse: matrix singular to machine\
+ precision, rcond = %g", rcond);
+		}
+	      else
+		atmp = a;
+
+	      FloatComplexMatrix result (atmp);
+
+	      btmp--;
+
+	      while (btmp > 0)
+		{
+		  if (btmp & 1)
+		    result = result * atmp;
+
+		  btmp >>= 1;
+
+		  if (btmp > 0)
+		    atmp = atmp * atmp;
+		}
+
+	      retval = result;
+	    }
+	}
+      else
+	{
+	  FloatEIG a_eig (a);
+
+	  if (! error_state)
+	    {
+	      FloatComplexColumnVector lambda (a_eig.eigenvalues ());
+	      FloatComplexMatrix Q (a_eig.eigenvectors ());
+
+	      for (octave_idx_type i = 0; i < nr; i++)
+		lambda(i) = std::pow (lambda(i), b);
+
+	      FloatComplexDiagMatrix D (lambda);
+
+	      retval = FloatComplexMatrix (Q * D * Q.inverse ());
+	    }
+	  else
+	    error ("xpow: matrix diagonalization failed");
+	}
+    }
+
+  return retval;
+}
+
+// -*- 12 -*-
+octave_value
+xpow (const FloatComplexMatrix& a, const FloatComplex& b)
+{
+  octave_value retval;
+
+  octave_idx_type nr = a.rows ();
+  octave_idx_type nc = a.cols ();
+
+  if (nr == 0 || nc == 0 || nr != nc)
+    error ("for A^b, A must be square");
+  else
+    {
+      FloatEIG a_eig (a);
+
+      if (! error_state)
+	{
+	  FloatComplexColumnVector lambda (a_eig.eigenvalues ());
+	  FloatComplexMatrix Q (a_eig.eigenvectors ());
+
+	  for (octave_idx_type i = 0; i < nr; i++)
+	    lambda(i) = std::pow (lambda(i), b);
+
+	  FloatComplexDiagMatrix D (lambda);
+
+	  retval = FloatComplexMatrix (Q * D * Q.inverse ());
+	}
+      else
+	error ("xpow: matrix diagonalization failed");
+    }
+
+  return retval;
+}
+
+// Safer pow functions that work elementwise for matrices.
+//
+//       op2 \ op1:   s   m   cs   cm
+//            +--   +---+---+----+----+
+//   scalar   |     | * | 3 |  * |  9 |
+//                  +---+---+----+----+
+//   matrix         | 1 | 4 |  7 | 10 |
+//                  +---+---+----+----+
+//   complex_scalar | * | 5 |  * | 11 |
+//                  +---+---+----+----+
+//   complex_matrix | 2 | 6 |  8 | 12 |
+//                  +---+---+----+----+
+//
+//   * -> not needed.
+
+// FIXME -- these functions need to be fixed so that things
+// like
+//
+//   a = -1; b = [ 0, 0.5, 1 ]; r = a .^ b
+//
+// and
+//
+//   a = -1; b = [ 0, 0.5, 1 ]; for i = 1:3, r(i) = a .^ b(i), end
+//
+// produce identical results.  Also, it would be nice if -1^0.5
+// produced a pure imaginary result instead of a complex number with a
+// small real part.  But perhaps that's really a problem with the math
+// library...
+
+// -*- 1 -*-
+octave_value
+elem_xpow (float a, const FloatMatrix& b)
+{
+  octave_value retval;
+
+  octave_idx_type nr = b.rows ();
+  octave_idx_type nc = b.cols ();
+
+  float d1, d2;
+
+  if (a < 0.0 && ! b.all_integers (d1, d2))
+    {
+      FloatComplex atmp (a);
+      FloatComplexMatrix result (nr, nc);
+
+      for (octave_idx_type j = 0; j < nc; j++)
+	for (octave_idx_type i = 0; i < nr; i++)
+	  {
+	    OCTAVE_QUIT;
+	    result (i, j) = std::pow (atmp, b (i, j));
+	  }
+
+      retval = result;
+    }
+  else
+    {
+      FloatMatrix result (nr, nc);
+
+      for (octave_idx_type j = 0; j < nc; j++)
+	for (octave_idx_type i = 0; i < nr; i++)
+	  {
+	    OCTAVE_QUIT;
+	    result (i, j) = std::pow (a, b (i, j));
+	  }
+
+      retval = result;
+    }
+
+  return retval;
+}
+
+// -*- 2 -*-
+octave_value
+elem_xpow (float a, const FloatComplexMatrix& b)
+{
+  octave_idx_type nr = b.rows ();
+  octave_idx_type nc = b.cols ();
+
+  FloatComplexMatrix result (nr, nc);
+  FloatComplex atmp (a);
+
+  for (octave_idx_type j = 0; j < nc; j++)
+    for (octave_idx_type i = 0; i < nr; i++)
+      {
+	OCTAVE_QUIT;
+	result (i, j) = std::pow (atmp, b (i, j));
+      }
+
+  return result;
+}
+
+// -*- 3 -*-
+octave_value
+elem_xpow (const FloatMatrix& a, float b)
+{
+  octave_value retval;
+
+  octave_idx_type nr = a.rows ();
+  octave_idx_type nc = a.cols ();
+
+  if (static_cast<int> (b) != b && a.any_element_is_negative ())
+    {
+      FloatComplexMatrix result (nr, nc);
+
+      for (octave_idx_type j = 0; j < nc; j++)
+	for (octave_idx_type i = 0; i < nr; i++)
+	  {
+	    OCTAVE_QUIT; 
+      
+	    FloatComplex atmp (a (i, j));
+
+	    result (i, j) = std::pow (atmp, b);
+	  }
+
+      retval = result;
+    }
+  else
+    {
+      FloatMatrix result (nr, nc);
+
+      for (octave_idx_type j = 0; j < nc; j++)
+	for (octave_idx_type i = 0; i < nr; i++)
+	  {
+	    OCTAVE_QUIT;
+	    result (i, j) = std::pow (a (i, j), b);
+	  }
+
+      retval = result;
+    }
+
+  return retval;
+}
+
+// -*- 4 -*-
+octave_value
+elem_xpow (const FloatMatrix& a, const FloatMatrix& b)
+{
+  octave_value retval;
+
+  octave_idx_type nr = a.rows ();
+  octave_idx_type nc = a.cols ();
+
+  octave_idx_type b_nr = b.rows ();
+  octave_idx_type b_nc = b.cols ();
+
+  if (nr != b_nr || nc != b_nc)
+    {
+      gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
+      return octave_value ();
+    }
+
+  int convert_to_complex = 0;
+  for (octave_idx_type j = 0; j < nc; j++)
+    for (octave_idx_type i = 0; i < nr; i++)
+      {
+	OCTAVE_QUIT;
+	float atmp = a (i, j);
+	float btmp = b (i, j);
+	if (atmp < 0.0 && static_cast<int> (btmp) != btmp)
+	  {
+	    convert_to_complex = 1;
+	    goto done;
+	  }
+      }
+
+done:
+
+  if (convert_to_complex)
+    {
+      FloatComplexMatrix complex_result (nr, nc);
+
+      for (octave_idx_type j = 0; j < nc; j++)
+	for (octave_idx_type i = 0; i < nr; i++)
+	  {
+	    OCTAVE_QUIT;
+	    FloatComplex atmp (a (i, j));
+	    FloatComplex btmp (b (i, j));
+	    complex_result (i, j) = std::pow (atmp, btmp);
+	  }
+
+      retval = complex_result;
+    }
+  else
+    {
+      FloatMatrix result (nr, nc);
+
+      for (octave_idx_type j = 0; j < nc; j++)
+	for (octave_idx_type i = 0; i < nr; i++)
+	  {
+	    OCTAVE_QUIT;
+	    result (i, j) = std::pow (a (i, j), b (i, j));
+	  }
+
+      retval = result;
+    }
+
+  return retval;
+}
+
+// -*- 5 -*-
+octave_value
+elem_xpow (const FloatMatrix& a, const FloatComplex& b)
+{
+  octave_idx_type nr = a.rows ();
+  octave_idx_type nc = a.cols ();
+
+  FloatComplexMatrix result (nr, nc);
+
+  for (octave_idx_type j = 0; j < nc; j++)
+    for (octave_idx_type i = 0; i < nr; i++)
+      {
+	OCTAVE_QUIT;
+	result (i, j) = std::pow (FloatComplex (a (i, j)), b);
+      }
+
+  return result;
+}
+
+// -*- 6 -*-
+octave_value
+elem_xpow (const FloatMatrix& a, const FloatComplexMatrix& b)
+{
+  octave_idx_type nr = a.rows ();
+  octave_idx_type nc = a.cols ();
+
+  octave_idx_type b_nr = b.rows ();
+  octave_idx_type b_nc = b.cols ();
+
+  if (nr != b_nr || nc != b_nc)
+    {
+      gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
+      return octave_value ();
+    }
+
+  FloatComplexMatrix result (nr, nc);
+
+  for (octave_idx_type j = 0; j < nc; j++)
+    for (octave_idx_type i = 0; i < nr; i++)
+      {
+	OCTAVE_QUIT;
+	result (i, j) = std::pow (FloatComplex (a (i, j)), b (i, j));
+      }
+
+  return result;
+}
+
+// -*- 7 -*-
+octave_value
+elem_xpow (const FloatComplex& a, const FloatMatrix& b)
+{
+  octave_idx_type nr = b.rows ();
+  octave_idx_type nc = b.cols ();
+
+  FloatComplexMatrix result (nr, nc);
+
+  for (octave_idx_type j = 0; j < nc; j++)
+    for (octave_idx_type i = 0; i < nr; i++)
+      {
+	OCTAVE_QUIT;
+	float btmp = b (i, j);
+	if (xisint (btmp))
+	  result (i, j) = std::pow (a, static_cast<int> (btmp));
+	else
+	  result (i, j) = std::pow (a, btmp);
+      }
+
+  return result;
+}
+
+// -*- 8 -*-
+octave_value
+elem_xpow (const FloatComplex& a, const FloatComplexMatrix& b)
+{
+  octave_idx_type nr = b.rows ();
+  octave_idx_type nc = b.cols ();
+
+  FloatComplexMatrix result (nr, nc);
+
+  for (octave_idx_type j = 0; j < nc; j++)
+    for (octave_idx_type i = 0; i < nr; i++)
+      {
+	OCTAVE_QUIT;
+	result (i, j) = std::pow (a, b (i, j));
+      }
+
+  return result;
+}
+
+// -*- 9 -*-
+octave_value
+elem_xpow (const FloatComplexMatrix& a, float b)
+{
+  octave_idx_type nr = a.rows ();
+  octave_idx_type nc = a.cols ();
+
+  FloatComplexMatrix result (nr, nc);
+
+  if (xisint (b))
+    {
+      for (octave_idx_type j = 0; j < nc; j++)
+	for (octave_idx_type i = 0; i < nr; i++)
+	  {
+	    OCTAVE_QUIT;
+	    result (i, j) = std::pow (a (i, j), static_cast<int> (b));
+	  }
+    }
+  else
+    {
+      for (octave_idx_type j = 0; j < nc; j++)
+	for (octave_idx_type i = 0; i < nr; i++)
+	  {
+	    OCTAVE_QUIT;
+	    result (i, j) = std::pow (a (i, j), b);
+	  }
+    }
+
+  return result;
+}
+
+// -*- 10 -*-
+octave_value
+elem_xpow (const FloatComplexMatrix& a, const FloatMatrix& b)
+{
+  octave_idx_type nr = a.rows ();
+  octave_idx_type nc = a.cols ();
+
+  octave_idx_type b_nr = b.rows ();
+  octave_idx_type b_nc = b.cols ();
+
+  if (nr != b_nr || nc != b_nc)
+    {
+      gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
+      return octave_value ();
+    }
+
+  FloatComplexMatrix result (nr, nc);
+
+  for (octave_idx_type j = 0; j < nc; j++)
+    for (octave_idx_type i = 0; i < nr; i++)
+      {
+	OCTAVE_QUIT;
+	float btmp = b (i, j);
+	if (xisint (btmp))
+	  result (i, j) = std::pow (a (i, j), static_cast<int> (btmp));
+	else
+	  result (i, j) = std::pow (a (i, j), btmp);
+      }
+
+  return result;
+}
+
+// -*- 11 -*-
+octave_value
+elem_xpow (const FloatComplexMatrix& a, const FloatComplex& b)
+{
+  octave_idx_type nr = a.rows ();
+  octave_idx_type nc = a.cols ();
+
+  FloatComplexMatrix result (nr, nc);
+
+  for (octave_idx_type j = 0; j < nc; j++)
+    for (octave_idx_type i = 0; i < nr; i++)
+      {
+	OCTAVE_QUIT;
+	result (i, j) = std::pow (a (i, j), b);
+      }
+
+  return result;
+}
+
+// -*- 12 -*-
+octave_value
+elem_xpow (const FloatComplexMatrix& a, const FloatComplexMatrix& b)
+{
+  octave_idx_type nr = a.rows ();
+  octave_idx_type nc = a.cols ();
+
+  octave_idx_type b_nr = b.rows ();
+  octave_idx_type b_nc = b.cols ();
+
+  if (nr != b_nr || nc != b_nc)
+    {
+      gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
+      return octave_value ();
+    }
+
+  FloatComplexMatrix result (nr, nc);
+
+  for (octave_idx_type j = 0; j < nc; j++)
+    for (octave_idx_type i = 0; i < nr; i++)
+      {
+	OCTAVE_QUIT;
+	result (i, j) = std::pow (a (i, j), b (i, j));
+      }
+
+  return result;
+}
+
+// Safer pow functions that work elementwise for N-d arrays.
+//
+//       op2 \ op1:   s   nd  cs   cnd
+//            +--   +---+---+----+----+
+//   scalar   |     | * | 3 |  * |  9 |
+//                  +---+---+----+----+
+//   N_d            | 1 | 4 |  7 | 10 |
+//                  +---+---+----+----+
+//   complex_scalar | * | 5 |  * | 11 |
+//                  +---+---+----+----+
+//   complex_N_d    | 2 | 6 |  8 | 12 |
+//                  +---+---+----+----+
+//
+//   * -> not needed.
+
+// FIXME -- these functions need to be fixed so that things
+// like
+//
+//   a = -1; b = [ 0, 0.5, 1 ]; r = a .^ b
+//
+// and
+//
+//   a = -1; b = [ 0, 0.5, 1 ]; for i = 1:3, r(i) = a .^ b(i), end
+//
+// produce identical results.  Also, it would be nice if -1^0.5
+// produced a pure imaginary result instead of a complex number with a
+// small real part.  But perhaps that's really a problem with the math
+// library...
+
+// -*- 1 -*-
+octave_value
+elem_xpow (float a, const FloatNDArray& b)
+{
+  octave_value retval;
+
+  float d1, d2;
+
+  if (a < 0.0 && ! b.all_integers (d1, d2))
+    {
+      FloatComplex atmp (a);
+      FloatComplexNDArray result (b.dims ());
+      for (octave_idx_type i = 0; i < b.length (); i++)
+	{
+	  OCTAVE_QUIT;
+	  result(i) = std::pow (atmp, b(i));
+	}
+
+      retval = result;
+    }
+  else
+    {
+      FloatNDArray result (b.dims ());
+      for (octave_idx_type i = 0; i < b.length (); i++)
+	{
+	  OCTAVE_QUIT;
+	  result (i) = std::pow (a, b(i));
+	}
+
+      retval = result;
+    }
+
+  return retval;
+}
+
+// -*- 2 -*-
+octave_value
+elem_xpow (float a, const FloatComplexNDArray& b)
+{
+  FloatComplexNDArray result (b.dims ());
+  FloatComplex atmp (a);
+
+  for (octave_idx_type i = 0; i < b.length (); i++)
+    {
+      OCTAVE_QUIT;
+      result(i) = std::pow (atmp, b(i));
+    }
+
+  return result;
+}
+
+// -*- 3 -*-
+octave_value
+elem_xpow (const FloatNDArray& a, float b)
+{
+  octave_value retval;
+
+  if (static_cast<int> (b) != b && a.any_element_is_negative ())
+    {
+      FloatComplexNDArray result (a.dims ());
+
+      for (octave_idx_type i = 0; i < a.length (); i++)
+	{
+	  OCTAVE_QUIT;
+
+	  FloatComplex atmp (a (i));
+
+	  result(i) = std::pow (atmp, b);
+	}
+
+      retval = result;
+    }
+  else
+    {
+      FloatNDArray result (a.dims ());
+
+      for (octave_idx_type i = 0; i < a.length (); i++)
+	{
+	  OCTAVE_QUIT;
+	  result(i) = std::pow (a(i), b);
+	}
+
+      retval = result;
+    }
+
+  return retval;
+}
+
+// -*- 4 -*-
+octave_value
+elem_xpow (const FloatNDArray& a, const FloatNDArray& b)
+{
+  octave_value retval;
+
+  dim_vector a_dims = a.dims ();
+  dim_vector b_dims = b.dims ();
+
+  if (a_dims != b_dims)
+    {
+      gripe_nonconformant ("operator .^", a_dims, b_dims);
+      return octave_value ();
+    }
+
+  int len = a.length ();
+
+  bool convert_to_complex = false;
+
+  for (octave_idx_type i = 0; i < len; i++)
+    {
+      OCTAVE_QUIT;
+      float atmp = a(i);
+      float btmp = b(i);
+      if (atmp < 0.0 && static_cast<int> (btmp) != btmp)
+	{
+	  convert_to_complex = true;
+	  goto done;
+	}
+    }
+
+done:
+
+  if (convert_to_complex)
+    {
+      FloatComplexNDArray complex_result (a_dims);
+
+      for (octave_idx_type i = 0; i < len; i++)
+	{
+	  OCTAVE_QUIT;
+	  FloatComplex atmp (a(i));
+	  FloatComplex btmp (b(i));
+	  complex_result(i) = std::pow (atmp, btmp);
+	}
+
+      retval = complex_result;
+    }
+  else
+    {
+      FloatNDArray result (a_dims);
+
+      for (octave_idx_type i = 0; i < len; i++)
+	{
+	  OCTAVE_QUIT;
+	  result(i) = std::pow (a(i), b(i));
+	}
+
+      retval = result;
+    }
+
+  return retval;
+}
+
+// -*- 5 -*-
+octave_value
+elem_xpow (const FloatNDArray& a, const FloatComplex& b)
+{
+  FloatComplexNDArray result (a.dims ());
+
+  for (octave_idx_type i = 0; i < a.length (); i++)
+    {
+      OCTAVE_QUIT;
+      result(i) = std::pow (FloatComplex (a(i)), b);
+    }
+
+  return result;
+}
+
+// -*- 6 -*-
+octave_value
+elem_xpow (const FloatNDArray& a, const FloatComplexNDArray& b)
+{
+  dim_vector a_dims = a.dims ();
+  dim_vector b_dims = b.dims ();
+
+  if (a_dims != b_dims)
+    {
+      gripe_nonconformant ("operator .^", a_dims, b_dims);
+      return octave_value ();
+    }
+
+  FloatComplexNDArray result (a_dims);
+
+  for (octave_idx_type i = 0; i < a.length (); i++)
+    {
+      OCTAVE_QUIT;
+      result(i) = std::pow (FloatComplex (a(i)), b(i));
+    }
+
+  return result;
+}
+
+// -*- 7 -*-
+octave_value
+elem_xpow (const FloatComplex& a, const FloatNDArray& b)
+{
+  FloatComplexNDArray result (b.dims ());
+
+  for (octave_idx_type i = 0; i < b.length (); i++)
+    {
+      OCTAVE_QUIT;
+      float btmp = b(i);
+      if (xisint (btmp))
+	result(i) = std::pow (a, static_cast<int> (btmp));
+      else
+	result(i) = std::pow (a, btmp);
+    }
+
+  return result;
+}
+
+// -*- 8 -*-
+octave_value
+elem_xpow (const FloatComplex& a, const FloatComplexNDArray& b)
+{
+  FloatComplexNDArray result (b.dims ());
+
+  for (octave_idx_type i = 0; i < b.length (); i++)
+    {
+      OCTAVE_QUIT;
+      result(i) = std::pow (a, b(i));
+    }
+
+  return result;
+}
+
+// -*- 9 -*-
+octave_value
+elem_xpow (const FloatComplexNDArray& a, float b)
+{
+  FloatComplexNDArray result (a.dims ());
+
+  if (xisint (b))
+    {
+      for (octave_idx_type i = 0; i < a.length (); i++)
+	{
+	  OCTAVE_QUIT;
+	  result(i) = std::pow (a(i), static_cast<int> (b));
+	}
+    }
+  else
+    {
+      for (octave_idx_type i = 0; i < a.length (); i++)
+	{
+	  OCTAVE_QUIT;
+	  result(i) = std::pow (a(i), b);
+	}
+    }
+
+  return result;
+}
+
+// -*- 10 -*-
+octave_value
+elem_xpow (const FloatComplexNDArray& a, const FloatNDArray& b)
+{
+  dim_vector a_dims = a.dims ();
+  dim_vector b_dims = b.dims ();
+
+  if (a_dims != b_dims)
+    {
+      gripe_nonconformant ("operator .^", a_dims, b_dims);
+      return octave_value ();
+    }
+
+  FloatComplexNDArray result (a_dims);
+
+  for (octave_idx_type i = 0; i < a.length (); i++)
+    {
+      OCTAVE_QUIT;
+      float btmp = b(i);
+      if (xisint (btmp))
+	result(i) = std::pow (a(i), static_cast<int> (btmp));
+      else
+	result(i) = std::pow (a(i), btmp);
+    }
+
+  return result;
+}
+
+// -*- 11 -*-
+octave_value
+elem_xpow (const FloatComplexNDArray& a, const FloatComplex& b)
+{
+  FloatComplexNDArray result (a.dims ());
+
+  for (octave_idx_type i = 0; i < a.length (); i++)
+    {
+      OCTAVE_QUIT;
+      result(i) = std::pow (a(i), b);
+    }
+
+  return result;
+}
+
+// -*- 12 -*-
+octave_value
+elem_xpow (const FloatComplexNDArray& a, const FloatComplexNDArray& b)
+{
+  dim_vector a_dims = a.dims ();
+  dim_vector b_dims = b.dims ();
+
+  if (a_dims != b_dims)
+    {
+      gripe_nonconformant ("operator .^", a_dims, b_dims);
+      return octave_value ();
+    }
+
+  FloatComplexNDArray result (a_dims);
+
+  for (octave_idx_type i = 0; i < a.length (); i++)
+    {
+      OCTAVE_QUIT;
+      result(i) = std::pow (a(i), b(i));
+    }
+
+  return result;
+}
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***