--- ./glucat/clifford_algebra_imp.h.org	2010-05-29 21:44:36.000000000 +0900
+++ ./glucat/clifford_algebra_imp.h	2011-08-10 15:40:54.000000000 +0900
@@ -444,7 +444,7 @@
   >
   inline
   const Multivector<Scalar_T,LO,HI>
-  pow(const Multivector<Scalar_T,LO,HI>& lhs, const RHS<Scalar_T,LO,HI>& rhs)
+  original_pow(const Multivector<Scalar_T,LO,HI>& lhs, const RHS<Scalar_T,LO,HI>& rhs)
   {
     typedef numeric_traits<Scalar_T> traits_t;
 
@@ -463,6 +463,38 @@
     return exp(log(lhs) * rhs);
   }
 
+  /// Multivector power of multivector
+  template
+  <
+    template<typename, const index_t, const index_t> class Multivector,
+    template<typename, const index_t, const index_t> class RHS,
+    typename Scalar_T, const index_t LO, const index_t HI
+  >
+  inline
+  const Multivector<Scalar_T,LO,HI>
+  simple_pow(const Multivector<Scalar_T,LO,HI>& lhs, const RHS<Scalar_T,LO,HI>& rhs)
+  {
+    return exp(rhs * log(lhs));
+  }
+
+  /// Multivector power of multivector
+  template
+  <
+    template<typename, const index_t, const index_t> class Multivector,
+    template<typename, const index_t, const index_t> class RHS,
+    typename Scalar_T, const index_t LO, const index_t HI
+  >
+  inline
+  const Multivector<Scalar_T,LO,HI>
+  pow(const Multivector<Scalar_T,LO,HI>& lhs, const RHS<Scalar_T,LO,HI>& rhs)
+  {
+#if defined(_GLUCAT_USE_SIMPLE_POW)
+    return simple_pow(lhs, rhs);
+#else
+    return original_pow(lhs, rhs);
+#endif
+  }
+
   /// Square root of -1 which commutes with all members of the frame of the given multivector
   template< template<typename, const index_t, const index_t> class Multivector,
             typename Scalar_T, const index_t LO, const index_t HI >
@@ -580,9 +612,17 @@
   {
     typedef numeric_traits<Scalar_T> traits_t;
     check_complex(val, i, prechecked);
+#if defined(_GLUCAT_USE_COMPATIBLE_COMPLEX_ACOSH)
+    if (val.isnan())
+      return traits_t::NaN();
+    Multivector<Scalar_T,LO,HI> s = sqrt(val*val - Scalar_T(1), i, true), r = log(val + s, i, true);
+    if (scalar(r) < Scalar_T(0)) r = log(val - s, i, true);
+    return scalar(r) + pure(r)*Scalar_T((abs(pure(val)) != Scalar_T(0) || !(traits_t::abs(scalar(val)) > Scalar_T(1))) ? 1 : -1);
+#else
     return val.isnan()
          ? traits_t::NaN()
          : log(val + sqrt(val*val - Scalar_T(1), i, true), i, true);
+#endif
   }
 
   /// Inverse hyperbolic cosine of multivector
@@ -637,12 +677,18 @@
     if (val.isnan())
       return traits_t::NaN();
 
+#if defined(_GLUCAT_USE_COMPATIBLE_COMPLEX_ACOS)
+    check_complex(val, i, prechecked);
+    Multivector<Scalar_T,LO,HI> r = -i * log(val + i * sqrt(Scalar_T(1) - val*val, i, true), i, true);
+    return ((abs(pure(val)) != Scalar_T(0) || scalar(val) < Scalar_T(1)) ? scalar(r) : Scalar_T(0)) + pure(r)*Scalar_T((abs(pure(val)) != Scalar_T(0) || scalar(val) > Scalar_T(0)) ? 1 : -1);
+#else
     const Scalar_T& realval = real(val);
     if (val == realval && traits_t::abs(realval) <= Scalar_T(1))
       return traits_t::acos(realval);
 
     check_complex(val, i, prechecked);
     return i * acosh(val, i, true);
+#endif
   }
 
   /// Inverse cosine of multivector
@@ -682,9 +728,16 @@
   {
     typedef numeric_traits<Scalar_T> traits_t;
     check_complex(val, i, prechecked);
+#if defined(_GLUCAT_USE_COMPATIBLE_COMPLEX_ASINH)
+    if (val.isnan())
+      return traits_t::NaN();
+    Multivector<Scalar_T,LO,HI> r = log(val + sqrt(val*val + Scalar_T(1), i, true), i, true);
+    return scalar(r)*Scalar_T((scalar(val) != Scalar_T(0) || scalar(i*val) > Scalar_T(1)) ? 1 : -1) + pure(r);
+#else
     return val.isnan()
          ? traits_t::NaN()
          : log(val + sqrt(val*val + Scalar_T(1), i, true), i, true);
+#endif
   }
 
   /// Inverse hyperbolic sine of multivector
