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 >
79 template<
typename Scalar_T, const index_t LO, const index_t HI >
80 template<
typename Other_Scalar_T >
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))
124 this->insert(
term_t(ist, crd));
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))
137 this->insert(
term_t(ist, crd));
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 >
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();
246 const Scalar_T eps = std::numeric_limits<Scalar_T>::epsilon();
248 #if defined(_GLUCAT_MAP_IS_HASH) 261 if ((abs_crd * abs_crd) > tol)
262 this->insert(
term_t(ist, crd));
267 template<
typename Scalar_T, const index_t LO, const index_t HI >
272 if (this->size() != rhs.size())
277 #if defined(_GLUCAT_MAP_IS_ORDERED) 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;
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 >
309 switch (this->size())
312 return scr == Scalar_T(0);
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 >
336 template<
typename Scalar_T, const index_t LO, const index_t HI >
343 rhs_it = rhs.begin();
351 template<
typename Scalar_T, const index_t LO, const index_t HI >
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 >
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 >
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();
417 const bool direct_mult = lhs_size * rhs_size <= double(algebra_dim)
418 #if defined(_GLUCAT_MAP_IS_HASH) 424 multivector_t result =
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) 446 else if (frm_count < Tune_P::mult_matrix_threshold)
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))
491 template<
typename Scalar_T, const index_t LO, const index_t HI >
496 {
return *
this = *
this * rhs; }
499 template<
typename Scalar_T, const index_t LO, const index_t HI >
505 if (lhs.empty() || rhs.empty())
514 multivector_t result;
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;
531 multivector_t result =
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 >
596 {
return *
this = *
this ^ rhs; }
599 template<
typename Scalar_T, const index_t LO, const index_t HI >
605 if (lhs.empty() || rhs.empty())
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;
626 multivector_t result =
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));
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 >
714 {
return *
this = *
this & rhs; }
717 template<
typename Scalar_T, const index_t LO, const index_t HI >
726 if (lhs.empty() || rhs.empty())
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;
781 multivector_t result =
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 >
850 {
return *
this = *
this % rhs; }
853 template<
typename Scalar_T, const index_t LO, const index_t HI >
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 >
921 if (rhs == Scalar_T(0))
922 return traits_t::NaN();
924 const index_set_t our_frame = lhs.
frame() | rhs.
frame();
929 template<
typename Scalar_T, const index_t LO, const index_t HI >
934 {
return *
this = *
this / rhs; }
937 template<
typename Scalar_T, const index_t LO, const index_t HI >
949 template<
typename Scalar_T, const index_t LO, const index_t HI >
954 {
return *
this = *
this | rhs; }
957 template<
typename Scalar_T, const index_t LO, const index_t HI >
968 template<
typename Scalar_T, const index_t LO, const index_t HI >
975 template<
typename Scalar_T, const index_t LO, const index_t HI >
982 throw error_t(
"outer_pow(int): negative exponent");
994 template<
typename Scalar_T, const index_t LO, const index_t 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 >
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 >
1033 if (this_it == this->end())
1036 return this_it->second;
1040 template<
typename Scalar_T, const index_t LO, const index_t HI >
1045 if ((grade < 0) || (grade > HI-LO))
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 >
1069 template<
typename Scalar_T, const index_t LO, const index_t HI >
1074 {
return *
this - this->
scalar(); }
1077 template<
typename Scalar_T, const index_t LO, const index_t HI >
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 >
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 >
1116 template<
typename Scalar_T, const index_t LO, const index_t HI >
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;
1138 template<
typename Scalar_T, const index_t LO, const index_t HI >
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 >
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 >
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 >
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 >
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 >
1280 random_generator_t& generator = random_generator_t::generator();
1283 (fill < Scalar_T(0))
1285 : (fill > Scalar_T(1))
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;
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 >
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 >
1335 for (
typename map_t::const_iterator
1336 val_it = val.begin();
1337 val_it != val.end();
1339 sorted_val.insert(*val_it);
1340 sorted_begin = sorted_val.begin();
1341 sorted_end = sorted_val.end();
1347 template<
typename Sorted_Map_T >
1356 : sorted_begin( val.begin() ),
1357 sorted_end( val.end() )
1364 template<
typename Scalar_T, const index_t LO, const index_t HI >
1366 operator<< (std::ostream& os, const framed_multi<Scalar_T,LO,HI>& val)
1370 else if (val.isnan())
1371 os << std::numeric_limits<Scalar_T>::quiet_NaN();
1377 typedef typename sorted_map_t::const_iterator sorted_iterator;
1378 sorted_map_t sorted_val;
1380 sorted_iterator sorted_it = sorted_val_range.
sorted_begin;
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 >
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);
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))
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 >
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 >
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 >
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 >
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 >
1659 throw error_t(
"centre_pm4_qp4(p,q): LO is too high to represent this value");
1660 if (this->
frame().max() > p-4)
1663 const index_set_t pm3210(index_pair_t(p-3,p),
true);
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) *
1684 term_t(term.first, term.second));
1687 result.insert(*this_it);
1697 template<
typename Scalar_T, const index_t LO, const index_t 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)
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);
1713 this_it = this->begin();
1714 this_it != this->end();
1718 if (ist.
min() < -q+4)
1728 result.insert(
term_t(term.first, term.second) *
1729 term_t(ist & ~qp0123, this_it->second));
1732 result.insert(*this_it);
1742 template<
typename Scalar_T, const index_t LO, const index_t 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");
1755 this_it = this->begin();
1756 this_it != this->end();
1765 if (n != 0 && ist[n])
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);
1797 template<
typename Scalar_T, const index_t LO, const index_t HI >
1806 const matrix_index_t dim = 1 << level;
1812 return matrix::unit<matrix_t>(1) * this->
scalar();
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);
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 >
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");
1886 const index_t level = (p + q)/2;
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>
1917 typedef std::pair<const index_set<LO,HI>, Scalar_T>
term_t;
1922 template<
typename Scalar_T, const index_t LO, const index_t HI >
1928 return traits_t::NaN();
1932 const Scalar_T realval = val.
scalar();
1935 if (realval < Scalar_T(0))
1945 template<
typename Scalar_T, const index_t LO, const index_t HI >
1951 return traits_t::NaN();
1953 const Scalar_T s =
scalar(val);
1957 const double size = val.size();
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);
1995 template<
typename Scalar_T, const index_t LO, const index_t HI >
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))
2017 #endif // _GLUCAT_FRAMED_MULTI_IMP_H multivector_t & centre_pm4_qp4(index_t &p, index_t &q)
Subalgebra isomorphism: R_{p,q} to R_{p-4,q+4}.
index_set< LO, HI > index_set_t
static Scalar_T 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)
Coordinate of product of terms.
multivector_t & centre_qp1_pm1(index_t &p, index_t &q)
Subalgebra isomorphism: R_{p,q} to R_{q+1,p-1}.
virtual index_t grade() const =0
Maximum of the grades of each term.
class var_term var_term_t
const Multivector< Scalar_T, LO, HI > sqrt(const Multivector< Scalar_T, LO, HI > &val, const Multivector< Scalar_T, LO, HI > &i, const bool prechecked=false)
Square root of multivector with specified complexifier.
index_t count_pos() const
Number of positive indices included in set.
matrix_multi_t::matrix_t matrix_t
Random number generator with single instance per Scalar_T.
Scalar_T abs(const Multivector< Scalar_T, LO, HI > &val)
Absolute value == sqrt(norm)
A framed_multi<Scalar_T,LO,HI> is a framed approximation to a multivector.
framed_multi multivector_t
friend const framed_multi_t operator*(const framed_multi_t &lhs, const framed_multi_t &rhs)
const Multivector< Scalar_T, LO, HI > log(const Multivector< Scalar_T, LO, HI > &val, const Multivector< Scalar_T, LO, HI > &i, const bool prechecked=false)
Natural logarithm of multivector with specified complexifier.
static const std::string classname()
Class name used in messages.
#define _GLUCAT_HASH_N(x)
virtual Scalar_T scalar() const =0
Scalar part.
map_t::const_iterator const_iterator
Sorted_Map_T sorted_map_t
virtual multivector_t & operator|=(const multivector_t &rhs)=0
Transformation via twisted adjoint action.
const RHS_T kron(const LHS_T &lhs, const RHS_T &rhs)
Kronecker tensor product of matrices - as per Matlab kron.
virtual multivector_t & operator*=(const Scalar_T &scr)=0
Product of multivector and scalar.
Sorted_Map_T sorted_map_t
virtual const vector_t vector_part() const =0
Vector part of multivector, as a vector_t with respect to frame()
static const index_t offset_to_super[]
Offsets between the current signature and that of the real superalgebra.
virtual multivector_t & operator/=(const Scalar_T &scr)=0
Quotient of multivector and scalar.
const matrix_t fast(const index_t level, const bool odd) const
Generalized FFT from framed_multi_t to matrix_t.
index_t max() const
Maximum member.
std::vector< Scalar_T > vector_t
_GLUCAT_CLIFFORD_ALGEBRA_OPERATIONS unsigned long nbr_terms() const
Number of terms.
virtual const multivector_t inv() const =0
Geometric multiplicative inverse.
virtual Scalar_T operator[](const index_set_t ist) const =0
Subscripting: map from index set to scalar coordinate.
sorted_iterator sorted_end
Extra traits which extend numeric limits.
virtual const multivector_t reverse() const =0
Reversion, eg. {1}*{2} -> {2}*{1}.
const Multivector< Scalar_T, LO, HI > clifford_exp(const Multivector< Scalar_T, LO, HI > &val)
Exponential of multivector.
const matrix_multi< Other_Scalar_T, LO, HI > fast_matrix_multi(const index_set_t frm) const
Use generalized FFT to construct a matrix_multi_t.
const Multivector< Scalar_T, LO, HI > pow(const Multivector< Scalar_T, LO, HI > &lhs, int rhs)
Integer power of multivector.
virtual const multivector_t even() const =0
Even part of multivector, sum of even grade terms.
const framed_multi_t fast_framed_multi() const
Use inverse generalized FFT to construct a framed_multi_t.
index_t count() const
Cardinality: Number of indices included in set.
static void check_complex(const Multivector< Scalar_T, LO, HI > &val, const Multivector< Scalar_T, LO, HI > &i, const bool prechecked=false)
Check that i is a valid complexifier for val.
sorted_iterator sorted_begin
std::pair< const multivector_t, const multivector_t > framed_pair_t
friend const framed_multi_t operator&(const framed_multi_t &lhs, const framed_multi_t &rhs)
A matrix_multi<Scalar_T,LO,HI> is a matrix approximation to a multivector.
virtual Scalar_T norm() const =0
Scalar_T norm == sum of norm of coordinates.
Sorted_Map_T::const_iterator sorted_iterator
index_t min() const
Minimum member.
bool is_contiguous() const
Determine if the index set is contiguous, ie. has no gaps.
virtual const multivector_t operator-() const =0
Unary -.
Sorted range for use with output.
multivector_t fold(const index_set_t frm) const
Subalgebra isomorphism: fold each term within the given frame.
virtual Scalar_T max_abs() const =0
Maximum of absolute values of components of multivector: multivector infinity norm.
const framed_multi< Scalar_T, LO, HI > exp(const framed_multi< Scalar_T, LO, HI > &val)
Exponential of multivector.
virtual const multivector_t operator()(index_t grade) const =0
Pure grade-vector part.
std::unordered_map< index_set_t, Scalar_T, index_set_hash< LO, HI > > map_t
Abstract exception class.
std::pair< const index_set_t, Scalar_T > term_t
std::pair< index_t, index_t > index_pair_t
virtual const multivector_t involute() const =0
Main involution, each {i} is replaced by -{i} in each term, eg. {1} -> -{1}.
multivector_t & centre_pp4_qm4(index_t &p, index_t &q)
Subalgebra isomorphism: R_{p,q} to R_{p+4,q-4}.
friend const framed_multi_t exp(const framed_multi_t &val)
static const precision_t function_precision
Precision used for exp, log and sqrt functions.
virtual const multivector_t pow(int m) const =0
*this to the m
virtual const multivector_t outer_pow(int m) const =0
Outer product power.
virtual void write(const std::string &msg="") const =0
Write formatted multivector to output.
multivector_t unfold(const index_set_t frm) const
Subalgebra isomorphism: unfold each term within the given frame.
const basis_matrix_t basis_element(const index_set< LO, HI > &ist) const
Create a basis element matrix within the current frame.
virtual const multivector_t odd() const =0
Odd part of multivector, sum of odd grade terms.
static double to_double(const Scalar_T &val)
Cast to double.
sorted_iterator sorted_begin
matrix_t::size_type matrix_index_t
virtual multivector_t & operator-=(const multivector_t &rhs)=0
Geometric difference.
friend const framed_multi_t operator^(const framed_multi_t &lhs, const framed_multi_t &rhs)
friend Scalar_T star(const framed_multi_t &lhs, const framed_multi_t &rhs)
friend const framed_multi_t operator%(const framed_multi_t &lhs, const framed_multi_t &rhs)
friend std::istream & operator>>(std::istream &s, multivector_t &val)
sorted_range(Sorted_Map_T &sorted_val, const Map_T &val)
matrix_t m_matrix
Matrix value representing the multivector within the folded frame.
error< multivector_t > error_t
Index set class based on std::bitset<> in Gnu standard C++ library.
virtual const multivector_t conj() const =0
Conjugation, reverse o involute == involute o reverse.
friend const framed_multi_t operator/(const framed_multi_t &lhs, const framed_multi_t &rhs)
virtual multivector_t & operator^=(const multivector_t &rhs)=0
Outer product.
const framed_pair_t divide(const index_set_t ist) const
Divide multivector into part divisible by index_set and remainder.
virtual const multivector_t truncated(const Scalar_T &limit=Scalar_T(DEFAULT_TRUNCATION)) const =0
Remove all terms with relative size smaller than limit.
virtual bool operator==(const multivector_t &val) const =0
Test for equality of multivectors.
#define _GLUCAT_HASH_SIZE_T(x)
sorted_range(Sorted_Map_T &sorted_val, const Sorted_Map_T &val)
matrix_multi< Scalar_T, LO, HI > matrix_multi_t
ublas::compressed_matrix< int, orientation_t > basis_matrix_t
friend const framed_multi_t operator|(const framed_multi_t &lhs, const framed_multi_t &rhs)
int index_t
Size of index_t should be enough to represent LO, HI.
unsigned long set_value_t
Size of set_value_t should be enough to contain index_set<LO,HI>
Sorted_Map_T::const_iterator sorted_iterator
virtual const multivector_t pure() const =0
Pure part.
framed_multi()
Default constructor.
LHS_T pos_mod(LHS_T lhs, RHS_T rhs)
Modulo function which works reliably for lhs < 0.
virtual multivector_t & operator&=(const multivector_t &rhs)=0
Inner product.
static const framed_multi_t random(const index_set_t frm, Scalar_T fill=Scalar_T(1))
Random multivector within a frame.
virtual Scalar_T quad() const =0
Scalar_T quadratic form == (rev(x)*x)(0)
index_t count_neg() const
Number of negative indices included in set.
virtual bool isnan() const =0
Check if a multivector contains any IEEE NaN values.
multivector_t & operator+=(const term_t &term)
Add a term, if non-zero.
virtual const index_set_t frame() const =0
Subalgebra generated by all generators of terms of given multivector.
std::map< index_set_t, Scalar_T, std::less< const index_set_t > > sorted_map_t
virtual multivector_t & operator%=(const multivector_t &rhs)=0
Contraction.
Matrix_T::size_type nnz(const Matrix_T &m)
Number of non-zeros.
sorted_iterator sorted_end