1 #ifndef _GLUCAT_FRAMED_MULTI_IMP_H
2 #define _GLUCAT_FRAMED_MULTI_IMP_H
39 #if defined(_GLUCAT_USE_BOOST_POOL_ALLOC)
41 #include <boost/pool/pool_alloc.hpp>
50 template<
typename Scalar_T, const index_t LO, const index_t HI >
54 {
return "framed_multi"; }
56 #if defined(_GLUCAT_MAP_IS_HASH)
57 #define _GLUCAT_HASH_N(x) (x)
58 #define _GLUCAT_HASH_SIZE_T(x) (typename multivector_t::hash_size_t)(x)
60 #define _GLUCAT_HASH_N(x)
61 #define _GLUCAT_HASH_SIZE_T(x)
65 template<
typename Scalar_T, const index_t LO, const index_t HI >
72 template<
typename Scalar_T, const index_t LO, const index_t HI >
73 framed_multi<Scalar_T,LO,HI>::
74 framed_multi(
const hash_size_t& hash_size)
79 template<
typename Scalar_T, const index_t LO, const index_t HI >
80 template<
typename Other_Scalar_T >
81 framed_multi<Scalar_T,LO,HI>::
86 typedef typename other_framed_multi_t::const_iterator other_const_iterator;
87 const other_const_iterator val_begin = val.begin();
88 const other_const_iterator val_end = val.end();
89 for (other_const_iterator val_it = val_begin; val_it != val_end; ++val_it)
94 template<
typename Scalar_T, const index_t LO, const index_t HI >
95 template<
typename Other_Scalar_T >
102 typedef typename other_framed_multi_t::const_iterator other_const_iterator;
103 const other_const_iterator val_begin = val.begin();
104 const other_const_iterator val_end = val.end();
105 for (other_const_iterator val_it = val_begin; val_it != val_end; ++val_it)
106 this->insert(
term_t(val_it->first,
static_cast<Scalar_T
>(val_it->second)));
110 template<
typename Scalar_T, const index_t LO, const index_t HI >
118 template<
typename Scalar_T, const index_t LO, const index_t HI >
123 if (crd != Scalar_T(0))
128 template<
typename Scalar_T, const index_t LO, const index_t HI >
134 if (!prechecked && (
ist | frm) != frm)
135 throw error_t(
"multivector_t(ist,crd,frm): cannot initialize with value outside of frame");
136 if (crd != Scalar_T(0))
141 template<
typename Scalar_T, const index_t LO, const index_t HI >
146 if (scr != Scalar_T(0))
151 template<
typename Scalar_T, const index_t LO, const index_t HI >
156 if (scr != Scalar_T(0))
161 template<
typename Scalar_T, const index_t LO, const index_t HI >
168 throw error_t(
"multivector_t(vec,frm): cannot initialize with vector not matching frame");
169 typename vector_t::const_iterator vec_it = vec.begin();
184 template<
typename Scalar_T, const index_t LO, const index_t HI >
189 std::istringstream ss(str);
192 throw error_t(
"multivector_t(str): could not parse string");
196 throw error_t(
"multivector_t(str): could not parse entire string");
200 template<
typename Scalar_T, const index_t LO, const index_t HI >
202 framed_multi(
const std::string& str,
const index_set_t frm,
const bool prechecked)
212 template<
typename Scalar_T, const index_t LO, const index_t HI >
213 template<
typename Other_Scalar_T >
218 if (val == Other_Scalar_T(0))
222 const matrix_index_t dim = val.
m_matrix.size1();
232 *
this = val.template fast_framed_multi<Scalar_T>();
240 const Scalar_T val_norm = traits_t::to_scalar_t(val.norm());
241 if (traits_t::isNaN_or_isInf(val_norm))
243 *
this = traits_t::NaN();
248 #if defined(_GLUCAT_MAP_IS_HASH)
261 if ((abs_crd * abs_crd) > tol)
267 template<
typename Scalar_T, const index_t LO, const index_t HI >
269 framed_multi<Scalar_T,LO,HI>::
270 operator== (
const multivector_t& rhs)
const
272 if (this->size() != rhs.size())
274 const const_iterator this_begin = this->begin();
275 const const_iterator this_end = this->end();
276 const const_iterator rhs_end = rhs.end();
277 #if defined(_GLUCAT_MAP_IS_ORDERED)
278 const_iterator this_it = this_begin;
279 const_iterator rhs_it = rhs.begin();
281 (this_it != this_end) && (rhs_it != rhs_end);
283 if (*this_it != *rhs_it)
285 return (this_it == this_end) && (rhs_it == rhs_end);
288 this_it = this_begin;
292 const const_iterator& rhs_it = rhs.find(this_it->first);
293 if (rhs_it == rhs_end)
295 else if (rhs_it->second != this_it->second)
303 template<
typename Scalar_T, const index_t LO, const index_t HI >
306 framed_multi<Scalar_T,LO,HI>::
307 operator== (
const Scalar_T& scr)
const
309 switch (this->size())
312 return scr == Scalar_T(0);
315 const const_iterator& this_it = this->begin();
316 return this_it->first == index_set_t() && this_it->second == scr;
325 template<
typename Scalar_T, const index_t LO, const index_t HI >
327 framed_multi<Scalar_T,LO,HI> &
331 *
this += term_t(index_set_t(),scr);
336 template<
typename Scalar_T, const index_t LO, const index_t HI >
338 framed_multi<Scalar_T,LO,HI> &
343 rhs_it = rhs.begin();
351 template<
typename Scalar_T, const index_t LO, const index_t HI >
353 framed_multi<Scalar_T,LO,HI> &
354 framed_multi<Scalar_T,LO,HI>::
355 operator-= (
const multivector_t& rhs)
358 rhs_it = rhs.begin();
361 *
this +=
term_t(rhs_it->first, -(rhs_it->second));
366 template<
typename Scalar_T, const index_t LO, const index_t HI >
371 {
return *
this * Scalar_T(-1); }
374 template<
typename Scalar_T, const index_t LO, const index_t HI >
375 framed_multi<Scalar_T,LO,HI> &
376 framed_multi<Scalar_T,LO,HI>::
377 operator*= (
const Scalar_T& scr)
379 typedef numeric_traits<Scalar_T> traits_t;
381 if (traits_t::isNaN_or_isInf(scr))
382 return *
this = traits_t::NaN();
383 if (scr == Scalar_T(0))
385 *
this = traits_t::NaN();
390 this_it = this->begin();
391 this_it != this->end();
393 this_it->second *= scr;
398 template<
typename Scalar_T, const index_t LO, const index_t HI >
399 const framed_multi<Scalar_T,LO,HI>
400 operator* (
const framed_multi<Scalar_T,LO,HI>& lhs,
const framed_multi<Scalar_T,LO,HI>& rhs)
402 typedef framed_multi<Scalar_T,LO,HI> multivector_t;
403 typedef numeric_traits<Scalar_T> traits_t;
404 typedef typename multivector_t::index_set_t index_set_t;
405 typedef typename multivector_t::term_t term_t;
406 typedef typename multivector_t::map_t map_t;
407 typedef typename map_t::const_iterator const_iterator;
409 if (lhs.isnan() || rhs.isnan())
410 return traits_t::NaN();
412 const double lhs_size = lhs.size();
413 const double rhs_size = rhs.size();
414 const index_set_t our_frame = lhs.frame() | rhs.frame();
415 const index_t frm_count = our_frame.count();
417 const bool direct_mult = lhs_size * rhs_size <= double(algebra_dim)
418 #if defined(_GLUCAT_MAP_IS_HASH)
424 multivector_t result =
425 multivector_t(
_GLUCAT_HASH_SIZE_T(
size_t(std::min(lhs_size * rhs_size,
double(algebra_dim)))));
426 const const_iterator lhs_begin = lhs.begin();
427 const const_iterator lhs_end = lhs.end();
428 const const_iterator rhs_begin = rhs.begin();
429 const const_iterator rhs_end = rhs.end();
436 const term_t& lhs_term = *lhs_it;
441 result += lhs_term * *rhs_it;
445 #if !defined(_GLUCAT_MAP_IS_HASH)
448 typedef std::vector<Scalar_T> array_t;
449 array_t result_array(algebra_dim, Scalar_T(0));
451 const const_iterator lhs_begin = lhs.begin();
452 const const_iterator lhs_end = lhs.end();
453 const const_iterator rhs_begin = rhs.begin();
454 const const_iterator rhs_end = rhs.end();
461 const term_t& lhs_term = *lhs_it;
467 const term_t& term = lhs_term * *rhs_it;
468 const set_value_t stv = term.first.value_of_fold(our_frame);
469 result_array[stv] += term.second;
472 multivector_t result;
477 if (result_array[stv] != Scalar_T(0))
478 result.insert(term_t(index_set_t(stv, our_frame,
true), result_array[stv]));
484 typedef typename multivector_t::matrix_multi_t matrix_multi_t;
485 return matrix_multi_t(lhs, our_frame,
true) *
486 matrix_multi_t(rhs, our_frame,
true);
491 template<
typename Scalar_T, const index_t LO, const index_t HI >
493 framed_multi<Scalar_T,LO,HI> &
494 framed_multi<Scalar_T,LO,HI>::
495 operator*= (
const multivector_t& rhs)
496 {
return *
this = *
this * rhs; }
499 template<
typename Scalar_T, const index_t LO, const index_t HI >
500 const framed_multi<Scalar_T,LO,HI>
501 operator^ (
const framed_multi<Scalar_T,LO,HI>& lhs,
const framed_multi<Scalar_T,LO,HI>& rhs)
505 if (lhs.empty() || rhs.empty())
508 typedef framed_multi<Scalar_T,LO,HI> multivector_t;
509 typedef typename multivector_t::index_set_t index_set_t;
510 typedef typename multivector_t::term_t term_t;
511 typedef typename multivector_t::map_t map_t;
512 typedef typename map_t::const_iterator const_iterator;
514 multivector_t result;
515 const index_set_t empty_set = index_set_t();
517 const const_iterator lhs_begin = lhs.begin();
518 const const_iterator lhs_end = lhs.end();
519 const const_iterator rhs_begin = rhs.begin();
520 const const_iterator rhs_end = rhs.end();
521 const double lhs_size = lhs.size();
522 const double rhs_size = rhs.size();
526 const index_set_t lhs_frame = lhs.frame();
527 const index_set_t rhs_frame = rhs.frame();
529 const index_set_t our_frame = lhs_frame | rhs_frame;
530 const set_value_t algebra_dim = 1 << our_frame.count();
531 multivector_t result =
532 multivector_t(
_GLUCAT_HASH_SIZE_T(
size_t(std::min(lhs_size * rhs_size,
double(algebra_dim)))));
535 result_stv != algebra_dim;
538 const index_set_t result_ist = index_set_t(result_stv, our_frame,
true);
539 const index_set_t lhs_result_frame = lhs_frame & result_ist;
540 const set_value_t lhs_result_dim = 1 << lhs_result_frame.count();
541 Scalar_T result_crd = Scalar_T(0);
544 lhs_stv != lhs_result_dim;
547 const index_set_t lhs_ist = index_set_t(lhs_stv, lhs_result_frame,
true);
548 const index_set_t rhs_ist = result_ist ^ lhs_ist;
549 if ((rhs_ist | rhs_frame) == rhs_frame)
551 const const_iterator lhs_it = lhs.find(lhs_ist);
552 if (lhs_it != lhs_end)
554 const const_iterator rhs_it = rhs.find(rhs_ist);
555 if (rhs_it != rhs_end)
560 if (result_crd != Scalar_T(0))
561 result.insert(term_t(result_ist, result_crd));
567 multivector_t result;
573 const term_t& rhs_term = *rhs_it;
574 const index_set_t rhs_ist = rhs_term.first;
580 const term_t& lhs_term = *lhs_it;
581 const index_set_t lhs_ist = lhs_term.first;
582 if ((lhs_ist & rhs_ist) == empty_set)
583 result += lhs_term * rhs_term;
591 template<
typename Scalar_T, const index_t LO, const index_t HI >
593 framed_multi<Scalar_T,LO,HI> &
594 framed_multi<Scalar_T,LO,HI>::
595 operator^= (
const multivector_t& rhs)
596 {
return *
this = *
this ^ rhs; }
599 template<
typename Scalar_T, const index_t LO, const index_t HI >
600 const framed_multi<Scalar_T,LO,HI>
601 operator& (
const framed_multi<Scalar_T,LO,HI>& lhs,
const framed_multi<Scalar_T,LO,HI>& rhs)
605 if (lhs.empty() || rhs.empty())
608 typedef framed_multi<Scalar_T,LO,HI> multivector_t;
609 typedef typename multivector_t::index_set_t index_set_t;
610 typedef typename multivector_t::term_t term_t;
611 typedef typename multivector_t::map_t map_t;
612 typedef typename map_t::const_iterator const_iterator;
614 const const_iterator lhs_end = lhs.end();
615 const const_iterator rhs_end = rhs.end();
616 const double lhs_size = lhs.size();
617 const double rhs_size = rhs.size();
621 const index_set_t lhs_frame = lhs.frame();
622 const index_set_t rhs_frame = rhs.frame();
624 const index_set_t our_frame = lhs_frame | rhs_frame;
625 const set_value_t algebra_dim = 1 << our_frame.count();
626 multivector_t result =
627 multivector_t(
_GLUCAT_HASH_SIZE_T(
size_t(std::min(lhs_size * rhs_size,
double(algebra_dim)))));
630 result_stv != algebra_dim;
633 const index_set_t result_ist = index_set_t(result_stv, our_frame,
true);
634 const index_set_t comp_frame = our_frame & ~result_ist;
635 const set_value_t comp_dim = 1 << comp_frame.count();
636 Scalar_T result_crd = Scalar_T(0);
639 comp_stv != comp_dim;
642 const index_set_t comp_ist = index_set_t(comp_stv, comp_frame,
true);
643 const index_set_t our_ist = result_ist ^ comp_ist;
644 if ((our_ist | lhs_frame) == lhs_frame)
646 const const_iterator lhs_it = lhs.find(our_ist);
647 if (lhs_it != lhs_end)
649 const const_iterator rhs_it = rhs.find(comp_ist);
650 if (rhs_it != rhs_end)
656 if ((our_ist | rhs_frame) == rhs_frame)
658 const const_iterator rhs_it = rhs.find(our_ist);
659 if (rhs_it != rhs_end)
661 const const_iterator lhs_it = lhs.find(comp_ist);
662 if (lhs_it != lhs_end)
668 if (result_crd != Scalar_T(0))
669 result.insert(term_t(result_ist, result_crd));
675 const index_set_t empty_set = index_set_t();
677 const const_iterator lhs_begin = lhs.begin();
678 const const_iterator rhs_begin = rhs.begin();
680 multivector_t result;
686 const term_t& lhs_term = *lhs_it;
687 const index_set_t lhs_ist = lhs_term.first;
688 if (lhs_ist != empty_set)
694 const term_t& rhs_term = *rhs_it;
695 const index_set_t rhs_ist = rhs_term.first;
696 if (rhs_ist != empty_set)
698 const index_set_t our_ist = lhs_ist | rhs_ist;
699 if ((lhs_ist == our_ist) || (rhs_ist == our_ist))
700 result += lhs_term * rhs_term;
709 template<
typename Scalar_T, const index_t LO, const index_t HI >
711 framed_multi<Scalar_T,LO,HI> &
712 framed_multi<Scalar_T,LO,HI>::
713 operator&= (
const multivector_t& rhs)
714 {
return *
this = *
this & rhs; }
717 template<
typename Scalar_T, const index_t LO, const index_t HI >
718 const framed_multi<Scalar_T,LO,HI>
719 operator% (
const framed_multi<Scalar_T,LO,HI>& lhs,
const framed_multi<Scalar_T,LO,HI>& rhs)
726 if (lhs.empty() || rhs.empty())
729 typedef framed_multi<Scalar_T,LO,HI> multivector_t;
730 typedef typename multivector_t::index_set_t index_set_t;
731 typedef typename multivector_t::term_t term_t;
732 typedef typename multivector_t::map_t map_t;
733 typedef typename map_t::const_iterator const_iterator;
735 #if defined(_GLUCAT_MAP_IS_ORDERED)
739 const const_iterator lhs_begin = lhs.begin();
741 typedef typename map_t::const_reverse_iterator const_reverse_iterator;
742 const const_reverse_iterator rhs_rbegin = rhs.rbegin();
743 const const_reverse_iterator rhs_rlower_bound =
744 static_cast<const_reverse_iterator
>(rhs.lower_bound(lhs_begin->first));
746 multivector_t result;
748 for (const_reverse_iterator
750 rhs_it != rhs_rlower_bound;
753 const term_t& rhs_term = *rhs_it;
754 const index_set_t rhs_ist = rhs_term.first;
755 const const_iterator lhs_upper_bound = lhs.upper_bound(rhs_ist);
758 lhs_it != lhs_upper_bound;
761 const term_t& lhs_term = *lhs_it;
762 const index_set_t lhs_ist = lhs_term.first;
763 if ((lhs_ist | rhs_ist) == rhs_ist)
764 result += lhs_term * rhs_term;
769 const const_iterator lhs_end = lhs.end();
770 const const_iterator rhs_end = rhs.end();
771 const double lhs_size = lhs.size();
772 const double rhs_size = rhs.size();
776 const index_set_t lhs_frame = lhs.frame();
777 const index_set_t rhs_frame = rhs.frame();
779 const index_set_t our_frame = lhs_frame | rhs_frame;
780 const set_value_t algebra_dim = 1 << our_frame.count();
781 multivector_t result =
782 multivector_t(
_GLUCAT_HASH_SIZE_T(
size_t(std::min(lhs_size * rhs_size,
double(algebra_dim)))));
785 result_stv != algebra_dim;
788 const index_set_t result_ist = index_set_t(result_stv, our_frame,
true);
789 const index_set_t comp_frame = lhs_frame & ~result_ist;
790 const set_value_t comp_dim = 1 << comp_frame.count();
791 Scalar_T result_crd = Scalar_T(0);
794 comp_stv != comp_dim;
797 const index_set_t comp_ist = index_set_t(comp_stv, comp_frame,
true);
798 const index_set_t rhs_ist = result_ist ^ comp_ist;
799 if ((rhs_ist | rhs_frame) == rhs_frame)
801 const const_iterator rhs_it = rhs.find(rhs_ist);
802 if (rhs_it != rhs_end)
804 const const_iterator lhs_it = lhs.find(comp_ist);
805 if (lhs_it != lhs_end)
810 if (result_crd != Scalar_T(0))
811 result.insert(term_t(result_ist, result_crd));
817 const const_iterator rhs_begin = rhs.begin();
818 const const_iterator lhs_begin = lhs.begin();
820 multivector_t result;
826 const term_t& rhs_term = *rhs_it;
827 const index_set_t rhs_ist = rhs_term.first;
833 const term_t& lhs_term = *lhs_it;
834 const index_set_t lhs_ist = lhs_term.first;
835 if ((lhs_ist | rhs_ist) == rhs_ist)
836 result += lhs_term * rhs_term;
845 template<
typename Scalar_T, const index_t LO, const index_t HI >
847 framed_multi<Scalar_T,LO,HI> &
848 framed_multi<Scalar_T,LO,HI>::
849 operator%= (
const multivector_t& rhs)
850 {
return *
this = *
this % rhs; }
853 template<
typename Scalar_T, const index_t LO, const index_t HI >
855 star(
const framed_multi<Scalar_T,LO,HI>& lhs,
const framed_multi<Scalar_T,LO,HI>& rhs)
857 typedef framed_multi<Scalar_T,LO,HI> multivector_t;
858 typedef typename multivector_t::map_t map_t;
859 typedef typename map_t::const_iterator const_iterator;
860 typedef typename multivector_t::index_set_t index_set_t;
862 Scalar_T result = Scalar_T(0);
863 const bool small_star_large = lhs.size() < rhs.size();
864 const multivector_t* smallp =
868 const multivector_t* largep =
874 small_it = smallp->begin();
875 small_it != smallp->end();
878 const index_set_t small_ist = small_it->first;
879 const Scalar_T large_crd = (*largep)[small_ist];
880 if (large_crd != Scalar_T(0))
881 result += small_ist.sign_of_square() * small_it->second * large_crd;
887 template<
typename Scalar_T, const index_t LO, const index_t HI >
894 if (traits_t::isNaN(scr))
895 return *
this = traits_t::NaN();
896 if (traits_t::isInf(scr))
898 *
this = traits_t::NaN();
903 this_it = this->begin();
904 this_it != this->end();
906 this_it->second /= scr;
911 template<
typename Scalar_T, const index_t LO, const index_t HI >
913 const framed_multi<Scalar_T,LO,HI>
914 operator/ (
const framed_multi<Scalar_T,LO,HI>& lhs,
const framed_multi<Scalar_T,LO,HI>& rhs)
916 typedef framed_multi<Scalar_T,LO,HI> multivector_t;
917 typedef numeric_traits<Scalar_T> traits_t;
918 typedef typename multivector_t::index_set_t index_set_t;
919 typedef typename multivector_t::matrix_multi_t matrix_multi_t;
921 if (rhs == Scalar_T(0))
922 return traits_t::NaN();
924 const index_set_t our_frame = lhs.frame() | rhs.frame();
925 return matrix_multi_t(lhs, our_frame,
true) / matrix_multi_t(rhs, our_frame,
true);
929 template<
typename Scalar_T, const index_t LO, const index_t HI >
931 framed_multi<Scalar_T,LO,HI> &
932 framed_multi<Scalar_T,LO,HI>::
933 operator/= (
const multivector_t& rhs)
934 {
return *
this = *
this / rhs; }
937 template<
typename Scalar_T, const index_t LO, const index_t HI >
939 const framed_multi<Scalar_T,LO,HI>
940 operator| (
const framed_multi<Scalar_T,LO,HI>& lhs,
const framed_multi<Scalar_T,LO,HI>& rhs)
942 typedef framed_multi<Scalar_T,LO,HI> multivector_t;
943 typedef typename multivector_t::matrix_multi_t matrix_multi_t;
945 return matrix_multi_t(rhs) * matrix_multi_t(lhs) / matrix_multi_t(rhs.involute());
949 template<
typename Scalar_T, const index_t LO, const index_t HI >
951 framed_multi<Scalar_T,LO,HI>&
952 framed_multi<Scalar_T,LO,HI>::
953 operator|= (
const multivector_t& rhs)
954 {
return *
this = *
this | rhs; }
957 template<
typename Scalar_T, const index_t LO, const index_t HI >
959 const framed_multi<Scalar_T,LO,HI>
963 matrix_multi_t result = matrix_multi_t(Scalar_T(1), this->frame());
964 return result /= matrix_multi_t(*
this);
968 template<
typename Scalar_T, const index_t LO, const index_t HI >
969 const framed_multi<Scalar_T,LO,HI>
975 template<
typename Scalar_T, const index_t LO, const index_t HI >
977 framed_multi<Scalar_T,LO,HI>
982 throw error_t(
"outer_pow(int): negative exponent");
983 multivector_t result = Scalar_T(1);
984 multivector_t a = *
this;
994 template<
typename Scalar_T, const index_t LO, const index_t HI >
996 const index_set<LO,HI>
997 framed_multi<Scalar_T,LO,HI>::
1002 this_it = this->begin();
1003 this_it != this->end();
1005 result |= this_it->first;
1010 template<
typename Scalar_T, const index_t LO, const index_t HI >
1013 framed_multi<Scalar_T,LO,HI>::
1018 this_it = this->begin();
1019 this_it != this->end();
1021 result = std::max( result, this_it->first.count() );
1026 template<
typename Scalar_T, const index_t LO, const index_t HI >
1029 framed_multi<Scalar_T,LO,HI>::
1030 operator[] (
const index_set_t
ist)
const
1032 const const_iterator& this_it = this->find(
ist);
1033 if (this_it == this->end())
1036 return this_it->second;
1040 template<
typename Scalar_T, const index_t LO, const index_t HI >
1041 const framed_multi<Scalar_T,LO,HI>
1042 framed_multi<Scalar_T,LO,HI>::
1043 operator() (
index_t grade)
const
1045 if ((grade < 0) || (grade > HI-LO))
1049 multivector_t result;
1051 this_it = this->begin();
1052 this_it != this->end();
1054 if (this_it->first.count() == grade)
1061 template<
typename Scalar_T, const index_t LO, const index_t HI >
1066 {
return (*
this)[index_set_t()]; }
1069 template<
typename Scalar_T, const index_t LO, const index_t HI >
1071 const framed_multi<Scalar_T,LO,HI>
1074 {
return *
this - this->
scalar(); }
1077 template<
typename Scalar_T, const index_t LO, const index_t HI >
1078 const framed_multi<Scalar_T,LO,HI>
1082 multivector_t result;
1084 this_it = this->begin();
1085 this_it != this->end();
1087 if ((this_it->first.count() % 2) == 0)
1088 result.insert(*this_it);
1093 template<
typename Scalar_T, const index_t LO, const index_t HI >
1094 const framed_multi<Scalar_T,LO,HI>
1098 multivector_t result;
1100 this_it = this->begin();
1101 this_it != this->end();
1103 if ((this_it->first.count() % 2) == 1)
1104 result.insert(*this_it);
1109 template<
typename Scalar_T, const index_t LO, const index_t HI >
1113 {
return this->
vector_part(this->frame(),
true); }
1116 template<
typename Scalar_T, const index_t LO, const index_t HI >
1119 vector_part(
const index_set_t frm,
const bool prechecked)
const
1121 if (!prechecked && (this->frame() | frm) != frm)
1122 throw error_t(
"vector_part(frm): value is outside of requested frame");
1124 result.reserve(frm.count());
1125 const index_t frm_end = frm.max()+1;
1133 result.push_back((*
this)[index_set_t(idx)]);
1138 template<
typename Scalar_T, const index_t LO, const index_t HI >
1139 const framed_multi<Scalar_T,LO,HI>
1143 multivector_t result = *
this;
1145 result_it = result.begin();
1146 result_it != result.end();
1149 if ((result_it->first.count() % 2) == 1)
1150 result_it->second *= Scalar_T(-1);
1156 template<
typename Scalar_T, const index_t LO, const index_t HI >
1157 const framed_multi<Scalar_T,LO,HI>
1161 multivector_t result = *
this;
1163 result_it = result.begin();
1164 result_it != result.end();
1168 switch (result_it->first.count() % 4)
1172 result_it->second *= Scalar_T(-1);
1181 template<
typename Scalar_T, const index_t LO, const index_t HI >
1182 const framed_multi<Scalar_T,LO,HI>
1186 multivector_t result = *
this;
1188 result_it = result.begin();
1189 result_it != result.end();
1193 switch (result_it->first.count() % 4)
1197 result_it->second *= Scalar_T(-1);
1206 template<
typename Scalar_T, const index_t LO, const index_t HI >
1213 Scalar_T result = Scalar_T(0);
1215 this_it = this->begin();
1216 this_it != this->end();
1219 const Scalar_T sign =
1220 (this_it->first.count_neg() % 2)
1223 result += sign * (this_it->second) * (this_it->second);
1229 template<
typename Scalar_T, const index_t LO, const index_t HI >
1234 typedef numeric_traits<Scalar_T> traits_t;
1236 Scalar_T result = Scalar_T(0);
1238 this_it = this->begin();
1239 this_it != this->end();
1243 result += abs_crd * abs_crd;
1249 template<
typename Scalar_T, const index_t LO, const index_t HI >
1254 typedef numeric_traits<Scalar_T> traits_t;
1256 Scalar_T result = Scalar_T(0);
1258 this_it = this->begin();
1259 this_it != this->end();
1263 if (abs_crd > result)
1270 template<
typename Scalar_T, const index_t LO, const index_t HI >
1271 const framed_multi<Scalar_T,LO,HI>
1273 random(
const index_set<LO,HI> frm, Scalar_T
fill)
1275 typedef framed_multi<Scalar_T,LO,HI> multivector_t;
1276 typedef typename multivector_t::index_set_t index_set_t;
1277 typedef typename multivector_t::term_t term_t;
1279 typedef random_generator<Scalar_T> random_generator_t;
1280 random_generator_t& generator = random_generator_t::generator();
1283 (
fill < Scalar_T(0))
1285 : (
fill > Scalar_T(1))
1289 typedef numeric_traits<Scalar_T> traits_t;
1290 const Scalar_T mean_abs =
traits_t::sqrt(Scalar_T(
double(algebra_dim)));
1291 multivector_t result;
1296 if (generator.uniform() <
fill)
1298 const Scalar_T& result_crd = generator.normal() / mean_abs;
1299 result.insert(term_t(index_set_t(stv, frm,
true), result_crd));
1305 template<
typename Scalar_T, const index_t LO, const index_t HI >
1309 write(
const std::string& msg)
const
1310 { std::cout << msg << std::endl <<
" " << (*this) << std::endl; }
1313 template<
typename Scalar_T, const index_t LO, const index_t HI >
1316 framed_multi<Scalar_T,LO,HI>::
1317 write(std::ofstream& ofile,
const std::string& msg)
const
1320 throw error_t(
"write(ofile,msg): cannot write to output file");
1321 ofile << msg << std::endl <<
" " << (*this) << std::endl;
1325 template<
typename Map_T,
typename Sorted_Map_T >
1329 typedef Map_T
map_t;
1333 sorted_range (Sorted_Map_T &sorted_val,
const Map_T& val)
1335 for (
typename map_t::const_iterator
1336 val_it = val.begin();
1337 val_it != val.end();
1339 sorted_val.insert(*val_it);
1347 template<
typename Sorted_Map_T >
1348 class sorted_range< Sorted_Map_T, Sorted_Map_T >
1351 typedef Sorted_Map_T
map_t;
1355 sorted_range (Sorted_Map_T &sorted_val,
const Sorted_Map_T& val)
1364 template<
typename Scalar_T, const index_t LO, const index_t HI >
1370 else if (val.isnan())
1371 os << std::numeric_limits<Scalar_T>::quiet_NaN();
1375 typedef typename multivector_t::map_t
map_t;
1376 typedef typename multivector_t::sorted_map_t
sorted_map_t;
1383 sorted_it != sorted_val_range.sorted_end;
1386 const Scalar_T& scr = sorted_it->second;
1396 template<
typename Scalar_T, const index_t LO, const index_t HI >
1398 operator<< (std::ostream& os,
const std::pair<
const index_set<LO,HI>, Scalar_T >& term)
1401 const bool use_double =
1402 (os.precision() <= std::numeric_limits<double>::digits10) ||
1403 (term.second == Scalar_T(second_as_double));
1404 if (term.first.count() == 0)
1406 os << second_as_double;
1409 else if (term.second == Scalar_T(-1))
1414 else if (term.second != Scalar_T(1))
1418 double tol =
std::pow(10.0,-os.precision());
1419 if ( std::fabs(second_as_double + 1.0) < tol )
1421 else if ( std::fabs(second_as_double - 1.0) >= tol )
1422 os << second_as_double;
1434 template<
typename Scalar_T, const index_t LO, const index_t HI >
1436 operator>> (std::istream& s, framed_multi<Scalar_T,LO,HI> & val)
1438 typedef framed_multi<Scalar_T,LO,HI> multivector_t;
1440 multivector_t local_val;
1443 bool negative =
false;
1444 bool expect_term =
true;
1447 if (s.good() && (c ==
int(
'+') || c ==
int(
'-')))
1449 negative = (c == int(
'-'));
1458 Scalar_T coordinate = Scalar_T(1);
1460 index_set<LO,HI>
ist;
1468 double coordinate_as_double;
1469 s >> coordinate_as_double;
1473 coordinate = Scalar_T(coordinate_as_double);
1485 if (s.good() && c ==
int(
'{'))
1495 expect_term =
false;
1496 if (coordinate != Scalar_T(0))
1503 typedef typename multivector_t::term_t term_t;
1504 local_val += term_t(
ist, coordinate);
1513 if (c ==
int(
'+') || c ==
int(
'-'))
1515 negative = (c == int(
'-'));
1532 s.clear(std::istream::failbit);
1542 template<
typename Scalar_T, const index_t LO, const index_t HI >
1546 {
return this->size(); }
1549 template<
typename Scalar_T, const index_t LO, const index_t HI >
1551 framed_multi<Scalar_T,LO,HI> &
1555 if (term.second != Scalar_T(0))
1557 const iterator& this_it = this->find(term.first);
1558 if (this_it == this->end())
1560 else if (this_it->second + term.second == Scalar_T(0))
1562 this->erase(this_it);
1564 this_it->second += term.second;
1570 template<
typename Scalar_T, const index_t LO, const index_t HI >
1577 if (std::numeric_limits<Scalar_T>::has_quiet_NaN)
1579 this_it = this->begin();
1580 this_it != this->end();
1582 if (traits_t::isNaN(this_it->second))
1588 template<
typename Scalar_T, const index_t LO, const index_t HI >
1599 multivector_t result;
1600 if (top != Scalar_T(0))
1602 this_it = this->begin();
1603 this_it != this->end();
1606 result.insert(term_t(this_it->first, this_it->second));
1611 template<
typename Scalar_T, const index_t LO, const index_t HI >
1612 framed_multi<Scalar_T,LO,HI>
1614 fold(
const index_set_t frm)
const
1616 if (frm.is_contiguous())
1620 multivector_t result;
1622 this_it = this->begin();
1623 this_it != this->end();
1625 result.insert(term_t(this_it->first.fold(frm), this_it->second));
1631 template<
typename Scalar_T, const index_t LO, const index_t HI >
1632 framed_multi<Scalar_T,LO,HI>
1634 unfold(
const index_set_t frm)
const
1636 if (frm.is_contiguous())
1640 multivector_t result;
1642 this_it = this->begin();
1643 this_it != this->end();
1645 result.insert(
term_t(this_it->first.unfold(frm,
true), this_it->second));
1652 template<
typename Scalar_T, const index_t LO, const index_t HI >
1653 framed_multi<Scalar_T,LO,HI>&
1659 throw error_t(
"centre_pm4_qp4(p,q): LO is too high to represent this value");
1660 if (this->frame().max() > p-4)
1662 typedef typename index_set_t::index_pair_t index_pair_t;
1664 const index_set_t qm4321(index_pair_t(-q-4,-q-1),
true);
1668 this_it = this->begin();
1669 this_it != this->end();
1673 if (
ist.max() > p-4)
1683 result.insert(
term_t(
ist & ~pm3210, this_it->second) *
1687 result.insert(*this_it);
1697 template<
typename Scalar_T, const index_t LO, const index_t HI >
1698 framed_multi<Scalar_T,LO,HI>&
1704 throw error_t(
"centre_pp4_qm4(p,q): HI is too low to represent this value");
1705 if (this->frame().min() < -q+4)
1707 typedef typename index_set_t::index_pair_t index_pair_t;
1708 const index_set_t qp0123(index_pair_t(-q,-q+3),
true);
1709 const index_set_t pp1234(index_pair_t(p+1,p+4),
true);
1710 const term_t& tpp1234 = term_t(pp1234, Scalar_T(1));
1711 multivector_t result;
1713 this_it = this->begin();
1714 this_it != this->end();
1717 index_set_t
ist = this_it->first;
1718 if (
ist.min() < -q+4)
1726 term *= term_t(index_set_t(n+p+1), Scalar_T(1)) * tpp1234;
1728 result.insert(term_t(term.first, term.second) *
1732 result.insert(*this_it);
1742 template<
typename Scalar_T, const index_t LO, const index_t HI >
1743 framed_multi<Scalar_T,LO,HI>&
1748 throw error_t(
"centre_qp1_pm1(p,q): HI is too low to represent this value");
1750 throw error_t(
"centre_qp1_pm1(p,q): LO is too high to represent this value");
1751 const index_set_t qp1 = index_set_t(q+1);
1752 const term_t& tqp1 = term_t(qp1, Scalar_T(1));
1753 multivector_t result;
1755 this_it = this->begin();
1756 this_it != this->end();
1759 const index_set_t
ist = this_it->first;
1760 var_term_t term = var_term_t(index_set_t(), this_it->second);
1765 if (n != 0 &&
ist[n])
1766 term *= term_t(index_set_t(-n) | qp1, Scalar_T(1));
1767 if (p != 0 &&
ist[p])
1769 result.insert(term_t(term.first, term.second));
1774 return *
this = result;
1778 template<
typename Scalar_T, const index_t LO, const index_t HI >
1786 this_it = this->begin();
1787 this_it != this->end();
1789 if ((this_it->first |
ist) == this_it->first)
1790 quo.insert(term_t(this_it->first ^
ist, this_it->second));
1792 rem.insert(*this_it);
1793 return framed_pair_t(quo, rem);
1797 template<
typename Scalar_T, const index_t LO, const index_t HI >
1805 typedef typename matrix_multi_t::matrix_index_t matrix_index_t;
1806 const matrix_index_t dim = 1 << level;
1807 matrix_t result(dim, dim);
1812 return matrix::unit<matrix_t>(1) * this->
scalar();
1814 typedef typename matrix_multi_t::basis_matrix_t basis_matrix_t;
1815 typedef typename basis_matrix_t::value_type basis_scalar_t;
1817 const basis_matrix_t& I = matrix::unit<basis_matrix_t>(2);
1818 basis_matrix_t J(2,2,2);
1820 J(0,1) = basis_scalar_t(-1);
1821 J(1,0) = basis_scalar_t( 1);
1822 basis_matrix_t K = J;
1823 K(0,1) = basis_scalar_t( 1);
1824 basis_matrix_t JK = I;
1825 JK(0,0) = basis_scalar_t(-1);
1827 const index_set_t ist_mn = index_set_t(-level);
1828 const index_set_t ist_pn = index_set_t(level);
1838 const framed_pair_t& pair_mn = this->divide(ist_mn);
1839 const framed_multi_t& quo_mn = pair_mn.first;
1840 const framed_multi_t& rem_mn = pair_mn.second;
1841 const framed_pair_t& pair_quo_mnpn = quo_mn.divide(ist_pn);
1842 const framed_multi_t& val_mnpn = pair_quo_mnpn.first;
1843 const framed_multi_t& val_mn = pair_quo_mnpn.second;
1844 const framed_pair_t& pair_rem_mnpn = rem_mn.divide(ist_pn);
1845 const framed_multi_t& val_pn = pair_rem_mnpn.first;
1846 const framed_multi_t& val_1 = pair_rem_mnpn.second;
1849 return -
kron(JK, val_1.fast (level-1, 1))
1850 +
kron(I, val_mnpn.fast(level-1, 1))
1851 +
kron(J, val_mn.fast (level-1, 0))
1852 +
kron(K, val_pn.fast (level-1, 0));
1854 return kron(I, val_1.fast (level-1, 0))
1855 +
kron(JK, val_mnpn.fast(level-1, 0))
1856 +
kron(K, val_mn.fast (level-1, 1))
1857 -
kron(J, val_pn.fast (level-1, 1));
1862 template<
typename Scalar_T, const index_t LO, const index_t HI >
1863 template<
typename Other_Scalar_T >
1864 const matrix_multi<Other_Scalar_T,LO,HI>
1869 multivector_t val = this->fold(frm);
1873 p += std::max(bott_offset,
index_t(0));
1874 q -= std::min(bott_offset,
index_t(0));
1876 throw error_t(
"fast_matrix_multi(frm): HI is too low to represent this value");
1878 throw error_t(
"fast_matrix_multi(frm): LO is too high to represent this value");
1881 val.centre_pm4_qp4(p, q);
1883 val.centre_pp4_qm4(p, q);
1885 val.centre_qp1_pm1(p, q);
1886 const index_t level = (p + q)/2;
1889 const multivector_t& ev_val = val.even();
1890 const multivector_t& od_val = val.odd();
1891 return matrix_multi<Other_Scalar_T,LO,HI>(ev_val.fast(level, 0) + od_val.fast(level, 1), frm);
1894 template<
typename Scalar_T, const index_t LO, const index_t HI >
1902 template<
typename Scalar_T, const index_t LO, const index_t HI >
1908 {
return lhs.first.sign_of_mult(rhs.first) * lhs.second * rhs.second; }
1911 template<
typename Scalar_T, const index_t LO, const index_t HI >
1913 const std::pair<const index_set<LO,HI>, Scalar_T>
1914 operator* (
const std::pair<
const index_set<LO,HI>, Scalar_T>& lhs,
1915 const std::pair<
const index_set<LO,HI>, Scalar_T>& rhs)
1917 typedef std::pair<const index_set<LO,HI>, Scalar_T> term_t;
1918 return term_t(lhs.first ^ rhs.first,
crd_of_mult(lhs, rhs));
1922 template<
typename Scalar_T, const index_t LO, const index_t HI >
1923 const framed_multi<Scalar_T,LO,HI>
1924 sqrt(
const framed_multi<Scalar_T,LO,HI>& val,
const framed_multi<Scalar_T,LO,HI>&
i,
bool prechecked)
1926 typedef numeric_traits<Scalar_T> traits_t;
1928 return traits_t::NaN();
1932 const Scalar_T realval = val.scalar();
1935 if (realval < Scalar_T(0))
1941 return sqrt(matrix_multi_t(val), matrix_multi_t(
i), prechecked);
1945 template<
typename Scalar_T, const index_t LO, const index_t HI >
1946 const framed_multi<Scalar_T,LO,HI>
1947 exp(
const framed_multi<Scalar_T,LO,HI>& val)
1949 typedef numeric_traits<Scalar_T> traits_t;
1951 return traits_t::NaN();
1953 const Scalar_T s =
scalar(val);
1957 const double size = val.size();
1958 const index_t frm_count = val.frame().count();
1967 typedef typename traits_t::demoted::type demoted_scalar_t;
1970 const demoted_multivector_t& demoted_val = demoted_multivector_t(val);
1976 typedef typename traits_t::promoted::type promoted_scalar_t;
1979 const promoted_multivector_t& promoted_val = promoted_multivector_t(val);
1989 typedef matrix_multi<Scalar_T,LO,HI> matrix_multi_t;
1990 return exp(matrix_multi_t(val));
1995 template<
typename Scalar_T, const index_t LO, const index_t HI >
1996 const framed_multi<Scalar_T,LO,HI>
1997 log(
const framed_multi<Scalar_T,LO,HI>& val,
const framed_multi<Scalar_T,LO,HI>&
i,
bool prechecked)
1999 typedef numeric_traits<Scalar_T> traits_t;
2000 if (val == Scalar_T(0) || val.isnan())
2001 return traits_t::NaN();
2005 const Scalar_T realval = val.scalar();
2008 if (realval < Scalar_T(0))
2014 return log(matrix_multi_t(val), matrix_multi_t(
i), prechecked);
2017 #endif // _GLUCAT_FRAMED_MULTI_IMP_H