@@ -739,12 +792,18 @@
     if (val.isnan())
       return traits_t::NaN();
 
+#if defined(_GLUCAT_USE_COMPATIBLE_COMPLEX_ASIN)
+    check_complex(val, i, prechecked);
+    Multivector<Scalar_T,LO,HI> r = -i * asinh(val * i, i, true);
+    return scalar(r) + pure(r)*Scalar_T((abs(pure(val)) != Scalar_T(0) || scalar(val) < Scalar_T(0)) ? 1 : -1);
+#else
     const Scalar_T& realval = real(val);
     if (val == realval && traits_t::abs(realval) <= Scalar_T(1))
       return traits_t::asin(realval);
 
     check_complex(val, i, prechecked);
     return -i * asinh(i * val, i, true);
+#endif
   }
 
   /// Inverse sine of multivector
@@ -784,11 +843,20 @@
   {
     typedef numeric_traits<Scalar_T> traits_t;
     check_complex(val, i, prechecked);
+#if defined(_GLUCAT_USE_COMPATIBLE_COMPLEX_ATANH)
+    if (val.isnan())
+      return traits_t::NaN();
+    Multivector<Scalar_T,LO,HI> r = (norm(val + Scalar_T(1)) > norm(val - Scalar_T(1))) ? 
+      (log(val + Scalar_T(1), i, true) - log(-val + Scalar_T(1), i, true)) / Scalar_T(2) :
+      log(Scalar_T(1)/(-val + Scalar_T(1))*(val + Scalar_T(1)), i, true) / Scalar_T(2);
+    return scalar(r) + pure(r)*Scalar_T((scalar(val) < Scalar_T(1) || abs(pure(val)) != Scalar_T(0)) ? 1 : -1);
+#else
     return val.isnan()
         ? traits_t::NaN()
         : (norm(val + Scalar_T(1)) > norm(val - Scalar_T(1)))
           ? (log(val + Scalar_T(1), i, true) - log(-val + Scalar_T(1), i, true)) / Scalar_T(2)
           : log((val + Scalar_T(1)) / (-val + Scalar_T(1)), i, true) / Scalar_T(2);
+#endif
   }
 
   /// Inverse hyperbolic tangent of multivector
@@ -839,12 +907,20 @@
     if (val.isnan())
       return traits_t::NaN();
 
+#if defined(_GLUCAT_USE_COMPATIBLE_COMPLEX_ATAN)
+    check_complex(val, i, prechecked);
+    Multivector<Scalar_T,LO,HI> val_i = val * i, r = (norm(val_i + Scalar_T(1)) > norm(val_i - Scalar_T(1))) ? 
+      -i * (log(val_i + Scalar_T(1), i, true) - log(-val_i + Scalar_T(1), i, true)) / Scalar_T(2) :
+      -i * log(Scalar_T(1)/(-val_i + Scalar_T(1))*(val_i + Scalar_T(1)), i, true) / Scalar_T(2);
+    return scalar(r)*Scalar_T((scalar(val) != Scalar_T(0) || scalar(i*val) < Scalar_T(1)) ? 1 : -1) + pure(r);
+#else
     const Scalar_T& s = scalar(val);
     if (val == s)
       return traits_t::atan(s);
 
     check_complex(val, i, prechecked);
     return -i * atanh(i * val, i, true);
+#endif
   }
 
   /// Inverse tangent of multivector
--- ./glucat/matrix_multi_imp.h.org	2010-05-29 22:58:05.000000000 +0900
+++ ./glucat/matrix_multi_imp.h	2011-08-09 18:59:58.000000000 +0900
@@ -1179,7 +1179,7 @@
   template< typename Scalar_T, const index_t LO, const index_t HI >
   static
   const matrix_multi<Scalar_T,LO,HI>
-  sqrt(const matrix_multi<Scalar_T,LO,HI>& val, const matrix_multi<Scalar_T,LO,HI>& i, bool prechecked)
+  original_sqrt(const matrix_multi<Scalar_T,LO,HI>& val, const matrix_multi<Scalar_T,LO,HI>& i, bool prechecked)
   {
     // Reference: [GW], Section 4.3, pp318-322
     // Reference: [GL], Section 11.3, p572-576
@@ -1236,6 +1236,11 @@
       return sqrt(-i * val, i, prechecked) * (i + Scalar_T(1)) / sqrt_2;
 #endif
 
+#if defined(_GLUCAT_USE_COMPATIBLE_COMPLEX_SQRT)
+    const Scalar_T max_norm = Scalar_T(1.0/9.0);
+    const Scalar_T scale = abs(val);
+    const Scalar_T sqrt_scale = traits_t::sqrt(traits_t::abs(scale));
+#else
     // Scale val towards abs(A) == 1 or towards A == 1 as appropriate
     const Scalar_T max_norm = Scalar_T(1.0/9.0);
     const Scalar_T scale =
@@ -1247,6 +1252,7 @@
     const Scalar_T sqrt_scale = traits_t::sqrt(traits_t::abs(scale));
     if (traits_t::isNaN_or_isInf(sqrt_scale))
       return traits_t::NaN();
+#endif
 
     multivector_t rescale = sqrt_scale;
     if (scale < Scalar_T(0))
@@ -1314,6 +1320,28 @@
 #endif
   }
 
