--- ./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 - pow(const Multivector& lhs, const RHS& rhs) + original_pow(const Multivector& lhs, const RHS& rhs) { typedef numeric_traits traits_t; @@ -463,6 +463,38 @@ return exp(log(lhs) * rhs); } + /// Multivector power of multivector + template + < + template class Multivector, + template class RHS, + typename Scalar_T, const index_t LO, const index_t HI + > + inline + const Multivector + simple_pow(const Multivector& lhs, const RHS& rhs) + { + return exp(rhs * log(lhs)); + } + + /// Multivector power of multivector + template + < + template class Multivector, + template class RHS, + typename Scalar_T, const index_t LO, const index_t HI + > + inline + const Multivector + pow(const Multivector& lhs, const RHS& 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 class Multivector, typename Scalar_T, const index_t LO, const index_t HI > @@ -580,9 +612,17 @@ { typedef numeric_traits traits_t; check_complex(val, i, prechecked); +#if defined(_GLUCAT_USE_COMPATIBLE_COMPLEX_ACOSH) + if (val.isnan()) + return traits_t::NaN(); + Multivector 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 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 traits_t; check_complex(val, i, prechecked); +#if defined(_GLUCAT_USE_COMPATIBLE_COMPLEX_ASINH) + if (val.isnan()) + return traits_t::NaN(); + Multivector 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 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 traits_t; check_complex(val, i, prechecked); +#if defined(_GLUCAT_USE_COMPATIBLE_COMPLEX_ATANH) + if (val.isnan()) + return traits_t::NaN(); + Multivector 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 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 - sqrt(const matrix_multi& val, const matrix_multi& i, bool prechecked) + original_sqrt(const matrix_multi& val, const matrix_multi& 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 + simple_sqrt(const matrix_multi& val, const matrix_multi& i, bool prechecked) + { + return pow(val, matrix_multi(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 + sqrt(const matrix_multi& val, const matrix_multi& 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 @@ -1417,7 +1445,12 @@ }; typedef numeric_traits 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 multivector_t; typedef numeric_traits 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 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;