1#ifndef _GLUCAT_MATRIX_MULTI_IMP_H
2#define _GLUCAT_MATRIX_MULTI_IMP_H
40# if defined(_GLUCAT_GCC_IGNORE_UNUSED_LOCAL_TYPEDEFS)
41# pragma GCC diagnostic push
42# pragma GCC diagnostic ignored "-Wunused-local-typedefs"
44# if defined(_GLUCAT_HAVE_BOOST_SERIALIZATION_ARRAY_WRAPPER_H)
45# include <boost/serialization/array_wrapper.hpp>
47#include <boost/numeric/ublas/matrix.hpp>
48#include <boost/numeric/ublas/matrix_expression.hpp>
49#include <boost/numeric/ublas/matrix_proxy.hpp>
50#include <boost/numeric/ublas/matrix_sparse.hpp>
51#include <boost/numeric/ublas/operation.hpp>
52#include <boost/numeric/ublas/operation_sparse.hpp>
53#include <boost/numeric/ublas/triangular.hpp>
54#include <boost/numeric/ublas/lu.hpp>
55#include <boost/numeric/ublas/io.hpp>
56# if defined(_GLUCAT_GCC_IGNORE_UNUSED_LOCAL_TYPEDEFS)
57# pragma GCC diagnostic pop
76 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
80 {
return "matrix_multi"; }
90 static const std::array<int, 8> offset_log2_dim = {0, 1, 0, 1, 1, 2, 1, 1};
92 return (p+q)/2 + offset_log2_dim[bott];
97 template<
typename Matrix_Index_T, const index_t LO, const index_t HI >
102 {
return 1 <<
offset_level(sub.count_pos(), sub.count_neg()); }
105 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
113 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
114 template<
typename Other_Scalar_T >
125 val_it2 = val_it1.begin();
126 val_it2 != val_it1.end();
132 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
133 template<
typename Other_Scalar_T >
143 this->
m_matrix.resize(dim, dim,
false);
150 val_it2 = val_it1.begin();
151 val_it2 != val_it1.end();
158 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
170 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
176 this->
m_matrix.resize(dim, dim,
false);
178 *
this +=
term_t(ist, crd);
182 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
187 if (!prechecked && (ist | frm) != frm)
188 throw error_t(
"multivector_t(ist,crd,frm): cannot initialize with value outside of frame");
190 this->
m_matrix.resize(dim, dim,
false);
192 *
this +=
term_t(ist, crd);
196 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
202 this->
m_matrix.resize(dim, dim,
false);
208 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
214 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
221 throw error_t(
"multivector_t(vec,frm): cannot initialize with vector not matching frame");
223 this->
m_matrix.resize(dim, dim,
false);
225 auto idx = frm.
min();
226 const auto frm_end = frm.
max()+1;
227 for (
auto& crd : vec)
232 idx != frm_end && !frm[idx];
239 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
245 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
251 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
252 template<
typename Other_Scalar_T >
257 if (val.size() >= Tune_P::fast_size_threshold)
260 *this = val.template fast_matrix_multi<Scalar_T>(this->m_frame);
266 this->
m_matrix.resize(dim, dim,
false);
269 for (
auto& val_term : val)
274 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
275 template<
typename Other_Scalar_T >
280 const auto our_frame = val.frame() | frm;
281 if (val.size() >= Tune_P::fast_size_threshold)
291 this->
m_matrix.resize(dim, dim,
false);
294 for (
auto& val_term : val)
299 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
300 template<
typename Matrix_T >
308 mtx_it1 = mtx.begin1();
309 mtx_it1 != mtx.end1();
312 mtx_it2 = mtx_it1.begin();
313 mtx_it2 != mtx_it1.end();
319 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
326 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
340 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
348 using framed_multi_t =
typename multivector_t::framed_multi_t;
350 index_set_t our_frame = lhs.m_frame | rhs.m_frame;
351 framed_multi_t framed_lhs;
352 framed_multi_t framed_rhs;
353 if ((lhs.m_frame != our_frame) || (rhs.m_frame != our_frame))
356 framed_lhs = framed_multi_t(lhs);
357 framed_rhs = framed_multi_t(rhs);
358 our_frame |= framed_lhs.frame() | framed_rhs.frame();
361 if (lhs.m_frame != our_frame)
362 lhs_reframed = multivector_t(framed_lhs, our_frame,
true);
363 if (rhs.m_frame != our_frame)
364 rhs_reframed = multivector_t(framed_rhs, our_frame,
true);
369 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
372 operator== (
const multivector_t& rhs)
const ->
bool
379 multivector_t lhs_reframed;
380 multivector_t rhs_reframed;
381 const index_set_t our_frame =
reframe(*
this, rhs, lhs_reframed, rhs_reframed);
382 const multivector_t& lhs_ref = (this->m_frame == our_frame)
385 const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
389 return ublas::norm_inf(lhs_ref.m_matrix - rhs_ref.m_matrix) == 0;
393 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
397 operator== (
const Scalar_T& scr)
const ->
bool
399 if (scr != Scalar_T(0))
400 return *
this == multivector_t(framed_multi_t(scr), this->m_frame,
true);
401 else if (ublas::norm_inf(this->m_matrix) != 0)
405 const matrix_index_t dim = this->m_matrix.size1();
406 return !(dim == 1 && this->
isnan());
411 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
419 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
423 operator+= (
const multivector_t& rhs) -> multivector_t&
427 return *
this *= Scalar_T(2);
430 multivector_t rhs_reframed;
431 const index_set_t our_frame =
reframe(*
this, rhs, *
this, rhs_reframed);
432 const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
436 noalias(this->m_matrix) += rhs_ref.m_matrix;
441 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
445 operator-= (
const Scalar_T& scr) -> multivector_t&
446 {
return *
this += term_t(index_set_t(), -scr); }
449 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
453 operator-= (
const multivector_t& rhs) -> multivector_t&
457 return *
this = Scalar_T(0);
460 multivector_t rhs_reframed;
461 const index_set_t our_frame =
reframe(*
this, rhs, *
this, rhs_reframed);
462 const multivector_t& rhs_ref = (
rhs.m_frame == our_frame)
466 noalias(this->m_matrix) -= rhs_ref.m_matrix;
471 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
475 operator- () const -> const multivector_t
476 {
return multivector_t(-(this->m_matrix), this->m_frame); }
479 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
483 operator*= (
const Scalar_T& scr) -> multivector_t&
487 if (traits_t::isNaN_or_isInf(scr) || this->
isnan())
488 return *
this = traits_t::NaN();
489 if (scr == Scalar_T(0))
492 this->m_matrix *= scr;
497 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
500 operator* (
const matrix_multi<Scalar_T,LO,HI,Tune_P>& lhs,
const matrix_multi<Scalar_T,LO,HI,Tune_P>& rhs) ->
const matrix_multi<Scalar_T,LO,HI,Tune_P>
503 using index_set_t =
typename multivector_t::index_set_t;
505 if (lhs.isnan() || rhs.isnan())
509 multivector_t lhs_reframed;
510 multivector_t rhs_reframed;
511 const index_set_t our_frame =
reframe(lhs, rhs, lhs_reframed, rhs_reframed);
512 const multivector_t& lhs_ref = (lhs.m_frame == our_frame)
515 const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
519 using matrix_t =
typename multivector_t::matrix_t;
520 using matrix_index_t =
typename matrix_t::size_type;
522 const matrix_index_t dim = lhs_ref.m_matrix.size1();
523 multivector_t result = multivector_t(matrix_t(dim, dim), our_frame);
524 result.m_matrix.clear();
525 ublas::axpy_prod(lhs_ref.m_matrix, rhs_ref.m_matrix, result.m_matrix,
true);
530 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
534 operator*= (
const multivector_t& rhs) -> multivector_t&
535 {
return *
this = *
this * rhs; }
538 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
541 operator^ (
const matrix_multi<Scalar_T,LO,HI,Tune_P>& lhs,
const matrix_multi<Scalar_T,LO,HI,Tune_P>& rhs) ->
const matrix_multi<Scalar_T,LO,HI,Tune_P>
544 using framed_multi_t =
typename multivector_t::framed_multi_t;
545 return framed_multi_t(lhs) ^ framed_multi_t(rhs);
549 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
553 operator^= (
const multivector_t& rhs) -> multivector_t&
554 {
return *
this = *
this ^ rhs; }
557 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
560 operator& (
const matrix_multi<Scalar_T,LO,HI,Tune_P>& lhs,
const matrix_multi<Scalar_T,LO,HI,Tune_P>& rhs) ->
const matrix_multi<Scalar_T,LO,HI,Tune_P>
563 using framed_multi_t =
typename multivector_t::framed_multi_t;
564 return framed_multi_t(lhs) & framed_multi_t(rhs);
568 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
572 operator&= (
const multivector_t& rhs) -> multivector_t&
573 {
return *
this = *
this & rhs; }
576 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
579 operator% (
const matrix_multi<Scalar_T,LO,HI,Tune_P>& lhs,
const matrix_multi<Scalar_T,LO,HI,Tune_P>& rhs) ->
const matrix_multi<Scalar_T,LO,HI,Tune_P>
582 using framed_multi_t =
typename multivector_t::framed_multi_t;
583 return framed_multi_t(lhs) % framed_multi_t(rhs);
587 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
591 operator%= (
const multivector_t& rhs) -> multivector_t&
592 {
return *
this = *
this % rhs; }
595 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
599 {
return (lhs * rhs).scalar(); }
602 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
606 operator/= (
const Scalar_T& scr) -> multivector_t&
607 {
return *
this *= Scalar_T(1)/scr; }
610 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
612 operator/ (
const matrix_multi<Scalar_T,LO,HI,Tune_P>& lhs,
const matrix_multi<Scalar_T,LO,HI,Tune_P>& rhs) ->
const matrix_multi<Scalar_T,LO,HI,Tune_P>
616 if (lhs.isnan() || rhs.isnan())
617 return traits_t::NaN();
619 if (rhs == Scalar_T(0))
620 return traits_t::NaN();
625 multivector_t lhs_reframed;
626 multivector_t rhs_reframed;
627 const auto our_frame =
reframe(lhs, rhs, lhs_reframed, rhs_reframed);
628 const auto& lhs_ref = (lhs.m_frame == our_frame)
631 const auto& rhs_ref = (rhs.m_frame == our_frame)
641 using matrix_t =
typename multivector_t::matrix_t;
642 using matrix_index_t =
typename matrix_t::size_type;
644 const auto& AT = matrix_t(ublas::trans(rhs_ref.m_matrix));
647 using permutation_t = ublas::permutation_matrix<matrix_index_t>;
649 auto pvector = permutation_t(AT.size1());
650 if (! ublas::lu_factorize(LU, pvector))
652 const auto& BT = matrix_t(ublas::trans(lhs_ref.m_matrix));
654 ublas::lu_substitute(LU, pvector, XT);
656 return traits_t::NaN();
661 if (Tune_P::div_max_steps > 0)
664 auto R = matrix_t(-BT);
665 ublas::axpy_prod(AT, XT, R,
false);
667 return traits_t::NaN();
669 auto nr = Scalar_T(ublas::norm_inf(R));
670 if ( nr != Scalar_T(0) && !traits_t::isNaN_or_isInf(nr) )
673 auto nrold = nr + Scalar_T(1);
676 step != Tune_P::div_max_steps &&
686 ublas::lu_substitute(LU, pvector, D);
690 ublas::axpy_prod(AT, XTnew, R,
false);
691 nr = ublas::norm_inf(R);
695 return multivector_t(ublas::trans(XT), our_frame);
699 return traits_t::NaN();
703 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
707 operator/= (
const multivector_t& rhs) -> multivector_t&
708 {
return *
this = *
this / rhs; }
711 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
714 operator| (
const matrix_multi<Scalar_T,LO,HI,Tune_P>& lhs,
const matrix_multi<Scalar_T,LO,HI,Tune_P>& rhs) ->
const matrix_multi<Scalar_T,LO,HI,Tune_P>
715 {
return rhs * lhs / rhs.involute(); }
718 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
722 operator|= (
const multivector_t& rhs) -> multivector_t&
723 {
return *
this = rhs * *
this / rhs.involute(); }
726 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
730 inv() const -> const multivector_t
731 {
return multivector_t(Scalar_T(1), this->m_frame) / *
this; }
734 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
738 pow(
int m)
const ->
const multivector_t
742 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
745 outer_pow(
int m)
const ->
const multivector_t
748 throw error_t(
"outer_pow(m): negative exponent");
749 framed_multi_t a = *
this;
750 return a.outer_pow(m);
754 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
759 {
return framed_multi_t(*this).grade(); }
762 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
766 frame() const -> const index_set_t
767 {
return this->m_frame; }
770 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
774 operator[] (
const index_set_t ist)
const -> Scalar_T
777 if ( (ist | this->m_frame) == this->m_frame)
784 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
790 if ((grade < 0) || (grade > HI-LO))
793 return (framed_multi_t(*
this))(grade);
797 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
801 scalar() const -> Scalar_T
803 const matrix_index_t dim = this->m_matrix.size1();
804 return matrix::trace(this->m_matrix) / Scalar_T(
double(dim) );
808 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
812 pure() const -> const multivector_t
813 {
return *
this - this->
scalar(); }
816 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
820 even() const -> const multivector_t
821 {
return framed_multi_t(*this).even(); }
824 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
828 odd() const -> const multivector_t
829 {
return framed_multi_t(*this).odd(); }
832 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
839 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
842 vector_part(
const index_set_t frm,
const bool prechecked)
const ->
const vector_t
844 if (!prechecked && (this->frame() | frm) != frm)
845 throw error_t(
"vector_part(frm): value is outside of requested frame");
848 if (this->frame() != frm)
849 return framed_multi_t(*this).vector_part(frm,
true);
851 const auto begin_index = frm.min();
852 const auto end_index = frm.max()+1;
867 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
871 involute() const -> const multivector_t
872 {
return framed_multi_t(*this).involute(); }
875 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
879 reverse() const -> const multivector_t
880 {
return framed_multi_t(*this).reverse(); }
883 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
887 conj() const -> const multivector_t
888 {
return framed_multi_t(*this).conj(); }
891 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
895 quad() const -> Scalar_T
898 return framed_multi_t(*this).quad();
902 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
906 norm() const -> Scalar_T
908 const matrix_index_t dim = this->m_matrix.size1();
913 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
918 {
return framed_multi_t(*this).max_abs(); }
921 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
930 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
934 write(
const std::string& msg)
const
935 { framed_multi_t(*this).write(msg); }
938 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
942 write(std::ofstream& ofile,
const std::string& msg)
const
945 throw error_t(
"write(ofile,msg): cannot write to output file");
946 framed_multi_t(*this).write(ofile, msg);
950 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
955 os << typename matrix_multi<Scalar_T,LO,HI,Tune_P>::framed_multi_t(val);
960 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
975 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
979 isinf() const ->
bool
981 if (std::numeric_limits<Scalar_T>::has_infinity)
988 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
992 isnan() const ->
bool
994 if (std::numeric_limits<Scalar_T>::has_quiet_NaN)
1001 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1005 truncated(
const Scalar_T& limit)
const ->
const multivector_t
1006 {
return framed_multi_t(*this).truncated(limit); }
1009 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1013 operator+= (
const term_t& term) -> multivector_t&
1015 if (term.second != Scalar_T(0))
1016 this->m_matrix.plus_assign(matrix_t(this->basis_element(term.first)) * term.second);
1021 template<
typename Multivector_T,
typename Matrix_T,
typename Basis_Matrix_T >
1026 using framed_multi_t = Multivector_T;
1028 using index_set_t =
typename framed_multi_t::index_set_t;
1029 using Scalar_T =
typename framed_multi_t::scalar_t;
1030 using matrix_t = Matrix_T;
1031 using basis_matrix_t = Basis_Matrix_T;
1032 using basis_scalar_t =
typename basis_matrix_t::value_type;
1036 return framed_multi_t(traits_t::to_scalar_t(X(0,0)));
1038 if (ublas::norm_inf(X) == 0)
1042 basis_matrix_t J(2,2,2);
1044 J(0,1) = basis_scalar_t(-1);
1045 J(1,0) = basis_scalar_t( 1);
1046 basis_matrix_t K = J;
1047 K(0,1) = basis_scalar_t( 1);
1048 basis_matrix_t JK = I;
1049 JK(0,0) = basis_scalar_t(-1);
1052 const index_set_t ist_mn = index_set_t(-level);
1053 const index_set_t ist_pn = index_set_t(level);
1054 const index_set_t ist_mnpn = ist_mn | ist_pn;
1057 using term_t =
typename framed_multi_t::term_t;
1058 const Scalar_T i_x = traits_t::to_scalar_t(signed_perm_nork(I, X)(0, 0));
1059 const Scalar_T j_x = traits_t::to_scalar_t(signed_perm_nork(J, X)(0, 0));
1060 const Scalar_T k_x = traits_t::to_scalar_t(signed_perm_nork(K, X)(0, 0));
1061 const Scalar_T jk_x = traits_t::to_scalar_t(signed_perm_nork(JK,X)(0, 0));
1064 result += term_t(ist_mn, j_x);
1065 result += term_t(ist_pn, k_x);
1066 return result += term_t(ist_mnpn, jk_x);
1070 const framed_multi_t& mn = framed_multi_t(ist_mn);
1071 const framed_multi_t& pn = framed_multi_t(ist_pn);
1072 const framed_multi_t& mnpn = framed_multi_t(ist_mnpn);
1074 (signed_perm_nork(I, X), level-1);
1076 (signed_perm_nork(J, X), level-1);
1078 (signed_perm_nork(K, X), level-1);
1080 (signed_perm_nork(JK,X), level-1);
1082 result = i_x.even() - jk_x.odd();
1083 result += (j_x.even() - k_x.odd()) * mn;
1084 result += (k_x.even() - j_x.odd()) * pn;
1085 return result += (jk_x.even() - i_x.odd()) * mnpn;
1090 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1103 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1104 template <
typename Other_Scalar_T>
1129 const index_t level = (p+q)/2;
1136 switch (
pos_mod(orig_p-orig_q, 8))
1146 if (orig_p-orig_q > 4)
1149 if (orig_p-orig_q < -3)
1158 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Matrix_T >
1160 public std::map< std::pair< const index_set<LO,HI>, const index_set<LO,HI> >,
1181 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1186 using index_set_pair_t = std::pair<const index_set_t, const index_set_t>;
1187 const auto& unfolded_pair = index_set_pair_t(ist, this->
m_frame);
1190 auto& basis_cache = basis_table_t::basis();
1192 const auto frame_count = this->
m_frame.count();
1193 const auto use_cache = frame_count <=
index_t(Tune_P::basis_max_count);
1197 const auto basis_it = basis_cache.find(unfolded_pair);
1198 if (basis_it != basis_cache.end())
1199 return *(basis_it->second);
1201 const auto folded_set = ist.fold(this->
m_frame);
1202 const auto folded_frame = this->
m_frame.fold();
1203 const auto& folded_pair = index_set_pair_t(folded_set, folded_frame);
1204 using basis_pair_t = std::pair<const index_set_pair_t, basis_matrix_t *>;
1207 const auto basis_it = basis_cache.find(folded_pair);
1208 if (basis_it != basis_cache.end())
1210 auto* result_ptr = basis_it->second;
1211 basis_cache.insert(basis_pair_t(unfolded_pair, result_ptr));
1215 const auto folded_max = folded_frame.max();
1216 const auto folded_min = folded_frame.min();
1217 const auto p = std::max(folded_max,
index_t(0));
1231 basis_cache.insert(basis_pair_t(folded_pair, result_ptr));
1232 basis_cache.insert(basis_pair_t(unfolded_pair, result_ptr));
1238 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P, const
size_t Size >
1243 const std::array<Scalar_T, Size>& numer,
1244 const std::array<Scalar_T, Size>& denom,
1255 return traits_t::NaN();
1258 const auto nbr_even_powers = Size/2 - 1;
1261 auto XX = std::vector<multivector_t>(nbr_even_powers);
1263 XX[1] = XX[0] * XX[0];
1266 k != nbr_even_powers;
1268 XX[k] = XX[k-2] * XX[1];
1271 auto N = multivector_t(numer[1]);
1274 k != nbr_even_powers;
1276 N += XX[k] * numer[2*k + 3];
1281 k != nbr_even_powers;
1283 N += XX[k] * numer[2*k + 2];
1284 auto D = multivector_t(denom[1]);
1287 k != nbr_even_powers;
1289 D += XX[k] * denom[2*k + 3];
1294 k != nbr_even_powers;
1296 D += XX[k] * denom[2*k + 2];
1301 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1308 const auto& invM =
inv(M);
1309 M = ((M + invM)/Scalar_T(2) + Scalar_T(1)) / Scalar_T(2);
1310 Y *= (invM + Scalar_T(1)) / Scalar_T(2);
1314 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1321 if (val == Scalar_T(0))
1324 static const auto sqrt_max_steps = Tune_P::db_sqrt_max_steps;
1330 step != sqrt_max_steps &&
norm(M - Scalar_T(1)) > norm_tol;
1341 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1348 if (val == Scalar_T(0))
1351 static const auto sqrt_max_steps = Tune_P::cr_sqrt_max_steps;
1352 auto Z = Scalar_T(2) * (Scalar_T(1) + val);
1353 auto Y = Scalar_T(1) - val;
1354 auto norm_Y =
norm(Y);
1357 step != sqrt_max_steps && norm_Y > norm_Y_tol;
1360 const auto old_norm_Y = norm_Y;
1363 if (Y.isnan() || (norm_Y > old_norm_Y * Scalar_T(2)))
1366 Z += Y * Scalar_T(2);
1368 return Z / Scalar_T(4);
1375 template<
typename Scalar_T >
1381 template<
typename Scalar_T >
1384 1.0, 27.0/4.0, 81.0/4.0, 2277.0/64.0,
1385 10395.0/256.0, 32319.0/1024.0, 8721.0/512.0, 26163.0/4096.0,
1386 53703.0/32768.0, 36465.0/131072.0, 3861.0/131072.0, 7371.0/4194304.0,
1387 819.0/16777216.0, 27.0/67108864.0
1392 template<
typename Scalar_T >
1398 template<
typename Scalar_T >
1401 1.0, 25.0/4.0, 69.0/4.0, 1771.0/64.0,
1402 7315.0/256.0, 20349.0/1024.0, 4845.0/512.0, 12597.0/4096.0,
1403 21879.0/32768.0, 12155.0/131072.0, 1001.0/131072.0, 1365.0/4194304.0,
1404 91.0/16777216.0, 1.0/67108864.0
1415 1.0, 19.0/4.0, 19.0/2.0, 665.0/64.0,
1416 1729.0/256.0, 2717.0/1024.0, 627.0/1024.0, 627.0/8192.0,
1417 285.0/65536.0, 19.0/262144.0
1427 1.0, 17.0/4.0, 15.0/2.0, 455.0/64.0,
1428 1001.0/256.0, 1287.0/1024.0, 231.0/1024.0, 165.0/8192.0,
1429 45.0/65536, 1.0/262144.0
1435 using array = std::array<long double, 18>;
1440 1.0L, 35.0L/4.0L, 35.0L, 5425.0L/64.0L,
1441 35525.0L/256.0L, 166257.0L/1024.0L, 143325.0L/1024.0L, 740025.0L/8192.0L,
1442 2877875.0L/65536.0L, 4206125.0L/262144.0L, 572033.0L/131072.0L, 1820105.0L/2097152.0L,
1443 1028755.0L/8388608.0L, 395675.0L/33554432.0L, 24225.0L/33554432.0L, 6783.0L/268435456.0L,
1444 1785.0L/4294967296.0L, 35.0L/17179869184.0L
1449 using array = std::array<long double, 18>;
1454 1.0L, 33.0L/4.0L, 31.0L, 4495.0L/64.0L,
1455 27405.0L/256.0L, 118755.0L/1024.0L, 94185.0L/1024.0L, 444015.0L/8192.0L,
1456 1562275.0L/65536.0L, 2042975.0L/262144.0L, 245157.0L/131072.0L, 676039.0L/2097152.0L,
1457 323323.0L/8388608.0L, 101745.0L/33554432.0L, 4845.0L/33554432.0L, 969.0L/268435456.0L,
1458 153.0L/4294967296.0L, 1.0L/17179869184.0L
1461#if defined(_GLUCAT_USE_QD)
1470 dd_real(
"1"), dd_real(
"43")/dd_real(
"4"),
1471 dd_real(
"215")/dd_real(
"4"), dd_real(
"10621")/dd_real(
"64"),
1472 dd_real(
"90687")/dd_real(
"256"), dd_real(
"567987")/dd_real(
"1024"),
1473 dd_real(
"168861")/dd_real(
"256"), dd_real(
"1246355")/dd_real(
"2048"),
1474 dd_real(
"7228859")/dd_real(
"16384"), dd_real(
"16583853")/dd_real(
"65536"),
1475 dd_real(
"7538115")/dd_real(
"65536"), dd_real(
"173376645")/dd_real(
"4194304"),
1476 dd_real(
"195747825")/dd_real(
"16777216"), dd_real(
"171655785")/dd_real(
"67108864"),
1477 dd_real(
"14375115")/dd_real(
"33554432"), dd_real(
"14375115")/dd_real(
"268435456"),
1478 dd_real(
"20764055")/dd_real(
"4294967296"), dd_real(
"5167525")/dd_real(
"17179869184"),
1479 dd_real(
"206701")/dd_real(
"17179869184"), dd_real(
"76153")/dd_real(
"274877906944"),
1480 dd_real(
"3311")/dd_real(
"1099511627776") , dd_real(
"43")/dd_real(
"4398046511104")
1490 dd_real(
"1"), dd_real(
"41")/dd_real(
"4"),
1491 dd_real(
"195")/dd_real(
"4"), dd_real(
"9139")/dd_real(
"64"),
1492 dd_real(
"73815")/dd_real(
"256"), dd_real(
"435897")/dd_real(
"1024"),
1493 dd_real(
"121737")/dd_real(
"256"), dd_real(
"840565")/dd_real(
"2048"),
1494 dd_real(
"4539051")/dd_real(
"16384"), dd_real(
"9641775")/dd_real(
"65536"),
1495 dd_real(
"4032015")/dd_real(
"65536"), dd_real(
"84672315")/dd_real(
"4194304"),
1496 dd_real(
"86493225")/dd_real(
"16777216"), dd_real(
"67863915")/dd_real(
"67108864"),
1497 dd_real(
"5014575")/dd_real(
"33554432"), dd_real(
"4345965")/dd_real(
"268435456"),
1498 dd_real(
"5311735")/dd_real(
"4294967296"), dd_real(
"1081575")/dd_real(
"17179869184"),
1499 dd_real(
"33649")/dd_real(
"17179869184"), dd_real(
"8855")/dd_real(
"274877906944"),
1500 dd_real(
"231")/dd_real(
"1099511627776"), dd_real(
"1")/dd_real(
"4398046511104")
1511 qd_real(
"1"), qd_real(
"67")/qd_real(
"4"),
1512 qd_real(
"134"), qd_real(
"43617")/qd_real(
"64"),
1513 qd_real(
"633485")/qd_real(
"256"), qd_real(
"6992857")/qd_real(
"1024"),
1514 qd_real(
"15246721")/qd_real(
"1024"), qd_real(
"215632197")/qd_real(
"8192"),
1515 qd_real(
"2518145487")/qd_real(
"65536"), qd_real(
"12301285425")/qd_real(
"262144"),
1516 qd_real(
"6344873535")/qd_real(
"131072"), qd_real(
"89075432355")/qd_real(
"2097152"),
1517 qd_real(
"267226297065")/qd_real(
"8388608"), qd_real(
"687479618945")/qd_real(
"33554432"),
1518 qd_real(
"379874182975")/qd_real(
"33554432"), qd_real(
"1443521895305")/qd_real(
"268435456"),
1519 qd_real(
"9425348845815")/qd_real(
"4294967296"), qd_real(
"13195488384141")/qd_real(
"17179869184"),
1520 qd_real(
"987417498133")/qd_real(
"4294967296"), qd_real(
"8055248011085")/qd_real(
"137438953472"),
1521 qd_real(
"6958363175533")/qd_real(
"549755813888"), qd_real(
"5056698705201")/qd_real(
"2199023255552"),
1522 qd_real(
"766166470485")/qd_real(
"2199023255552"), qd_real(
"766166470485")/qd_real(
"17592186044416"),
1523 qd_real(
"623623871325")/qd_real(
"140737488355328"), qd_real(
"203123203803")/qd_real(
"562949953421312"),
1524 qd_real(
"6478601247")/qd_real(
"281474976710656"), qd_real(
"5038912081")/qd_real(
"4503599627370496"),
1525 qd_real(
"719844583")/qd_real(
"18014398509481984"), qd_real(
"71853815")/qd_real(
"72057594037927936"),
1526 qd_real(
"1165197")/qd_real(
"72057594037927936"), qd_real(
"87703")/qd_real(
"576460752303423488"),
1527 qd_real(
"12529")/qd_real(
"18446744073709551616"), qd_real(
"67")/qd_real(
"73786976294838206464")
1537 qd_real(
"1"), qd_real(
"65")/qd_real(
"4"),
1538 qd_real(
"126"), qd_real(
"39711")/qd_real(
"64"),
1539 qd_real(
"557845")/qd_real(
"256"), qd_real(
"5949147")/qd_real(
"1024"),
1540 qd_real(
"12515965")/qd_real(
"1024"), qd_real(
"170574723")/qd_real(
"8192"),
1541 qd_real(
"1916797311")/qd_real(
"65536"), qd_real(
"8996462475")/qd_real(
"262144"),
1542 qd_real(
"4450881435")/qd_real(
"131072"), qd_real(
"59826782925")/qd_real(
"2097152"),
1543 qd_real(
"171503444385")/qd_real(
"8388608"), qd_real(
"420696483235")/qd_real(
"33554432"),
1544 qd_real(
"221120793075")/qd_real(
"33554432"), qd_real(
"797168807855")/qd_real(
"268435456"),
1545qd_real(
"4923689695575")/qd_real(
"4294967296"), qd_real(
"6499270398159")/qd_real(
"17179869184"),
1546 qd_real(
"456864812569")/qd_real(
"4294967296"), qd_real(
"3486599885395")/qd_real(
"137438953472"),
1547qd_real(
"2804116503573")/qd_real(
"549755813888"), qd_real(
"1886827875075")/qd_real(
"2199023255552"),
1548 qd_real(
"263012370465")/qd_real(
"2199023255552"), qd_real(
"240141729555")/qd_real(
"17592186044416"),
1549 qd_real(
"176848560525")/qd_real(
"140737488355328"), qd_real(
"51538723353")/qd_real(
"562949953421312"),
1550 qd_real(
"1450433115")/qd_real(
"281474976710656"), qd_real(
"977699359")/qd_real(
"4503599627370496"),
1551 qd_real(
"118183439")/qd_real(
"18014398509481984"), qd_real(
"9652005")/qd_real(
"72057594037927936"),
1552 qd_real(
"121737")/qd_real(
"72057594037927936"), qd_real(
"6545")/qd_real(
"576460752303423488"),
1553 qd_real(
"561")/qd_real(
"18446744073709551616"), qd_real(
"1")/qd_real(
"73786976294838206464")
1561 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1574 return traits_t::NaN();
1576 const auto scr_val = val.
scalar();
1579 if (scr_val < Scalar_T(0))
1580 return i * traits_t::sqrt(-scr_val);
1582 return traits_t::sqrt(scr_val);
1587 (scr_val != Scalar_T(0) &&
norm(val/scr_val - Scalar_T(1)) < Scalar_T(1))
1589 : (scr_val < Scalar_T(0))
1592 const auto sqrt_scale = traits_t::sqrt(traits_t::abs(scale));
1593 if (traits_t::isNaN_or_isInf(sqrt_scale))
1594 return traits_t::NaN();
1597 auto rescale = multivector_t(sqrt_scale);
1598 if (scale < Scalar_T(0))
1599 rescale = i * sqrt_scale;
1601 const auto& unitval = val / scale;
1602 static const auto max_norm = Scalar_T(1.0/4.0);
1603 auto use_approx_sqrt =
true;
1604 auto use_cr_sqrt =
false;
1605 auto scaled_result = multivector_t();
1606#if defined(_GLUCAT_USE_EIGENVALUES)
1607 static const auto sqrt_2 = traits_t::sqrt(Scalar_T(2));
1613 (genus.m_is_singular)
1616 switch (genus.m_eig_case)
1618 case matrix::neg_real_eigs:
1619 scaled_result =
matrix_sqrt(-i * unitval, i, next_level) * (i + Scalar_T(1)) / sqrt_2;
1620 use_approx_sqrt =
false;
1622 case matrix::both_eigs:
1624 const auto safe_arg = genus.m_safe_arg;
1625 scaled_result =
matrix_sqrt(
exp(i*safe_arg) * unitval, i, next_level) *
exp(-i*safe_arg / Scalar_T(2));
1627 use_approx_sqrt =
false;
1632 use_cr_sqrt = genus.m_is_singular;
1635 if (use_approx_sqrt)
1638 (
norm(unitval - Scalar_T(1)) < max_norm)
1642 unitval - Scalar_T(1))
1648 return (scaled_result.isnan() ||
1651 : scaled_result * rescale;
1655 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1657 sqrt(
const matrix_multi<Scalar_T,LO,HI,Tune_P>& val,
const matrix_multi<Scalar_T,LO,HI,Tune_P>& i,
bool prechecked) ->
const matrix_multi<Scalar_T,LO,HI,Tune_P>
1666 return traits_t::NaN();
1670 switch (Tune_P::function_precision)
1672 case precision_demoted:
1674 using demoted_scalar_t =
typename traits_t::demoted::type;
1677 const auto& demoted_val = demoted_multivector_t(val);
1678 const auto& demoted_i = demoted_multivector_t(i);
1683 case precision_promoted:
1685 using promoted_scalar_t =
typename traits_t::promoted::type;
1688 const auto& promoted_val = promoted_multivector_t(val);
1689 const auto& promoted_i = promoted_multivector_t(i);
1703 template<
typename Scalar_T >
1709 template<
typename Scalar_T >
1712 0.0, 1.0, 6.0, 4741.0/300.0,
1713 1441.0/60.0, 107091.0/4600.0, 8638.0/575.0, 263111.0/40250.0,
1714 153081.0/80500.0, 395243.0/1101240.0, 28549.0/688275.0, 605453.0/228813200.0,
1715 785633.0/10296594000.0, 1145993.0/1873980108000.0
1720 template<
typename Scalar_T >
1726 template<
typename Scalar_T >
1729 1.0, 13.0/2.0, 468.0/25.0, 1573.0/50.0,
1730 1573.0/46.0, 11583.0/460.0, 10296.0/805.0, 2574.0/575.0,
1731 11583.0/10925.0, 143.0/874.0, 572.0/37145.0, 117.0/148580.0,
1732 13.0/742900.0, 1.0/10400600.0
1743 0.0, 1.0, 4.0, 1337.0/204.0,
1744 385.0/68.0, 1879.0/680.0, 193.0/255.0, 197.0/1820.0,
1745 419.0/61880.0, 7129.0/61261200.0
1755 1.0, 9.0/2.0, 144.0/17.0, 147.0/17.0,
1756 441.0/85.0, 63.0/34.0, 84.0/221.0, 9.0/221.0,
1757 9.0/4862.0, 1.0/48620.0
1763 using array = std::array<long double, 18>;
1768 0.0L, 1.0L, 8.0L, 3835.0L/132.0L,
1769 8365.0L/132.0L, 11363807.0L/122760.0L, 162981.0L/1705.0L, 9036157.0L/125860.0L,
1770 18009875.0L/453096.0L, 44211925.0L/2718576.0L, 4149566.0L/849555.0L, 16973929.0L/16020180.0L,
1771 172459.0L/1068012.0L, 116317061.0L/7025382936.0L, 19679783.0L/18441630207.0L, 23763863.0L/614721006900.0L,
1772 50747.0L/79318839600.0L, 42142223.0L/14295951736466400.0L
1777 using array = std::array<long double, 18>;
1782 1.0L, 17.0L/2.0L, 1088.0L/33.0L, 850.0L/11.0L,
1783 41650.0L/341.0L, 140777.0L/1023.0L, 1126216.0L/9889.0L, 63206.0L/899.0L,
1784 790075.0L/24273.0L, 60775.0L/5394.0L, 38896.0L/13485.0L, 21658.0L/40455.0L,
1785 21658.0L/310155.0L, 4165.0L/682341.0L, 680.0L/2047023.0L, 34.0L/3411705.0L,
1786 17.0L/129644790.0L, 1.0L/2333606220
1788#if defined(_GLUCAT_USE_QD)
1797 dd_real(
"0"), dd_real(
"1"),
1798 dd_real(
"10"), dd_real(
"22781")/dd_real(
"492"),
1799 dd_real(
"21603")/dd_real(
"164"), dd_real(
"5492649")/dd_real(
"21320"),
1800 dd_real(
"978724")/dd_real(
"2665"), dd_real(
"4191605")/dd_real(
"10619"),
1801 dd_real(
"12874933")/dd_real(
"39442"), dd_real(
"11473457")/dd_real(
"54612"),
1802 dd_real(
"2406734")/dd_real(
"22755"), dd_real(
"166770367")/dd_real(
"4004880"),
1803 dd_real(
"30653165")/dd_real(
"2402928"), dd_real(
"647746389")/dd_real(
"215195552"),
1804 dd_real(
"25346331")/dd_real(
"47074027"), dd_real(
"278270613")/dd_real(
"3900419380"),
1805 dd_real(
"105689791")/dd_real(
"15601677520"), dd_real(
"606046475")/dd_real(
"1379188292768"),
1806 dd_real(
"969715")/dd_real(
"53502994116"), dd_real(
"11098301")/dd_real(
"26204577562592"),
1807 dd_real(
"118999")/dd_real(
"26204577562592"), dd_real(
"18858053")/dd_real(
"1392249205900512960")
1817 dd_real(
"1"), dd_real(
"21")/dd_real(
"2"),
1818 dd_real(
"2100")/dd_real(
"41"), dd_real(
"12635")/dd_real(
"82"),
1819 dd_real(
"341145")/dd_real(
"1066"), dd_real(
"1037799")/dd_real(
"2132"),
1820 dd_real(
"11069856")/dd_real(
"19721"), dd_real(
"9883800")/dd_real(
"19721"),
1821 dd_real(
"6918660")/dd_real(
"19721"), dd_real(
"293930")/dd_real(
"1517"),
1822 dd_real(
"1410864")/dd_real(
"16687"), dd_real(
"88179")/dd_real(
"3034"),
1823 dd_real(
"734825")/dd_real(
"94054"), dd_real(
"305235")/dd_real(
"188108"),
1824 dd_real(
"348840")/dd_real(
"1363783"), dd_real(
"40698")/dd_real(
"1363783"),
1825 dd_real(
"6783")/dd_real(
"2727566"), dd_real(
"9975")/dd_real(
"70916716"),
1826 dd_real(
"266")/dd_real(
"53187537"), dd_real(
"7")/dd_real(
"70916716"),
1827 dd_real(
"7")/dd_real(
"8155422340"), dd_real(
"1")/dd_real(
"538257874440")
1838 qd_real(
"0"), qd_real(
"1"),
1839 qd_real(
"16"), qd_real(
"95201")/qd_real(
"780"),
1840 qd_real(
"30721")/qd_real(
"52"), qd_real(
"7416257")/qd_real(
"3640"),
1841 qd_real(
"1039099")/qd_real(
"195"), qd_real(
"6097772319")/qd_real(
"555100"),
1842 qd_real(
"1564058073")/qd_real(
"85400"), qd_real(
"30404640205")/qd_real(
"1209264"),
1843 qd_real(
"725351278")/qd_real(
"25193"), qd_real(
"4092322670789")/qd_real(
"147429436"),
1844 qd_real(
"4559713849589")/qd_real(
"201040140"), qd_real(
"5049361751189")/qd_real(
"320023080"),
1845 qd_real(
"74979677195")/qd_real(
"8000577"), qd_real(
"16569850691873")/qd_real(
"3481514244"),
1846 qd_real(
"1065906022369")/qd_real(
"515779888"), qd_real(
"335956770855841")/qd_real(
"438412904800"),
1847 qd_real(
"1462444287585964")/qd_real(
"6041877844275"), qd_real(
"397242326339851")/qd_real(
"6122436215532"),
1848 qd_real(
"64211291334131")/qd_real(
"4373168725380"), qd_real(
"142322343550859")/qd_real(
"51080680851480"),
1849 qd_real(
"154355972958659")/qd_real(
"351179680853925"), qd_real(
"167483568676259")/qd_real(
"2937139148960100"),
1850 qd_real(
"4230788929433")/qd_real(
"704913395750424"), qd_real(
"197968763176019")/qd_real(
"392923948371995600"),
1851 qd_real(
"10537522306718")/qd_real(
"319250708052246425"), qd_real(
"236648286272519")/qd_real(
"144249197475035425500"),
1852 qd_real(
"260715545088119")/qd_real(
"4375558990076074573500"), qd_real(
"289596255666839")/qd_real(
"192874640282553367199880"),
1853 qd_real(
"8802625510547")/qd_real(
"361639950529787563499775"), qd_real(
"373831661521439")/qd_real(
"1659204093030665341336967700"),
1854 qd_real(
"446033437968239")/qd_real(
"464577146048586295574350956000"), qd_real(
"53676090078349")/qd_real(
"47386868896955802148583797512000")
1864 qd_real(
"1"), qd_real(
"33")/qd_real(
"2"),
1865 qd_real(
"8448")/qd_real(
"65"), qd_real(
"42284")/qd_real(
"65"),
1866 qd_real(
"211420")/qd_real(
"91"), qd_real(
"573562")/qd_real(
"91"),
1867 qd_real(
"32119472")/qd_real(
"2379"), qd_real(
"92917044")/qd_real(
"3965"),
1868 qd_real(
"603960786")/qd_real(
"17995"), qd_real(
"144626625")/qd_real(
"3599"),
1869 qd_real(
"2776831200")/qd_real(
"68381"), qd_real(
"16692542100")/qd_real(
"478667"),
1870 qd_real(
"12241197540")/qd_real(
"478667"), qd_real(
"1098569010")/qd_real(
"68381"),
1871 qd_real(
"31387686000")/qd_real(
"3624193"), qd_real(
"9939433900")/qd_real(
"2479711"),
1872 qd_real(
"67091178825")/qd_real(
"42155087"), qd_real(
"2683647153")/qd_real(
"4959422"),
1873 qd_real(
"19083713088")/qd_real(
"121505839"), qd_real(
"4708152900")/qd_real(
"121505839"),
1874 qd_real(
"941630580")/qd_real(
"116546417"), qd_real(
"88704330")/qd_real(
"62755763"),
1875 qd_real(
"12902448")/qd_real(
"62755763"), qd_real(
"1542684")/qd_real(
"62755763"),
1876 qd_real(
"6427850")/qd_real(
"2698497809"), qd_real(
"3471039")/qd_real(
"18889484663"),
1877 qd_real(
"8544096")/qd_real(
"774468871183"), qd_real(
"39556")/qd_real(
"79027435835"),
1878 qd_real(
"118668")/qd_real(
"7191496660985"), qd_real(
"10230")/qd_real(
"27327687311743"),
1879 qd_real(
"5456")/qd_real(
"1011124430534491"), qd_real(
"44")/qd_real(
"1011124430534491"),
1880 qd_real(
"11")/qd_real(
"70778710137414370"), qd_real(
"1")/qd_real(
"7219428434016265740")
1887 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1898 if (val == Scalar_T(0) || val.
isnan())
1899 return traits_t::NaN();
1907 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1915 if (val == Scalar_T(0) || val.
isnan())
1916 return traits_t::NaN();
1918 using limits_t = std::numeric_limits<Scalar_T>;
1919 static const auto epsilon = limits_t::epsilon();
1920 static const auto max_inner_norm = traits_t::pow(
epsilon, 2);
1921 static const auto max_outer_norm = Scalar_T(6.0/limits_t::digits);
1923 auto E = multivector_t(Scalar_T(0));
1925 auto pow_2_outer_step = Scalar_T(1);
1926 auto pow_4_outer_step = Scalar_T(1);
1928 for (outer_step = 0, norm_Y_1 =
norm(Y - Scalar_T(1));
1929 outer_step != Tune_P::log_max_outer_steps && norm_Y_1 * pow_2_outer_step > max_outer_norm;
1930 ++outer_step, norm_Y_1 =
norm(Y - Scalar_T(1)))
1932 if (Y == Scalar_T(0) || Y.isnan())
1933 return traits_t::NaN();
1939 inner_step != Tune_P::log_max_inner_steps &&
1940 norm(M - Scalar_T(1)) * pow_4_outer_step > max_inner_norm;
1944 E += (M - Scalar_T(1)) * pow_2_outer_step;
1945 pow_2_outer_step *= Scalar_T(2);
1946 pow_4_outer_step *= Scalar_T(4);
1948 if (outer_step == Tune_P::log_max_outer_steps && norm_Y_1 * pow_2_outer_step > max_outer_norm)
1949 return traits_t::NaN();
1951 return pade_log(Y) * pow_2_outer_step - E;
1955 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1965 if (val == Scalar_T(0) || val.
isnan())
1966 return traits_t::NaN();
1968 static const auto pi = traits_t::pi();
1969 const auto scr_val = val.
scalar();
1972 if (scr_val < Scalar_T(0))
1973 return i * pi + traits_t::log(-scr_val);
1975 return traits_t::log(scr_val);
1979 const auto max_norm = Scalar_T(1.0/9.0);
1981 (scr_val != Scalar_T(0) &&
norm(val/scr_val - Scalar_T(1)) < max_norm)
1983 : (scr_val < Scalar_T(0))
1986 if (scale == Scalar_T(0))
1987 return traits_t::NaN();
1990 const auto log_scale = traits_t::log(traits_t::abs(scale));
1991 auto rescale = multivector_t(log_scale);
1992 if (scale < Scalar_T(0))
1993 rescale = i * pi + log_scale;
1994 const auto unitval = val/scale;
1995 if (
inv(unitval).isnan())
1996 return traits_t::NaN();
1998#if defined(_GLUCAT_USE_EIGENVALUES)
1999 auto scaled_result = multivector_t();
2004 switch (genus.m_eig_case)
2006 case matrix::neg_real_eigs:
2007 scaled_result =
matrix_log(-i * unitval, i, level + 1) + i * pi/Scalar_T(2);
2009 case matrix::both_eigs:
2011 const Scalar_T safe_arg = genus.m_safe_arg;
2012 scaled_result =
matrix_log(
exp(i*safe_arg) * unitval, i, level + 1) - i * safe_arg;
2025 return (scaled_result.isnan())
2027 : scaled_result + rescale;
2031 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
2033 log(
const matrix_multi<Scalar_T,LO,HI,Tune_P>& val,
const matrix_multi<Scalar_T,LO,HI,Tune_P>& i,
bool prechecked) ->
const matrix_multi<Scalar_T,LO,HI,Tune_P>
2037 if (val == Scalar_T(0) || val.
isnan())
2038 return traits_t::NaN();
2042 switch (Tune_P::function_precision)
2044 case precision_demoted:
2046 using demoted_scalar_t =
typename traits_t::demoted::type;
2049 const auto& demoted_val = demoted_multivector_t(val);
2050 const auto& demoted_i = demoted_multivector_t(i);
2052 return matrix_log(demoted_val, demoted_i, 0);
2055 case precision_promoted:
2057 using promoted_scalar_t =
typename traits_t::promoted::type;
2060 const auto& promoted_val = promoted_multivector_t(val);
2061 const auto& promoted_i = promoted_multivector_t(i);
2063 return matrix_log(promoted_val, promoted_i, 0);
2072 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
2078 return traits_t::NaN();
2080 const auto scr_val = val.
scalar();
2082 return traits_t::exp(scr_val);
2084 switch (Tune_P::function_precision)
2086 case precision_demoted:
2088 using demoted_scalar_t =
typename traits_t::demoted::type;
2091 const auto& demoted_val = demoted_multivector_t(val);
2095 case precision_promoted:
2097 using promoted_scalar_t =
typename traits_t::promoted::type;
2100 const auto& promoted_val = promoted_multivector_t(val);
Table of basis elements used as a cache by basis_element()
auto operator=(const basis_table &) -> basis_table &=delete
static auto basis() -> basis_table &
Single instance of basis table.
friend class friend_for_private_destructor
basis_table(const basis_table &)=delete
virtual void write(const std::string &msg="") const=0
virtual auto involute() const -> const multivector_t=0
virtual auto outer_pow(int m) const -> const multivector_t=0
virtual auto isinf() const -> bool=0
virtual auto quad() const -> double=0
virtual auto scalar() const -> double=0
virtual auto operator()(index_t grade) const -> const multivector_t=0
virtual auto conj() const -> const multivector_t=0
virtual auto operator-=(const multivector_t &rhs) -> multivector_t &=0
virtual auto grade() const -> index_t=0
virtual auto even() const -> const multivector_t=0
virtual auto frame() const -> const index_set_t=0
virtual auto operator/=(const double &scr) -> multivector_t &=0
virtual auto max_abs() const -> double=0
virtual auto odd() const -> const multivector_t=0
virtual auto vector_part() const -> const vector_t=0
virtual auto reverse() const -> const multivector_t=0
virtual auto operator%=(const multivector_t &rhs) -> multivector_t &=0
virtual auto operator-() const -> const multivector_t=0
virtual auto truncated(const Scalar_T &limit=default_truncation) const -> const multivector_t=0
Remove all terms with relative size smaller than limit.
virtual auto operator|=(const multivector_t &rhs) -> multivector_t &=0
virtual auto pure() const -> const multivector_t=0
virtual auto operator==(const multivector_t &val) const -> bool=0
virtual auto operator^=(const multivector_t &rhs) -> multivector_t &=0
virtual auto isnan() const -> bool=0
virtual auto operator*=(const double &scr) -> multivector_t &=0
virtual auto norm() const -> double=0
virtual auto pow(int m) const -> const multivector_t=0
virtual auto operator&=(const multivector_t &rhs) -> multivector_t &=0
virtual auto operator[](const index_set_t ist) const -> double=0
virtual auto inv() const -> const multivector_t=0
A framed_multi<Scalar_T,LO,HI,Tune_P> is a framed approximation to a multivector.
auto centre_pm4_qp4(index_t &p, index_t &q) -> multivector_t &
Subalgebra isomorphism: R_{p,q} to R_{p-4,q+4}.
static auto random(const index_set_t frm, Scalar_T fill=Scalar_T(1)) -> const multivector_t
Random multivector within a frame.
auto centre_qp1_pm1(index_t &p, index_t &q) -> multivector_t &
Subalgebra isomorphism: R_{p,q} to R_{q+1,p-1}.
auto unfold(const index_set_t frm) const -> multivector_t
Subalgebra isomorphism: unfold each term within the given frame.
auto centre_pp4_qm4(index_t &p, index_t &q) -> multivector_t &
Subalgebra isomorphism: R_{p,q} to R_{p+4,q-4}.
static auto generator() -> generator_table< Matrix_T > &
Single instance of generator table.
Abstract exception class.
Index set class based on std::bitset<> in Gnu standard C++ library.
auto count() const -> index_t
Cardinality: Number of indices included in set.
auto min() const -> index_t
Minimum member.
auto max() const -> index_t
Maximum member.
A matrix_multi<Scalar_T,LO,HI,Tune_P> is a matrix approximation to a multivector.
_GLUCAT_CLIFFORD_ALGEBRA_OPERATIONS auto operator=(const multivector_t &rhs) -> multivector_t &
Assignment operator.
typename matrix_t::size_type matrix_index_t
friend class framed_multi
static auto classname() -> const std::string
Class name used in messages.
index_set< LO, HI > index_set_t
auto basis_element(const index_set< LO, HI > &ist) const -> const basis_matrix_t
Create a basis element matrix within the current frame.
std::pair< const index_set_t, Scalar_T > term_t
auto operator+=(const term_t &rhs) -> multivector_t &
Add a term, if non-zero.
matrix_t m_matrix
Matrix value representing the multivector within the folded frame.
auto fast_framed_multi() const -> const framed_multi< Other_Scalar_T, LO, HI, Tune_P >
Use inverse generalized FFT to construct a framed_multi_t.
auto fast_matrix_multi(const index_set_t frm) const -> const matrix_multi_t
Use generalized FFT to construct a matrix_multi_t.
ublas::compressed_matrix< int, orientation_t > basis_matrix_t
index_set_t m_frame
Index set representing the frame for the subalgebra which contains the multivector.
friend class matrix_multi
ublas::matrix< Scalar_T, orientation_t > matrix_t
matrix_multi multivector_t
error< multivector_t > error_t
std::vector< Scalar_T > vector_t
static auto random(const index_set_t frm, Scalar_T fill=Scalar_T(1)) -> const matrix_multi_t
Random multivector within a frame.
framed_multi< Scalar_T, LO, HI, Tune_P > framed_multi_t
Extra traits which extend numeric limits.
static auto NaN() -> Scalar_T
Smart NaN.
static auto to_scalar_t(const Other_Scalar_T &val) -> Scalar_T
Cast to Scalar_T.
static const std::array< index_t, 8 > offset_to_super
Offsets between the current signature and that of the real superalgebra.
auto inner(const LHS_T &lhs, const RHS_T &rhs) -> Scalar_T
Inner product: sum(x(i,j)*y(i,j))/x.nrows()
auto signed_perm_nork(const LHS_T &lhs, const RHS_T &rhs) -> const RHS_T
Left inverse of Kronecker product where lhs is a signed permutation matrix.
auto trace(const Matrix_T &val) -> typename Matrix_T::value_type
Matrix trace.
auto isinf(const Matrix_T &m) -> bool
Infinite.
auto norm_frob2(const Matrix_T &val) -> typename Matrix_T::value_type
Square of Frobenius norm.
auto classify_eigenvalues(const Matrix_T &val) -> eig_genus< Matrix_T >
Classify the eigenvalues of a matrix.
auto mono_prod(const ublas::matrix_expression< LHS_T > &lhs, const ublas::matrix_expression< RHS_T > &rhs) -> const typename RHS_T::expression_type
Product of monomial matrices.
auto unit(const typename Matrix_T::size_type n) -> const Matrix_T
Unit matrix - as per Matlab eye.
auto isnan(const Matrix_T &m) -> bool
Not a Number.
auto operator<<(std::ostream &os, const framed_multi< Scalar_T, LO, HI, Tune_P > &val) -> std::ostream &
Write multivector to output.
auto operator|(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const RHS< Scalar_T, LO, HI, Tune_P > &rhs) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Transformation via twisted adjoint action.
auto operator*(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const Scalar_T &scr) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Product of multivector and scalar.
auto operator&(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const RHS< Scalar_T, LO, HI, Tune_P > &rhs) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Inner product.
auto exp(const framed_multi< Scalar_T, LO, HI, Tune_P > &val) -> const framed_multi< Scalar_T, LO, HI, Tune_P >
Exponential of multivector.
auto scalar(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> Scalar_T
Scalar part.
auto approx_equal(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const RHS< Scalar_T, LO, HI, Tune_P > &rhs, const Scalar_T threshold, const Scalar_T tolerance) -> bool
Test for approximate equality of multivectors.
static auto pade_approx(const std::array< Scalar_T, Size > &numer, const std::array< Scalar_T, Size > &denom, const matrix_multi< Scalar_T, LO, HI, Tune_P > &X) -> const matrix_multi< Scalar_T, LO, HI, Tune_P >
Pade' approximation.
static auto fast(const Matrix_T &X, index_t level) -> Multivector_T
Inverse generalized Fast Fourier Transform.
static void check_complex(const Multivector< Scalar_T, LO, HI, Tune_P > &val, const Multivector< Scalar_T, LO, HI, Tune_P > &i, const bool prechecked=false)
Check that i is a valid complexifier for val.
static auto cr_sqrt(const matrix_multi< Scalar_T, LO, HI, Tune_P > &val, Scalar_T norm_Y_tol=std::pow(std::numeric_limits< Scalar_T >::epsilon(), 1)) -> const matrix_multi< Scalar_T, LO, HI, Tune_P >
Cyclic reduction square root iteration.
static auto db_sqrt(const matrix_multi< Scalar_T, LO, HI, Tune_P > &val, Scalar_T norm_tol=std::pow(std::numeric_limits< Scalar_T >::epsilon(), 4)) -> const matrix_multi< Scalar_T, LO, HI, Tune_P >
Product form of Denman-Beavers square root iteration.
auto operator%(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const RHS< Scalar_T, LO, HI, Tune_P > &rhs) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Left contraction.
static void db_step(matrix_multi< Scalar_T, LO, HI, Tune_P > &M, matrix_multi< Scalar_T, LO, HI, Tune_P > &Y)
Single step of product form of Denman-Beavers square root iteration.
auto reframe(const matrix_multi< Scalar_T, LO, HI, Tune_P > &lhs, const matrix_multi< Scalar_T, LO, HI, Tune_P > &rhs, matrix_multi< Scalar_T, LO, HI, Tune_P > &lhs_reframed, matrix_multi< Scalar_T, LO, HI, Tune_P > &rhs_reframed) -> const index_set< LO, HI >
Find a common frame for operands of a binary operator.
auto pos_mod(LHS_T lhs, RHS_T rhs) -> LHS_T
Modulo function which works reliably for lhs < 0.
auto pow(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, int rhs) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Integer power of multivector.
auto abs(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> Scalar_T
Absolute value == sqrt(norm)
auto inv(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Geometric multiplicative inverse.
static auto pade_log(const matrix_multi< Scalar_T, LO, HI, Tune_P > &val) -> const matrix_multi< Scalar_T, LO, HI, Tune_P >
Pade' approximation of log.
auto operator>>(std::istream &s, framed_multi< Scalar_T, LO, HI, Tune_P > &val) -> std::istream &
Read multivector from input.
int index_t
Size of index_t should be enough to represent LO, HI.
auto clifford_exp(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Exponential of multivector.
static auto cascade_log(const matrix_multi< Scalar_T, LO, HI, Tune_P > &val) -> const matrix_multi< Scalar_T, LO, HI, Tune_P >
Incomplete square root cascade and Pade' approximation of log.
auto operator/(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const Scalar_T &scr) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Quotient of multivector and scalar.
auto log(const Multivector< Scalar_T, LO, HI, Tune_P > &val, const Multivector< Scalar_T, LO, HI, Tune_P > &i, const bool prechecked=false) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Natural logarithm of multivector with specified complexifier.
auto norm(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> Scalar_T
Scalar_T norm == sum of norm of coordinates.
auto matrix_sqrt(const matrix_multi< Scalar_T, LO, HI, Tune_P > &val, const matrix_multi< Scalar_T, LO, HI, Tune_P > &i, const index_t level) -> const matrix_multi< Scalar_T, LO, HI, Tune_P >
Square root of multivector with specified complexifier.
auto star(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const RHS< Scalar_T, LO, HI, Tune_P > &rhs) -> Scalar_T
Hestenes scalar product.
auto operator^(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const RHS< Scalar_T, LO, HI, Tune_P > &rhs) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Outer product.
auto matrix_log(const matrix_multi< Scalar_T, LO, HI, Tune_P > &val, const matrix_multi< Scalar_T, LO, HI, Tune_P > &i, const index_t level) -> const matrix_multi< Scalar_T, LO, HI, Tune_P >
Natural logarithm of multivector with specified complexifier.
auto sqrt(const Multivector< Scalar_T, LO, HI, Tune_P > &val, const Multivector< Scalar_T, LO, HI, Tune_P > &i, const bool prechecked=false) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Square root of multivector with specified complexifier.
static auto folded_dim(const index_set< LO, HI > &sub) -> Matrix_Index_T
Determine the matrix dimension of the fold of a subalegbra.
auto vector_part(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const std::vector< Scalar_T >
Vector part of multivector, as a vector_t with respect to frame()
auto offset_level(const index_t p, const index_t q) -> index_t
Determine the log2 dim corresponding to signature p, q.
std::array< dd_real, 22 > array
std::array< float, 10 > array
std::array< long double, 18 > array
std::array< qd_real, 34 > array
Coefficients of denominator polynomials of Pade approximations produced by Pade1(log(1+x),...
std::array< Scalar_T, 14 > array
std::array< dd_real, 22 > array
std::array< float, 10 > array
std::array< long double, 18 > array
std::array< qd_real, 34 > array
Coefficients of numerator polynomials of Pade approximations produced by Pade1(log(1+x),...
std::array< Scalar_T, 14 > array
std::array< dd_real, 22 > array
std::array< float, 10 > array
std::array< long double, 18 > array
std::array< qd_real, 34 > array
Coefficients of denominator polynomials of Pade approximations produced by Pade1(sqrt(1+x),...
std::array< Scalar_T, 14 > array
std::array< dd_real, 22 > array
std::array< float, 10 > array
std::array< long double, 18 > array
std::array< qd_real, 34 > array
Coefficients of numerator polynomials of Pade approximations produced by Pade1(sqrt(1+x),...
std::array< Scalar_T, 14 > array