1#ifndef _GLUCAT_FRAMED_MULTI_IMP_H
2#define _GLUCAT_FRAMED_MULTI_IMP_H
48 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
52 {
return "framed_multi"; }
54#define _GLUCAT_HASH_N(x) (x)
55#define _GLUCAT_HASH_SIZE_T(x) (typename multivector_t::hash_size_t)(x)
58 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
65 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
72 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
73 template<
typename Other_Scalar_T >
78 for (
auto& val_term : val)
83 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
84 template<
typename Other_Scalar_T >
90 if (!prechecked && (val.
frame() | frm) != frm)
91 throw error_t(
"multivector_t(val,frm): cannot initialize with value outside of frame");
92 for (
auto& val_term : val)
97 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
103 if (!prechecked && (val.
frame() | frm) != frm)
104 throw error_t(
"multivector_t(val,frm): cannot initialize with value outside of frame");
105 for (
auto& val_term : val)
106 this->insert(val_term);
110 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
115 if (crd != Scalar_T(0))
116 this->insert(
term_t(ist, crd));
120 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
126 if (!prechecked && (ist | frm) != frm)
127 throw error_t(
"multivector_t(ist,crd,frm): cannot initialize with value outside of frame");
128 if (crd != Scalar_T(0))
129 this->insert(
term_t(ist, crd));
133 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
138 if (scr != Scalar_T(0))
143 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
148 if (scr != Scalar_T(0))
153 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
160 throw error_t(
"multivector_t(vec,frm): cannot initialize with vector not matching frame");
161 auto idx = frm.
min();
162 const auto frm_end = frm.
max()+1;
163 for (
auto& crd : vec)
168 idx != frm_end && !frm[idx];
175 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
180 std::istringstream ss(str);
183 throw error_t(
"multivector_t(str): could not parse string");
187 throw error_t(
"multivector_t(str): could not parse entire string");
191 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
203 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
204 template<
typename Other_Scalar_T >
209 if (val == Other_Scalar_T(0))
212 const auto dim = val.
m_matrix.size1();
219 if (dim >= Tune_P::inv_fast_dim_threshold)
222 *
this = (val.template fast_framed_multi<Scalar_T>()).
truncated();
228 const auto val_norm = traits_t::to_scalar_t(val.
norm());
229 if (traits_t::isNaN_or_isInf(val_norm))
231 *
this = traits_t::NaN();
234 const auto frm = val.
frame();
235 const auto algebra_dim =
set_value_t(1) << frm.count();
246 if (crd != Scalar_T(0))
247 result.insert(
term_t(ist, crd));
249 *
this = result.truncated();
253 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
256 operator== (
const multivector_t& rhs)
const ->
bool
258 if (this->size() != rhs.size())
260 const auto rhs_end = rhs.end();
261 for (
auto& this_term : *
this)
263 const const_iterator& rhs_it = rhs.find(this_term.first);
264 if (rhs_it == rhs_end || rhs_it->second != this_term.second)
271 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
275 operator== (
const Scalar_T& scr)
const ->
bool
277 switch (this->size())
280 return scr == Scalar_T(0);
283 const auto& this_it = this->begin();
284 return this_it->first == index_set_t() && this_it->second == scr;
292 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
303 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
307 operator+= (
const multivector_t& rhs) -> multivector_t&
309 for (
auto& rhs_term : rhs)
315 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
319 operator-= (
const Scalar_T& scr) -> multivector_t&
321 *
this += term_t(index_set_t(), -scr);
326 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
330 operator-= (
const multivector_t& rhs) -> multivector_t&
332 for (
auto& rhs_term : rhs)
333 *
this += term_t(rhs_term.first, -(rhs_term.second));
338 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
342 operator- () const -> const multivector_t
345 for (
auto& result_term : result)
346 result_term.second *= Scalar_T(-1);
351 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
354 operator*= (
const Scalar_T& scr) -> multivector_t&
356 using traits_t = numeric_traits<Scalar_T>;
358 if (traits_t::isNaN_or_isInf(scr))
359 return *
this = traits_t::NaN();
360 if (scr == Scalar_T(0))
362 *
this = traits_t::NaN();
366 for (
auto& this_term : *
this)
367 this_term.second *= scr;
372 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
374 operator* (
const framed_multi<Scalar_T,LO,HI,Tune_P>& lhs,
const framed_multi<Scalar_T,LO,HI,Tune_P>& rhs) ->
const framed_multi<Scalar_T,LO,HI,Tune_P>
379 if (lhs.isnan() || rhs.isnan())
380 return traits_t::NaN();
382 const double lhs_size = lhs.size();
383 const double rhs_size = rhs.size();
384 const auto our_frame = lhs.frame() | rhs.frame();
385 const auto frm_count = our_frame.count();
386 const auto algebra_dim =
set_value_t(1) << frm_count;
387 const auto direct_mult = lhs_size * rhs_size <= double(algebra_dim);
390 auto result = multivector_t(
392 for (
auto& lhs_term : lhs)
393 for (
auto& rhs_term : rhs)
394 result += lhs_term * rhs_term;
399 using matrix_multi_t =
typename multivector_t::matrix_multi_t;
400 return matrix_multi_t(lhs, our_frame,
true) *
401 matrix_multi_t(rhs, our_frame,
true);
406 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
410 operator*= (
const multivector_t& rhs) -> multivector_t&
411 {
return *
this = *
this * rhs; }
414 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
416 operator^ (
const framed_multi<Scalar_T,LO,HI,Tune_P>& lhs,
const framed_multi<Scalar_T,LO,HI,Tune_P>& rhs) ->
const framed_multi<Scalar_T,LO,HI,Tune_P>
420 if (lhs.empty() || rhs.empty())
424 using index_set_t =
typename multivector_t::index_set_t;
425 using term_t =
typename multivector_t::term_t;
427 const auto empty_set = index_set_t();
429 const double lhs_size = lhs.size();
430 const double rhs_size = rhs.size();
431 const auto lhs_frame = lhs.frame();
432 const auto rhs_frame = rhs.frame();
433 const auto our_frame = lhs_frame | rhs_frame;
434 const auto algebra_dim =
set_value_t(1) << our_frame.count();
435 auto result = multivector_t(
437 const auto lhs_end = lhs.end();
438 const auto rhs_end = rhs.end();
440 if (lhs_size * rhs_size >
double(Tune_P::products_size_threshold))
444 result_stv != algebra_dim;
447 const auto result_ist = index_set_t(result_stv, our_frame,
true);
448 const auto lhs_result_frame = lhs_frame & result_ist;
449 const auto lhs_result_dim =
set_value_t(1) << lhs_result_frame.count();
450 auto result_crd = Scalar_T(0);
453 lhs_stv != lhs_result_dim;
456 const auto lhs_ist = index_set_t(lhs_stv, lhs_result_frame,
true);
457 const auto rhs_ist = result_ist ^ lhs_ist;
458 if ((rhs_ist | rhs_frame) == rhs_frame)
460 const auto lhs_it = lhs.find(lhs_ist);
461 if (lhs_it != lhs_end)
463 const auto rhs_it = rhs.find(rhs_ist);
464 if (rhs_it != rhs_end)
469 if (result_crd != Scalar_T(0))
470 result.insert(term_t(result_ist, result_crd));
476 for (
auto& lhs_term : lhs)
477 for (
auto& rhs_term : rhs)
478 if ((lhs_term.first & rhs_term.first) == empty_set)
479 result += lhs_term * rhs_term;
485 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
489 operator^= (
const multivector_t& rhs) -> multivector_t&
490 {
return *
this = *
this ^ rhs; }
493 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
495 operator& (
const framed_multi<Scalar_T,LO,HI,Tune_P>& lhs,
const framed_multi<Scalar_T,LO,HI,Tune_P>& rhs) ->
const framed_multi<Scalar_T,LO,HI,Tune_P>
499 if (lhs.empty() || rhs.empty())
503 using index_set_t =
typename multivector_t::index_set_t;
504 using term_t =
typename multivector_t::term_t;
506 const auto lhs_end = lhs.end();
507 const auto rhs_end = rhs.end();
508 const double lhs_size = lhs.size();
509 const double rhs_size = rhs.size();
511 const auto lhs_frame = lhs.frame();
512 const auto rhs_frame = rhs.frame();
514 const auto our_frame = lhs_frame | rhs_frame;
515 const auto algebra_dim =
set_value_t(1) << our_frame.count();
516 auto result = multivector_t(
518 if (lhs_size * rhs_size >
double(Tune_P::products_size_threshold))
522 result_stv != algebra_dim;
525 const auto result_ist = index_set_t(result_stv, our_frame,
true);
526 const auto comp_frame = our_frame & ~result_ist;
527 const auto comp_dim =
set_value_t(1) << comp_frame.count();
528 auto result_crd = Scalar_T(0);
531 comp_stv != comp_dim;
534 const auto comp_ist = index_set_t(comp_stv, comp_frame,
true);
535 const auto our_ist = result_ist ^ comp_ist;
536 if ((our_ist | lhs_frame) == lhs_frame)
538 const auto lhs_it = lhs.find(our_ist);
539 if (lhs_it != lhs_end)
541 const auto rhs_it = rhs.find(comp_ist);
542 if (rhs_it != rhs_end)
548 if ((our_ist | rhs_frame) == rhs_frame)
550 const auto rhs_it = rhs.find(our_ist);
551 if (rhs_it != rhs_end)
553 const auto lhs_it = lhs.find(comp_ist);
554 if (lhs_it != lhs_end)
560 if (result_crd != Scalar_T(0))
561 result.insert(term_t(result_ist, result_crd));
566 const auto empty_set = index_set_t();
567 for (
auto& lhs_term : lhs)
569 const auto lhs_ist = lhs_term.first;
570 if (lhs_ist != empty_set)
571 for (
auto& rhs_term : rhs)
573 const auto rhs_ist = rhs_term.first;
574 if (rhs_ist != empty_set)
576 const auto our_ist = lhs_ist | rhs_ist;
577 if ((lhs_ist == our_ist) || (rhs_ist == our_ist))
578 result += lhs_term * rhs_term;
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 >
597 operator% (
const framed_multi<Scalar_T,LO,HI,Tune_P>& lhs,
const framed_multi<Scalar_T,LO,HI,Tune_P>& rhs) ->
const framed_multi<Scalar_T,LO,HI,Tune_P>
604 if (lhs.empty() || rhs.empty())
608 using index_set_t =
typename multivector_t::index_set_t;
609 using term_t =
typename multivector_t::term_t;
610 using map_t =
typename multivector_t::map_t;
612 const auto lhs_end = lhs.end();
613 const auto rhs_end = rhs.end();
614 const double lhs_size = lhs.size();
615 const double rhs_size = rhs.size();
616 const auto lhs_frame = lhs.frame();
617 const auto rhs_frame = rhs.frame();
619 const auto our_frame = lhs_frame | rhs_frame;
620 const auto algebra_dim =
set_value_t(1) << our_frame.count();
621 auto result = multivector_t(
624 if (lhs_size * rhs_size >
double(Tune_P::products_size_threshold))
628 result_stv != algebra_dim;
631 const auto result_ist = index_set_t(result_stv, our_frame,
true);
632 const auto comp_frame = lhs_frame & ~result_ist;
633 const auto comp_dim =
set_value_t(1) << comp_frame.count();
634 auto result_crd = Scalar_T(0);
637 comp_stv != comp_dim;
640 const auto comp_ist = index_set_t(comp_stv, comp_frame,
true);
641 const auto rhs_ist = result_ist ^ comp_ist;
642 if ((rhs_ist | rhs_frame) == rhs_frame)
644 const auto rhs_it = rhs.find(rhs_ist);
645 if (rhs_it != rhs_end)
647 const auto lhs_it = lhs.find(comp_ist);
648 if (lhs_it != lhs_end)
653 if (result_crd != Scalar_T(0))
654 result.insert(term_t(result_ist, result_crd));
659 for (
auto& rhs_term : rhs)
661 const auto rhs_ist = rhs_term.first;
662 for (
auto& lhs_term : lhs)
664 const index_set_t lhs_ist = lhs_term.first;
665 if ((lhs_ist | rhs_ist) == rhs_ist)
666 result += lhs_term * rhs_term;
674 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
678 operator%= (
const multivector_t& rhs) -> multivector_t&
679 {
return *
this = *
this % rhs; }
682 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
688 auto result = Scalar_T(0);
689 const auto small_star_large = lhs.size() < rhs.size();
699 for (
auto& small_term : *smallp)
701 const auto small_ist = small_term.first;
702 const auto large_crd = (*largep)[small_ist];
703 if (large_crd != Scalar_T(0))
704 result += small_ist.sign_of_square() * small_term.second * large_crd;
710 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
713 operator/= (
const Scalar_T& scr) -> multivector_t&
715 using traits_t = numeric_traits<Scalar_T>;
717 if (traits_t::isNaN(scr))
718 return *
this = traits_t::NaN();
719 if (traits_t::isInf(scr))
721 *
this = traits_t::NaN();
725 for (
auto& this_term : *
this)
726 this_term.second /= scr;
731 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
734 operator/ (
const framed_multi<Scalar_T,LO,HI,Tune_P>& lhs,
const framed_multi<Scalar_T,LO,HI,Tune_P>& rhs) ->
const framed_multi<Scalar_T,LO,HI,Tune_P>
738 using index_set_t =
typename multivector_t::index_set_t;
739 using matrix_multi_t =
typename multivector_t::matrix_multi_t;
741 if (rhs == Scalar_T(0))
742 return traits_t::NaN();
744 const auto our_frame = lhs.frame() | rhs.frame();
745 return matrix_multi_t(lhs, our_frame,
true) / matrix_multi_t(rhs, our_frame,
true);
749 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
753 operator/= (
const multivector_t& rhs) -> multivector_t&
754 {
return *
this = *
this / rhs; }
757 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
760 operator| (
const framed_multi<Scalar_T,LO,HI,Tune_P>& lhs,
const framed_multi<Scalar_T,LO,HI,Tune_P>& rhs) ->
const framed_multi<Scalar_T,LO,HI,Tune_P>
763 using matrix_multi_t =
typename multivector_t::matrix_multi_t;
765 return matrix_multi_t(rhs) * matrix_multi_t(lhs) / matrix_multi_t(rhs.involute());
769 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
773 operator|= (
const multivector_t& rhs) -> multivector_t&
774 {
return *
this = *
this | rhs; }
777 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
781 inv() const -> const multivector_t
783 auto result = matrix_multi_t(Scalar_T(1), this->frame());
784 return result /= matrix_multi_t(*
this);
788 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
791 pow(
int m)
const ->
const multivector_t
795 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
798 outer_pow(
int m)
const ->
const multivector_t
801 throw error_t(
"outer_pow(int): negative exponent");
802 auto result = multivector_t(Scalar_T(1));
813 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
817 frame() const -> const index_set_t
819 auto result = index_set_t();
820 for (
auto& this_term : *
this)
821 result |= this_term.first;
826 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
833 for (
auto& this_term : *
this)
834 result = std::max(result, this_term.first.count());
839 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
843 operator[] (
const index_set_t ist)
const -> Scalar_T
845 const auto& this_it = this->find(ist);
846 if (this_it == this->end())
849 return this_it->second;
853 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
858 if ((grade < 0) || (grade > HI-LO))
862 auto result = multivector_t();
863 for (
auto& this_term : *
this)
864 if (this_term.first.count() == grade)
871 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
875 scalar() const -> Scalar_T
876 {
return (*
this)[index_set_t()]; }
879 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
883 pure() const -> const multivector_t
884 {
return *
this - this->
scalar(); }
887 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
890 even() const -> const multivector_t
892 auto result = multivector_t();
893 for (
auto& this_term : *
this)
894 if ((this_term.first.count() % 2) == 0)
895 result.insert(this_term);
900 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
903 odd() const -> const multivector_t
905 auto result = multivector_t();
906 for (
auto& this_term : *
this)
907 if ((this_term.first.count() % 2) == 1)
908 result.insert(this_term);
913 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
920 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
923 vector_part(
const index_set_t frm,
const bool prechecked)
const ->
const vector_t
925 if (!prechecked && (this->frame() | frm) != frm)
926 throw error_t(
"vector_part(frm): value is outside of requested frame");
927 auto result = vector_t();
928 result.reserve(frm.count());
929 const auto frm_end = frm.max()+1;
937 result.push_back((*
this)[index_set_t(idx)]);
942 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
945 involute() const -> const multivector_t
948 for (
auto& result_term : result)
950 if ((result_term.first.count() % 2) == 1)
951 result_term.second *= Scalar_T(-1);
957 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
960 reverse() const -> const multivector_t
963 for (
auto& result_term : result)
966 switch (result_term.first.count() % 4)
970 result_term.second *= Scalar_T(-1);
979 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
982 conj() const -> const multivector_t
985 for (
auto& result_term : result)
988 switch (result_term.first.count() % 4)
992 result_term.second *= Scalar_T(-1);
1001 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1004 quad() const -> Scalar_T
1008 auto result = Scalar_T(0);
1009 for (
auto& this_term : *
this)
1012 (this_term.first.count_neg() % 2)
1015 result += sign * (this_term.second) * (this_term.second);
1021 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1024 norm() const -> Scalar_T
1026 using traits_t = numeric_traits<Scalar_T>;
1028 auto result = Scalar_T(0);
1029 for (
auto& this_term : *
this)
1031 const auto abs_crd = traits_t::abs(this_term.second);
1032 result += abs_crd * abs_crd;
1038 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1043 using traits_t = numeric_traits<Scalar_T>;
1045 auto result = Scalar_T(0);
1046 for (
auto& this_term : *
this)
1048 const auto abs_crd = traits_t::abs(this_term.second);
1049 if (abs_crd > result)
1056 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1066 auto& generator = random_generator_t::generator();
1069 (fill < Scalar_T(0))
1071 : (fill > Scalar_T(1))
1074 const auto algebra_dim =
set_value_t(1) << frm.count();
1076 const auto mean_abs = traits_t::sqrt(Scalar_T(
double(algebra_dim)));
1082 if (generator.uniform() < fill)
1084 const auto& result_crd = generator.normal() / mean_abs;
1091 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1095 write(
const std::string& msg)
const
1096 { std::cout << msg << std::endl <<
" " << (*this) << std::endl; }
1099 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1103 write(std::ofstream& ofile,
const std::string& msg)
const
1106 throw error_t(
"write(ofile,msg): cannot write to output file");
1107 ofile << msg << std::endl <<
" " << (*this) << std::endl;
1111 template<
typename Map_T,
typename Sorted_Map_T >
1121 for (
auto& val_term : val)
1122 sorted_val.insert(val_term);
1130 template<
typename Sorted_Map_T >
1147 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1151 using limits_t = std::numeric_limits<Scalar_T>;
1154 else if (val.isnan())
1155 os << limits_t::quiet_NaN();
1156 else if (val.isinf())
1158 const Scalar_T& inf = limits_t::infinity();
1159 os << (
scalar(val) < 0.0 ? -inf : inf);
1165 Scalar_T truncation;
1166 switch (os.flags() & std::ios::floatfield)
1168 case std::ios_base::scientific:
1169 truncation = Scalar_T(1) / traits_t::pow(Scalar_T(10),
int(os.precision()) + 1);
1171 case std::ios_base::fixed:
1172 truncation = Scalar_T(1) / (traits_t::pow(Scalar_T(10),
int(os.precision())) * val.max_abs());
1174 case std::ios_base::fixed | std::ios_base::scientific:
1175 truncation = multivector_t::default_truncation;
1178 truncation = Scalar_T(1) / traits_t::pow(Scalar_T(10),
int(os.precision()));
1181 auto truncated_val = val.truncated(truncation);
1182 if (truncated_val.empty())
1186 using map_t =
typename multivector_t::map_t;
1187 using sorted_map_t =
typename multivector_t::sorted_map_t;
1188 using sorted_iterator =
typename sorted_map_t::const_iterator;
1189 auto sorted_val = sorted_map_t();
1191 auto sorted_it = sorted_val_range.sorted_begin;
1194 sorted_it != sorted_val_range.sorted_end;
1197 const Scalar_T& scr = sorted_it->second;
1208 template<
typename Scalar_T, const index_t LO, const index_t HI >
1213 const auto use_double =
1214 (os.precision() <= std::numeric_limits<double>::digits10) ||
1215 (term.second == Scalar_T(second_as_double));
1216 if (term.first.count() == 0)
1218 os << second_as_double;
1221 else if (term.second == Scalar_T(-1))
1226 else if (term.second != Scalar_T(1))
1230 auto tol = std::pow(10.0,-os.precision());
1231 if ( std::fabs(second_as_double + 1.0) < tol )
1233 else if ( std::fabs(second_as_double - 1.0) >= tol )
1234 os << second_as_double;
1246 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1252 auto local_val = multivector_t();
1255 auto negative =
false;
1256 auto expect_term =
true;
1259 if (s.good() && (c ==
int(
'+') || c ==
int(
'-')))
1261 negative = (c == int(
'-'));
1270 auto coordinate = Scalar_T(1);
1280 double coordinate_as_double;
1281 s >> coordinate_as_double;
1285 coordinate = Scalar_T(coordinate_as_double);
1297 if (s.good() && c ==
int(
'{'))
1307 expect_term =
false;
1308 if (coordinate != Scalar_T(0))
1315 using term_t =
typename multivector_t::term_t;
1316 local_val += term_t(ist, coordinate);
1325 if (c ==
int(
'+') || c ==
int(
'-'))
1327 negative = (c == int(
'-'));
1344 s.clear(std::istream::failbit);
1354 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1358 {
return this->size(); }
1361 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1365 operator+= (
const term_t& term) -> multivector_t&
1367 if (term.second != Scalar_T(0))
1369 const auto& this_it = this->find(term.first);
1370 if (this_it == this->end())
1372 else if (this_it->second + term.second == Scalar_T(0))
1374 this->erase(this_it);
1376 this_it->second += term.second;
1382 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1385 isinf() const ->
bool
1387 using traits_t = numeric_traits<Scalar_T>;
1389 if (std::numeric_limits<Scalar_T>::has_infinity)
1390 for (
auto& this_term : *
this)
1391 if (traits_t::isInf(this_term.second))
1397 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1400 isnan() const ->
bool
1402 using traits_t = numeric_traits<Scalar_T>;
1404 if (std::numeric_limits<Scalar_T>::has_quiet_NaN)
1405 for (
auto& this_term : *
this)
1406 if (traits_t::isNaN(this_term.second))
1412 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1415 truncated(
const Scalar_T& limit)
const ->
const multivector_t
1417 using traits_t = numeric_traits<Scalar_T>;
1421 const auto truncation = traits_t::abs(limit);
1423 auto result = multivector_t();
1424 if (top != Scalar_T(0))
1425 for (
auto& this_term : *
this)
1426 if (traits_t::abs(this_term.second) > top * truncation)
1427 result.insert(this_term);
1432 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1437 if (frm.is_contiguous())
1442 for (
auto& this_term : *
this)
1443 result.insert(
term_t(this_term.first.fold(frm), this_term.second));
1449 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1454 if (frm.is_contiguous())
1459 for (
auto& this_term : *
this)
1460 result.insert(
term_t(this_term.first.unfold(frm), this_term.second));
1467 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1474 throw error_t(
"centre_pm4_qp4(p,q): LO is too high to represent this value");
1475 if (this->frame().max() > p-4)
1478 const auto pm3210 =
index_set_t(index_pair_t(p-3,p),
true);
1479 const auto qm4321 =
index_set_t(index_pair_t(-q-4,-q-1),
true);
1480 const auto& tqm4321 =
term_t(qm4321, Scalar_T(1));
1482 for (
auto& this_term : *
this)
1484 const auto ist = this_term.first;
1485 if (ist.
max() > p-4)
1495 result.insert(
term_t(ist & ~pm3210, this_term.second) *
1499 result.insert(this_term);
1509 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1516 throw error_t(
"centre_pp4_qm4(p,q): HI is too low to represent this value");
1517 if (this->frame().min() < -q+4)
1520 const auto qp0123 =
index_set_t(index_pair_t(-q,-q+3),
true);
1521 const auto pp1234 =
index_set_t(index_pair_t(p+1,p+4),
true);
1522 const auto& tpp1234 =
term_t(pp1234, Scalar_T(1));
1524 for (
auto& this_term : *
this)
1527 if (ist.
min() < -q+4)
1538 term_t(ist & ~qp0123, this_term.second));
1541 result.insert(this_term);
1551 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1557 throw error_t(
"centre_qp1_pm1(p,q): HI is too low to represent this value");
1559 throw error_t(
"centre_qp1_pm1(p,q): LO is too high to represent this value");
1561 const auto& tqp1 =
term_t(qp1, Scalar_T(1));
1563 for (
auto& this_term : *
this)
1565 const auto ist = this_term.first;
1571 if (n != 0 && ist[n])
1573 if (p != 0 && ist[p])
1580 return *
this = result;
1584 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1591 for (
auto& this_term : *
this)
1592 if ((this_term.first | ist) == this_term.first)
1593 quo.insert(
term_t(this_term.first ^ ist, this_term.second));
1595 rem.insert(this_term);
1600 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1609 const auto dim = matrix_index_t(1) << level;
1615 return matrix::unit<matrix_t>(1) * this->
scalar();
1618 using basis_scalar_t =
typename basis_matrix_t::value_type;
1620 const auto& I = matrix::unit<basis_matrix_t>(2);
1621 auto J = basis_matrix_t(2,2,2);
1623 J(0,1) = basis_scalar_t(-1);
1624 J(1,0) = basis_scalar_t( 1);
1626 K(0,1) = basis_scalar_t( 1);
1628 JK(0,0) = basis_scalar_t(-1);
1641 const auto& pair_mn = this->divide(ist_mn);
1642 const auto& quo_mn = pair_mn.first;
1643 const auto& rem_mn = pair_mn.second;
1644 const auto& pair_quo_mnpn = quo_mn.divide(ist_pn);
1645 const auto& val_mnpn = pair_quo_mnpn.first;
1646 const auto& val_mn = pair_quo_mnpn.second;
1647 const auto& pair_rem_mnpn = rem_mn.divide(ist_pn);
1648 const auto& val_pn = pair_rem_mnpn.first;
1649 const auto& val_1 = pair_rem_mnpn.second;
1652 return - kron(JK, val_1.fast (level-1, 1))
1653 + kron(I, val_mnpn.fast(level-1, 1))
1654 + kron(J, val_mn.fast (level-1, 0))
1655 + kron(K, val_pn.fast (level-1, 0));
1657 return kron(I, val_1.fast (level-1, 0))
1658 + kron(JK, val_mnpn.fast(level-1, 0))
1659 + kron(K, val_mn.fast (level-1, 1))
1660 - kron(J, val_pn.fast (level-1, 1));
1665 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1666 template<
typename Other_Scalar_T >
1672 auto val = this->fold(frm);
1673 auto p = frm.count_pos();
1674 auto q = frm.count_neg();
1676 p += std::max(bott_offset,
index_t(0));
1677 q -= std::min(bott_offset,
index_t(0));
1679 throw error_t(
"fast_matrix_multi(frm): HI is too low to represent this value");
1681 throw error_t(
"fast_matrix_multi(frm): LO is too high to represent this value");
1684 val.centre_pm4_qp4(p, q);
1686 val.centre_pp4_qm4(p, q);
1688 val.centre_qp1_pm1(p, q);
1689 const index_t level = (p + q)/2;
1692 const auto& ev_val = val.even();
1693 const auto& od_val = val.odd();
1697 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1705 template<
typename Scalar_T, const index_t LO, const index_t HI >
1711 {
return lhs.first.sign_of_mult(rhs.first) * lhs.second * rhs.second; }
1714 template<
typename Scalar_T, const index_t LO, const index_t HI >
1718 const std::pair<
const index_set<LO,HI>, Scalar_T>& rhs) ->
const std::pair<const index_set<LO,HI>, Scalar_T>
1720 using term_t = std::pair<const index_set<LO,HI>, Scalar_T>;
1721 return term_t(lhs.first ^ rhs.first,
crd_of_mult(lhs, rhs));
1725 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1727 sqrt(
const framed_multi<Scalar_T,LO,HI,Tune_P>& val,
const framed_multi<Scalar_T,LO,HI,Tune_P>& i,
bool prechecked) ->
const framed_multi<Scalar_T,LO,HI,Tune_P>
1731 return traits_t::NaN();
1735 const auto realval = val.scalar();
1738 if (realval < Scalar_T(0))
1739 return i * traits_t::sqrt(-realval);
1741 return traits_t::sqrt(realval);
1744 return sqrt(matrix_multi_t(val), matrix_multi_t(i), prechecked);
1748 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1754 return traits_t::NaN();
1756 const auto s =
scalar(val);
1758 return traits_t::exp(s);
1760 const double size = val.size();
1761 const auto frm_count = val.frame().count();
1762 const auto algebra_dim =
set_value_t(1) << frm_count;
1764 if( (size * size <=
double(algebra_dim)) || (frm_count < Tune_P::mult_matrix_threshold))
1766 switch (Tune_P::function_precision)
1768 case precision_demoted:
1770 using demoted_scalar_t =
typename traits_t::demoted::type;
1773 const auto& demoted_val = demoted_multivector_t(val);
1777 case precision_promoted:
1779 using promoted_scalar_t =
typename traits_t::promoted::type;
1782 const auto& promoted_val = promoted_multivector_t(val);
1793 return exp(matrix_multi_t(val));
1798 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Tune_P >
1800 log(
const framed_multi<Scalar_T,LO,HI,Tune_P>& val,
const framed_multi<Scalar_T,LO,HI,Tune_P>& i,
bool prechecked) ->
const framed_multi<Scalar_T,LO,HI,Tune_P>
1803 if (val == Scalar_T(0) || val.isnan())
1804 return traits_t::NaN();
1808 const auto realval = val.scalar();
1811 if (realval < Scalar_T(0))
1812 return i * traits_t::pi() + traits_t::log(-realval);
1814 return traits_t::log(realval);
1817 return log(matrix_multi_t(val), matrix_multi_t(i), prechecked);
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.
Specific exception class.
A framed_multi<Scalar_T,LO,HI,Tune_P> is a framed approximation to a multivector.
auto divide(const index_set_t ist) const -> const framed_pair_t
Divide multivector into part divisible by index_set and remainder.
auto fast(const index_t level, const bool odd) const -> const matrix_t
Generalized FFT from multivector_t to matrix_t.
auto operator+=(const term_t &term) -> multivector_t &
Add a term, if non-zero.
friend class framed_multi
framed_multi multivector_t
std::vector< Scalar_T > vector_t
auto fast_matrix_multi(const index_set_t frm) const -> const matrix_multi< Other_Scalar_T, LO, HI, Tune_P >
Use generalized FFT to construct a matrix_multi_t.
index_set< LO, HI > index_set_t
std::pair< const multivector_t, const multivector_t > framed_pair_t
auto centre_pm4_qp4(index_t &p, index_t &q) -> multivector_t &
Subalgebra isomorphism: R_{p,q} to R_{p-4,q+4}.
std::unordered_map< index_set_t, Scalar_T, index_set_hash< LO, HI > > map_t
typename matrix_multi_t::matrix_t matrix_t
auto fold(const index_set_t frm) const -> multivector_t
Subalgebra isomorphism: fold each term within the given frame.
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 fast_framed_multi() const -> const framed_multi_t
Use inverse generalized FFT to construct a framed_multi_t.
auto unfold(const index_set_t frm) const -> multivector_t
Subalgebra isomorphism: unfold each term within the given frame.
static auto classname() -> const std::string
Class name used in messages.
auto centre_pp4_qm4(index_t &p, index_t &q) -> multivector_t &
Subalgebra isomorphism: R_{p,q} to R_{p+4,q-4}.
std::pair< const index_set_t, Scalar_T > term_t
error< multivector_t > error_t
_GLUCAT_CLIFFORD_ALGEBRA_OPERATIONS auto nbr_terms() const -> unsigned long
Number of terms.
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.
std::pair< index_t, index_t > index_pair_t
auto max() const -> index_t
Maximum member.
A matrix_multi<Scalar_T,LO,HI,Tune_P> is a matrix approximation to a multivector.
typename matrix_t::size_type matrix_index_t
auto basis_element(const index_set< LO, HI > &ist) const -> const basis_matrix_t
Create a basis element matrix within the current frame.
matrix_t m_matrix
Matrix value representing the multivector within the folded frame.
ublas::compressed_matrix< int, orientation_t > basis_matrix_t
Extra traits which extend numeric limits.
static auto to_double(const Scalar_T &val) -> double
Cast to double.
Random number generator with single instance per Scalar_T.
sorted_iterator sorted_begin
Sorted_Map_T sorted_map_t
sorted_range(Sorted_Map_T &sorted_val, const Sorted_Map_T &val)
sorted_iterator sorted_end
typename Sorted_Map_T::const_iterator sorted_iterator
Sorted range for use with output.
sorted_iterator sorted_begin
typename Sorted_Map_T::const_iterator sorted_iterator
sorted_iterator sorted_end
sorted_range(Sorted_Map_T &sorted_val, const Map_T &val)
Sorted_Map_T sorted_map_t
#define _GLUCAT_HASH_N(x)
#define _GLUCAT_HASH_SIZE_T(x)
static const std::array< index_t, 8 > offset_to_super
Offsets between the current signature and that of the real superalgebra.
auto isinf(const Matrix_T &m) -> bool
Infinite.
auto nnz(const Matrix_T &m) -> typename Matrix_T::size_type
Number of non-zeros.
auto kron(const LHS_T &lhs, const RHS_T &rhs) -> const RHS_T
Kronecker tensor product of matrices - as per Matlab kron.
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.
int index_t
Size of index_t should be enough to represent LO, HI.
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 crd_of_mult(const std::pair< const index_set< LO, HI >, Scalar_T > &lhs, const std::pair< const index_set< LO, HI >, Scalar_T > &rhs) -> Scalar_T
Coordinate of product of terms.
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.
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 operator>>(std::istream &s, framed_multi< Scalar_T, LO, HI, Tune_P > &val) -> std::istream &
Read multivector from input.
unsigned long set_value_t
Size of set_value_t should be enough to contain index_set<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.
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 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 max_abs(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> Scalar_T
Maximum of absolute values of components of multivector: multivector infinity norm.
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.
auto odd(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Odd part.
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()