+  /// Square root of multivector with specified complexifier
+  template< typename Scalar_T, const index_t LO, const index_t HI >
+  static
+  const matrix_multi<Scalar_T,LO,HI>
+  simple_sqrt(const matrix_multi<Scalar_T,LO,HI>& val, const matrix_multi<Scalar_T,LO,HI>& i, bool prechecked)
+  {
+    return pow(val, matrix_multi<Scalar_T,LO,HI>(0.5));
+  }
+
+  /// Square root of multivector with specified complexifier
+  template< typename Scalar_T, const index_t LO, const index_t HI >
+  inline
+  const matrix_multi<Scalar_T,LO,HI>
+  sqrt(const matrix_multi<Scalar_T,LO,HI>& val, const matrix_multi<Scalar_T,LO,HI>& i, bool prechecked)
+  {
+#if defined(_GLUCAT_USE_SIMPLE_SQRT)
+    return simple_sqrt(val, i, prechecked);
+#else
+    return original_sqrt(val, i, prechecked);
+#endif
+  }
+
   /// Exponential of multivector
   template< typename Scalar_T, const index_t LO, const index_t HI >
   const matrix_multi<Scalar_T,LO,HI>
@@ -1417,7 +1445,12 @@
     };
     typedef numeric_traits<Scalar_T> traits_t;
 
-    if (val == Scalar_T(0) || val.isnan())
+    if (
+#if defined(_GLUCAT_USE_COMPATIBLE_COMPLEX_LOG)
+#else
+	val == Scalar_T(0) ||
+#endif
+	val.isnan())
       return traits_t::NaN();
     else
       return pade_approx(a, b, -val + Scalar_T(1));
@@ -1433,7 +1466,12 @@
     typedef matrix_multi<Scalar_T,LO,HI> multivector_t;
     typedef numeric_traits<Scalar_T> traits_t;
 
-    if (val == Scalar_T(0) || val.isnan())
+    if (
+#if defined(_GLUCAT_USE_COMPATIBLE_COMPLEX_LOG)
+#else
+	val == Scalar_T(0) ||
+#endif
+	val.isnan())
       return traits_t::NaN();
 
     multivector_t Y = val;
@@ -1449,7 +1487,12 @@
         outer_step != Tune_P::log_max_outer_steps && norm_Y_1 * pow_2_outer_step > max_outer_norm;
         ++outer_step,    norm_Y_1 = norm(Y - Scalar_T(1)))
     {
-      if (Y == Scalar_T(0) || Y.isnan())
+      if (
+#if defined(_GLUCAT_USE_COMPATIBLE_COMPLEX_LOG)
+#else
+	  Y == Scalar_T(0) ||
+#endif
+	  Y.isnan())
         return traits_t::NaN();
 
       // Incomplete product form of Denman-Beavers square root iteration
@@ -1481,7 +1524,12 @@
     //std::cout << "Log called" << std::endl;
     typedef numeric_traits<Scalar_T> traits_t;
 
-    if (val == Scalar_T(0) || val.isnan())
+    if (
+#if defined(_GLUCAT_USE_COMPATIBLE_COMPLEX_LOG)
+#else
+	val == Scalar_T(0) ||
+#endif
+	val.isnan())
       return traits_t::NaN();
 
     static const Scalar_T pi = traits_t::pi();
@@ -1506,6 +1554,9 @@
       return log(-i * val, i, prechecked) + i * pi/Scalar_T(2);
     }
 #endif
+#if defined(_GLUCAT_USE_COMPATIBLE_COMPLEX_LOG)
+    const Scalar_T scale = abs(val);
+#else
     // Scale val towards abs(A) == 1 or towards A == 1 as appropriate
     const Scalar_T max_norm = Scalar_T(1.0/9.0);
     const Scalar_T scale =
@@ -1516,6 +1567,7 @@
           : abs(val);
     if (scale == Scalar_T(0))
       return traits_t::NaN();
+#endif
 
     const Scalar_T log_scale = traits_t::log(traits_t::abs(scale));
     multivector_t rescale = log_scale;
