glucat 0.13.0
framed_multi_imp.h
Go to the documentation of this file.
1#ifndef _GLUCAT_FRAMED_MULTI_IMP_H
2#define _GLUCAT_FRAMED_MULTI_IMP_H
3/***************************************************************************
4 GluCat : Generic library of universal Clifford algebra templates
5 framed_multi_imp.h : Implement the coordinate map representation of a
6 Clifford algebra element
7 -------------------
8 begin : Sun 2001-12-09
9 copyright : (C) 2001-2021 by Paul C. Leopardi
10 ***************************************************************************
11
12 This library is free software: you can redistribute it and/or modify
13 it under the terms of the GNU Lesser General Public License as published
14 by the Free Software Foundation, either version 3 of the License, or
15 (at your option) any later version.
16
17 This library is distributed in the hope that it will be useful,
18 but WITHOUT ANY WARRANTY; without even the implied warranty of
19 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 GNU Lesser General Public License for more details.
21
22 You should have received a copy of the GNU Lesser General Public License
23 along with this library. If not, see <http://www.gnu.org/licenses/>.
24
25 ***************************************************************************
26 This library is based on a prototype written by Arvind Raja and was
27 licensed under the LGPL with permission of the author. See Arvind Raja,
28 "Object-oriented implementations of Clifford algebras in C++: a prototype",
29 in Ablamowicz, Lounesto and Parra (eds.)
30 "Clifford algebras with numeric and symbolic computations", Birkhauser, 1996.
31 ***************************************************************************
32 See also Arvind Raja's original header comments in glucat.h
33 ***************************************************************************/
34
35#include "glucat/framed_multi.h"
36
37#include "glucat/scalar.h"
38#include "glucat/random.h"
39#include "glucat/generation.h"
40#include "glucat/matrix.h"
41
42#include <sstream>
43#include <fstream>
44
45namespace glucat
46{
48 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
49 auto
51 classname() -> const std::string
52 { return "framed_multi"; }
53
54#define _GLUCAT_HASH_N(x) (x)
55#define _GLUCAT_HASH_SIZE_T(x) (typename multivector_t::hash_size_t)(x)
56
58 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
63
65 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
70
72 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
73 template< typename Other_Scalar_T >
76 : map_t(_GLUCAT_HASH_N(val.size()))
77 {
78 for (auto& val_term : val)
79 this->insert(term_t(val_term.first, numeric_traits<Scalar_T>::to_scalar_t(val_term.second)));
80 }
81
83 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
84 template< typename Other_Scalar_T >
87 const index_set_t frm, const bool prechecked)
88 : map_t(_GLUCAT_HASH_N(val.size()))
89 {
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)
93 this->insert(term_t(val_term.first, numeric_traits<Scalar_T>::to_scalar_t(val_term.second)));
94 }
95
97 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
100 const index_set_t frm, const bool prechecked)
101 : map_t(_GLUCAT_HASH_N(val.size()))
102 {
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);
107 }
108
110 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
112 framed_multi(const index_set_t ist, const Scalar_T& crd)
114 {
115 if (crd != Scalar_T(0))
116 this->insert(term_t(ist, crd));
117 }
118
120 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
122 framed_multi(const index_set_t ist, const Scalar_T& crd,
123 const index_set_t frm, const bool prechecked)
125 {
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));
130 }
131
133 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
135 framed_multi(const Scalar_T& scr, const index_set_t frm)
137 {
138 if (scr != Scalar_T(0))
139 this->insert(term_t(index_set_t(), scr));
140 }
141
143 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
145 framed_multi(const int scr, const index_set_t frm)
147 {
148 if (scr != Scalar_T(0))
149 this->insert(term_t(index_set_t(), Scalar_T(scr)));
150 }
151
153 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
155 framed_multi(const vector_t& vec,
156 const index_set_t frm, const bool prechecked)
157 : map_t(_GLUCAT_HASH_N(vec.size()))
158 {
159 if (!prechecked && index_t(vec.size()) != frm.count())
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)
164 {
165 *this += term_t(index_set_t(idx), crd);
166 for (
167 ++idx;
168 idx != frm_end && !frm[idx];
169 ++idx)
170 ;
172 }
173
175 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
177 framed_multi(const std::string& str)
180 std::istringstream ss(str);
181 ss >> *this;
182 if (!ss)
183 throw error_t("multivector_t(str): could not parse string");
184 // Peek to see if the end of the string has been reached.
185 ss.peek();
186 if (!ss.eof())
187 throw error_t("multivector_t(str): could not parse entire string");
188 }
191 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
193 framed_multi(const std::string& str, const index_set_t frm, const bool prechecked)
195 {
196 if (prechecked)
197 *this = multivector_t(str);
198 else
199 *this = multivector_t(multivector_t(str), frm, false);
200 }
203 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
204 template< typename Other_Scalar_T >
208 {
209 if (val == Other_Scalar_T(0))
210 return;
211
212 const auto dim = val.m_matrix.size1();
213 using traits_t = numeric_traits<Scalar_T>;
214 if (dim == 1)
215 {
216 this->insert(term_t(index_set_t(), traits_t::to_scalar_t(val.m_matrix(0, 0))));
217 return;
218 }
219 if (dim >= Tune_P::inv_fast_dim_threshold)
220 try
221 {
222 *this = (val.template fast_framed_multi<Scalar_T>()).truncated();
223 return;
224 }
225 catch (const glucat_error& e)
226 { }
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();
232 return;
233 }
234 const auto frm = val.frame();
235 const auto algebra_dim = set_value_t(1) << frm.count();
236 auto result = multivector_t(
237 _GLUCAT_HASH_SIZE_T(std::min<size_t>(algebra_dim, matrix::nnz(val.m_matrix))));
238 for (auto
239 stv = set_value_t(0);
240 stv != algebra_dim;
241 stv++)
242 {
243 const auto ist = index_set_t(stv, frm, true);
244 const auto crd =
245 traits_t::to_scalar_t(matrix::inner<Other_Scalar_T>(val.basis_element(ist), val.m_matrix));
246 if (crd != Scalar_T(0))
247 result.insert(term_t(ist, crd));
248 }
249 *this = result.truncated();
250 }
251
253 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
254 auto
256 operator== (const multivector_t& rhs) const -> bool
257 {
258 if (this->size() != rhs.size())
259 return false;
260 const auto rhs_end = rhs.end();
261 for (auto& this_term : *this)
262 {
263 const const_iterator& rhs_it = rhs.find(this_term.first);
264 if (rhs_it == rhs_end || rhs_it->second != this_term.second)
265 return false;
267 return true;
269
270 /// Test for equality of multivector and scalar
271 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
272 inline
273 auto
275 operator== (const Scalar_T& scr) const -> bool
277 switch (this->size())
278 {
279 case 0:
280 return scr == Scalar_T(0);
281 case 1:
282 {
283 const auto& this_it = this->begin();
284 return this_it->first == index_set_t() && this_it->second == scr;
285 }
286 default:
287 return false;
288 }
289 }
290
292 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
293 inline
294 auto
296 operator+= (const Scalar_T& scr) -> multivector_t&
297 {
298 *this += term_t(index_set_t(), scr);
299 return *this;
300 }
301
303 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
304 inline
305 auto
307 operator+= (const multivector_t& rhs) -> multivector_t&
308 { // simply add terms
309 for (auto& rhs_term : rhs)
310 *this += rhs_term;
311 return *this;
312 }
313
315 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
316 inline
317 auto
319 operator-= (const Scalar_T& scr) -> multivector_t&
320 {
321 *this += term_t(index_set_t(), -scr);
322 return *this;
323 }
324
326 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
327 inline
328 auto
330 operator-= (const multivector_t& rhs) -> multivector_t&
331 {
332 for (auto& rhs_term : rhs)
333 *this += term_t(rhs_term.first, -(rhs_term.second));
334 return *this;
335 }
336
338 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
339 inline
340 auto
342 operator- () const -> const multivector_t
343 { // multiply coordinates of all terms by -1
344 auto result = *this;
345 for (auto& result_term : result)
346 result_term.second *= Scalar_T(-1);
347 return result;
348 }
349
351 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
352 auto
354 operator*= (const Scalar_T& scr) -> multivector_t&
355 { // multiply coordinates of all terms by scalar
356 using traits_t = numeric_traits<Scalar_T>;
357
358 if (traits_t::isNaN_or_isInf(scr))
359 return *this = traits_t::NaN();
360 if (scr == Scalar_T(0))
361 if (this->isnan())
362 *this = traits_t::NaN();
363 else
364 this->clear();
365 else
366 for (auto& this_term : *this)
367 this_term.second *= scr;
368 return *this;
369 }
370
372 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
373 auto
375 {
376 using multivector_t = framed_multi<Scalar_T,LO,HI,Tune_P>;
377 using traits_t = numeric_traits<Scalar_T>;
378
379 if (lhs.isnan() || rhs.isnan())
380 return traits_t::NaN();
381
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);
388 if (direct_mult)
389 { // If we have a sparse multiply, store the result directly
390 auto result = multivector_t(
391 _GLUCAT_HASH_SIZE_T(size_t(std::min(lhs_size * rhs_size, double(algebra_dim)))));
392 for (auto& lhs_term : lhs)
393 for (auto& rhs_term : rhs)
394 result += lhs_term * rhs_term;
395 return result;
396 }
397 else
398 { // Past a certain threshold, the matrix algorithm is fastest
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);
402 }
403 }
404
406 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
407 inline
408 auto
410 operator*= (const multivector_t& rhs) -> multivector_t&
411 { return *this = *this * rhs; }
412
414 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
415 auto
417 { // Arvind Raja's original reference:
418 // "old clical, outerproduct(p,q:pterm):pterm in file compmod.pas"
419
420 if (lhs.empty() || rhs.empty())
421 return Scalar_T(0);
422
423 using multivector_t = framed_multi<Scalar_T,LO,HI,Tune_P>;
424 using index_set_t = typename multivector_t::index_set_t;
425 using term_t = typename multivector_t::term_t;
426
427 const auto empty_set = index_set_t();
428
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(
436 _GLUCAT_HASH_SIZE_T(size_t(std::min(lhs_size * rhs_size, double(algebra_dim)))));
437 const auto lhs_end = lhs.end();
438 const auto rhs_end = rhs.end();
439
440 if (lhs_size * rhs_size > double(Tune_P::products_size_threshold))
441 {
442 for (auto
443 result_stv = set_value_t(0);
444 result_stv != algebra_dim;
445 ++result_stv)
446 {
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);
451 for (auto
452 lhs_stv = set_value_t(0);
453 lhs_stv != lhs_result_dim;
454 ++lhs_stv)
455 {
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)
459 {
460 const auto lhs_it = lhs.find(lhs_ist);
461 if (lhs_it != lhs_end)
462 {
463 const auto rhs_it = rhs.find(rhs_ist);
464 if (rhs_it != rhs_end)
465 result_crd += crd_of_mult(*lhs_it, *rhs_it);
466 }
467 }
468 }
469 if (result_crd != Scalar_T(0))
470 result.insert(term_t(result_ist, result_crd));
471 }
472 return result;
473 }
474 else
475 {
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;
480 return result;
481 }
482 }
483
485 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
486 inline
487 auto
489 operator^= (const multivector_t& rhs) -> multivector_t&
490 { return *this = *this ^ rhs; }
491
493 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
494 auto
496 { // Arvind Raja's original reference:
497 // "old clical, innerproduct(p,q:pterm):pterm in file compmod.pas"
498
499 if (lhs.empty() || rhs.empty())
500 return Scalar_T(0);
501
502 using multivector_t = framed_multi<Scalar_T,LO,HI,Tune_P>;
503 using index_set_t = typename multivector_t::index_set_t;
504 using term_t = typename multivector_t::term_t;
505
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();
510
511 const auto lhs_frame = lhs.frame();
512 const auto rhs_frame = rhs.frame();
513
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(
517 _GLUCAT_HASH_SIZE_T(size_t(std::min(lhs_size * rhs_size, double(algebra_dim)))));
518 if (lhs_size * rhs_size > double(Tune_P::products_size_threshold))
519 {
520 for (auto
521 result_stv = set_value_t(0);
522 result_stv != algebra_dim;
523 ++result_stv)
524 {
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);
529 for (auto
530 comp_stv = set_value_t(1);
531 comp_stv != comp_dim;
532 ++comp_stv)
533 {
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)
537 {
538 const auto lhs_it = lhs.find(our_ist);
539 if (lhs_it != lhs_end)
540 {
541 const auto rhs_it = rhs.find(comp_ist);
542 if (rhs_it != rhs_end)
543 result_crd += crd_of_mult(*lhs_it, *rhs_it);
544 }
545 }
546 if (result_stv != 0)
547 {
548 if ((our_ist | rhs_frame) == rhs_frame)
549 {
550 const auto rhs_it = rhs.find(our_ist);
551 if (rhs_it != rhs_end)
552 {
553 const auto lhs_it = lhs.find(comp_ist);
554 if (lhs_it != lhs_end)
555 result_crd += crd_of_mult(*lhs_it, *rhs_it);
556 }
557 }
558 }
559 }
560 if (result_crd != Scalar_T(0))
561 result.insert(term_t(result_ist, result_crd));
562 }
563 }
564 else
565 {
566 const auto empty_set = index_set_t();
567 for (auto& lhs_term : lhs)
568 {
569 const auto lhs_ist = lhs_term.first;
570 if (lhs_ist != empty_set)
571 for (auto& rhs_term : rhs)
572 {
573 const auto rhs_ist = rhs_term.first;
574 if (rhs_ist != empty_set)
575 {
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;
579 }
580 }
581 }
582 }
583 return result;
584 }
585
587 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
588 inline
589 auto
591 operator&= (const multivector_t& rhs) -> multivector_t&
592 { return *this = *this & rhs; }
593
595 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
596 auto
598 {
599 // Reference: Leo Dorst, "Honing geometric algebra for its use in the computer sciences",
600 // in Geometric Computing with Clifford Algebras, ed. G. Sommer,
601 // Springer 2001, Chapter 6, pp. 127-152.
602 // http://staff.science.uva.nl/~leo/clifford/index.html
603
604 if (lhs.empty() || rhs.empty())
605 return Scalar_T(0);
606
607 using multivector_t = framed_multi<Scalar_T,LO,HI,Tune_P>;
608 using index_set_t = typename multivector_t::index_set_t;
609 using term_t = typename multivector_t::term_t;
610
611 const auto lhs_end = lhs.end();
612 const auto rhs_end = rhs.end();
613 const double lhs_size = lhs.size();
614 const double rhs_size = rhs.size();
615 const auto lhs_frame = lhs.frame();
616 const auto rhs_frame = rhs.frame();
617
618 const auto our_frame = lhs_frame | rhs_frame;
619 const auto algebra_dim = set_value_t(1) << our_frame.count();
620 auto result = multivector_t(
621 _GLUCAT_HASH_SIZE_T(size_t(std::min(lhs_size * rhs_size, double(algebra_dim)))));
622
623 if (lhs_size * rhs_size > double(Tune_P::products_size_threshold))
624 {
625 for (auto
626 result_stv = set_value_t(0);
627 result_stv != algebra_dim;
628 ++result_stv)
629 {
630 const auto result_ist = index_set_t(result_stv, our_frame, true);
631 const auto comp_frame = lhs_frame & ~result_ist;
632 const auto comp_dim = set_value_t(1) << comp_frame.count();
633 auto result_crd = Scalar_T(0);
634 for (auto
635 comp_stv = set_value_t(0);
636 comp_stv != comp_dim;
637 ++comp_stv)
638 {
639 const auto comp_ist = index_set_t(comp_stv, comp_frame, true);
640 const auto rhs_ist = result_ist ^ comp_ist;
641 if ((rhs_ist | rhs_frame) == rhs_frame)
642 {
643 const auto rhs_it = rhs.find(rhs_ist);
644 if (rhs_it != rhs_end)
645 {
646 const auto lhs_it = lhs.find(comp_ist);
647 if (lhs_it != lhs_end)
648 result_crd += crd_of_mult(*lhs_it, *rhs_it);
649 }
650 }
651 }
652 if (result_crd != Scalar_T(0))
653 result.insert(term_t(result_ist, result_crd));
654 }
655 }
656 else
657 {
658 for (auto& rhs_term : rhs)
659 {
660 const auto rhs_ist = rhs_term.first;
661 for (auto& lhs_term : lhs)
662 {
663 const index_set_t lhs_ist = lhs_term.first;
664 if ((lhs_ist | rhs_ist) == rhs_ist)
665 result += lhs_term * rhs_term;
666 }
667 }
668 }
669 return result;
670 }
671
673 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
674 inline
675 auto
677 operator%= (const multivector_t& rhs) -> multivector_t&
678 { return *this = *this % rhs; }
679
681 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
682 auto
684 {
685 auto result = Scalar_T(0);
686 const auto small_star_large = lhs.size() < rhs.size();
687 const auto* smallp =
688 small_star_large
689 ? &lhs
690 : &rhs;
691 const auto* largep =
692 small_star_large
693 ? &rhs
694 : &lhs;
695
696 for (auto& small_term : *smallp)
697 {
698 const auto small_ist = small_term.first;
699 const auto large_crd = (*largep)[small_ist];
700 if (large_crd != Scalar_T(0))
701 result += small_ist.sign_of_square() * small_term.second * large_crd;
702 }
703 return result;
704 }
705
707 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
708 auto
710 operator/= (const Scalar_T& scr) -> multivector_t&
711 { // Divide coordinates of all terms by scr
712 using traits_t = numeric_traits<Scalar_T>;
713
714 if (traits_t::isNaN(scr))
715 return *this = traits_t::NaN();
716 if (traits_t::isInf(scr))
717 if (this->isnan())
718 *this = traits_t::NaN();
719 else
720 this->clear();
721 else
722 for (auto& this_term : *this)
723 this_term.second /= scr;
724 return *this;
725 }
726
728 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
729 inline
730 auto
732 {
733 using multivector_t = framed_multi<Scalar_T,LO,HI,Tune_P>;
734 using traits_t = numeric_traits<Scalar_T>;
735 using matrix_multi_t = typename multivector_t::matrix_multi_t;
736
737 if (rhs == Scalar_T(0))
738 return traits_t::NaN();
739
740 const auto our_frame = lhs.frame() | rhs.frame();
741 return matrix_multi_t(lhs, our_frame, true) / matrix_multi_t(rhs, our_frame, true);
742 }
743
745 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
746 inline
747 auto
749 operator/= (const multivector_t& rhs) -> multivector_t&
750 { return *this = *this / rhs; }
751
753 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
754 inline
755 auto
757 {
758 using multivector_t = framed_multi<Scalar_T,LO,HI,Tune_P>;
759 using matrix_multi_t = typename multivector_t::matrix_multi_t;
760
761 return matrix_multi_t(rhs) * matrix_multi_t(lhs) / matrix_multi_t(rhs.involute());
762 }
763
765 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
766 inline
767 auto
769 operator|= (const multivector_t& rhs) -> multivector_t&
770 { return *this = *this | rhs; }
771
773 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
774 inline
775 auto
777 inv() const -> const multivector_t
778 {
779 auto result = matrix_multi_t(Scalar_T(1), this->frame());
780 return result /= matrix_multi_t(*this);
781 }
782
784 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
785 auto
787 pow(int m) const -> const multivector_t
788 { return glucat::pow(*this, m); }
789
791 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
792 auto
794 outer_pow(int m) const -> const multivector_t
795 {
796 if (m < 0)
797 throw error_t("outer_pow(int): negative exponent");
798 auto result = multivector_t(Scalar_T(1));
799 auto a = *this;
800 for (;
801 m != 0;
802 m >>= 1, a = a ^ a)
803 if (m & 1)
804 result ^= a;
805 return result;
806 }
807
809 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
810 inline
811 auto
813 frame() const -> const index_set_t
814 {
815 auto result = index_set_t();
816 for (auto& this_term : *this)
817 result |= this_term.first;
818 return result;
819 }
820
822 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
823 inline
824 auto
826 grade() const -> index_t
827 {
828 auto result = index_t(0);
829 for (auto& this_term : *this)
830 result = std::max(result, this_term.first.count());
831 return result;
832 }
833
835 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
836 inline
837 auto
839 operator[] (const index_set_t ist) const -> Scalar_T
840 {
841 const auto& this_it = this->find(ist);
842 if (this_it == this->end())
843 return Scalar_T(0);
844 else
845 return this_it->second;
846 }
847
849 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
850 auto
852 operator() (index_t grade) const -> const multivector_t
853 {
854 if ((grade < 0) || (grade > HI-LO))
855 return Scalar_T(0);
856 else
857 {
858 auto result = multivector_t();
859 for (auto& this_term : *this)
860 if (this_term.first.count() == grade)
861 result += this_term;
862 return result;
863 }
864 }
865
867 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
868 inline
869 auto
871 scalar() const -> Scalar_T
872 { return (*this)[index_set_t()]; }
873
875 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
876 inline
877 auto
879 pure() const -> const multivector_t
880 { return *this - this->scalar(); }
881
883 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
884 auto
886 even() const -> const multivector_t
887 { // even part of x, sum of the pure(count) with even count
888 auto result = multivector_t();
889 for (auto& this_term : *this)
890 if ((this_term.first.count() % 2) == 0)
891 result.insert(this_term);
892 return result;
893 }
894
896 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
897 auto
899 odd() const -> const multivector_t
900 { // even part of x, sum of the pure(count) with even count
901 auto result = multivector_t();
902 for (auto& this_term : *this)
903 if ((this_term.first.count() % 2) == 1)
904 result.insert(this_term);
905 return result;
906 }
907
909 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
910 auto
912 vector_part() const -> const vector_t
913 { return this->vector_part(this->frame(), true); }
914
916 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
917 auto
919 vector_part(const index_set_t frm, const bool prechecked) const -> const vector_t
920 {
921 if (!prechecked && (this->frame() | frm) != frm)
922 throw error_t("vector_part(frm): value is outside of requested frame");
923 auto result = vector_t();
924 result.reserve(frm.count());
925 const auto frm_end = frm.max()+1;
926 for (auto
927 idx = frm.min();
928 idx != frm_end;
929 ++idx)
930 // Frame may contain indices which do not correspond to a grade 1 term but
931 // frame cannot omit any index corresponding to a grade 1 term
932 if (frm[idx])
933 result.push_back((*this)[index_set_t(idx)]);
934 return result;
935 }
936
938 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
939 auto
941 involute() const -> const multivector_t
942 {
943 auto result = *this;
944 for (auto& result_term : result)
945 { // for a k-vector u, involute(u) == (-1)^k * u
946 if ((result_term.first.count() % 2) == 1)
947 result_term.second *= Scalar_T(-1);
948 }
949 return result;
950 }
951
953 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
954 auto
956 reverse() const -> const multivector_t
957 {
958 auto result = *this;
959 for (auto& result_term : result)
960 // For a k-vector u, reverse(u) = { -u, k == 2,3 (mod 4)
961 // { u, k == 0,1 (mod 4)
962 switch (result_term.first.count() % 4)
963 {
964 case 2:
965 case 3:
966 result_term.second *= Scalar_T(-1);
967 break;
968 default:
969 break;
970 }
971 return result;
972 }
973
975 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
976 auto
978 conj() const -> const multivector_t
979 {
980 auto result = *this;
981 for (auto& result_term : result)
982 // For a k-vector u, conj(u) = { -u, k == 1,2 (mod 4)
983 // { u, k == 0,3 (mod 4)
984 switch (result_term.first.count() % 4)
985 {
986 case 1:
987 case 2:
988 result_term.second *= Scalar_T(-1);
989 break;
990 default:
991 break;
992 }
993 return result;
994 }
995
997 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
998 auto
1000 quad() const -> Scalar_T
1001 {
1002 // scalar(conj(x)*x) = 2*quad(even(x)) - quad(x)
1003 // ref: old clical: quadfunction(p:pter):pterm in file compmod.pas
1004 auto result = Scalar_T(0);
1005 for (auto& this_term : *this)
1006 {
1007 const auto sign =
1008 (this_term.first.count_neg() % 2)
1009 ? -Scalar_T(1)
1010 : Scalar_T(1);
1011 result += sign * (this_term.second) * (this_term.second);
1012 }
1013 return result;
1014 }
1015
1017 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1018 auto
1020 norm() const -> Scalar_T
1021 {
1022 using traits_t = numeric_traits<Scalar_T>;
1023
1024 auto result = Scalar_T(0);
1025 for (auto& this_term : *this)
1026 {
1027 const auto abs_crd = traits_t::abs(this_term.second);
1028 result += abs_crd * abs_crd;
1029 }
1030 return result;
1031 }
1032
1034 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1035 auto
1037 max_abs() const -> Scalar_T
1038 {
1039 using traits_t = numeric_traits<Scalar_T>;
1040
1041 auto result = Scalar_T(0);
1042 for (auto& this_term : *this)
1043 {
1044 const auto abs_crd = traits_t::abs(this_term.second);
1045 if (abs_crd > result)
1046 result = abs_crd;
1047 }
1048 return result;
1049 }
1050
1052 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1053 auto
1055 random(const index_set_t frm, Scalar_T fill) -> const multivector_t
1056 {
1058 using index_set_t = typename multivector_t::index_set_t;
1059 using term_t = typename multivector_t::term_t;
1060
1061 using random_generator_t = random_generator<Scalar_T>;
1062 auto& generator = random_generator_t::generator();
1063
1064 fill =
1065 (fill < Scalar_T(0))
1066 ? Scalar_T(0)
1067 : (fill > Scalar_T(1))
1068 ? Scalar_T(1)
1069 : fill;
1070 const auto algebra_dim = set_value_t(1) << frm.count();
1071 using traits_t = numeric_traits<Scalar_T>;
1072 const auto mean_abs = traits_t::sqrt(Scalar_T(double(algebra_dim)));
1073 auto result = multivector_t();
1074 for (auto
1075 stv = set_value_t(0);
1076 stv != algebra_dim;
1077 ++stv)
1078 if (generator.uniform() < fill)
1079 {
1080 const auto& result_crd = generator.normal() / mean_abs;
1081 result.insert(term_t(index_set_t(stv, frm, true), result_crd));
1082 }
1083 return result;
1084 }
1085
1087 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1088 inline
1089 void
1091 write(const std::string& msg) const
1092 { std::cout << msg << std::endl << " " << (*this) << std::endl; }
1093
1095 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1096 inline
1097 void
1099 write(std::ofstream& ofile, const std::string& msg) const
1100 {
1101 if (!ofile)
1102 throw error_t("write(ofile,msg): cannot write to output file");
1103 ofile << msg << std::endl << " " << (*this) << std::endl;
1104 }
1105
1107 template< typename Map_T,typename Sorted_Map_T >
1109 {
1110 public:
1111 using map_t = Map_T;
1112 using sorted_map_t = Sorted_Map_T;
1113 using sorted_iterator = typename Sorted_Map_T::const_iterator;
1114
1115 sorted_range (Sorted_Map_T &sorted_val, const Map_T& val)
1116 {
1117 for (auto& val_term : val)
1118 sorted_val.insert(val_term);
1119 sorted_begin = sorted_val.begin();
1120 sorted_end = sorted_val.end();
1121 }
1124 };
1125
1126 template< typename Sorted_Map_T >
1127 class sorted_range< Sorted_Map_T, Sorted_Map_T >
1128 {
1129 public:
1130 using map_t = Sorted_Map_T;
1131 using sorted_map_t = Sorted_Map_T;
1132 using sorted_iterator = typename Sorted_Map_T::const_iterator;
1133
1134 sorted_range (Sorted_Map_T &sorted_val, const Sorted_Map_T& val)
1135 : sorted_begin( val.begin() ),
1136 sorted_end( val.end() )
1137 { }
1140 };
1141
1143 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1144 auto
1145 operator<< (std::ostream& os, const framed_multi<Scalar_T,LO,HI,Tune_P>& val) -> std::ostream&
1146 {
1147 using limits_t = std::numeric_limits<Scalar_T>;
1148 if (val.empty())
1149 os << 0;
1150 else if (val.isnan())
1151 os << limits_t::quiet_NaN();
1152 else if (val.isinf())
1153 {
1154 const Scalar_T& inf = limits_t::infinity();
1155 os << (scalar(val) < 0.0 ? -inf : inf);
1156 }
1157 else
1158 {
1159 using traits_t = numeric_traits<Scalar_T>;
1160 using multivector_t = framed_multi<Scalar_T,LO,HI,Tune_P>;
1161 Scalar_T truncation;
1162 switch (os.flags() & std::ios::floatfield)
1163 {
1164 case std::ios_base::scientific:
1165 truncation = Scalar_T(1) / traits_t::pow(Scalar_T(10), int(os.precision()) + 1);
1166 break;
1167 case std::ios_base::fixed:
1168 truncation = Scalar_T(1) / (traits_t::pow(Scalar_T(10), int(os.precision())) * val.max_abs());
1169 break;
1170 case std::ios_base::fixed | std::ios_base::scientific:
1171 truncation = multivector_t::default_truncation;
1172 break;
1173 default:
1174 truncation = Scalar_T(1) / traits_t::pow(Scalar_T(10), int(os.precision()));
1175 break;
1176 }
1177 auto truncated_val = val.truncated(truncation);
1178 if (truncated_val.empty())
1179 os << 0;
1180 else
1181 {
1182 using map_t = typename multivector_t::map_t;
1183 using sorted_map_t = typename multivector_t::sorted_map_t;
1184 auto sorted_val = sorted_map_t();
1185 const auto sorted_val_range = sorted_range< map_t, sorted_map_t >(sorted_val, truncated_val);
1186 auto sorted_it = sorted_val_range.sorted_begin;
1187 os << *sorted_it;
1188 for (++sorted_it;
1189 sorted_it != sorted_val_range.sorted_end;
1190 ++sorted_it)
1191 {
1192 const Scalar_T& scr = sorted_it->second;
1193 if (scr >= 0.0)
1194 os << '+';
1195 os << *sorted_it;
1196 }
1197 }
1198 }
1199 return os;
1200 }
1201
1203 template< typename Scalar_T, const index_t LO, const index_t HI >
1204 auto
1205 operator<< (std::ostream& os, const std::pair< const index_set<LO,HI>, Scalar_T >& term) -> std::ostream&
1206 {
1207 const auto second_as_double = numeric_traits<Scalar_T>::to_double(term.second);
1208 const auto use_double =
1209 (os.precision() <= std::numeric_limits<double>::digits10) ||
1210 (term.second == Scalar_T(second_as_double));
1211 if (term.first.count() == 0)
1212 if (use_double)
1213 os << second_as_double;
1214 else
1215 os << term.second;
1216 else if (term.second == Scalar_T(-1))
1217 {
1218 os << '-';
1219 os << term.first;
1220 }
1221 else if (term.second != Scalar_T(1))
1222 {
1223 if (use_double)
1224 {
1225 auto tol = std::pow(10.0,-os.precision());
1226 if ( std::fabs(second_as_double + 1.0) < tol )
1227 os << '-';
1228 else if ( std::fabs(second_as_double - 1.0) >= tol )
1229 os << second_as_double;
1230 }
1231 else
1232 os << term.second;
1233 os << term.first;
1234 }
1235 else
1236 os << term.first;
1237 return os;
1238 }
1239
1241 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1242 auto
1243 operator>> (std::istream& s, framed_multi<Scalar_T,LO,HI,Tune_P> & val) -> std::istream&
1244 { // Input looks like 1.0-2.0{1,2}+3.2{3,4}.
1245 using multivector_t = framed_multi<Scalar_T,LO,HI,Tune_P>;
1246 // Parsing variables.
1247 auto local_val = multivector_t();
1248 auto c = 0;
1249 // Parsing control variables.
1250 auto negative = false;
1251 auto expect_term = true;
1252 // The multivector may begin with '+' or '-'. Check for this.
1253 c = s.peek();
1254 if (s.good() && (c == int('+') || c == int('-')))
1255 { // A '-' here negates the following term.
1256 negative = (c == int('-'));
1257 // Consume the '+' or '-'.
1258 s.get();
1259 }
1260 while (s.good())
1261 { // Parse a term.
1262 // A term consists of an optional scalar, followed by an optional index set.
1263 // At least one of the two must be present.
1264 // Default coordinate is Scalar_T(1).
1265 auto coordinate = Scalar_T(1);
1266 // Default index set is empty.
1267 auto ist = index_set<LO,HI>();
1268 // First, check for an opening brace.
1269 c = s.peek();
1270 if (s.good())
1271 { // If the character is not an opening brace,
1272 // a coordinate value is expected here.
1273 if (c != int('{'))
1274 { // Try to read a coordinate value.
1275 double coordinate_as_double;
1276 s >> coordinate_as_double;
1277 // Reading the coordinate may have resulted in an end of file condition.
1278 // This is not a failure.
1279 if (s)
1280 coordinate = Scalar_T(coordinate_as_double);
1281 }
1282 }
1283 else
1284 { // End of file here ends parsing while a term may still be expected.
1285 break;
1286 }
1287 // Coordinate is now Scalar_T(1) or a Scalar_T value.
1288 // Parse an optional index set.
1289 if (s.good())
1290 {
1291 c = s.peek();
1292 if (s.good() && c == int('{'))
1293 { // Try to read index set.
1294 s >> ist;
1295 }
1296 }
1297 // Reading the term may have resulted in an end of file condition.
1298 // This is not a failure.
1299 if (s)
1300 {
1301 // Immediately after parsing a term, another term is not expected.
1302 expect_term = false;
1303 if (coordinate != Scalar_T(0))
1304 {
1305 // Add the term to the local multivector.
1306 coordinate =
1307 negative
1308 ? -coordinate
1309 : coordinate;
1310 using term_t = typename multivector_t::term_t;
1311 local_val += term_t(ist, coordinate);
1312 }
1313 }
1314 // Check if anything follows the current term.
1315 if (s.good())
1316 {
1317 c = s.peek();
1318 if (s.good())
1319 { // Only '+' and '-' are valid here.
1320 if (c == int('+') || c == int('-'))
1321 { // A '-' here negates the following term.
1322 negative = (c == int('-'));
1323 // Consume the '+' or '-'.
1324 s.get();
1325 // Immediately after '+' or '-',
1326 // expect another term.
1327 expect_term = true;
1328 }
1329 else
1330 { // Any other character here is a not failure,
1331 // but still ends the parsing of the multivector.
1332 break;
1333 }
1334 }
1335 }
1336 }
1337 // If a term is still expected, this is a failure.
1338 if (expect_term)
1339 s.clear(std::istream::failbit);
1340 // End of file is not a failure.
1341 if (s)
1342 { // The multivector has been successfully parsed.
1343 val = local_val;
1344 }
1345 return s;
1346 }
1347
1349 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1350 auto
1352 nbr_terms () const -> unsigned long
1353 { return this->size(); }
1354
1356 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1357 inline
1358 auto
1360 operator+= (const term_t& term) -> multivector_t&
1361 { // Do not insert terms with 0 coordinate
1362 if (term.second != Scalar_T(0))
1363 {
1364 const auto& this_it = this->find(term.first);
1365 if (this_it == this->end())
1366 this->insert(term);
1367 else if (this_it->second + term.second == Scalar_T(0))
1368 // Erase term if resulting coordinate is 0
1369 this->erase(this_it);
1370 else
1371 this_it->second += term.second;
1372 }
1373 return *this;
1374 }
1375
1377 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1378 auto
1380 isinf() const -> bool
1381 {
1382 using traits_t = numeric_traits<Scalar_T>;
1383
1384 if (std::numeric_limits<Scalar_T>::has_infinity)
1385 for (auto& this_term : *this)
1386 if (traits_t::isInf(this_term.second))
1387 return true;
1388 return false;
1389 }
1390
1392 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1393 auto
1395 isnan() const -> bool
1396 {
1397 using traits_t = numeric_traits<Scalar_T>;
1398
1399 if (std::numeric_limits<Scalar_T>::has_quiet_NaN)
1400 for (auto& this_term : *this)
1401 if (traits_t::isNaN(this_term.second))
1402 return true;
1403 return false;
1404 }
1405
1407 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1408 auto
1410 truncated(const Scalar_T& limit) const -> const multivector_t
1411 {
1412 using traits_t = numeric_traits<Scalar_T>;
1413
1414 if (this->isnan() || this->isinf())
1415 return *this;
1416 const auto truncation = traits_t::abs(limit);
1417 const auto top = max_abs();
1418 auto result = multivector_t();
1419 if (top != Scalar_T(0))
1420 for (auto& this_term : *this)
1421 if (traits_t::abs(this_term.second) > top * truncation)
1422 result.insert(this_term);
1423 return result;
1424 }
1425
1427 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1428 auto
1430 fold(const index_set_t frm) const -> multivector_t
1431 {
1432 if (frm.is_contiguous())
1433 return *this;
1434 else
1435 {
1436 auto result = multivector_t();
1437 for (auto& this_term : *this)
1438 result.insert(term_t(this_term.first.fold(frm), this_term.second));
1439 return result;
1440 }
1441 }
1442
1444 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1445 auto
1447 unfold(const index_set_t frm) const -> multivector_t
1448 {
1449 if (frm.is_contiguous())
1450 return *this;
1451 else
1452 {
1453 auto result = multivector_t();
1454 for (auto& this_term : *this)
1455 result.insert(term_t(this_term.first.unfold(frm), this_term.second));
1456 return result;
1457 }
1458 }
1459
1461 // Reference: [L] 16.4 Periodicity of 8, p216
1462 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1463 auto
1466 {
1467 // We add 4 to q by subtracting 4 from p
1468 if (q+4 > -LO)
1469 throw error_t("centre_pm4_qp4(p,q): LO is too high to represent this value");
1470 if (this->frame().max() > p-4)
1471 {
1472 using index_pair_t = typename index_set_t::index_pair_t;
1473 const auto pm3210 = index_set_t(index_pair_t(p-3,p), true);
1474 const auto qm4321 = index_set_t(index_pair_t(-q-4,-q-1), true);
1475 const auto& tqm4321 = term_t(qm4321, Scalar_T(1));
1476 auto result = multivector_t();
1477 for (auto& this_term : *this)
1478 {
1479 const auto ist = this_term.first;
1480 if (ist.max() > p-4)
1481 {
1482 auto var_term = var_term_t();
1483 for (auto
1484 n = index_t(0);
1485 n != index_t(4);
1486 ++n)
1487 if (ist[n+p-3])
1488 var_term *= term_t(index_set_t(n-q-4), Scalar_T(1)) * tqm4321;
1489 // Mask out {p-3}..{p}
1490 result.insert(term_t(ist & ~pm3210, this_term.second) *
1491 term_t(var_term.first, var_term.second));
1492 }
1493 else
1494 result.insert(this_term);
1495 }
1496 *this = result;
1497 }
1498 p -=4; q += 4;
1499 return *this;
1500 }
1501
1503 // Reference: [L] 16.4 Periodicity of 8, p216
1504 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1505 auto
1508 {
1509 // We add 4 to p by subtracting 4 from q
1510 if (p+4 > HI)
1511 throw error_t("centre_pp4_qm4(p,q): HI is too low to represent this value");
1512 if (this->frame().min() < -q+4)
1513 {
1514 using index_pair_t = typename index_set_t::index_pair_t;
1515 const auto qp0123 = index_set_t(index_pair_t(-q,-q+3), true);
1516 const auto pp1234 = index_set_t(index_pair_t(p+1,p+4), true);
1517 const auto& tpp1234 = term_t(pp1234, Scalar_T(1));
1518 auto result = multivector_t();
1519 for (auto& this_term : *this)
1520 {
1521 index_set_t ist = this_term.first;
1522 if (ist.min() < -q+4)
1523 {
1524 auto var_term = var_term_t();
1525 for (auto
1526 n = index_t(0);
1527 n != index_t(4);
1528 ++n)
1529 if (ist[n-q])
1530 var_term *= term_t(index_set_t(n+p+1), Scalar_T(1)) * tpp1234;
1531 // Mask out {-q}..{-q+3}
1532 result.insert(term_t(var_term.first, var_term.second) *
1533 term_t(ist & ~qp0123, this_term.second));
1534 }
1535 else
1536 result.insert(this_term);
1537 }
1538 *this = result;
1539 }
1540 p +=4; q -= 4;
1541 return *this;
1542 }
1543
1545 // Reference: [P] Proposition 15.20, p 131
1546 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1547 auto
1550 {
1551 if (q+1 > HI)
1552 throw error_t("centre_qp1_pm1(p,q): HI is too low to represent this value");
1553 if (p-1 > -LO)
1554 throw error_t("centre_qp1_pm1(p,q): LO is too high to represent this value");
1555 const auto qp1 = index_set_t(q+1);
1556 const auto& tqp1 = term_t(qp1, Scalar_T(1));
1557 auto result = multivector_t();
1558 for (auto& this_term : *this)
1559 {
1560 const auto ist = this_term.first;
1561 auto var_term = var_term_t(index_set_t(), this_term.second);
1562 for (auto
1563 n = -q;
1564 n != p;
1565 ++n)
1566 if (n != 0 && ist[n])
1567 var_term *= term_t(index_set_t(-n) | qp1, Scalar_T(1));
1568 if (p != 0 && ist[p])
1569 var_term *= tqp1;
1570 result.insert(term_t(var_term.first, var_term.second));
1571 }
1572 index_t orig_p = p;
1573 p = q+1;
1574 q = orig_p-1;
1575 return *this = result;
1576 }
1577
1579 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1580 auto
1582 divide(const index_set_t ist) const -> const framed_pair_t
1583 {
1584 auto quo = multivector_t();
1585 auto rem = multivector_t();
1586 for (auto& this_term : *this)
1587 if ((this_term.first | ist) == this_term.first)
1588 quo.insert(term_t(this_term.first ^ ist, this_term.second));
1589 else
1590 rem.insert(this_term);
1591 return framed_pair_t(quo, rem);
1592 }
1593
1595 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1596 auto
1598 fast(const index_t level, const bool odd) const -> const matrix_t
1599 {
1600 // Assume val is already folded and centred
1601 if (this->empty())
1602 {
1603 using matrix_index_t = typename matrix_multi_t::matrix_index_t;
1604 const auto dim = matrix_index_t(1) << level;
1605 auto result =matrix_t(dim, dim);
1606 result.clear();
1607 return result;
1608 }
1609 if (level == 0)
1610 return matrix::unit<matrix_t>(1) * this->scalar();
1611
1612 using basis_matrix_t = typename matrix_multi_t::basis_matrix_t;
1613 using basis_scalar_t = typename basis_matrix_t::value_type;
1614
1615 const auto& I = matrix::unit<basis_matrix_t>(2);
1616 auto J = basis_matrix_t(2,2,2);
1617 J.clear();
1618 J(0,1) = basis_scalar_t(-1);
1619 J(1,0) = basis_scalar_t( 1);
1620 auto K = J;
1621 K(0,1) = basis_scalar_t( 1);
1622 auto JK = I;
1623 JK(0,0) = basis_scalar_t(-1);
1624
1625 const auto ist_mn = index_set_t(-level);
1626 const auto ist_pn = index_set_t(level);
1627 if (level == 1)
1628 {
1629 if (odd)
1630 return matrix_t(J) * (*this)[ist_mn] + matrix_t(K) * (*this)[ist_pn];
1631 else
1632 return matrix_t(I) * this->scalar() + matrix_t(JK) * (*this)[ist_mn ^ ist_pn];
1633 }
1634 else
1635 {
1636 const auto& pair_mn = this->divide(ist_mn);
1637 const auto& quo_mn = pair_mn.first;
1638 const auto& rem_mn = pair_mn.second;
1639 const auto& pair_quo_mnpn = quo_mn.divide(ist_pn);
1640 const auto& val_mnpn = pair_quo_mnpn.first;
1641 const auto& val_mn = pair_quo_mnpn.second;
1642 const auto& pair_rem_mnpn = rem_mn.divide(ist_pn);
1643 const auto& val_pn = pair_rem_mnpn.first;
1644 const auto& val_1 = pair_rem_mnpn.second;
1645 using matrix::kron;
1646 if (odd)
1647 return - kron(JK, val_1.fast (level-1, 1))
1648 + kron(I, val_mnpn.fast(level-1, 1))
1649 + kron(J, val_mn.fast (level-1, 0))
1650 + kron(K, val_pn.fast (level-1, 0));
1651 else
1652 return kron(I, val_1.fast (level-1, 0))
1653 + kron(JK, val_mnpn.fast(level-1, 0))
1654 + kron(K, val_mn.fast (level-1, 1))
1655 - kron(J, val_pn.fast (level-1, 1));
1656 }
1657 }
1658
1660 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1661 template< typename Other_Scalar_T >
1662 auto
1665 {
1666 // Fold val
1667 auto val = this->fold(frm);
1668 auto p = frm.count_pos();
1669 auto q = frm.count_neg();
1670 const auto bott_offset = gen::offset_to_super[pos_mod(p - q, 8)];
1671 p += std::max(bott_offset,index_t(0));
1672 q -= std::min(bott_offset,index_t(0));
1673 if (p > HI)
1674 throw error_t("fast_matrix_multi(frm): HI is too low to represent this value");
1675 if (q > -LO)
1676 throw error_t("fast_matrix_multi(frm): LO is too high to represent this value");
1677 // Centre val
1678 while (p - q > 4)
1679 val.centre_pm4_qp4(p, q);
1680 while (p - q < -3)
1681 val.centre_pp4_qm4(p, q);
1682 if (p - q > 1)
1683 val.centre_qp1_pm1(p, q);
1684 const index_t level = (p + q)/2;
1685
1686 // Do the fast transform
1687 const auto& ev_val = val.even();
1688 const auto& od_val = val.odd();
1689 return matrix_multi<Other_Scalar_T,LO,HI,Tune_P>(ev_val.fast(level, 0) + od_val.fast(level, 1), frm);
1690 }
1691
1692 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1693 inline
1694 auto
1698
1700 template< typename Scalar_T, const index_t LO, const index_t HI >
1701 inline
1702 static
1703 auto
1704 crd_of_mult(const std::pair<const index_set<LO,HI>, Scalar_T>& lhs,
1705 const std::pair<const index_set<LO,HI>, Scalar_T>& rhs) -> Scalar_T
1706 { return lhs.first.sign_of_mult(rhs.first) * lhs.second * rhs.second; }
1707
1709 template< typename Scalar_T, const index_t LO, const index_t HI >
1710 inline
1711 auto
1712 operator* (const std::pair<const index_set<LO,HI>, Scalar_T>& lhs,
1713 const std::pair<const index_set<LO,HI>, Scalar_T>& rhs) -> const std::pair<const index_set<LO,HI>, Scalar_T>
1714 {
1715 using term_t = std::pair<const index_set<LO,HI>, Scalar_T>;
1716 return term_t(lhs.first ^ rhs.first, crd_of_mult(lhs, rhs));
1717 }
1718
1720 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1721 auto
1723 {
1724 using traits_t = numeric_traits<Scalar_T>;
1725 if (val.isnan())
1726 return traits_t::NaN();
1727
1728 check_complex(val, i, prechecked);
1729
1730 const auto realval = val.scalar();
1731 if (val == realval)
1732 {
1733 if (realval < Scalar_T(0))
1734 return i * traits_t::sqrt(-realval);
1735 else
1736 return traits_t::sqrt(realval);
1737 }
1738 using matrix_multi_t = typename framed_multi<Scalar_T,LO,HI,Tune_P>::matrix_multi_t;
1739 return sqrt(matrix_multi_t(val), matrix_multi_t(i), prechecked);
1740 }
1741
1743 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1744 auto
1746 {
1747 using traits_t = numeric_traits<Scalar_T>;
1748 if (val.isnan())
1749 return traits_t::NaN();
1750
1751 const auto s = scalar(val);
1752 if (val == s)
1753 return traits_t::exp(s);
1754
1755 const double size = val.size();
1756 const auto frm_count = val.frame().count();
1757 const auto algebra_dim = set_value_t(1) << frm_count;
1758
1759 if( (size * size <= double(algebra_dim)) || (frm_count < Tune_P::mult_matrix_threshold))
1760 {
1761 switch (Tune_P::function_precision)
1762 {
1763 case precision_demoted:
1764 {
1765 using demoted_scalar_t = typename traits_t::demoted::type;
1766 using demoted_multivector_t = framed_multi<demoted_scalar_t,LO,HI,Tune_P>;
1767
1768 const auto& demoted_val = demoted_multivector_t(val);
1769 return clifford_exp(demoted_val);
1770 }
1771 break;
1772 case precision_promoted:
1773 {
1774 using promoted_scalar_t = typename traits_t::promoted::type;
1775 using promoted_multivector_t = framed_multi<promoted_scalar_t,LO,HI,Tune_P>;
1776
1777 const auto& promoted_val = promoted_multivector_t(val);
1778 return clifford_exp(promoted_val);
1779 }
1780 break;
1781 default:
1782 return clifford_exp(val);
1783 }
1784 }
1785 else
1786 {
1787 using matrix_multi_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
1788 return exp(matrix_multi_t(val));
1789 }
1790 }
1791
1793 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1794 auto
1796 {
1797 using traits_t = numeric_traits<Scalar_T>;
1798 if (val == Scalar_T(0) || val.isnan())
1799 return traits_t::NaN();
1800
1801 check_complex(val, i, prechecked);
1802
1803 const auto realval = val.scalar();
1804 if (val == realval)
1805 {
1806 if (realval < Scalar_T(0))
1807 return i * traits_t::pi() + traits_t::log(-realval);
1808 else
1809 return traits_t::log(realval);
1810 }
1811 using matrix_multi_t = typename framed_multi<Scalar_T,LO,HI,Tune_P>::matrix_multi_t;
1812 return log(matrix_multi_t(val), matrix_multi_t(i), prechecked);
1813 }
1814}
1815#endif // _GLUCAT_FRAMED_MULTI_IMP_H
virtual auto isinf() const -> bool=0
Check if a multivector contains any infinite values.
virtual auto frame() const -> const index_set_t=0
Subalgebra generated by all generators of terms of given multivector.
virtual auto isnan() const -> bool=0
Check if a multivector contains any IEEE NaN values.
virtual auto norm() const -> Scalar_T=0
Scalar_T norm == sum of norm of coordinates.
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
typename map_t::const_iterator const_iterator
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.
std::unordered_map< index_set_t, Scalar_T, index_set_hash< LO, HI > > map_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}.
matrix_multi< Scalar_T, LO, HI, Tune_P > matrix_multi_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.
class var_term var_term_t
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.
Definition errors.h:42
Index set class based on std::bitset<> in Gnu standard C++ library.
Definition index_set.h:75
auto count() const -> index_t
Cardinality: Number of indices included in set.
auto min() const -> index_t
Minimum member.
auto is_contiguous() const -> bool
Determine if the index set is contiguous, ie. has no gaps.
auto count_pos() const -> index_t
Number of positive indices included in set.
std::pair< index_t, index_t > index_pair_t
Definition index_set.h:85
auto count_neg() const -> index_t
Number of negative indices included in set.
auto max() const -> index_t
Maximum member.
A matrix_multi<Scalar_T,LO,HI,Tune_P> is a matrix approximation to a multivector.
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.
Definition scalar.h:48
static auto to_scalar_t(const Other_Scalar_T &val) -> Scalar_T
Cast to Scalar_T.
Definition scalar.h:141
static auto to_double(const Scalar_T &val) -> double
Cast to double.
Definition scalar.h:133
Random number generator with single instance per Scalar_T.
Definition random.h:43
sorted_range(Sorted_Map_T &sorted_val, const Sorted_Map_T &val)
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)
#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.
Definition generation.h:86
auto inner(const LHS_T &lhs, const RHS_T &rhs) -> Scalar_T
Inner product: sum(x(i,j)*y(i,j))/x.nrows()
Definition matrix_imp.h:368
auto isinf(const Matrix_T &m) -> bool
Infinite.
Definition matrix_imp.h:270
auto nnz(const Matrix_T &m) -> typename Matrix_T::size_type
Number of non-zeros.
Definition matrix_imp.h:253
auto kron(const LHS_T &lhs, const RHS_T &rhs) -> const RHS_T
Kronecker tensor product of matrices - as per Matlab kron.
Definition matrix_imp.h:78
auto unit(const typename Matrix_T::size_type n) -> const Matrix_T
Unit matrix - as per Matlab eye.
Definition matrix_imp.h:305
auto isnan(const Matrix_T &m) -> bool
Not a Number.
Definition matrix_imp.h:287
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.
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.
Definition global.h:117
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.
unsigned long set_value_t
Size of set_value_t should be enough to contain index_set<LO,HI>
Definition global.h:79
auto operator>>(std::istream &s, framed_multi< Scalar_T, LO, HI, Tune_P > &val) -> std::istream &
Read multivector from input.
int index_t
Size of index_t should be enough to represent LO, HI.
Definition global.h:77
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 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()