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 >
117 : m_frame( val.m_frame ), m_matrix( val.m_matrix.size1(), val.m_matrix.size2() )
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 >
175 const auto dim = folded_dim<matrix_index_t>(this->
m_frame);
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 >
201 const auto dim = folded_dim<matrix_index_t>(frm);
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");
222 const auto dim = folded_dim<matrix_index_t>(frm);
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 >
255 : m_frame( val.frame() )
257 if (val.size() >= Tune_P::fast_size_threshold)
260 *
this = val.template fast_matrix_multi<Scalar_T>(this->
m_frame);
265 const auto dim = folded_dim<matrix_index_t>(this->
m_frame);
266 this->
m_matrix.resize(dim, dim,
false);
270 for (
auto& val_term : val)
275 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
276 template<
typename Other_Scalar_T >
281 const auto our_frame = val.frame() | frm;
282 if (val.size() >= Tune_P::fast_size_threshold)
285 *
this = val.template fast_matrix_multi<Scalar_T>(our_frame);
290 this->m_frame = our_frame;
291 const auto dim = folded_dim<matrix_index_t>(our_frame);
292 this->m_matrix.resize(dim, dim,
false);
293 this->m_matrix.clear();
296 for (
auto& val_term : val)
301 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
302 template<
typename Matrix_T >
305 : m_frame( frm ), m_matrix( mtx.size1(), mtx.size2() )
310 mtx_it1 = mtx.begin1();
311 mtx_it1 != mtx.end1();
314 mtx_it2 = mtx_it1.begin();
315 mtx_it2 != mtx_it1.end();
321 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
324 : m_frame( frm ), m_matrix( mtx )
328 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
336 this->m_frame = rhs.m_frame;
337 this->m_matrix = rhs.m_matrix;
342 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
350 using framed_multi_t =
typename multivector_t::framed_multi_t;
352 index_set_t our_frame = lhs.m_frame | rhs.m_frame;
353 framed_multi_t framed_lhs;
354 framed_multi_t framed_rhs;
355 if ((lhs.m_frame != our_frame) || (rhs.m_frame != our_frame))
358 framed_lhs = framed_multi_t(lhs);
359 framed_rhs = framed_multi_t(rhs);
360 our_frame |= framed_lhs.frame() | framed_rhs.frame();
363 if (lhs.m_frame != our_frame)
364 lhs_reframed = multivector_t(framed_lhs, our_frame,
true);
365 if (rhs.m_frame != our_frame)
366 rhs_reframed = multivector_t(framed_rhs, our_frame,
true);
371 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
374 operator== (
const multivector_t& rhs)
const ->
bool
381 multivector_t lhs_reframed;
382 multivector_t rhs_reframed;
383 const index_set_t our_frame =
reframe(*
this, rhs, lhs_reframed, rhs_reframed);
384 const multivector_t& lhs_ref = (this->m_frame == our_frame)
387 const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
391 return ublas::norm_inf(lhs_ref.m_matrix - rhs_ref.m_matrix) == 0;
395 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
399 operator== (
const Scalar_T& scr)
const ->
bool
401 if (scr != Scalar_T(0))
402 return *
this == multivector_t(framed_multi_t(scr), this->m_frame,
true);
403 else if (ublas::norm_inf(this->m_matrix) != 0)
407 const matrix_index_t dim = this->m_matrix.size1();
408 return !(dim == 1 && this->
isnan());
413 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
421 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
425 operator+= (
const multivector_t& rhs) -> multivector_t&
429 return *
this *= Scalar_T(2);
432 multivector_t rhs_reframed;
433 const index_set_t our_frame =
reframe(*
this, rhs, *
this, rhs_reframed);
434 const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
438 noalias(this->m_matrix) += rhs_ref.m_matrix;
443 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
447 operator-= (
const Scalar_T& scr) -> multivector_t&
448 {
return *
this += term_t(index_set_t(), -scr); }
451 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
455 operator-= (
const multivector_t& rhs) -> multivector_t&
459 return *
this = Scalar_T(0);
462 multivector_t rhs_reframed;
463 const index_set_t our_frame =
reframe(*
this, rhs, *
this, rhs_reframed);
464 const multivector_t& rhs_ref = (
rhs.m_frame == our_frame)
468 noalias(this->m_matrix) -= rhs_ref.m_matrix;
473 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
477 operator- () const -> const multivector_t
478 {
return multivector_t(-(this->m_matrix), this->m_frame); }
481 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
485 operator*= (
const Scalar_T& scr) -> multivector_t&
488 using traits_t = numeric_traits<Scalar_T>;
489 if (traits_t::isNaN_or_isInf(scr) || this->
isnan())
490 return *
this = traits_t::NaN();
491 if (scr == Scalar_T(0))
494 this->m_matrix *= scr;
499 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
502 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>
505 using index_set_t =
typename multivector_t::index_set_t;
507 if (lhs.isnan() || rhs.isnan())
511 multivector_t lhs_reframed;
512 multivector_t rhs_reframed;
513 const index_set_t our_frame =
reframe(lhs, rhs, lhs_reframed, rhs_reframed);
514 const multivector_t& lhs_ref = (lhs.m_frame == our_frame)
517 const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
521 using matrix_t =
typename multivector_t::matrix_t;
522 using matrix_index_t =
typename matrix_t::size_type;
524 const matrix_index_t dim = lhs_ref.m_matrix.size1();
525 multivector_t result = multivector_t(matrix_t(dim, dim), our_frame);
526 result.m_matrix.clear();
527 ublas::axpy_prod(lhs_ref.m_matrix, rhs_ref.m_matrix, result.m_matrix,
true);
532 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
536 operator*= (
const multivector_t& rhs) -> multivector_t&
537 {
return *
this = *
this * rhs; }
540 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
543 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>
546 using framed_multi_t =
typename multivector_t::framed_multi_t;
547 return framed_multi_t(lhs) ^ framed_multi_t(rhs);
551 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
555 operator^= (
const multivector_t& rhs) -> multivector_t&
556 {
return *
this = *
this ^ rhs; }
559 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
562 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>
565 using framed_multi_t =
typename multivector_t::framed_multi_t;
566 return framed_multi_t(lhs) & framed_multi_t(rhs);
570 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
574 operator&= (
const multivector_t& rhs) -> multivector_t&
575 {
return *
this = *
this & rhs; }
578 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
581 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>
584 using framed_multi_t =
typename multivector_t::framed_multi_t;
585 return framed_multi_t(lhs) % framed_multi_t(rhs);
589 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
593 operator%= (
const multivector_t& rhs) -> multivector_t&
594 {
return *
this = *
this % rhs; }
597 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
601 {
return (lhs * rhs).scalar(); }
604 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
608 operator/= (
const Scalar_T& scr) -> multivector_t&
609 {
return *
this *= Scalar_T(1)/scr; }
612 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
614 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>
618 if (lhs.isnan() || rhs.isnan())
619 return traits_t::NaN();
621 if (rhs == Scalar_T(0))
622 return traits_t::NaN();
625 using index_set_t =
typename multivector_t::index_set_t;
628 multivector_t lhs_reframed;
629 multivector_t rhs_reframed;
630 const auto our_frame =
reframe(lhs, rhs, lhs_reframed, rhs_reframed);
631 const auto& lhs_ref = (lhs.m_frame == our_frame)
634 const auto& rhs_ref = (rhs.m_frame == our_frame)
644 using matrix_t =
typename multivector_t::matrix_t;
645 using matrix_index_t =
typename matrix_t::size_type;
647 const auto& AT = matrix_t(ublas::trans(rhs_ref.m_matrix));
650 using permutation_t = ublas::permutation_matrix<matrix_index_t>;
652 auto pvector = permutation_t(AT.size1());
653 if (! ublas::lu_factorize(LU, pvector))
655 const auto& BT = matrix_t(ublas::trans(lhs_ref.m_matrix));
657 ublas::lu_substitute(LU, pvector, XT);
659 return traits_t::NaN();
664 if (Tune_P::div_max_steps > 0)
667 auto R = matrix_t(-BT);
668 ublas::axpy_prod(AT, XT, R,
false);
670 return traits_t::NaN();
672 auto nr = Scalar_T(ublas::norm_inf(R));
673 if ( nr != Scalar_T(0) && !traits_t::isNaN_or_isInf(nr) )
676 auto nrold = nr + Scalar_T(1);
679 step != Tune_P::div_max_steps &&
689 ublas::lu_substitute(LU, pvector, D);
693 ublas::axpy_prod(AT, XTnew, R,
false);
694 nr = ublas::norm_inf(R);
698 return multivector_t(ublas::trans(XT), our_frame);
702 return traits_t::NaN();
706 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
710 operator/= (
const multivector_t& rhs) -> multivector_t&
711 {
return *
this = *
this / rhs; }
714 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
717 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>
718 {
return rhs * lhs / rhs.involute(); }
721 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
725 operator|= (
const multivector_t& rhs) -> multivector_t&
726 {
return *
this = rhs * *
this / rhs.involute(); }
729 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
733 inv() const -> const multivector_t
734 {
return multivector_t(Scalar_T(1), this->m_frame) / *
this; }
737 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
741 pow(
int m)
const ->
const multivector_t
745 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
748 outer_pow(
int m)
const ->
const multivector_t
751 throw error_t(
"outer_pow(m): negative exponent");
752 framed_multi_t a = *
this;
753 return a.outer_pow(m);
757 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
762 {
return framed_multi_t(*this).grade(); }
765 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
769 frame() const -> const index_set_t
770 {
return this->m_frame; }
773 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
777 operator[] (
const index_set_t ist)
const -> Scalar_T
780 if ( (ist | this->m_frame) == this->m_frame)
781 return matrix::inner<Scalar_T>(this->basis_element(ist), this->m_matrix);
787 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
793 if ((grade < 0) || (grade > HI-LO))
796 return (framed_multi_t(*
this))(grade);
800 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
804 scalar() const -> Scalar_T
806 const matrix_index_t dim = this->m_matrix.size1();
807 return matrix::trace(this->m_matrix) / Scalar_T(
double(dim) );
811 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
815 pure() const -> const multivector_t
816 {
return *
this - this->
scalar(); }
819 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
823 even() const -> const multivector_t
824 {
return framed_multi_t(*this).even(); }
827 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
831 odd() const -> const multivector_t
832 {
return framed_multi_t(*this).odd(); }
835 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
842 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
845 vector_part(
const index_set_t frm,
const bool prechecked)
const ->
const vector_t
847 if (!prechecked && (this->frame() | frm) != frm)
848 throw error_t(
"vector_part(frm): value is outside of requested frame");
851 if (this->frame() != frm)
852 return framed_multi_t(*this).vector_part(frm,
true);
854 const auto begin_index = frm.min();
855 const auto end_index = frm.max()+1;
864 matrix::inner<Scalar_T>(this->basis_element(index_set_t(idx)),
870 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
874 involute() const -> const multivector_t
875 {
return framed_multi_t(*this).involute(); }
878 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
882 reverse() const -> const multivector_t
883 {
return framed_multi_t(*this).reverse(); }
886 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
890 conj() const -> const multivector_t
891 {
return framed_multi_t(*this).conj(); }
894 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
898 quad() const -> Scalar_T
901 return framed_multi_t(*this).quad();
905 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
909 norm() const -> Scalar_T
911 const matrix_index_t dim = this->m_matrix.size1();
916 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
921 {
return framed_multi_t(*this).max_abs(); }
924 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
933 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
937 write(
const std::string& msg)
const
938 { framed_multi_t(*this).write(msg); }
941 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
945 write(std::ofstream& ofile,
const std::string& msg)
const
948 throw error_t(
"write(ofile,msg): cannot write to output file");
949 framed_multi_t(*this).write(ofile, msg);
953 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
958 os << typename matrix_multi<Scalar_T,LO,HI,Tune_P>::framed_multi_t(val);
963 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
978 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
982 isinf() const ->
bool
984 if (std::numeric_limits<Scalar_T>::has_infinity)
991 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
995 isnan() const ->
bool
997 if (std::numeric_limits<Scalar_T>::has_quiet_NaN)
1004 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1008 truncated(
const Scalar_T& limit)
const ->
const multivector_t
1009 {
return framed_multi_t(*this).truncated(limit); }
1012 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1016 operator+= (
const term_t& term) -> multivector_t&
1018 if (term.second != Scalar_T(0))
1019 this->m_matrix.plus_assign(matrix_t(this->basis_element(term.first)) * term.second);
1024 template<
typename Multivector_T,
typename Matrix_T,
typename Basis_Matrix_T >
1029 using framed_multi_t = Multivector_T;
1031 using index_set_t =
typename framed_multi_t::index_set_t;
1032 using Scalar_T =
typename framed_multi_t::scalar_t;
1033 using matrix_t = Matrix_T;
1034 using basis_matrix_t = Basis_Matrix_T;
1035 using basis_scalar_t =
typename basis_matrix_t::value_type;
1039 return framed_multi_t(traits_t::to_scalar_t(X(0,0)));
1041 if (ublas::norm_inf(X) == 0)
1044 const basis_matrix_t& I = matrix::unit<basis_matrix_t>(2);
1045 basis_matrix_t J(2,2,2);
1047 J(0,1) = basis_scalar_t(-1);
1048 J(1,0) = basis_scalar_t( 1);
1049 basis_matrix_t K = J;
1050 K(0,1) = basis_scalar_t( 1);
1051 basis_matrix_t JK = I;
1052 JK(0,0) = basis_scalar_t(-1);
1055 const index_set_t ist_mn = index_set_t(-level);
1056 const index_set_t ist_pn = index_set_t(level);
1057 const index_set_t ist_mnpn = ist_mn | ist_pn;
1060 using term_t =
typename framed_multi_t::term_t;
1061 const Scalar_T i_x = traits_t::to_scalar_t(signed_perm_nork(I, X)(0, 0));
1062 const Scalar_T j_x = traits_t::to_scalar_t(signed_perm_nork(J, X)(0, 0));
1063 const Scalar_T k_x = traits_t::to_scalar_t(signed_perm_nork(K, X)(0, 0));
1064 const Scalar_T jk_x = traits_t::to_scalar_t(signed_perm_nork(JK,X)(0, 0));
1067 result += term_t(ist_mn, j_x);
1068 result += term_t(ist_pn, k_x);
1069 return result += term_t(ist_mnpn, jk_x);
1073 const framed_multi_t& mn = framed_multi_t(ist_mn);
1074 const framed_multi_t& pn = framed_multi_t(ist_pn);
1075 const framed_multi_t& mnpn = framed_multi_t(ist_mnpn);
1076 const framed_multi_t& i_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1077 (signed_perm_nork(I, X), level-1);
1078 const framed_multi_t& j_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1079 (signed_perm_nork(J, X), level-1);
1080 const framed_multi_t& k_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1081 (signed_perm_nork(K, X), level-1);
1082 const framed_multi_t& jk_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1083 (signed_perm_nork(JK,X), level-1);
1085 result = i_x.even() - jk_x.odd();
1086 result += (j_x.even() - k_x.odd()) * mn;
1087 result += (k_x.even() - j_x.odd()) * pn;
1088 return result += (jk_x.even() - i_x.odd()) * mnpn;
1093 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1099 if (this->m_frame == frm)
1102 return (this->
template fast_framed_multi<Scalar_T>()).template fast_matrix_multi<Scalar_T>(frm);
1106 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1107 template <
typename Other_Scalar_T>
1113 index_t p = this->m_frame.count_pos();
1114 index_t q = this->m_frame.count_neg();
1132 const index_t level = (p+q)/2;
1136 framed_multi_t val = fast<framed_multi_t, matrix_t, basis_matrix_t>(this->m_matrix, level);
1139 switch (
pos_mod(orig_p-orig_q, 8))
1149 if (orig_p-orig_q > 4)
1152 if (orig_p-orig_q < -3)
1157 return val.
unfold(this->m_frame);
1161 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Matrix_T >
1163 public std::map< std::pair< const index_set<LO,HI>, const index_set<LO,HI> >,
1184 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1189 using index_set_pair_t = std::pair<const index_set_t, const index_set_t>;
1190 const auto& unfolded_pair = index_set_pair_t(ist, this->m_frame);
1193 auto& basis_cache = basis_table_t::basis();
1195 const auto frame_count = this->m_frame.count();
1196 const auto use_cache = frame_count <=
index_t(Tune_P::basis_max_count);
1200 const auto basis_it = basis_cache.find(unfolded_pair);
1201 if (basis_it != basis_cache.end())
1202 return *(basis_it->second);
1204 const auto folded_set = ist.
fold(this->m_frame);
1205 const auto folded_frame = this->m_frame.fold();
1206 const auto& folded_pair = index_set_pair_t(folded_set, folded_frame);
1207 using basis_pair_t = std::pair<const index_set_pair_t, basis_matrix_t *>;
1210 const auto basis_it = basis_cache.find(folded_pair);
1211 if (basis_it != basis_cache.end())
1213 auto* result_ptr = basis_it->second;
1214 basis_cache.insert(basis_pair_t(unfolded_pair, result_ptr));
1218 const auto folded_max = folded_frame.max();
1219 const auto folded_min = folded_frame.min();
1220 const auto p = std::max(folded_max,
index_t(0));
1224 auto result = matrix::unit<basis_matrix_t>(dim);
1234 basis_cache.insert(basis_pair_t(folded_pair, result_ptr));
1235 basis_cache.insert(basis_pair_t(unfolded_pair, result_ptr));
1241 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P, const
size_t Size >
1246 const std::array<Scalar_T, Size>& numer,
1247 const std::array<Scalar_T, Size>& denom,
1258 return traits_t::NaN();
1261 const auto nbr_even_powers = Size/2 - 1;
1264 auto XX = std::vector<multivector_t>(nbr_even_powers);
1266 XX[1] = XX[0] * XX[0];
1269 k != nbr_even_powers;
1271 XX[k] = XX[k-2] * XX[1];
1274 auto N = multivector_t(numer[1]);
1277 k != nbr_even_powers;
1279 N += XX[k] * numer[2*k + 3];
1284 k != nbr_even_powers;
1286 N += XX[k] * numer[2*k + 2];
1287 auto D = multivector_t(denom[1]);
1290 k != nbr_even_powers;
1292 D += XX[k] * denom[2*k + 3];
1297 k != nbr_even_powers;
1299 D += XX[k] * denom[2*k + 2];
1304 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1311 const auto& invM =
inv(M);
1312 M = ((M + invM)/Scalar_T(2) + Scalar_T(1)) / Scalar_T(2);
1313 Y *= (invM + Scalar_T(1)) / Scalar_T(2);
1317 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1326 if (val == Scalar_T(0))
1329 static const auto sqrt_max_steps = Tune_P::db_sqrt_max_steps;
1335 step != sqrt_max_steps &&
norm(M - Scalar_T(1)) > norm_tol;
1346 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1355 if (val == Scalar_T(0))
1358 static const auto sqrt_max_steps = Tune_P::cr_sqrt_max_steps;
1359 auto Z = Scalar_T(2) * (Scalar_T(1) + val);
1360 auto Y = Scalar_T(1) - val;
1362 auto norm_Y =
norm(Y);
1365 step != sqrt_max_steps && norm_Y > norm_Y_tol;
1368 const auto old_norm_Y = norm_Y;
1371 if (Y.isnan() || (norm_Y > old_norm_Y * Scalar_T(2)))
1374 Z += Y * Scalar_T(2);
1376 return Z / Scalar_T(4);
1383 template<
typename Scalar_T >
1389 template<
typename Scalar_T >
1392 1.0, 27.0/4.0, 81.0/4.0, 2277.0/64.0,
1393 10395.0/256.0, 32319.0/1024.0, 8721.0/512.0, 26163.0/4096.0,
1394 53703.0/32768.0, 36465.0/131072.0, 3861.0/131072.0, 7371.0/4194304.0,
1395 819.0/16777216.0, 27.0/67108864.0
1400 template<
typename Scalar_T >
1406 template<
typename Scalar_T >
1409 1.0, 25.0/4.0, 69.0/4.0, 1771.0/64.0,
1410 7315.0/256.0, 20349.0/1024.0, 4845.0/512.0, 12597.0/4096.0,
1411 21879.0/32768.0, 12155.0/131072.0, 1001.0/131072.0, 1365.0/4194304.0,
1412 91.0/16777216.0, 1.0/67108864.0
1423 1.0, 19.0/4.0, 19.0/2.0, 665.0/64.0,
1424 1729.0/256.0, 2717.0/1024.0, 627.0/1024.0, 627.0/8192.0,
1425 285.0/65536.0, 19.0/262144.0
1435 1.0, 17.0/4.0, 15.0/2.0, 455.0/64.0,
1436 1001.0/256.0, 1287.0/1024.0, 231.0/1024.0, 165.0/8192.0,
1437 45.0/65536, 1.0/262144.0
1443 using array = std::array<long double, 18>;
1448 1.0L, 35.0L/4.0L, 35.0L, 5425.0L/64.0L,
1449 35525.0L/256.0L, 166257.0L/1024.0L, 143325.0L/1024.0L, 740025.0L/8192.0L,
1450 2877875.0L/65536.0L, 4206125.0L/262144.0L, 572033.0L/131072.0L, 1820105.0L/2097152.0L,
1451 1028755.0L/8388608.0L, 395675.0L/33554432.0L, 24225.0L/33554432.0L, 6783.0L/268435456.0L,
1452 1785.0L/4294967296.0L, 35.0L/17179869184.0L
1457 using array = std::array<long double, 18>;
1462 1.0L, 33.0L/4.0L, 31.0L, 4495.0L/64.0L,
1463 27405.0L/256.0L, 118755.0L/1024.0L, 94185.0L/1024.0L, 444015.0L/8192.0L,
1464 1562275.0L/65536.0L, 2042975.0L/262144.0L, 245157.0L/131072.0L, 676039.0L/2097152.0L,
1465 323323.0L/8388608.0L, 101745.0L/33554432.0L, 4845.0L/33554432.0L, 969.0L/268435456.0L,
1466 153.0L/4294967296.0L, 1.0L/17179869184.0L
1469#if defined(_GLUCAT_USE_QD)
1478 dd_real(
"1"), dd_real(
"43")/dd_real(
"4"),
1479 dd_real(
"215")/dd_real(
"4"), dd_real(
"10621")/dd_real(
"64"),
1480 dd_real(
"90687")/dd_real(
"256"), dd_real(
"567987")/dd_real(
"1024"),
1481 dd_real(
"168861")/dd_real(
"256"), dd_real(
"1246355")/dd_real(
"2048"),
1482 dd_real(
"7228859")/dd_real(
"16384"), dd_real(
"16583853")/dd_real(
"65536"),
1483 dd_real(
"7538115")/dd_real(
"65536"), dd_real(
"173376645")/dd_real(
"4194304"),
1484 dd_real(
"195747825")/dd_real(
"16777216"), dd_real(
"171655785")/dd_real(
"67108864"),
1485 dd_real(
"14375115")/dd_real(
"33554432"), dd_real(
"14375115")/dd_real(
"268435456"),
1486 dd_real(
"20764055")/dd_real(
"4294967296"), dd_real(
"5167525")/dd_real(
"17179869184"),
1487 dd_real(
"206701")/dd_real(
"17179869184"), dd_real(
"76153")/dd_real(
"274877906944"),
1488 dd_real(
"3311")/dd_real(
"1099511627776") , dd_real(
"43")/dd_real(
"4398046511104")
1498 dd_real(
"1"), dd_real(
"41")/dd_real(
"4"),
1499 dd_real(
"195")/dd_real(
"4"), dd_real(
"9139")/dd_real(
"64"),
1500 dd_real(
"73815")/dd_real(
"256"), dd_real(
"435897")/dd_real(
"1024"),
1501 dd_real(
"121737")/dd_real(
"256"), dd_real(
"840565")/dd_real(
"2048"),
1502 dd_real(
"4539051")/dd_real(
"16384"), dd_real(
"9641775")/dd_real(
"65536"),
1503 dd_real(
"4032015")/dd_real(
"65536"), dd_real(
"84672315")/dd_real(
"4194304"),
1504 dd_real(
"86493225")/dd_real(
"16777216"), dd_real(
"67863915")/dd_real(
"67108864"),
1505 dd_real(
"5014575")/dd_real(
"33554432"), dd_real(
"4345965")/dd_real(
"268435456"),
1506 dd_real(
"5311735")/dd_real(
"4294967296"), dd_real(
"1081575")/dd_real(
"17179869184"),
1507 dd_real(
"33649")/dd_real(
"17179869184"), dd_real(
"8855")/dd_real(
"274877906944"),
1508 dd_real(
"231")/dd_real(
"1099511627776"), dd_real(
"1")/dd_real(
"4398046511104")
1519 qd_real(
"1"), qd_real(
"67")/qd_real(
"4"),
1520 qd_real(
"134"), qd_real(
"43617")/qd_real(
"64"),
1521 qd_real(
"633485")/qd_real(
"256"), qd_real(
"6992857")/qd_real(
"1024"),
1522 qd_real(
"15246721")/qd_real(
"1024"), qd_real(
"215632197")/qd_real(
"8192"),
1523 qd_real(
"2518145487")/qd_real(
"65536"), qd_real(
"12301285425")/qd_real(
"262144"),
1524 qd_real(
"6344873535")/qd_real(
"131072"), qd_real(
"89075432355")/qd_real(
"2097152"),
1525 qd_real(
"267226297065")/qd_real(
"8388608"), qd_real(
"687479618945")/qd_real(
"33554432"),
1526 qd_real(
"379874182975")/qd_real(
"33554432"), qd_real(
"1443521895305")/qd_real(
"268435456"),
1527 qd_real(
"9425348845815")/qd_real(
"4294967296"), qd_real(
"13195488384141")/qd_real(
"17179869184"),
1528 qd_real(
"987417498133")/qd_real(
"4294967296"), qd_real(
"8055248011085")/qd_real(
"137438953472"),
1529 qd_real(
"6958363175533")/qd_real(
"549755813888"), qd_real(
"5056698705201")/qd_real(
"2199023255552"),
1530 qd_real(
"766166470485")/qd_real(
"2199023255552"), qd_real(
"766166470485")/qd_real(
"17592186044416"),
1531 qd_real(
"623623871325")/qd_real(
"140737488355328"), qd_real(
"203123203803")/qd_real(
"562949953421312"),
1532 qd_real(
"6478601247")/qd_real(
"281474976710656"), qd_real(
"5038912081")/qd_real(
"4503599627370496"),
1533 qd_real(
"719844583")/qd_real(
"18014398509481984"), qd_real(
"71853815")/qd_real(
"72057594037927936"),
1534 qd_real(
"1165197")/qd_real(
"72057594037927936"), qd_real(
"87703")/qd_real(
"576460752303423488"),
1535 qd_real(
"12529")/qd_real(
"18446744073709551616"), qd_real(
"67")/qd_real(
"73786976294838206464")
1545 qd_real(
"1"), qd_real(
"65")/qd_real(
"4"),
1546 qd_real(
"126"), qd_real(
"39711")/qd_real(
"64"),
1547 qd_real(
"557845")/qd_real(
"256"), qd_real(
"5949147")/qd_real(
"1024"),
1548 qd_real(
"12515965")/qd_real(
"1024"), qd_real(
"170574723")/qd_real(
"8192"),
1549 qd_real(
"1916797311")/qd_real(
"65536"), qd_real(
"8996462475")/qd_real(
"262144"),
1550 qd_real(
"4450881435")/qd_real(
"131072"), qd_real(
"59826782925")/qd_real(
"2097152"),
1551 qd_real(
"171503444385")/qd_real(
"8388608"), qd_real(
"420696483235")/qd_real(
"33554432"),
1552 qd_real(
"221120793075")/qd_real(
"33554432"), qd_real(
"797168807855")/qd_real(
"268435456"),
1553qd_real(
"4923689695575")/qd_real(
"4294967296"), qd_real(
"6499270398159")/qd_real(
"17179869184"),
1554 qd_real(
"456864812569")/qd_real(
"4294967296"), qd_real(
"3486599885395")/qd_real(
"137438953472"),
1555qd_real(
"2804116503573")/qd_real(
"549755813888"), qd_real(
"1886827875075")/qd_real(
"2199023255552"),
1556 qd_real(
"263012370465")/qd_real(
"2199023255552"), qd_real(
"240141729555")/qd_real(
"17592186044416"),
1557 qd_real(
"176848560525")/qd_real(
"140737488355328"), qd_real(
"51538723353")/qd_real(
"562949953421312"),
1558 qd_real(
"1450433115")/qd_real(
"281474976710656"), qd_real(
"977699359")/qd_real(
"4503599627370496"),
1559 qd_real(
"118183439")/qd_real(
"18014398509481984"), qd_real(
"9652005")/qd_real(
"72057594037927936"),
1560 qd_real(
"121737")/qd_real(
"72057594037927936"), qd_real(
"6545")/qd_real(
"576460752303423488"),
1561 qd_real(
"561")/qd_real(
"18446744073709551616"), qd_real(
"1")/qd_real(
"73786976294838206464")
1569 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1582 return traits_t::NaN();
1584 const auto scr_val = val.
scalar();
1587 if (scr_val < Scalar_T(0))
1588 return i * traits_t::sqrt(-scr_val);
1590 return traits_t::sqrt(scr_val);
1595 (scr_val != Scalar_T(0) &&
norm(val/scr_val - Scalar_T(1)) < Scalar_T(1))
1597 : (scr_val < Scalar_T(0))
1600 const auto sqrt_scale = traits_t::sqrt(traits_t::abs(scale));
1601 if (traits_t::isNaN_or_isInf(sqrt_scale))
1602 return traits_t::NaN();
1605 auto rescale = multivector_t(sqrt_scale);
1606 if (scale < Scalar_T(0))
1607 rescale = i * sqrt_scale;
1609 const auto& unitval = val / scale;
1610 static const auto max_norm = Scalar_T(1.0/4.0);
1611 auto use_approx_sqrt =
true;
1612 auto use_cr_sqrt =
false;
1613 auto scaled_result = multivector_t();
1614#if defined(_GLUCAT_USE_EIGENVALUES)
1615 static const auto sqrt_2 = traits_t::sqrt(Scalar_T(2));
1618 using matrix_t =
typename multivector_t::matrix_t;
1623 (genus.m_is_singular)
1626 switch (genus.m_eig_case)
1628 case matrix::neg_real_eigs:
1629 scaled_result =
matrix_sqrt(-i * unitval, i, next_level) * (i + Scalar_T(1)) / sqrt_2;
1630 use_approx_sqrt =
false;
1632 case matrix::both_eigs:
1634 const auto safe_arg = genus.m_safe_arg;
1635 scaled_result =
matrix_sqrt(
exp(i*safe_arg) * unitval, i, next_level) *
exp(-i*safe_arg / Scalar_T(2));
1637 use_approx_sqrt =
false;
1642 use_cr_sqrt = genus.m_is_singular;
1645 if (use_approx_sqrt)
1648 (
norm(unitval - Scalar_T(1)) < max_norm)
1652 unitval - Scalar_T(1))
1658 return (scaled_result.isnan() ||
1661 : scaled_result * rescale;
1665 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1667 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>
1676 return traits_t::NaN();
1680 switch (Tune_P::function_precision)
1682 case precision_demoted:
1684 using demoted_scalar_t =
typename traits_t::demoted::type;
1687 const auto& demoted_val = demoted_multivector_t(val);
1688 const auto& demoted_i = demoted_multivector_t(i);
1693 case precision_promoted:
1695 using promoted_scalar_t =
typename traits_t::promoted::type;
1698 const auto& promoted_val = promoted_multivector_t(val);
1699 const auto& promoted_i = promoted_multivector_t(i);
1713 template<
typename Scalar_T >
1719 template<
typename Scalar_T >
1722 0.0, 1.0, 6.0, 4741.0/300.0,
1723 1441.0/60.0, 107091.0/4600.0, 8638.0/575.0, 263111.0/40250.0,
1724 153081.0/80500.0, 395243.0/1101240.0, 28549.0/688275.0, 605453.0/228813200.0,
1725 785633.0/10296594000.0, 1145993.0/1873980108000.0
1730 template<
typename Scalar_T >
1736 template<
typename Scalar_T >
1739 1.0, 13.0/2.0, 468.0/25.0, 1573.0/50.0,
1740 1573.0/46.0, 11583.0/460.0, 10296.0/805.0, 2574.0/575.0,
1741 11583.0/10925.0, 143.0/874.0, 572.0/37145.0, 117.0/148580.0,
1742 13.0/742900.0, 1.0/10400600.0
1753 0.0, 1.0, 4.0, 1337.0/204.0,
1754 385.0/68.0, 1879.0/680.0, 193.0/255.0, 197.0/1820.0,
1755 419.0/61880.0, 7129.0/61261200.0
1765 1.0, 9.0/2.0, 144.0/17.0, 147.0/17.0,
1766 441.0/85.0, 63.0/34.0, 84.0/221.0, 9.0/221.0,
1767 9.0/4862.0, 1.0/48620.0
1773 using array = std::array<long double, 18>;
1778 0.0L, 1.0L, 8.0L, 3835.0L/132.0L,
1779 8365.0L/132.0L, 11363807.0L/122760.0L, 162981.0L/1705.0L, 9036157.0L/125860.0L,
1780 18009875.0L/453096.0L, 44211925.0L/2718576.0L, 4149566.0L/849555.0L, 16973929.0L/16020180.0L,
1781 172459.0L/1068012.0L, 116317061.0L/7025382936.0L, 19679783.0L/18441630207.0L, 23763863.0L/614721006900.0L,
1782 50747.0L/79318839600.0L, 42142223.0L/14295951736466400.0L
1787 using array = std::array<long double, 18>;
1792 1.0L, 17.0L/2.0L, 1088.0L/33.0L, 850.0L/11.0L,
1793 41650.0L/341.0L, 140777.0L/1023.0L, 1126216.0L/9889.0L, 63206.0L/899.0L,
1794 790075.0L/24273.0L, 60775.0L/5394.0L, 38896.0L/13485.0L, 21658.0L/40455.0L,
1795 21658.0L/310155.0L, 4165.0L/682341.0L, 680.0L/2047023.0L, 34.0L/3411705.0L,
1796 17.0L/129644790.0L, 1.0L/2333606220
1798#if defined(_GLUCAT_USE_QD)
1807 dd_real(
"0"), dd_real(
"1"),
1808 dd_real(
"10"), dd_real(
"22781")/dd_real(
"492"),
1809 dd_real(
"21603")/dd_real(
"164"), dd_real(
"5492649")/dd_real(
"21320"),
1810 dd_real(
"978724")/dd_real(
"2665"), dd_real(
"4191605")/dd_real(
"10619"),
1811 dd_real(
"12874933")/dd_real(
"39442"), dd_real(
"11473457")/dd_real(
"54612"),
1812 dd_real(
"2406734")/dd_real(
"22755"), dd_real(
"166770367")/dd_real(
"4004880"),
1813 dd_real(
"30653165")/dd_real(
"2402928"), dd_real(
"647746389")/dd_real(
"215195552"),
1814 dd_real(
"25346331")/dd_real(
"47074027"), dd_real(
"278270613")/dd_real(
"3900419380"),
1815 dd_real(
"105689791")/dd_real(
"15601677520"), dd_real(
"606046475")/dd_real(
"1379188292768"),
1816 dd_real(
"969715")/dd_real(
"53502994116"), dd_real(
"11098301")/dd_real(
"26204577562592"),
1817 dd_real(
"118999")/dd_real(
"26204577562592"), dd_real(
"18858053")/dd_real(
"1392249205900512960")
1827 dd_real(
"1"), dd_real(
"21")/dd_real(
"2"),
1828 dd_real(
"2100")/dd_real(
"41"), dd_real(
"12635")/dd_real(
"82"),
1829 dd_real(
"341145")/dd_real(
"1066"), dd_real(
"1037799")/dd_real(
"2132"),
1830 dd_real(
"11069856")/dd_real(
"19721"), dd_real(
"9883800")/dd_real(
"19721"),
1831 dd_real(
"6918660")/dd_real(
"19721"), dd_real(
"293930")/dd_real(
"1517"),
1832 dd_real(
"1410864")/dd_real(
"16687"), dd_real(
"88179")/dd_real(
"3034"),
1833 dd_real(
"734825")/dd_real(
"94054"), dd_real(
"305235")/dd_real(
"188108"),
1834 dd_real(
"348840")/dd_real(
"1363783"), dd_real(
"40698")/dd_real(
"1363783"),
1835 dd_real(
"6783")/dd_real(
"2727566"), dd_real(
"9975")/dd_real(
"70916716"),
1836 dd_real(
"266")/dd_real(
"53187537"), dd_real(
"7")/dd_real(
"70916716"),
1837 dd_real(
"7")/dd_real(
"8155422340"), dd_real(
"1")/dd_real(
"538257874440")
1848 qd_real(
"0"), qd_real(
"1"),
1849 qd_real(
"16"), qd_real(
"95201")/qd_real(
"780"),
1850 qd_real(
"30721")/qd_real(
"52"), qd_real(
"7416257")/qd_real(
"3640"),
1851 qd_real(
"1039099")/qd_real(
"195"), qd_real(
"6097772319")/qd_real(
"555100"),
1852 qd_real(
"1564058073")/qd_real(
"85400"), qd_real(
"30404640205")/qd_real(
"1209264"),
1853 qd_real(
"725351278")/qd_real(
"25193"), qd_real(
"4092322670789")/qd_real(
"147429436"),
1854 qd_real(
"4559713849589")/qd_real(
"201040140"), qd_real(
"5049361751189")/qd_real(
"320023080"),
1855 qd_real(
"74979677195")/qd_real(
"8000577"), qd_real(
"16569850691873")/qd_real(
"3481514244"),
1856 qd_real(
"1065906022369")/qd_real(
"515779888"), qd_real(
"335956770855841")/qd_real(
"438412904800"),
1857 qd_real(
"1462444287585964")/qd_real(
"6041877844275"), qd_real(
"397242326339851")/qd_real(
"6122436215532"),
1858 qd_real(
"64211291334131")/qd_real(
"4373168725380"), qd_real(
"142322343550859")/qd_real(
"51080680851480"),
1859 qd_real(
"154355972958659")/qd_real(
"351179680853925"), qd_real(
"167483568676259")/qd_real(
"2937139148960100"),
1860 qd_real(
"4230788929433")/qd_real(
"704913395750424"), qd_real(
"197968763176019")/qd_real(
"392923948371995600"),
1861 qd_real(
"10537522306718")/qd_real(
"319250708052246425"), qd_real(
"236648286272519")/qd_real(
"144249197475035425500"),
1862 qd_real(
"260715545088119")/qd_real(
"4375558990076074573500"), qd_real(
"289596255666839")/qd_real(
"192874640282553367199880"),
1863 qd_real(
"8802625510547")/qd_real(
"361639950529787563499775"), qd_real(
"373831661521439")/qd_real(
"1659204093030665341336967700"),
1864 qd_real(
"446033437968239")/qd_real(
"464577146048586295574350956000"), qd_real(
"53676090078349")/qd_real(
"47386868896955802148583797512000")
1874 qd_real(
"1"), qd_real(
"33")/qd_real(
"2"),
1875 qd_real(
"8448")/qd_real(
"65"), qd_real(
"42284")/qd_real(
"65"),
1876 qd_real(
"211420")/qd_real(
"91"), qd_real(
"573562")/qd_real(
"91"),
1877 qd_real(
"32119472")/qd_real(
"2379"), qd_real(
"92917044")/qd_real(
"3965"),
1878 qd_real(
"603960786")/qd_real(
"17995"), qd_real(
"144626625")/qd_real(
"3599"),
1879 qd_real(
"2776831200")/qd_real(
"68381"), qd_real(
"16692542100")/qd_real(
"478667"),
1880 qd_real(
"12241197540")/qd_real(
"478667"), qd_real(
"1098569010")/qd_real(
"68381"),
1881 qd_real(
"31387686000")/qd_real(
"3624193"), qd_real(
"9939433900")/qd_real(
"2479711"),
1882 qd_real(
"67091178825")/qd_real(
"42155087"), qd_real(
"2683647153")/qd_real(
"4959422"),
1883 qd_real(
"19083713088")/qd_real(
"121505839"), qd_real(
"4708152900")/qd_real(
"121505839"),
1884 qd_real(
"941630580")/qd_real(
"116546417"), qd_real(
"88704330")/qd_real(
"62755763"),
1885 qd_real(
"12902448")/qd_real(
"62755763"), qd_real(
"1542684")/qd_real(
"62755763"),
1886 qd_real(
"6427850")/qd_real(
"2698497809"), qd_real(
"3471039")/qd_real(
"18889484663"),
1887 qd_real(
"8544096")/qd_real(
"774468871183"), qd_real(
"39556")/qd_real(
"79027435835"),
1888 qd_real(
"118668")/qd_real(
"7191496660985"), qd_real(
"10230")/qd_real(
"27327687311743"),
1889 qd_real(
"5456")/qd_real(
"1011124430534491"), qd_real(
"44")/qd_real(
"1011124430534491"),
1890 qd_real(
"11")/qd_real(
"70778710137414370"), qd_real(
"1")/qd_real(
"7219428434016265740")
1897 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1908 if (val == Scalar_T(0) || val.
isnan())
1909 return traits_t::NaN();
1917 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1925 if (val == Scalar_T(0) || val.
isnan())
1926 return traits_t::NaN();
1928 using limits_t = std::numeric_limits<Scalar_T>;
1929 static const auto epsilon = limits_t::epsilon();
1930 static const auto max_inner_norm = traits_t::pow(
epsilon, 2);
1931 static const auto max_outer_norm = Scalar_T(6.0/limits_t::digits);
1933 auto E = multivector_t(Scalar_T(0));
1935 auto pow_2_outer_step = Scalar_T(1);
1936 auto pow_4_outer_step = Scalar_T(1);
1938 for (outer_step = 0, norm_Y_1 =
norm(Y - Scalar_T(1));
1939 outer_step != Tune_P::log_max_outer_steps && norm_Y_1 * pow_2_outer_step > max_outer_norm;
1940 ++outer_step, norm_Y_1 =
norm(Y - Scalar_T(1)))
1942 if (Y == Scalar_T(0) || Y.isnan())
1943 return traits_t::NaN();
1949 inner_step != Tune_P::log_max_inner_steps &&
1950 norm(M - Scalar_T(1)) * pow_4_outer_step > max_inner_norm;
1954 E += (M - Scalar_T(1)) * pow_2_outer_step;
1955 pow_2_outer_step *= Scalar_T(2);
1956 pow_4_outer_step *= Scalar_T(4);
1958 if (outer_step == Tune_P::log_max_outer_steps && norm_Y_1 * pow_2_outer_step > max_outer_norm)
1959 return traits_t::NaN();
1961 return pade_log(Y) * pow_2_outer_step - E;
1965 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1975 if (val == Scalar_T(0) || val.
isnan())
1976 return traits_t::NaN();
1978 static const auto pi = traits_t::pi();
1979 const auto scr_val = val.
scalar();
1982 if (scr_val < Scalar_T(0))
1983 return i * pi + traits_t::log(-scr_val);
1985 return traits_t::log(scr_val);
1989 const auto max_norm = Scalar_T(1.0/9.0);
1991 (scr_val != Scalar_T(0) &&
norm(val/scr_val - Scalar_T(1)) < max_norm)
1993 : (scr_val < Scalar_T(0))
1996 if (scale == Scalar_T(0))
1997 return traits_t::NaN();
2000 const auto log_scale = traits_t::log(traits_t::abs(scale));
2001 auto rescale = multivector_t(log_scale);
2002 if (scale < Scalar_T(0))
2003 rescale = i * pi + log_scale;
2004 const auto unitval = val/scale;
2005 if (
inv(unitval).isnan())
2006 return traits_t::NaN();
2008#if defined(_GLUCAT_USE_EIGENVALUES)
2009 auto scaled_result = multivector_t();
2012 using matrix_t =
typename multivector_t::matrix_t;
2016 switch (genus.m_eig_case)
2018 case matrix::neg_real_eigs:
2019 scaled_result =
matrix_log(-i * unitval, i, level + 1) + i * pi/Scalar_T(2);
2021 case matrix::both_eigs:
2023 const Scalar_T safe_arg = genus.m_safe_arg;
2024 scaled_result =
matrix_log(
exp(i*safe_arg) * unitval, i, level + 1) - i * safe_arg;
2037 return (scaled_result.isnan())
2039 : scaled_result + rescale;
2043 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
2045 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>
2049 if (val == Scalar_T(0) || val.
isnan())
2050 return traits_t::NaN();
2054 switch (Tune_P::function_precision)
2056 case precision_demoted:
2058 using demoted_scalar_t =
typename traits_t::demoted::type;
2061 const auto& demoted_val = demoted_multivector_t(val);
2062 const auto& demoted_i = demoted_multivector_t(i);
2064 return matrix_log(demoted_val, demoted_i, 0);
2067 case precision_promoted:
2069 using promoted_scalar_t =
typename traits_t::promoted::type;
2072 const auto& promoted_val = promoted_multivector_t(val);
2073 const auto& promoted_i = promoted_multivector_t(i);
2075 return matrix_log(promoted_val, promoted_i, 0);
2084 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
2090 return traits_t::NaN();
2092 const auto scr_val = val.
scalar();
2094 return traits_t::exp(scr_val);
2096 switch (Tune_P::function_precision)
2098 case precision_demoted:
2100 using demoted_scalar_t =
typename traits_t::demoted::type;
2103 const auto& demoted_val = demoted_multivector_t(val);
2107 case precision_promoted:
2109 using promoted_scalar_t =
typename traits_t::promoted::type;
2112 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
Write formatted multivector to output.
virtual auto involute() const -> const multivector_t=0
Main involution, each {i} is replaced by -{i} in each term, eg. {1} -> -{1}.
virtual auto outer_pow(int m) const -> const multivector_t=0
Outer product power.
virtual auto isinf() const -> bool=0
Check if a multivector contains any infinite values.
virtual auto quad() const -> Scalar_T=0
Scalar_T quadratic form == (rev(x)*x)(0)
virtual auto scalar() const -> Scalar_T=0
Scalar part.
virtual auto operator()(index_t grade) const -> const multivector_t=0
Pure grade-vector part.
virtual auto conj() const -> const multivector_t=0
Conjugation, reverse o involute == involute o reverse.
virtual auto operator-=(const multivector_t &rhs) -> multivector_t &=0
Geometric difference.
virtual auto grade() const -> index_t=0
Maximum of the grades of each term.
virtual auto even() const -> const multivector_t=0
Even part of multivector, sum of even grade terms.
virtual auto frame() const -> const index_set_t=0
Subalgebra generated by all generators of terms of given multivector.
virtual auto operator/=(const Scalar_T &scr) -> multivector_t &=0
Quotient of multivector and scalar.
virtual auto max_abs() const -> Scalar_T=0
Maximum of absolute values of components of multivector: multivector infinity norm.
virtual auto odd() const -> const multivector_t=0
Odd part of multivector, sum of odd grade terms.
virtual auto vector_part() const -> const vector_t=0
Vector part of multivector, as a vector_t with respect to frame()
virtual auto reverse() const -> const multivector_t=0
Reversion, eg. {1}*{2} -> {2}*{1}.
virtual auto operator%=(const multivector_t &rhs) -> multivector_t &=0
Contraction.
virtual auto operator-() const -> const multivector_t=0
Unary -.
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
Transformation via twisted adjoint action.
virtual auto pure() const -> const multivector_t=0
Pure part.
virtual auto operator==(const multivector_t &val) const -> bool=0
Test for equality of multivectors.
virtual auto operator^=(const multivector_t &rhs) -> multivector_t &=0
Outer product.
virtual auto isnan() const -> bool=0
Check if a multivector contains any IEEE NaN values.
virtual auto operator*=(const Scalar_T &scr) -> multivector_t &=0
Product of multivector and scalar.
virtual auto norm() const -> Scalar_T=0
Scalar_T norm == sum of norm of coordinates.
virtual auto pow(int m) const -> const multivector_t=0
*this to the m
virtual auto operator&=(const multivector_t &rhs) -> multivector_t &=0
Inner product.
virtual auto operator[](const index_set_t ist) const -> Scalar_T=0
Subscripting: map from index set to scalar coordinate.
virtual auto inv() const -> const multivector_t=0
Geometric multiplicative inverse.
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.
auto fold() const -> const index_set_t
Fold this index set within itself as a frame.
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
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 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 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.
int index_t
Size of index_t should be enough to represent LO, HI.
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.
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