1 #ifndef _GLUCAT_MATRIX_MULTI_IMP_H
2 #define _GLUCAT_MATRIX_MULTI_IMP_H
39 # if defined(_GLUCAT_GCC_IGNORE_UNUSED_LOCAL_TYPEDEFS)
40 # pragma GCC diagnostic push
41 # pragma GCC diagnostic ignored "-Wunused-local-typedefs"
43 # if defined(_GLUCAT_HAVE_BOOST_SERIALIZATION_ARRAY_WRAPPER_H)
44 # include <boost/serialization/array_wrapper.hpp>
46 #include <boost/numeric/ublas/matrix.hpp>
47 #include <boost/numeric/ublas/matrix_expression.hpp>
48 #include <boost/numeric/ublas/matrix_proxy.hpp>
49 #include <boost/numeric/ublas/matrix_sparse.hpp>
50 #include <boost/numeric/ublas/operation.hpp>
51 #include <boost/numeric/ublas/operation_sparse.hpp>
52 #include <boost/numeric/ublas/triangular.hpp>
53 #include <boost/numeric/ublas/lu.hpp>
54 # if defined(_GLUCAT_GCC_IGNORE_UNUSED_LOCAL_TYPEDEFS)
55 # pragma GCC diagnostic pop
69 template<
typename Scalar_T, const index_t LO, const index_t HI >
73 {
return "matrix_multi"; }
83 static const int offset_log2_dim[] = {0, 1, 0, 1, 1, 2, 1, 1};
85 return (p+q)/2 + offset_log2_dim[bott];
90 template<
typename Matrix_Index_T, const index_t LO, const index_t HI >
95 {
return 1 <<
offset_level(sub.count_pos(), sub.count_neg()); }
98 template<
typename Scalar_T, const index_t LO, const index_t HI >
106 template<
typename Scalar_T, const index_t LO, const index_t HI >
107 template<
typename Other_Scalar_T >
110 : m_frame( val.m_frame ), m_matrix( val.m_matrix.size1(), val.m_matrix.size2() )
112 this->m_matrix.clear();
114 typedef typename other_matrix_t::const_iterator1 other_const_iterator1;
115 typedef typename other_matrix_t::const_iterator2 other_const_iterator2;
116 for (other_const_iterator1
120 for (other_const_iterator2
121 val_it2 = val_it1.begin();
122 val_it2 != val_it1.end();
128 template<
typename Scalar_T, const index_t LO, const index_t HI >
129 template<
typename Other_Scalar_T >
139 this->
m_matrix.resize(dim, dim,
false);
142 typedef typename other_matrix_t::const_iterator1 other_const_iterator1;
143 typedef typename other_matrix_t::const_iterator2 other_const_iterator2;
144 for (other_const_iterator1
148 for (other_const_iterator2
149 val_it2 = val_it1.begin();
150 val_it2 != val_it1.end();
157 template<
typename Scalar_T, const index_t LO, const index_t HI >
169 template<
typename Scalar_T, const index_t LO, const index_t HI >
175 this->
m_matrix.resize(dim, dim,
false);
177 *
this +=
term_t(ist, crd);
181 template<
typename Scalar_T, const index_t LO, const index_t HI >
183 matrix_multi(
const index_set_t
ist,
const Scalar_T& crd,
const index_set_t frm,
const bool prechecked)
186 if (!prechecked && (
ist | frm) != frm)
187 throw error_t(
"multivector_t(ist,crd,frm): cannot initialize with value outside of frame");
189 this->
m_matrix.resize(dim, dim,
false);
191 *
this +=
term_t(ist, crd);
195 template<
typename Scalar_T, const index_t LO, const index_t HI >
197 matrix_multi(
const Scalar_T& scr,
const index_set_t frm)
201 this->m_matrix.resize(dim, dim,
false);
202 this->m_matrix.clear();
207 template<
typename Scalar_T, const index_t LO, const index_t HI >
210 { *
this = multivector_t(Scalar_T(scr), frm); }
213 template<
typename Scalar_T, const index_t LO, const index_t HI >
220 throw error_t(
"multivector_t(vec,frm): cannot initialize with vector not matching frame");
222 this->
m_matrix.resize(dim, dim,
false);
224 typename vector_t::const_iterator vec_it = vec.begin();
240 template<
typename Scalar_T, const index_t LO, const index_t HI >
246 template<
typename Scalar_T, const index_t LO, const index_t HI >
248 matrix_multi(
const std::string& str,
const index_set_t frm,
const bool prechecked)
249 { *
this = multivector_t(framed_multi_t(str), frm, prechecked); }
252 template<
typename Scalar_T, const index_t LO, const index_t HI >
253 template<
typename Other_Scalar_T >
255 matrix_multi(
const framed_multi<Other_Scalar_T,LO,HI>& val)
256 : m_frame( val.frame() )
261 *
this = val.template fast_matrix_multi<Scalar_T>(this->
m_frame);
264 catch (
const glucat_error&
e)
267 this->
m_matrix.resize(dim, dim,
false);
272 val_it = val.begin();
279 template<
typename Scalar_T, const index_t LO, const index_t HI >
280 template<
typename Other_Scalar_T >
288 *
this = val.template fast_matrix_multi<Scalar_T>(our_frame);
293 this->m_frame = our_frame;
295 this->m_matrix.resize(dim, dim,
false);
296 this->m_matrix.clear();
300 val_it = val.begin();
307 template<
typename Scalar_T, const index_t LO, const index_t HI >
308 template<
typename Matrix_T >
311 : m_frame( frm ), m_matrix( mtx.size1(), mtx.size2() )
315 typedef typename Matrix_T::const_iterator1 const_iterator1;
316 typedef typename Matrix_T::const_iterator2 const_iterator2;
318 mtx_it1 = mtx.begin1();
319 mtx_it1 != mtx.end1();
322 mtx_it2 = mtx_it1.begin();
323 mtx_it2 != mtx_it1.end();
329 template<
typename Scalar_T, const index_t LO, const index_t HI >
331 matrix_multi(
const matrix_t& mtx,
const index_set_t frm)
332 : m_frame( frm ), m_matrix( mtx )
336 template<
typename Scalar_T, const index_t LO, const index_t HI >
337 matrix_multi<Scalar_T,LO,HI>&
350 template<
typename Scalar_T, const index_t LO, const index_t HI >
358 typedef typename multivector_t::framed_multi_t framed_multi_t;
368 our_frame |= framed_lhs.frame() | framed_rhs.frame();
372 lhs_reframed = multivector_t(framed_lhs, our_frame,
true);
374 rhs_reframed = multivector_t(framed_rhs, our_frame,
true);
379 template<
typename Scalar_T, const index_t LO, const index_t HI >
381 matrix_multi<Scalar_T,LO,HI>::
382 operator== (
const multivector_t& rhs)
const
389 multivector_t lhs_reframed;
390 multivector_t rhs_reframed;
391 const index_set_t our_frame =
reframe(*
this, rhs, lhs_reframed, rhs_reframed);
392 const multivector_t& lhs_ref = (this->m_frame == our_frame)
395 const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
399 #if defined(_GLUCAT_USE_DENSE_MATRICES)
400 return ublas::norm_inf(lhs_ref.m_matrix - rhs_ref.m_matrix) == 0;
402 typedef typename matrix_t::const_iterator1 const_iterator1;
403 typedef typename matrix_t::const_iterator2 const_iterator2;
407 it1 = lhs_ref.m_matrix.begin1();
408 it1 != lhs_ref.m_matrix.end1();
415 return ublas::norm_inf(lhs_ref.m_matrix - rhs_ref.m_matrix) == 0;
417 it1 = rhs_ref.m_matrix.begin1();
418 it1 != rhs_ref.m_matrix.end1();
425 return ublas::norm_inf(lhs_ref.m_matrix - rhs_ref.m_matrix) == 0;
428 const_iterator1 this_it1 = lhs_ref.m_matrix.begin1();
429 const_iterator1 it1 = rhs_ref.m_matrix.begin1();
431 (this_it1 != lhs_ref.m_matrix.end1()) && (it1 != rhs_ref.m_matrix.end1());
434 if ( this_it1.index1() != it1.index1() )
436 const_iterator2 this_it2 = this_it1.begin();
437 const_iterator2 it2 = it1.begin();
439 (this_it2 != this_it1.end()) && (it2 != it1.end());
441 if ( (this_it2.index2() != it2.index2()) || (*this_it2 != *it2) )
443 if ( (this_it2 != this_it1.end()) || (it2 != it1.end()) )
446 return (this_it1 == lhs_ref.m_matrix.end1()) && (it1 == rhs_ref.m_matrix.end1());
451 template<
typename Scalar_T, const index_t LO, const index_t HI >
454 matrix_multi<Scalar_T,LO,HI>::
455 operator== (
const Scalar_T& scr)
const
457 if (scr != Scalar_T(0))
458 return *
this == multivector_t(framed_multi_t(scr), this->m_frame,
true);
459 else if (ublas::norm_inf(this->m_matrix) != 0)
463 const matrix_index_t dim = this->m_matrix.size1();
464 return !(dim == 1 && this->
isnan());
469 template<
typename Scalar_T, const index_t LO, const index_t HI >
471 matrix_multi<Scalar_T,LO,HI>&
474 {
return *
this += term_t(index_set_t(), scr); }
477 template<
typename Scalar_T, const index_t LO, const index_t HI >
479 matrix_multi<Scalar_T,LO,HI>&
485 return *
this *= Scalar_T(2);
488 multivector_t rhs_reframed;
489 const index_set_t our_frame =
reframe(*
this, rhs, *
this, rhs_reframed);
490 const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
494 noalias(this->m_matrix) += rhs_ref.m_matrix;
499 template<
typename Scalar_T, const index_t LO, const index_t HI >
507 return *
this = Scalar_T(0);
516 noalias(this->m_matrix) -= rhs_ref.
m_matrix;
521 template<
typename Scalar_T, const index_t LO, const index_t HI >
526 {
return multivector_t(-(this->m_matrix), this->m_frame); }
529 template<
typename Scalar_T, const index_t LO, const index_t HI >
531 matrix_multi<Scalar_T,LO,HI>&
532 matrix_multi<Scalar_T,LO,HI>::
533 operator*= (
const Scalar_T& scr)
536 typedef numeric_traits<Scalar_T> traits_t;
537 if (traits_t::isNaN_or_isInf(scr) || this->
isnan())
538 return *
this = traits_t::NaN();
539 if (scr == Scalar_T(0))
542 this->m_matrix *= scr;
547 template<
typename Scalar_T, const index_t LO, const index_t HI >
549 const matrix_multi<Scalar_T,LO,HI>
550 operator* (
const matrix_multi<Scalar_T,LO,HI>& lhs,
const matrix_multi<Scalar_T,LO,HI>& rhs)
552 typedef matrix_multi<Scalar_T,LO,HI> multivector_t;
553 typedef typename multivector_t::index_set_t index_set_t;
555 #if defined(_GLUCAT_CHECK_ISNAN)
556 if (lhs.isnan() || rhs.isnan())
561 multivector_t lhs_reframed;
562 multivector_t rhs_reframed;
563 const index_set_t our_frame =
reframe(lhs, rhs, lhs_reframed, rhs_reframed);
564 const multivector_t& lhs_ref = (lhs.m_frame == our_frame)
567 const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
571 typedef typename multivector_t::matrix_t matrix_t;
572 #if defined(_GLUCAT_USE_DENSE_MATRICES)
573 typedef typename matrix_t::size_type matrix_index_t;
575 const matrix_index_t dim = lhs_ref.m_matrix.size1();
576 multivector_t result = multivector_t(matrix_t(dim, dim), our_frame);
577 result.m_matrix.clear();
578 ublas::axpy_prod(lhs_ref.m_matrix, rhs_ref.m_matrix, result.m_matrix,
true);
581 typedef typename matrix_t::expression_type expression_t;
584 multivector_t(ublas::sparse_prod<expression_t>(lhs_ref.m_matrix, rhs_ref.m_matrix), our_frame);
589 template<
typename Scalar_T, const index_t LO, const index_t HI >
591 matrix_multi<Scalar_T,LO,HI>&
592 matrix_multi<Scalar_T,LO,HI>::
593 operator*= (
const multivector_t& rhs)
594 {
return *
this = *
this * rhs; }
597 template<
typename Scalar_T, const index_t LO, const index_t HI >
599 const matrix_multi<Scalar_T,LO,HI>
600 operator^ (
const matrix_multi<Scalar_T,LO,HI>& lhs,
const matrix_multi<Scalar_T,LO,HI>& rhs)
602 typedef matrix_multi<Scalar_T,LO,HI> multivector_t;
603 typedef typename multivector_t::framed_multi_t framed_multi_t;
604 return framed_multi_t(lhs) ^ framed_multi_t(rhs);
608 template<
typename Scalar_T, const index_t LO, const index_t HI >
610 matrix_multi<Scalar_T,LO,HI>&
611 matrix_multi<Scalar_T,LO,HI>::
612 operator^= (
const multivector_t& rhs)
613 {
return *
this = *
this ^ rhs; }
616 template<
typename Scalar_T, const index_t LO, const index_t HI >
618 const matrix_multi<Scalar_T,LO,HI>
619 operator& (
const matrix_multi<Scalar_T,LO,HI>& lhs,
const matrix_multi<Scalar_T,LO,HI>& rhs)
621 typedef matrix_multi<Scalar_T,LO,HI> multivector_t;
622 typedef typename multivector_t::framed_multi_t framed_multi_t;
623 return framed_multi_t(lhs) & framed_multi_t(rhs);
627 template<
typename Scalar_T, const index_t LO, const index_t HI >
632 {
return *
this = *
this & rhs; }
635 template<
typename Scalar_T, const index_t LO, const index_t HI >
637 const matrix_multi<Scalar_T,LO,HI>
638 operator% (
const matrix_multi<Scalar_T,LO,HI>& lhs,
const matrix_multi<Scalar_T,LO,HI>& rhs)
640 typedef matrix_multi<Scalar_T,LO,HI> multivector_t;
641 typedef typename multivector_t::framed_multi_t framed_multi_t;
642 return framed_multi_t(lhs) % framed_multi_t(rhs);
646 template<
typename Scalar_T, const index_t LO, const index_t HI >
651 {
return *
this = *
this % rhs; }
654 template<
typename Scalar_T, const index_t LO, const index_t HI >
657 star(
const matrix_multi<Scalar_T,LO,HI>& lhs,
const matrix_multi<Scalar_T,LO,HI>& rhs)
658 {
return (lhs * rhs).scalar(); }
661 template<
typename Scalar_T, const index_t LO, const index_t HI >
663 matrix_multi<Scalar_T,LO,HI>&
664 matrix_multi<Scalar_T,LO,HI>::
665 operator/= (
const Scalar_T& scr)
666 {
return *
this *= Scalar_T(1)/scr; }
669 template<
typename Scalar_T, const index_t LO, const index_t HI >
675 #if defined(_GLUCAT_CHECK_ISNAN)
676 if (lhs.isnan() || rhs.isnan())
677 return traits_t::NaN();
680 if (rhs == Scalar_T(0))
681 return traits_t::NaN();
684 typedef typename multivector_t::index_set_t index_set_t;
687 multivector_t lhs_reframed;
688 multivector_t rhs_reframed;
689 const index_set_t our_frame =
reframe(lhs, rhs, lhs_reframed, rhs_reframed);
690 const multivector_t& lhs_ref = (lhs.
m_frame == our_frame)
693 const multivector_t& rhs_ref = (rhs.
m_frame == our_frame)
703 typedef typename multivector_t::matrix_t matrix_t;
704 typedef typename matrix_t::size_type matrix_index_t;
706 const matrix_t& AT = ublas::trans(rhs_ref.m_matrix);
709 typedef ublas::permutation_matrix<matrix_index_t> permutation_t;
711 permutation_t pvector(AT.size1());
712 if (! ublas::lu_factorize(LU, pvector))
714 const matrix_t& BT = ublas::trans(lhs_ref.m_matrix);
716 ublas::lu_substitute(LU, pvector, XT);
717 #if defined(_GLUCAT_CHECK_ISNAN)
719 return traits_t::NaN();
729 ublas::axpy_prod(AT, XT, R,
false);
730 #if defined(_GLUCAT_CHECK_ISNAN)
732 return traits_t::NaN();
735 Scalar_T nr = ublas::norm_inf(R);
736 if ( nr != Scalar_T(0) && !traits_t::isNaN_or_isInf(nr) )
739 Scalar_T nrold = nr + Scalar_T(1);
752 ublas::lu_substitute(LU, pvector, D);
756 ublas::axpy_prod(AT, XTnew, R,
false);
757 nr = ublas::norm_inf(R);
761 return multivector_t(ublas::trans(XT), our_frame);
765 return traits_t::NaN();
769 template<
typename Scalar_T, const index_t LO, const index_t HI >
771 matrix_multi<Scalar_T,LO,HI>&
772 matrix_multi<Scalar_T,LO,HI>::
773 operator/= (
const multivector_t& rhs)
774 {
return *
this = *
this / rhs; }
777 template<
typename Scalar_T, const index_t LO, const index_t HI >
779 const matrix_multi<Scalar_T,LO,HI>
780 operator| (
const matrix_multi<Scalar_T,LO,HI>& lhs,
const matrix_multi<Scalar_T,LO,HI>& rhs)
781 {
return rhs * lhs / rhs.involute(); }
784 template<
typename Scalar_T, const index_t LO, const index_t HI >
786 matrix_multi<Scalar_T,LO,HI>&
787 matrix_multi<Scalar_T,LO,HI>::
788 operator|= (
const multivector_t& rhs)
789 {
return *
this = rhs * *
this / rhs.involute(); }
792 template<
typename Scalar_T, const index_t LO, const index_t HI >
794 const matrix_multi<Scalar_T,LO,HI>
797 {
return multivector_t(Scalar_T(1), this->m_frame) / *
this; }
800 template<
typename Scalar_T, const index_t LO, const index_t HI >
802 const matrix_multi<Scalar_T,LO,HI>
808 template<
typename Scalar_T, const index_t LO, const index_t HI >
814 throw error_t(
"outer_pow(m): negative exponent");
815 framed_multi_t result = Scalar_T(1);
816 framed_multi_t a = *
this;
829 template<
typename Scalar_T, const index_t LO, const index_t HI >
832 matrix_multi<Scalar_T,LO,HI>::
834 {
return framed_multi_t(*this).grade(); }
837 template<
typename Scalar_T, const index_t LO, const index_t HI >
839 const index_set<LO,HI>
840 matrix_multi<Scalar_T,LO,HI>::
842 {
return this->m_frame; }
845 template<
typename Scalar_T, const index_t LO, const index_t HI >
848 matrix_multi<Scalar_T,LO,HI>::
849 operator[] (
const index_set_t
ist)
const
852 if ( (
ist | this->m_frame) == this->m_frame)
853 return matrix::inner<Scalar_T>(this->basis_element(
ist), this->m_matrix);
859 template<
typename Scalar_T, const index_t LO, const index_t HI >
861 const matrix_multi<Scalar_T,LO,HI>
862 matrix_multi<Scalar_T,LO,HI>::
863 operator() (
index_t grade)
const
865 if ((grade < 0) || (grade > HI-LO))
868 return (framed_multi_t(*
this))(grade);
872 template<
typename Scalar_T, const index_t LO, const index_t HI >
878 const matrix_index_t dim = this->m_matrix.size1();
879 return matrix::trace(this->m_matrix) / Scalar_T(
double(dim) );
883 template<
typename Scalar_T, const index_t LO, const index_t HI >
885 const matrix_multi<Scalar_T,LO,HI>
888 {
return *
this - this->
scalar(); }
891 template<
typename Scalar_T, const index_t LO, const index_t HI >
893 const matrix_multi<Scalar_T,LO,HI>
896 {
return framed_multi_t(*this).even(); }
899 template<
typename Scalar_T, const index_t LO, const index_t HI >
901 const matrix_multi<Scalar_T,LO,HI>
904 {
return framed_multi_t(*this).odd(); }
907 template<
typename Scalar_T, const index_t LO, const index_t HI >
914 template<
typename Scalar_T, const index_t LO, const index_t HI >
917 vector_part(
const index_set_t frm,
const bool prechecked)
const
919 if (!prechecked && (this->frame() | frm) != frm)
920 throw error_t(
"vector_part(frm): value is outside of requested frame");
923 if (this->frame() != frm)
924 return framed_multi_t(*this).vector_part(frm,
true);
926 const index_t begin_index = frm.min();
927 const index_t end_index = frm.max()+1;
936 matrix::inner<Scalar_T>(this->basis_element(index_set_t(idx)),
942 template<
typename Scalar_T, const index_t LO, const index_t HI >
944 const matrix_multi<Scalar_T,LO,HI>
947 {
return framed_multi_t(*this).involute(); }
950 template<
typename Scalar_T, const index_t LO, const index_t HI >
952 const matrix_multi<Scalar_T,LO,HI>
955 {
return framed_multi_t(*this).reverse(); }
958 template<
typename Scalar_T, const index_t LO, const index_t HI >
960 const matrix_multi<Scalar_T,LO,HI>
963 {
return framed_multi_t(*this).conj(); }
966 template<
typename Scalar_T, const index_t LO, const index_t HI >
973 return framed_multi_t(*this).quad();
977 template<
typename Scalar_T, const index_t LO, const index_t HI >
983 const matrix_index_t dim = this->m_matrix.size1();
988 template<
typename Scalar_T, const index_t LO, const index_t HI >
993 {
return framed_multi_t(*this).max_abs(); }
996 template<
typename Scalar_T, const index_t LO, const index_t HI >
997 const matrix_multi<Scalar_T,LO,HI>
999 random(
const index_set<LO,HI> frm, Scalar_T
fill)
1005 template<
typename Scalar_T, const index_t LO, const index_t HI >
1008 matrix_multi<Scalar_T,LO,HI>::
1009 write(
const std::string& msg)
const
1010 { framed_multi_t(*this).write(msg); }
1013 template<
typename Scalar_T, const index_t LO, const index_t HI >
1016 matrix_multi<Scalar_T,LO,HI>::
1017 write(std::ofstream& ofile,
const std::string& msg)
const
1020 throw error_t(
"write(ofile,msg): cannot write to output file");
1021 framed_multi_t(*this).write(ofile, msg);
1025 template<
typename Scalar_T, const index_t LO, const index_t HI >
1030 os << typename matrix_multi<Scalar_T,LO,HI>::framed_multi_t(val);
1035 template<
typename Scalar_T, const index_t LO, const index_t HI >
1050 template<
typename Scalar_T, const index_t LO, const index_t HI >
1056 if (std::numeric_limits<Scalar_T>::has_quiet_NaN)
1063 template<
typename Scalar_T, const index_t LO, const index_t HI >
1065 const matrix_multi<Scalar_T,LO,HI>
1066 matrix_multi<Scalar_T,LO,HI>::
1067 truncated(
const Scalar_T& limit)
const
1068 {
return framed_multi_t(*this).truncated(limit); }
1071 template<
typename Scalar_T, const index_t LO, const index_t HI >
1073 matrix_multi<Scalar_T,LO,HI>&
1077 if (term.second != Scalar_T(0))
1078 this->m_matrix.plus_assign(matrix_t(this->basis_element(term.first)) * term.second);
1083 template<
typename Multivector_T,
typename Matrix_T,
typename Basis_Matrix_T >
1088 typedef Multivector_T framed_multi_t;
1090 typedef typename framed_multi_t::index_set_t index_set_t;
1092 typedef Matrix_T matrix_t;
1093 typedef Basis_Matrix_T basis_matrix_t;
1094 typedef typename basis_matrix_t::value_type basis_scalar_t;
1095 typedef numeric_traits<Scalar_T> traits_t;
1098 return framed_multi_t(traits_t::to_scalar_t(X(0,0)));
1100 if (ublas::norm_inf(X) == 0)
1103 const basis_matrix_t& I = matrix::unit<basis_matrix_t>(2);
1104 basis_matrix_t J(2,2,2);
1106 J(0,1) = basis_scalar_t(-1);
1107 J(1,0) = basis_scalar_t( 1);
1108 basis_matrix_t K = J;
1109 K(0,1) = basis_scalar_t( 1);
1110 basis_matrix_t JK = I;
1111 JK(0,0) = basis_scalar_t(-1);
1114 const index_set_t ist_mn = index_set_t(-level);
1115 const index_set_t ist_pn = index_set_t(level);
1116 const index_set_t ist_mnpn = ist_mn | ist_pn;
1119 typedef typename framed_multi_t::term_t term_t;
1126 result += term_t(ist_mn, j_x);
1127 result += term_t(ist_pn, k_x);
1128 return result += term_t(ist_mnpn, jk_x);
1132 const framed_multi_t& mn = framed_multi_t(ist_mn);
1133 const framed_multi_t& pn = framed_multi_t(ist_pn);
1134 const framed_multi_t& mnpn = framed_multi_t(ist_mnpn);
1135 const framed_multi_t& i_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1137 const framed_multi_t& j_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1139 const framed_multi_t& k_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1141 const framed_multi_t& jk_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1144 result = i_x.even() - jk_x.odd();
1145 result += (j_x.even() - k_x.odd()) * mn;
1146 result += (k_x.even() - j_x.odd()) * pn;
1147 return result += (jk_x.even() - i_x.odd()) * mnpn;
1152 template<
typename Scalar_T, const index_t LO, const index_t HI >
1154 const matrix_multi<Scalar_T,LO,HI>
1158 if (this->m_frame == frm)
1161 return (this->
template fast_framed_multi<Scalar_T>()).template fast_matrix_multi<Scalar_T>(frm);
1165 template<
typename Scalar_T, const index_t LO, const index_t HI >
1166 template <
typename Other_Scalar_T>
1167 const framed_multi<Other_Scalar_T,LO,HI>
1172 index_t p = this->m_frame.count_pos();
1173 index_t q = this->m_frame.count_neg();
1191 const index_t level = (p+q)/2;
1194 typedef framed_multi<Other_Scalar_T,LO,HI> framed_multi_t;
1195 framed_multi_t val = fast<framed_multi_t, matrix_t, basis_matrix_t>(this->m_matrix, level);
1198 switch (
pos_mod(orig_p-orig_q, 8))
1203 val.centre_qp1_pm1(p, q);
1208 if (orig_p-orig_q > 4)
1210 val.centre_pp4_qm4(p, q);
1211 if (orig_p-orig_q < -3)
1213 val.centre_pm4_qp4(p, q);
1216 return val.unfold(this->m_frame);
1220 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Matrix_T >
1222 public std::map< std::pair< const index_set<LO,HI>, const index_set<LO,HI> >,
1243 template<
typename Scalar_T, const index_t LO, const index_t HI >
1248 typedef std::pair<const index_set_t, const index_set_t> index_set_pair_t;
1249 const index_set_pair_t& unfolded_pair = index_set_pair_t(
ist, this->m_frame);
1252 typedef typename basis_table_t::const_iterator basis_const_iterator_t;
1253 basis_table_t& basis_cache = basis_table_t::basis();
1255 const index_t frame_count = this->m_frame.count();
1260 const basis_const_iterator_t basis_it = basis_cache.find(unfolded_pair);
1261 if (basis_it != basis_cache.end())
1262 return *(basis_it->second);
1264 const index_set_t folded_set =
ist.fold(this->m_frame);
1265 const index_set_t folded_frame = this->m_frame.fold();
1266 const index_set_pair_t& folded_pair = index_set_pair_t(folded_set, folded_frame);
1267 typedef std::pair<const index_set_pair_t, basis_matrix_t*> basis_pair_t;
1270 const basis_const_iterator_t basis_it = basis_cache.find(folded_pair);
1271 if (basis_it != basis_cache.end())
1273 basis_matrix_t* result_ptr = basis_it->second;
1274 basis_cache.insert(basis_pair_t(unfolded_pair, result_ptr));
1278 const index_t folded_max = folded_frame.max();
1279 const index_t folded_min = folded_frame.min();
1284 basis_matrix_t result = matrix::unit<basis_matrix_t>(dim);
1293 basis_matrix_t* result_ptr =
new basis_matrix_t(result);
1294 basis_cache.insert(basis_pair_t(folded_pair, result_ptr));
1295 basis_cache.insert(basis_pair_t(unfolded_pair, result_ptr));
1301 template<
typename Scalar_T, const index_t LO, const index_t HI >
1304 const matrix_multi<Scalar_T,LO,HI>
1305 pade_approx(
const int array_size,
const Scalar_T a[],
const Scalar_T b[],
const matrix_multi<Scalar_T,LO,HI>& X)
1311 typedef matrix_multi<Scalar_T,LO,HI> multivector_t;
1312 typedef numeric_traits<Scalar_T> traits_t;
1315 return traits_t::NaN();
1318 const int nbr_even_powers = array_size/2 - 1;
1321 std::vector<multivector_t> XX(nbr_even_powers);
1323 XX[1] = XX[0] * XX[0];
1326 k != nbr_even_powers;
1328 XX[k] = XX[k-2] * XX[1];
1331 multivector_t N = a[1];
1334 k != nbr_even_powers;
1336 N += XX[k] * a[2*k + 3];
1341 k != nbr_even_powers;
1343 N += XX[k] * a[2*k + 2];
1344 multivector_t D = b[1];
1347 k != nbr_even_powers;
1349 D += XX[k] * b[2*k + 3];
1354 k != nbr_even_powers;
1356 D += XX[k] * b[2*k + 2];
1361 template<
typename Scalar_T, const index_t LO, const index_t HI >
1365 db_step(matrix_multi<Scalar_T,LO,HI>& M, matrix_multi<Scalar_T,LO,HI>& Y)
1368 typedef matrix_multi<Scalar_T,LO,HI> multivector_t;
1369 const multivector_t& invM =
inv(M);
1370 M = ((M + invM)/Scalar_T(2) + Scalar_T(1)) / Scalar_T(2);
1371 Y *= (invM + Scalar_T(1)) / Scalar_T(2);
1375 template<
typename Scalar_T, const index_t LO, const index_t HI >
1377 const matrix_multi<Scalar_T,LO,HI>
1378 db_sqrt(
const matrix_multi<Scalar_T,LO,HI>& val)
1381 typedef matrix_multi<Scalar_T,LO,HI> multivector_t;
1383 if (val == Scalar_T(0))
1386 typedef std::numeric_limits<Scalar_T> limits_t;
1388 static const Scalar_T tol2 = tol * tol;
1390 multivector_t M = val;
1391 multivector_t Y = val;
1392 Scalar_T norm_M_1 =
norm(M - Scalar_T(1));
1393 typedef numeric_traits<Scalar_T> traits_t;
1396 step != sqrt_max_steps && norm_M_1 > tol2;
1400 return traits_t::NaN();
1402 norm_M_1 =
norm(M - Scalar_T(1));
1404 if (norm_M_1 > tol2)
1405 return traits_t::NaN();
1414 template<
typename Scalar_T >
1417 static const int array_size = 14;
1418 static const Scalar_T array[array_size];
1423 template<
typename Scalar_T >
1426 static const int array_size = 14;
1427 static const Scalar_T array[array_size];
1430 template<
typename Scalar_T >
1431 const Scalar_T pade_sqrt_a<Scalar_T>::array[pade_sqrt_a<Scalar_T>::array_size] =
1433 1.0, 27.0/4.0, 81.0/4.0, 2277.0/64.0,
1434 10395.0/256.0, 32319.0/1024.0, 8721.0/512.0, 26163.0/4096.0,
1435 53703.0/32768.0, 36465.0/131072.0, 3861.0/131072.0, 7371.0/4194304.0,
1436 819.0/16777216.0, 27.0/67108864.0
1438 template<
typename Scalar_T >
1439 const Scalar_T pade_sqrt_b<Scalar_T>::array[pade_sqrt_b<Scalar_T>::array_size] =
1441 1.0, 25.0/4.0, 69.0/4.0, 1771.0/64.0,
1442 7315.0/256.0, 20349.0/1024.0, 4845.0/512.0, 12597.0/4096.0,
1443 21879.0/32768.0, 12155.0/131072.0, 1001.0/131072.0, 1365.0/4194304.0,
1444 91.0/16777216.0, 1.0/67108864.0
1448 struct pade_sqrt_a<float>
1450 static const int array_size = 10;
1451 static const float array[array_size];
1454 struct pade_sqrt_b<float>
1456 static const int array_size = 10;
1457 static const float array[array_size];
1459 const float pade_sqrt_a<float>::array[pade_sqrt_a<float>::array_size] =
1461 1.0, 19.0/4.0, 19.0/2.0, 665.0/64.0,
1462 1729.0/256.0, 2717.0/1024.0, 627.0/1024.0, 627.0/8192.0,
1463 285.0/65536.0, 19.0/262144.0
1465 const float pade_sqrt_b<float>::array[pade_sqrt_a<float>::array_size] =
1467 1.0, 17.0/4.0, 15.0/2.0, 455.0/64.0,
1468 1001.0/256.0, 1287.0/1024.0, 231.0/1024.0, 165.0/8192.0,
1469 45.0/65536, 1.0/262144.0
1473 struct pade_sqrt_a<long double>
1475 static const int array_size = 18;
1476 static const long double array[array_size];
1479 struct pade_sqrt_b<long double>
1481 static const int array_size = 18;
1482 static const long double array[array_size];
1484 const long double pade_sqrt_a<long double>::array[pade_sqrt_a<long double>::array_size] =
1486 1.0L, 35.0L/4.0L, 35.0L, 5425.0L/64.0L,
1487 35525.0L/256.0L, 166257.0L/1024.0L, 143325.0L/1024.0L, 740025.0L/8192.0L,
1488 2877875.0L/65536.0L, 4206125.0L/262144.0L, 572033.0L/131072.0L, 1820105.0L/2097152.0L,
1489 1028755.0L/8388608.0L, 395675.0L/33554432.0L, 24225.0L/33554432.0L, 6783.0L/268435456.0L,
1490 1785.0L/4294967296.0L, 35.0L/17179869184.0L
1492 const long double pade_sqrt_b<long double>::array[pade_sqrt_a<long double>::array_size] =
1494 1.0L, 33.0L/4.0L, 31.0L, 4495.0L/64.0L,
1495 27405.0L/256.0L, 118755.0L/1024.0L, 94185.0L/1024.0L, 444015.0L/8192.0L,
1496 1562275.0L/65536.0L, 2042975.0L/262144.0L, 245157.0L/131072.0L, 676039.0L/2097152.0L,
1497 323323.0L/8388608.0L, 101745.0L/33554432.0L, 4845.0L/33554432.0L, 969.0L/268435456.0L,
1498 153.0L/4294967296.0L, 1.0L/17179869184.0L
1501 #if defined(_GLUCAT_USE_QD)
1503 struct pade_sqrt_a<dd_real>
1505 static const int array_size = 22;
1506 static const dd_real array[array_size];
1509 struct pade_sqrt_b<dd_real>
1511 static const int array_size = 22;
1512 static const dd_real array[array_size];
1514 const dd_real pade_sqrt_a<dd_real>::array[pade_sqrt_a<dd_real>::array_size] =
1516 dd_real(
"1"), dd_real(
"43")/dd_real(
"4"),
1517 dd_real(
"215")/dd_real(
"4"), dd_real(
"10621")/dd_real(
"64"),
1518 dd_real(
"90687")/dd_real(
"256"), dd_real(
"567987")/dd_real(
"1024"),
1519 dd_real(
"168861")/dd_real(
"256"), dd_real(
"1246355")/dd_real(
"2048"),
1520 dd_real(
"7228859")/dd_real(
"16384"), dd_real(
"16583853")/dd_real(
"65536"),
1521 dd_real(
"7538115")/dd_real(
"65536"), dd_real(
"173376645")/dd_real(
"4194304"),
1522 dd_real(
"195747825")/dd_real(
"16777216"), dd_real(
"171655785")/dd_real(
"67108864"),
1523 dd_real(
"14375115")/dd_real(
"33554432"), dd_real(
"14375115")/dd_real(
"268435456"),
1524 dd_real(
"20764055")/dd_real(
"4294967296"), dd_real(
"5167525")/dd_real(
"17179869184"),
1525 dd_real(
"206701")/dd_real(
"17179869184"), dd_real(
"76153")/dd_real(
"274877906944"),
1526 dd_real(
"3311")/dd_real(
"1099511627776") , dd_real(
"43")/dd_real(
"4398046511104")
1528 const dd_real pade_sqrt_b<dd_real>::array[pade_sqrt_a<dd_real>::array_size] =
1530 dd_real(
"1"), dd_real(
"41")/dd_real(
"4"),
1531 dd_real(
"195")/dd_real(
"4"), dd_real(
"9139")/dd_real(
"64"),
1532 dd_real(
"73815")/dd_real(
"256"), dd_real(
"435897")/dd_real(
"1024"),
1533 dd_real(
"121737")/dd_real(
"256"), dd_real(
"840565")/dd_real(
"2048"),
1534 dd_real(
"4539051")/dd_real(
"16384"), dd_real(
"9641775")/dd_real(
"65536"),
1535 dd_real(
"4032015")/dd_real(
"65536"), dd_real(
"84672315")/dd_real(
"4194304"),
1536 dd_real(
"86493225")/dd_real(
"16777216"), dd_real(
"67863915")/dd_real(
"67108864"),
1537 dd_real(
"5014575")/dd_real(
"33554432"), dd_real(
"4345965")/dd_real(
"268435456"),
1538 dd_real(
"5311735")/dd_real(
"4294967296"), dd_real(
"1081575")/dd_real(
"17179869184"),
1539 dd_real(
"33649")/dd_real(
"17179869184"), dd_real(
"8855")/dd_real(
"274877906944"),
1540 dd_real(
"231")/dd_real(
"1099511627776"), dd_real(
"1")/dd_real(
"4398046511104")
1544 struct pade_sqrt_a<qd_real>
1546 static const int array_size = 34;
1547 static const qd_real array[array_size];
1550 struct pade_sqrt_b<qd_real>
1552 static const int array_size = 34;
1553 static const qd_real array[array_size];
1555 const qd_real pade_sqrt_a<qd_real>::array[pade_sqrt_a<qd_real>::array_size] =
1557 qd_real(
"1"), qd_real(
"67")/qd_real(
"4"),
1558 qd_real(
"134"), qd_real(
"43617")/qd_real(
"64"),
1559 qd_real(
"633485")/qd_real(
"256"), qd_real(
"6992857")/qd_real(
"1024"),
1560 qd_real(
"15246721")/qd_real(
"1024"), qd_real(
"215632197")/qd_real(
"8192"),
1561 qd_real(
"2518145487")/qd_real(
"65536"), qd_real(
"12301285425")/qd_real(
"262144"),
1562 qd_real(
"6344873535")/qd_real(
"131072"), qd_real(
"89075432355")/qd_real(
"2097152"),
1563 qd_real(
"267226297065")/qd_real(
"8388608"), qd_real(
"687479618945")/qd_real(
"33554432"),
1564 qd_real(
"379874182975")/qd_real(
"33554432"), qd_real(
"1443521895305")/qd_real(
"268435456"),
1565 qd_real(
"9425348845815")/qd_real(
"4294967296"), qd_real(
"13195488384141")/qd_real(
"17179869184"),
1566 qd_real(
"987417498133")/qd_real(
"4294967296"), qd_real(
"8055248011085")/qd_real(
"137438953472"),
1567 qd_real(
"6958363175533")/qd_real(
"549755813888"), qd_real(
"5056698705201")/qd_real(
"2199023255552"),
1568 qd_real(
"766166470485")/qd_real(
"2199023255552"), qd_real(
"766166470485")/qd_real(
"17592186044416"),
1569 qd_real(
"623623871325")/qd_real(
"140737488355328"), qd_real(
"203123203803")/qd_real(
"562949953421312"),
1570 qd_real(
"6478601247")/qd_real(
"281474976710656"), qd_real(
"5038912081")/qd_real(
"4503599627370496"),
1571 qd_real(
"719844583")/qd_real(
"18014398509481984"), qd_real(
"71853815")/qd_real(
"72057594037927936"),
1572 qd_real(
"1165197")/qd_real(
"72057594037927936"), qd_real(
"87703")/qd_real(
"576460752303423488"),
1573 qd_real(
"12529")/qd_real(
"18446744073709551616"), qd_real(
"67")/qd_real(
"73786976294838206464")
1575 const qd_real pade_sqrt_b<qd_real>::array[pade_sqrt_a<qd_real>::array_size] =
1577 qd_real(
"1"), qd_real(
"65")/qd_real(
"4"),
1578 qd_real(
"126"), qd_real(
"39711")/qd_real(
"64"),
1579 qd_real(
"557845")/qd_real(
"256"), qd_real(
"5949147")/qd_real(
"1024"),
1580 qd_real(
"12515965")/qd_real(
"1024"), qd_real(
"170574723")/qd_real(
"8192"),
1581 qd_real(
"1916797311")/qd_real(
"65536"), qd_real(
"8996462475")/qd_real(
"262144"),
1582 qd_real(
"4450881435")/qd_real(
"131072"), qd_real(
"59826782925")/qd_real(
"2097152"),
1583 qd_real(
"171503444385")/qd_real(
"8388608"), qd_real(
"420696483235")/qd_real(
"33554432"),
1584 qd_real(
"221120793075")/qd_real(
"33554432"), qd_real(
"797168807855")/qd_real(
"268435456"),
1585 qd_real(
"4923689695575")/qd_real(
"4294967296"), qd_real(
"6499270398159")/qd_real(
"17179869184"),
1586 qd_real(
"456864812569")/qd_real(
"4294967296"), qd_real(
"3486599885395")/qd_real(
"137438953472"),
1587 qd_real(
"2804116503573")/qd_real(
"549755813888"), qd_real(
"1886827875075")/qd_real(
"2199023255552"),
1588 qd_real(
"263012370465")/qd_real(
"2199023255552"), qd_real(
"240141729555")/qd_real(
"17592186044416"),
1589 qd_real(
"176848560525")/qd_real(
"140737488355328"), qd_real(
"51538723353")/qd_real(
"562949953421312"),
1590 qd_real(
"1450433115")/qd_real(
"281474976710656"), qd_real(
"977699359")/qd_real(
"4503599627370496"),
1591 qd_real(
"118183439")/qd_real(
"18014398509481984"), qd_real(
"9652005")/qd_real(
"72057594037927936"),
1592 qd_real(
"121737")/qd_real(
"72057594037927936"), qd_real(
"6545")/qd_real(
"576460752303423488"),
1593 qd_real(
"561")/qd_real(
"18446744073709551616"), qd_real(
"1")/qd_real(
"73786976294838206464")
1601 template<
typename Scalar_T, const index_t LO, const index_t HI >
1602 const matrix_multi<Scalar_T,LO,HI>
1612 return traits_t::NaN();
1620 typedef typename traits_t::demoted::type demoted_scalar_t;
1623 const demoted_multivector_t& demoted_val = demoted_multivector_t(val);
1624 const demoted_multivector_t& demoted_i = demoted_multivector_t(
i);
1631 typedef typename traits_t::promoted::type promoted_scalar_t;
1634 const promoted_multivector_t& promoted_val = promoted_multivector_t(val);
1635 const promoted_multivector_t& promoted_i = promoted_multivector_t(
i);
1646 template<
typename Scalar_T, const index_t LO, const index_t HI >
1647 const matrix_multi<Scalar_T,LO,HI>
1657 return traits_t::NaN();
1661 const Scalar_T realval = val.scalar();
1664 if (realval < Scalar_T(0))
1672 #if !defined(_GLUCAT_USE_EIGENVALUES)
1673 const multivector_t val2 = val*val;
1674 const Scalar_T real_val2 = val2.scalar();
1675 if (val2 == real_val2 && real_val2 > Scalar_T(0))
1680 const Scalar_T scale =
1681 (realval != Scalar_T(0) &&
norm(val/realval - Scalar_T(1)) < Scalar_T(1))
1683 : (realval < Scalar_T(0))
1687 if (traits_t::isNaN_or_isInf(sqrt_scale))
1688 return traits_t::NaN();
1691 multivector_t rescale = sqrt_scale;
1692 if (scale < Scalar_T(0))
1693 rescale =
i * sqrt_scale;
1695 const multivector_t& unitval = val / scale;
1696 const Scalar_T max_norm = Scalar_T(1.0/4.0);
1698 #if defined(_GLUCAT_USE_EIGENVALUES)
1699 multivector_t scaled_result;
1700 typedef typename multivector_t::matrix_t matrix_t;
1707 scaled_result =
matrix_sqrt(-
i * unitval,
i) * (
i + Scalar_T(1)) / sqrt_2;
1717 (
norm(unitval - Scalar_T(1)) < max_norm)
1720 pade_sqrt_a<Scalar_T>::array,
1721 pade_sqrt_b<Scalar_T>::array,
1722 unitval - Scalar_T(1))
1727 if (scaled_result.isnan())
1728 return traits_t::NaN();
1730 return scaled_result * rescale;
1732 const multivector_t& scaled_result =
1733 (
norm(unitval - Scalar_T(1)) < max_norm)
1736 pade_sqrt_a<Scalar_T>::array,
1737 pade_sqrt_b<Scalar_T>::array,
1738 unitval - Scalar_T(1))
1741 if (scaled_result.isnan())
1743 if (
inv(unitval).isnan())
1744 return traits_t::NaN();
1746 const multivector_t& mi_unitval = -
i * unitval;
1748 const multivector_t& scaled_mi_result =
1749 (
norm(mi_unitval - Scalar_T(1)) < max_norm)
1752 pade_sqrt_a<Scalar_T>::array,
1753 pade_sqrt_b<Scalar_T>::array,
1754 mi_unitval - Scalar_T(1))
1757 if (scaled_mi_result.isnan())
1758 return traits_t::NaN();
1760 return scaled_mi_result * rescale * (
i + Scalar_T(1)) / sqrt_2;
1763 return scaled_result * rescale;
1771 template<
typename Scalar_T >
1774 static const int array_size = 14;
1775 static const Scalar_T array[array_size];
1780 template<
typename Scalar_T >
1783 static const int array_size = 14;
1784 static const Scalar_T array[array_size];
1786 template<
typename Scalar_T >
1787 const Scalar_T pade_log_a<Scalar_T>::array[pade_log_a<Scalar_T>::array_size] =
1789 0.0, 1.0, 6.0, 4741.0/300.0,
1790 1441.0/60.0, 107091.0/4600.0, 8638.0/575.0, 263111.0/40250.0,
1791 153081.0/80500.0, 395243.0/1101240.0, 28549.0/688275.0, 605453.0/228813200.0,
1792 785633.0/10296594000.0, 1145993.0/1873980108000.0
1794 template<
typename Scalar_T >
1795 const Scalar_T pade_log_b<Scalar_T>::array[pade_log_b<Scalar_T>::array_size] =
1797 1.0, 13.0/2.0, 468.0/25.0, 1573.0/50.0,
1798 1573.0/46.0, 11583.0/460.0, 10296.0/805.0, 2574.0/575.0,
1799 11583.0/10925.0, 143.0/874.0, 572.0/37145.0, 117.0/148580.0,
1800 13.0/742900.0, 1.0/10400600.0
1804 struct pade_log_a<float>
1806 static const int array_size = 10;
1807 static const float array[array_size];
1810 struct pade_log_b<float>
1812 static const int array_size = 10;
1813 static const float array[array_size];
1815 const float pade_log_a<float>::array[pade_log_a<float>::array_size] =
1817 0.0, 1.0, 4.0, 1337.0/204.0,
1818 385.0/68.0, 1879.0/680.0, 193.0/255.0, 197.0/1820.0,
1819 419.0/61880.0, 7129.0/61261200.0
1821 const float pade_log_b<float>::array[pade_log_a<float>::array_size] =
1823 1.0, 9.0/2.0, 144.0/17.0, 147.0/17.0,
1824 441.0/85.0, 63.0/34.0, 84.0/221.0, 9.0/221.0,
1825 9.0/4862.0, 1.0/48620.0
1829 struct pade_log_a<long double>
1831 static const int array_size = 18;
1832 static const long double array[array_size];
1835 struct pade_log_b<long double>
1837 static const int array_size = 18;
1838 static const long double array[array_size];
1840 const long double pade_log_a<long double>::array[pade_log_a<long double>::array_size] =
1842 0.0L, 1.0L, 8.0L, 3835.0L/132.0L,
1843 8365.0L/132.0L, 11363807.0L/122760.0L, 162981.0L/1705.0L, 9036157.0L/125860.0L,
1844 18009875.0L/453096.0L, 44211925.0L/2718576.0L, 4149566.0L/849555.0L, 16973929.0L/16020180.0L,
1845 172459.0L/1068012.0L, 116317061.0L/7025382936.0L, 19679783.0L/18441630207.0L, 23763863.0L/614721006900.0L,
1846 50747.0L/79318839600.0L, 42142223.0L/14295951736466400.0L
1848 const long double pade_log_b<long double>::array[pade_log_a<long double>::array_size] =
1850 1.0L, 17.0L/2.0L, 1088.0L/33.0L, 850.0L/11.0L,
1851 41650.0L/341.0L, 140777.0L/1023.0L, 1126216.0L/9889.0L, 63206.0L/899.0L,
1852 790075.0L/24273.0L, 60775.0L/5394.0L, 38896.0L/13485.0L, 21658.0L/40455.0L,
1853 21658.0L/310155.0L, 4165.0L/682341.0L, 680.0L/2047023.0L, 34.0L/3411705.0L,
1854 17.0L/129644790.0L, 1.0L/2333606220
1856 #if defined(_GLUCAT_USE_QD)
1858 struct pade_log_a<dd_real>
1860 static const int array_size = 22;
1861 static const dd_real array[array_size];
1864 struct pade_log_b<dd_real>
1866 static const int array_size = 22;
1867 static const dd_real array[array_size];
1869 const dd_real pade_log_a<dd_real>::array[pade_log_a<dd_real>::array_size] =
1871 dd_real(
"0"), dd_real(
"1"),
1872 dd_real(
"10"), dd_real(
"22781")/dd_real(
"492"),
1873 dd_real(
"21603")/dd_real(
"164"), dd_real(
"5492649")/dd_real(
"21320"),
1874 dd_real(
"978724")/dd_real(
"2665"), dd_real(
"4191605")/dd_real(
"10619"),
1875 dd_real(
"12874933")/dd_real(
"39442"), dd_real(
"11473457")/dd_real(
"54612"),
1876 dd_real(
"2406734")/dd_real(
"22755"), dd_real(
"166770367")/dd_real(
"4004880"),
1877 dd_real(
"30653165")/dd_real(
"2402928"), dd_real(
"647746389")/dd_real(
"215195552"),
1878 dd_real(
"25346331")/dd_real(
"47074027"), dd_real(
"278270613")/dd_real(
"3900419380"),
1879 dd_real(
"105689791")/dd_real(
"15601677520"), dd_real(
"606046475")/dd_real(
"1379188292768"),
1880 dd_real(
"969715")/dd_real(
"53502994116"), dd_real(
"11098301")/dd_real(
"26204577562592"),
1881 dd_real(
"118999")/dd_real(
"26204577562592"), dd_real(
"18858053")/dd_real(
"1392249205900512960")
1883 const dd_real pade_log_b<dd_real>::array[pade_log_a<dd_real>::array_size] =
1885 dd_real(
"1"), dd_real(
"21")/dd_real(
"2"),
1886 dd_real(
"2100")/dd_real(
"41"), dd_real(
"12635")/dd_real(
"82"),
1887 dd_real(
"341145")/dd_real(
"1066"), dd_real(
"1037799")/dd_real(
"2132"),
1888 dd_real(
"11069856")/dd_real(
"19721"), dd_real(
"9883800")/dd_real(
"19721"),
1889 dd_real(
"6918660")/dd_real(
"19721"), dd_real(
"293930")/dd_real(
"1517"),
1890 dd_real(
"1410864")/dd_real(
"16687"), dd_real(
"88179")/dd_real(
"3034"),
1891 dd_real(
"734825")/dd_real(
"94054"), dd_real(
"305235")/dd_real(
"188108"),
1892 dd_real(
"348840")/dd_real(
"1363783"), dd_real(
"40698")/dd_real(
"1363783"),
1893 dd_real(
"6783")/dd_real(
"2727566"), dd_real(
"9975")/dd_real(
"70916716"),
1894 dd_real(
"266")/dd_real(
"53187537"), dd_real(
"7")/dd_real(
"70916716"),
1895 dd_real(
"7")/dd_real(
"8155422340"), dd_real(
"1")/dd_real(
"538257874440")
1899 struct pade_log_a<qd_real>
1901 static const int array_size = 34;
1902 static const qd_real array[array_size];
1905 struct pade_log_b<qd_real>
1907 static const int array_size = 34;
1908 static const qd_real array[array_size];
1910 const qd_real pade_log_a<qd_real>::array[pade_log_a<qd_real>::array_size] =
1912 qd_real(
"0"), qd_real(
"1"),
1913 qd_real(
"16"), qd_real(
"95201")/qd_real(
"780"),
1914 qd_real(
"30721")/qd_real(
"52"), qd_real(
"7416257")/qd_real(
"3640"),
1915 qd_real(
"1039099")/qd_real(
"195"), qd_real(
"6097772319")/qd_real(
"555100"),
1916 qd_real(
"1564058073")/qd_real(
"85400"), qd_real(
"30404640205")/qd_real(
"1209264"),
1917 qd_real(
"725351278")/qd_real(
"25193"), qd_real(
"4092322670789")/qd_real(
"147429436"),
1918 qd_real(
"4559713849589")/qd_real(
"201040140"), qd_real(
"5049361751189")/qd_real(
"320023080"),
1919 qd_real(
"74979677195")/qd_real(
"8000577"), qd_real(
"16569850691873")/qd_real(
"3481514244"),
1920 qd_real(
"1065906022369")/qd_real(
"515779888"), qd_real(
"335956770855841")/qd_real(
"438412904800"),
1921 qd_real(
"1462444287585964")/qd_real(
"6041877844275"), qd_real(
"397242326339851")/qd_real(
"6122436215532"),
1922 qd_real(
"64211291334131")/qd_real(
"4373168725380"), qd_real(
"142322343550859")/qd_real(
"51080680851480"),
1923 qd_real(
"154355972958659")/qd_real(
"351179680853925"), qd_real(
"167483568676259")/qd_real(
"2937139148960100"),
1924 qd_real(
"4230788929433")/qd_real(
"704913395750424"), qd_real(
"197968763176019")/qd_real(
"392923948371995600"),
1925 qd_real(
"10537522306718")/qd_real(
"319250708052246425"), qd_real(
"236648286272519")/qd_real(
"144249197475035425500"),
1926 qd_real(
"260715545088119")/qd_real(
"4375558990076074573500"), qd_real(
"289596255666839")/qd_real(
"192874640282553367199880"),
1927 qd_real(
"8802625510547")/qd_real(
"361639950529787563499775"), qd_real(
"373831661521439")/qd_real(
"1659204093030665341336967700"),
1928 qd_real(
"446033437968239")/qd_real(
"464577146048586295574350956000"), qd_real(
"53676090078349")/qd_real(
"47386868896955802148583797512000")
1930 const qd_real pade_log_b<qd_real>::array[pade_log_a<qd_real>::array_size] =
1932 qd_real(
"1"), qd_real(
"33")/qd_real(
"2"),
1933 qd_real(
"8448")/qd_real(
"65"), qd_real(
"42284")/qd_real(
"65"),
1934 qd_real(
"211420")/qd_real(
"91"), qd_real(
"573562")/qd_real(
"91"),
1935 qd_real(
"32119472")/qd_real(
"2379"), qd_real(
"92917044")/qd_real(
"3965"),
1936 qd_real(
"603960786")/qd_real(
"17995"), qd_real(
"144626625")/qd_real(
"3599"),
1937 qd_real(
"2776831200")/qd_real(
"68381"), qd_real(
"16692542100")/qd_real(
"478667"),
1938 qd_real(
"12241197540")/qd_real(
"478667"), qd_real(
"1098569010")/qd_real(
"68381"),
1939 qd_real(
"31387686000")/qd_real(
"3624193"), qd_real(
"9939433900")/qd_real(
"2479711"),
1940 qd_real(
"67091178825")/qd_real(
"42155087"), qd_real(
"2683647153")/qd_real(
"4959422"),
1941 qd_real(
"19083713088")/qd_real(
"121505839"), qd_real(
"4708152900")/qd_real(
"121505839"),
1942 qd_real(
"941630580")/qd_real(
"116546417"), qd_real(
"88704330")/qd_real(
"62755763"),
1943 qd_real(
"12902448")/qd_real(
"62755763"), qd_real(
"1542684")/qd_real(
"62755763"),
1944 qd_real(
"6427850")/qd_real(
"2698497809"), qd_real(
"3471039")/qd_real(
"18889484663"),
1945 qd_real(
"8544096")/qd_real(
"774468871183"), qd_real(
"39556")/qd_real(
"79027435835"),
1946 qd_real(
"118668")/qd_real(
"7191496660985"), qd_real(
"10230")/qd_real(
"27327687311743"),
1947 qd_real(
"5456")/qd_real(
"1011124430534491"), qd_real(
"44")/qd_real(
"1011124430534491"),
1948 qd_real(
"11")/qd_real(
"70778710137414370"), qd_real(
"1")/qd_real(
"7219428434016265740")
1955 template<
typename Scalar_T, const index_t LO, const index_t HI >
1957 const matrix_multi<Scalar_T,LO,HI>
1966 if (val == Scalar_T(0) || val.isnan())
1967 return traits_t::NaN();
1969 return pade_approx(pade_log_a<Scalar_T>::array_size,
1970 pade_log_a<Scalar_T>::array,
1971 pade_log_b<Scalar_T>::array,
1976 template<
typename Scalar_T, const index_t LO, const index_t HI >
1978 const matrix_multi<Scalar_T,LO,HI>
1984 if (val == Scalar_T(0) || val.isnan())
1985 return traits_t::NaN();
1987 typedef std::numeric_limits<Scalar_T> limits_t;
1990 static const Scalar_T max_outer_norm = Scalar_T(6.0/limits_t::digits);
1991 multivector_t Y = val;
1992 multivector_t E = Scalar_T(0);
1995 Scalar_T pow_2_outer_step = Scalar_T(1);
1996 Scalar_T pow_4_outer_step = Scalar_T(1);
1997 for (outer_step = 0, norm_Y_1 =
norm(Y - Scalar_T(1));
1999 ++outer_step, norm_Y_1 =
norm(Y - Scalar_T(1)))
2001 if (Y == Scalar_T(0) || Y.isnan())
2002 return traits_t::NaN();
2005 multivector_t M = Y;
2009 norm(M - Scalar_T(1)) * pow_4_outer_step > max_inner_norm;
2013 E += (M - Scalar_T(1)) * pow_2_outer_step;
2014 pow_2_outer_step *= Scalar_T(2);
2015 pow_4_outer_step *= Scalar_T(4);
2018 return traits_t::NaN();
2020 return pade_log(Y) * pow_2_outer_step - E;
2024 template<
typename Scalar_T, const index_t LO, const index_t HI >
2025 const matrix_multi<Scalar_T,LO,HI>
2030 if (val == Scalar_T(0) || val.isnan())
2031 return traits_t::NaN();
2039 typedef typename traits_t::demoted::type demoted_scalar_t;
2042 const demoted_multivector_t& demoted_val = demoted_multivector_t(val);
2043 const demoted_multivector_t& demoted_i = demoted_multivector_t(
i);
2050 typedef typename traits_t::promoted::type promoted_scalar_t;
2053 const promoted_multivector_t& promoted_val = promoted_multivector_t(val);
2054 const promoted_multivector_t& promoted_i = promoted_multivector_t(
i);
2065 template<
typename Scalar_T, const index_t LO, const index_t HI >
2066 const matrix_multi<Scalar_T,LO,HI>
2073 if (val == Scalar_T(0) || val.isnan())
2074 return traits_t::NaN();
2077 const Scalar_T realval = val.scalar();
2080 if (realval < Scalar_T(0))
2086 #if !defined(_GLUCAT_USE_EIGENVALUES)
2087 const multivector_t val2 = val*val;
2088 const Scalar_T real_val2 = val2.scalar();
2089 if (val2 == real_val2 && real_val2 > 0)
2093 const Scalar_T max_norm = Scalar_T(1.0/9.0);
2094 const Scalar_T scale =
2095 (realval != Scalar_T(0) &&
norm(val/realval - Scalar_T(1)) < max_norm)
2097 : (realval < Scalar_T(0))
2100 if (scale == Scalar_T(0))
2101 return traits_t::NaN();
2104 multivector_t rescale = log_scale;
2105 if (scale < Scalar_T(0))
2106 rescale =
i *
pi + log_scale;
2107 const multivector_t unitval = val/scale;
2109 return traits_t::NaN();
2110 #if defined(_GLUCAT_USE_EIGENVALUES)
2111 multivector_t scaled_result;
2112 typedef typename multivector_t::matrix_t matrix_t;
2132 multivector_t scaled_result =
cascade_log(unitval);
2134 if (scaled_result.isnan())
2135 return traits_t::NaN();
2137 return scaled_result + rescale;
2141 template<
typename Scalar_T, const index_t LO, const index_t HI >
2142 const matrix_multi<Scalar_T,LO,HI>
2147 return traits_t::NaN();
2149 const Scalar_T s =
scalar(val);
2157 typedef typename traits_t::demoted::type demoted_scalar_t;
2160 const demoted_multivector_t& demoted_val = demoted_multivector_t(val);
2166 typedef typename traits_t::promoted::type promoted_scalar_t;
2169 const promoted_multivector_t& promoted_val = promoted_multivector_t(val);
2178 #endif // _GLUCAT_MATRIX_MULTI_IMP_H