glucat 0.12.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 ;
171 }
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)
179 {
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 }
189
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 }
201
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 { }
227
228 const auto val_norm = traits_t::to_scalar_t(val.norm());
229 if (traits_t::isNaN_or_isInf(val_norm))
230 {
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;
266 }
267 return true;
268 }
269
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
276 {
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 using map_t = typename multivector_t::map_t;
611
612 const auto lhs_end = lhs.end();
613 const auto rhs_end = rhs.end();
614 const double lhs_size = lhs.size();
615 const double rhs_size = rhs.size();
616 const auto lhs_frame = lhs.frame();
617 const auto rhs_frame = rhs.frame();
618
619 const auto our_frame = lhs_frame | rhs_frame;
620 const auto algebra_dim = set_value_t(1) << our_frame.count();
621 auto result = multivector_t(
622 _GLUCAT_HASH_SIZE_T(size_t(std::min(lhs_size * rhs_size, double(algebra_dim)))));
623
624 if (lhs_size * rhs_size > double(Tune_P::products_size_threshold))
625 {
626 for (auto
627 result_stv = set_value_t(0);
628 result_stv != algebra_dim;
629 ++result_stv)
630 {
631 const auto result_ist = index_set_t(result_stv, our_frame, true);
632 const auto comp_frame = lhs_frame & ~result_ist;
633 const auto comp_dim = set_value_t(1) << comp_frame.count();
634 auto result_crd = Scalar_T(0);
635 for (auto
636 comp_stv = set_value_t(0);
637 comp_stv != comp_dim;
638 ++comp_stv)
639 {
640 const auto comp_ist = index_set_t(comp_stv, comp_frame, true);
641 const auto rhs_ist = result_ist ^ comp_ist;
642 if ((rhs_ist | rhs_frame) == rhs_frame)
643 {
644 const auto rhs_it = rhs.find(rhs_ist);
645 if (rhs_it != rhs_end)
646 {
647 const auto lhs_it = lhs.find(comp_ist);
648 if (lhs_it != lhs_end)
649 result_crd += crd_of_mult(*lhs_it, *rhs_it);
650 }
651 }
652 }
653 if (result_crd != Scalar_T(0))
654 result.insert(term_t(result_ist, result_crd));
655 }
656 }
657 else
658 {
659 for (auto& rhs_term : rhs)
660 {
661 const auto rhs_ist = rhs_term.first;
662 for (auto& lhs_term : lhs)
663 {
664 const index_set_t lhs_ist = lhs_term.first;
665 if ((lhs_ist | rhs_ist) == rhs_ist)
666 result += lhs_term * rhs_term;
667 }
668 }
669 }
670 return result;
671 }
672
674 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
675 inline
676 auto
678 operator%= (const multivector_t& rhs) -> multivector_t&
679 { return *this = *this % rhs; }
680
682 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
683 auto
685 {
686 using multivector_t = framed_multi<Scalar_T,LO,HI,Tune_P>;
687
688 auto result = Scalar_T(0);
689 const auto small_star_large = lhs.size() < rhs.size();
690 const auto* smallp =
691 small_star_large
692 ? &lhs
693 : &rhs;
694 const auto* largep =
695 small_star_large
696 ? &rhs
697 : &lhs;
698
699 for (auto& small_term : *smallp)
700 {
701 const auto small_ist = small_term.first;
702 const auto large_crd = (*largep)[small_ist];
703 if (large_crd != Scalar_T(0))
704 result += small_ist.sign_of_square() * small_term.second * large_crd;
705 }
706 return result;
707 }
708
710 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
711 auto
713 operator/= (const Scalar_T& scr) -> multivector_t&
714 { // Divide coordinates of all terms by scr
715 using traits_t = numeric_traits<Scalar_T>;
716
717 if (traits_t::isNaN(scr))
718 return *this = traits_t::NaN();
719 if (traits_t::isInf(scr))
720 if (this->isnan())
721 *this = traits_t::NaN();
722 else
723 this->clear();
724 else
725 for (auto& this_term : *this)
726 this_term.second /= scr;
727 return *this;
728 }
729
731 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
732 inline
733 auto
735 {
736 using multivector_t = framed_multi<Scalar_T,LO,HI,Tune_P>;
737 using traits_t = numeric_traits<Scalar_T>;
738 using index_set_t = typename multivector_t::index_set_t;
739 using matrix_multi_t = typename multivector_t::matrix_multi_t;
740
741 if (rhs == Scalar_T(0))
742 return traits_t::NaN();
743
744 const auto our_frame = lhs.frame() | rhs.frame();
745 return matrix_multi_t(lhs, our_frame, true) / matrix_multi_t(rhs, our_frame, true);
746 }
747
749 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
750 inline
751 auto
753 operator/= (const multivector_t& rhs) -> multivector_t&
754 { return *this = *this / rhs; }
755
757 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
758 inline
759 auto
761 {
762 using multivector_t = framed_multi<Scalar_T,LO,HI,Tune_P>;
763 using matrix_multi_t = typename multivector_t::matrix_multi_t;
764
765 return matrix_multi_t(rhs) * matrix_multi_t(lhs) / matrix_multi_t(rhs.involute());
766 }
767
769 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
770 inline
771 auto
773 operator|= (const multivector_t& rhs) -> multivector_t&
774 { return *this = *this | rhs; }
775
777 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
778 inline
779 auto
781 inv() const -> const multivector_t
782 {
783 auto result = matrix_multi_t(Scalar_T(1), this->frame());
784 return result /= matrix_multi_t(*this);
785 }
786
788 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
789 auto
791 pow(int m) const -> const multivector_t
792 { return glucat::pow(*this, m); }
793
795 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
796 auto
798 outer_pow(int m) const -> const multivector_t
799 {
800 if (m < 0)
801 throw error_t("outer_pow(int): negative exponent");
802 auto result = multivector_t(Scalar_T(1));
803 auto a = *this;
804 for (;
805 m != 0;
806 m >>= 1, a = a ^ a)
807 if (m & 1)
808 result ^= a;
809 return result;
810 }
811
813 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
814 inline
815 auto
817 frame() const -> const index_set_t
818 {
819 auto result = index_set_t();
820 for (auto& this_term : *this)
821 result |= this_term.first;
822 return result;
823 }
824
826 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
827 inline
828 auto
830 grade() const -> index_t
831 {
832 auto result = index_t(0);
833 for (auto& this_term : *this)
834 result = std::max(result, this_term.first.count());
835 return result;
836 }
837
839 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
840 inline
841 auto
843 operator[] (const index_set_t ist) const -> Scalar_T
844 {
845 const auto& this_it = this->find(ist);
846 if (this_it == this->end())
847 return Scalar_T(0);
848 else
849 return this_it->second;
850 }
851
853 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
854 auto
856 operator() (index_t grade) const -> const multivector_t
857 {
858 if ((grade < 0) || (grade > HI-LO))
859 return Scalar_T(0);
860 else
861 {
862 auto result = multivector_t();
863 for (auto& this_term : *this)
864 if (this_term.first.count() == grade)
865 result += this_term;
866 return result;
867 }
868 }
869
871 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
872 inline
873 auto
875 scalar() const -> Scalar_T
876 { return (*this)[index_set_t()]; }
877
879 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
880 inline
881 auto
883 pure() const -> const multivector_t
884 { return *this - this->scalar(); }
885
887 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
888 auto
890 even() const -> const multivector_t
891 { // even part of x, sum of the pure(count) with even count
892 auto result = multivector_t();
893 for (auto& this_term : *this)
894 if ((this_term.first.count() % 2) == 0)
895 result.insert(this_term);
896 return result;
897 }
898
900 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
901 auto
903 odd() const -> const multivector_t
904 { // even part of x, sum of the pure(count) with even count
905 auto result = multivector_t();
906 for (auto& this_term : *this)
907 if ((this_term.first.count() % 2) == 1)
908 result.insert(this_term);
909 return result;
910 }
911
913 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
914 auto
916 vector_part() const -> const vector_t
917 { return this->vector_part(this->frame(), true); }
918
920 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
921 auto
923 vector_part(const index_set_t frm, const bool prechecked) const -> const vector_t
924 {
925 if (!prechecked && (this->frame() | frm) != frm)
926 throw error_t("vector_part(frm): value is outside of requested frame");
927 auto result = vector_t();
928 result.reserve(frm.count());
929 const auto frm_end = frm.max()+1;
930 for (auto
931 idx = frm.min();
932 idx != frm_end;
933 ++idx)
934 // Frame may contain indices which do not correspond to a grade 1 term but
935 // frame cannot omit any index corresponding to a grade 1 term
936 if (frm[idx])
937 result.push_back((*this)[index_set_t(idx)]);
938 return result;
939 }
940
942 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
943 auto
945 involute() const -> const multivector_t
946 {
947 auto result = *this;
948 for (auto& result_term : result)
949 { // for a k-vector u, involute(u) == (-1)^k * u
950 if ((result_term.first.count() % 2) == 1)
951 result_term.second *= Scalar_T(-1);
952 }
953 return result;
954 }
955
957 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
958 auto
960 reverse() const -> const multivector_t
961 {
962 auto result = *this;
963 for (auto& result_term : result)
964 // For a k-vector u, reverse(u) = { -u, k == 2,3 (mod 4)
965 // { u, k == 0,1 (mod 4)
966 switch (result_term.first.count() % 4)
967 {
968 case 2:
969 case 3:
970 result_term.second *= Scalar_T(-1);
971 break;
972 default:
973 break;
974 }
975 return result;
976 }
977
979 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
980 auto
982 conj() const -> const multivector_t
983 {
984 auto result = *this;
985 for (auto& result_term : result)
986 // For a k-vector u, conj(u) = { -u, k == 1,2 (mod 4)
987 // { u, k == 0,3 (mod 4)
988 switch (result_term.first.count() % 4)
989 {
990 case 1:
991 case 2:
992 result_term.second *= Scalar_T(-1);
993 break;
994 default:
995 break;
996 }
997 return result;
998 }
999
1001 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1002 auto
1004 quad() const -> Scalar_T
1005 {
1006 // scalar(conj(x)*x) = 2*quad(even(x)) - quad(x)
1007 // ref: old clical: quadfunction(p:pter):pterm in file compmod.pas
1008 auto result = Scalar_T(0);
1009 for (auto& this_term : *this)
1010 {
1011 const auto sign =
1012 (this_term.first.count_neg() % 2)
1013 ? -Scalar_T(1)
1014 : Scalar_T(1);
1015 result += sign * (this_term.second) * (this_term.second);
1016 }
1017 return result;
1018 }
1019
1021 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1022 auto
1024 norm() const -> Scalar_T
1025 {
1026 using traits_t = numeric_traits<Scalar_T>;
1027
1028 auto result = Scalar_T(0);
1029 for (auto& this_term : *this)
1030 {
1031 const auto abs_crd = traits_t::abs(this_term.second);
1032 result += abs_crd * abs_crd;
1033 }
1034 return result;
1035 }
1036
1038 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1039 auto
1041 max_abs() const -> Scalar_T
1042 {
1043 using traits_t = numeric_traits<Scalar_T>;
1044
1045 auto result = Scalar_T(0);
1046 for (auto& this_term : *this)
1047 {
1048 const auto abs_crd = traits_t::abs(this_term.second);
1049 if (abs_crd > result)
1050 result = abs_crd;
1051 }
1052 return result;
1053 }
1054
1056 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1057 auto
1059 random(const index_set_t frm, Scalar_T fill) -> const multivector_t
1060 {
1062 using index_set_t = typename multivector_t::index_set_t;
1063 using term_t = typename multivector_t::term_t;
1064
1065 using random_generator_t = random_generator<Scalar_T>;
1066 auto& generator = random_generator_t::generator();
1067
1068 fill =
1069 (fill < Scalar_T(0))
1070 ? Scalar_T(0)
1071 : (fill > Scalar_T(1))
1072 ? Scalar_T(1)
1073 : fill;
1074 const auto algebra_dim = set_value_t(1) << frm.count();
1075 using traits_t = numeric_traits<Scalar_T>;
1076 const auto mean_abs = traits_t::sqrt(Scalar_T(double(algebra_dim)));
1077 auto result = multivector_t();
1078 for (auto
1079 stv = set_value_t(0);
1080 stv != algebra_dim;
1081 ++stv)
1082 if (generator.uniform() < fill)
1083 {
1084 const auto& result_crd = generator.normal() / mean_abs;
1085 result.insert(term_t(index_set_t(stv, frm, true), result_crd));
1086 }
1087 return result;
1088 }
1089
1091 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1092 inline
1093 void
1095 write(const std::string& msg) const
1096 { std::cout << msg << std::endl << " " << (*this) << std::endl; }
1097
1099 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1100 inline
1101 void
1103 write(std::ofstream& ofile, const std::string& msg) const
1104 {
1105 if (!ofile)
1106 throw error_t("write(ofile,msg): cannot write to output file");
1107 ofile << msg << std::endl << " " << (*this) << std::endl;
1108 }
1109
1111 template< typename Map_T,typename Sorted_Map_T >
1113 {
1114 public:
1115 using map_t = Map_T;
1116 using sorted_map_t = Sorted_Map_T;
1117 using sorted_iterator = typename Sorted_Map_T::const_iterator;
1118
1119 sorted_range (Sorted_Map_T &sorted_val, const Map_T& val)
1120 {
1121 for (auto& val_term : val)
1122 sorted_val.insert(val_term);
1123 sorted_begin = sorted_val.begin();
1124 sorted_end = sorted_val.end();
1125 }
1128 };
1129
1130 template< typename Sorted_Map_T >
1131 class sorted_range< Sorted_Map_T, Sorted_Map_T >
1132 {
1133 public:
1134 using map_t = Sorted_Map_T;
1135 using sorted_map_t = Sorted_Map_T;
1136 using sorted_iterator = typename Sorted_Map_T::const_iterator;
1137
1138 sorted_range (Sorted_Map_T &sorted_val, const Sorted_Map_T& val)
1139 : sorted_begin( val.begin() ),
1140 sorted_end( val.end() )
1141 { }
1144 };
1145
1147 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1148 auto
1149 operator<< (std::ostream& os, const framed_multi<Scalar_T,LO,HI,Tune_P>& val) -> std::ostream&
1150 {
1151 using limits_t = std::numeric_limits<Scalar_T>;
1152 if (val.empty())
1153 os << 0;
1154 else if (val.isnan())
1155 os << limits_t::quiet_NaN();
1156 else if (val.isinf())
1157 {
1158 const Scalar_T& inf = limits_t::infinity();
1159 os << (scalar(val) < 0.0 ? -inf : inf);
1160 }
1161 else
1162 {
1163 using traits_t = numeric_traits<Scalar_T>;
1164 using multivector_t = framed_multi<Scalar_T,LO,HI,Tune_P>;
1165 Scalar_T truncation;
1166 switch (os.flags() & std::ios::floatfield)
1167 {
1168 case std::ios_base::scientific:
1169 truncation = Scalar_T(1) / traits_t::pow(Scalar_T(10), int(os.precision()) + 1);
1170 break;
1171 case std::ios_base::fixed:
1172 truncation = Scalar_T(1) / (traits_t::pow(Scalar_T(10), int(os.precision())) * val.max_abs());
1173 break;
1174 case std::ios_base::fixed | std::ios_base::scientific:
1175 truncation = multivector_t::default_truncation;
1176 break;
1177 default:
1178 truncation = Scalar_T(1) / traits_t::pow(Scalar_T(10), int(os.precision()));
1179 break;
1180 }
1181 auto truncated_val = val.truncated(truncation);
1182 if (truncated_val.empty())
1183 os << 0;
1184 else
1185 {
1186 using map_t = typename multivector_t::map_t;
1187 using sorted_map_t = typename multivector_t::sorted_map_t;
1188 using sorted_iterator = typename sorted_map_t::const_iterator;
1189 auto sorted_val = sorted_map_t();
1190 const auto sorted_val_range = sorted_range< map_t, sorted_map_t >(sorted_val, truncated_val);
1191 auto sorted_it = sorted_val_range.sorted_begin;
1192 os << *sorted_it;
1193 for (++sorted_it;
1194 sorted_it != sorted_val_range.sorted_end;
1195 ++sorted_it)
1196 {
1197 const Scalar_T& scr = sorted_it->second;
1198 if (scr >= 0.0)
1199 os << '+';
1200 os << *sorted_it;
1201 }
1202 }
1203 }
1204 return os;
1205 }
1206
1208 template< typename Scalar_T, const index_t LO, const index_t HI >
1209 auto
1210 operator<< (std::ostream& os, const std::pair< const index_set<LO,HI>, Scalar_T >& term) -> std::ostream&
1211 {
1212 const auto second_as_double = numeric_traits<Scalar_T>::to_double(term.second);
1213 const auto use_double =
1214 (os.precision() <= std::numeric_limits<double>::digits10) ||
1215 (term.second == Scalar_T(second_as_double));
1216 if (term.first.count() == 0)
1217 if (use_double)
1218 os << second_as_double;
1219 else
1220 os << term.second;
1221 else if (term.second == Scalar_T(-1))
1222 {
1223 os << '-';
1224 os << term.first;
1225 }
1226 else if (term.second != Scalar_T(1))
1227 {
1228 if (use_double)
1229 {
1230 auto tol = std::pow(10.0,-os.precision());
1231 if ( std::fabs(second_as_double + 1.0) < tol )
1232 os << '-';
1233 else if ( std::fabs(second_as_double - 1.0) >= tol )
1234 os << second_as_double;
1235 }
1236 else
1237 os << term.second;
1238 os << term.first;
1239 }
1240 else
1241 os << term.first;
1242 return os;
1243 }
1244
1246 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1247 auto
1248 operator>> (std::istream& s, framed_multi<Scalar_T,LO,HI,Tune_P> & val) -> std::istream&
1249 { // Input looks like 1.0-2.0{1,2}+3.2{3,4}.
1250 using multivector_t = framed_multi<Scalar_T,LO,HI,Tune_P>;
1251 // Parsing variables.
1252 auto local_val = multivector_t();
1253 auto c = 0;
1254 // Parsing control variables.
1255 auto negative = false;
1256 auto expect_term = true;
1257 // The multivector may begin with '+' or '-'. Check for this.
1258 c = s.peek();
1259 if (s.good() && (c == int('+') || c == int('-')))
1260 { // A '-' here negates the following term.
1261 negative = (c == int('-'));
1262 // Consume the '+' or '-'.
1263 s.get();
1264 }
1265 while (s.good())
1266 { // Parse a term.
1267 // A term consists of an optional scalar, followed by an optional index set.
1268 // At least one of the two must be present.
1269 // Default coordinate is Scalar_T(1).
1270 auto coordinate = Scalar_T(1);
1271 // Default index set is empty.
1272 auto ist = index_set<LO,HI>();
1273 // First, check for an opening brace.
1274 c = s.peek();
1275 if (s.good())
1276 { // If the character is not an opening brace,
1277 // a coordinate value is expected here.
1278 if (c != int('{'))
1279 { // Try to read a coordinate value.
1280 double coordinate_as_double;
1281 s >> coordinate_as_double;
1282 // Reading the coordinate may have resulted in an end of file condition.
1283 // This is not a failure.
1284 if (s)
1285 coordinate = Scalar_T(coordinate_as_double);
1286 }
1287 }
1288 else
1289 { // End of file here ends parsing while a term may still be expected.
1290 break;
1291 }
1292 // Coordinate is now Scalar_T(1) or a Scalar_T value.
1293 // Parse an optional index set.
1294 if (s.good())
1295 {
1296 c = s.peek();
1297 if (s.good() && c == int('{'))
1298 { // Try to read index set.
1299 s >> ist;
1300 }
1301 }
1302 // Reading the term may have resulted in an end of file condition.
1303 // This is not a failure.
1304 if (s)
1305 {
1306 // Immediately after parsing a term, another term is not expected.
1307 expect_term = false;
1308 if (coordinate != Scalar_T(0))
1309 {
1310 // Add the term to the local multivector.
1311 coordinate =
1312 negative
1313 ? -coordinate
1314 : coordinate;
1315 using term_t = typename multivector_t::term_t;
1316 local_val += term_t(ist, coordinate);
1317 }
1318 }
1319 // Check if anything follows the current term.
1320 if (s.good())
1321 {
1322 c = s.peek();
1323 if (s.good())
1324 { // Only '+' and '-' are valid here.
1325 if (c == int('+') || c == int('-'))
1326 { // A '-' here negates the following term.
1327 negative = (c == int('-'));
1328 // Consume the '+' or '-'.
1329 s.get();
1330 // Immediately after '+' or '-',
1331 // expect another term.
1332 expect_term = true;
1333 }
1334 else
1335 { // Any other character here is a not failure,
1336 // but still ends the parsing of the multivector.
1337 break;
1338 }
1339 }
1340 }
1341 }
1342 // If a term is still expected, this is a failure.
1343 if (expect_term)
1344 s.clear(std::istream::failbit);
1345 // End of file is not a failure.
1346 if (s)
1347 { // The multivector has been successfully parsed.
1348 val = local_val;
1349 }
1350 return s;
1351 }
1352
1354 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1355 auto
1357 nbr_terms () const -> unsigned long
1358 { return this->size(); }
1359
1361 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1362 inline
1363 auto
1365 operator+= (const term_t& term) -> multivector_t&
1366 { // Do not insert terms with 0 coordinate
1367 if (term.second != Scalar_T(0))
1368 {
1369 const auto& this_it = this->find(term.first);
1370 if (this_it == this->end())
1371 this->insert(term);
1372 else if (this_it->second + term.second == Scalar_T(0))
1373 // Erase term if resulting coordinate is 0
1374 this->erase(this_it);
1375 else
1376 this_it->second += term.second;
1377 }
1378 return *this;
1379 }
1380
1382 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1383 auto
1385 isinf() const -> bool
1386 {
1387 using traits_t = numeric_traits<Scalar_T>;
1388
1389 if (std::numeric_limits<Scalar_T>::has_infinity)
1390 for (auto& this_term : *this)
1391 if (traits_t::isInf(this_term.second))
1392 return true;
1393 return false;
1394 }
1395
1397 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1398 auto
1400 isnan() const -> bool
1401 {
1402 using traits_t = numeric_traits<Scalar_T>;
1403
1404 if (std::numeric_limits<Scalar_T>::has_quiet_NaN)
1405 for (auto& this_term : *this)
1406 if (traits_t::isNaN(this_term.second))
1407 return true;
1408 return false;
1409 }
1410
1412 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1413 auto
1415 truncated(const Scalar_T& limit) const -> const multivector_t
1416 {
1417 using traits_t = numeric_traits<Scalar_T>;
1418
1419 if (this->isnan() || this->isinf())
1420 return *this;
1421 const auto truncation = traits_t::abs(limit);
1422 const auto top = max_abs();
1423 auto result = multivector_t();
1424 if (top != Scalar_T(0))
1425 for (auto& this_term : *this)
1426 if (traits_t::abs(this_term.second) > top * truncation)
1427 result.insert(this_term);
1428 return result;
1429 }
1430
1432 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1433 auto
1435 fold(const index_set_t frm) const -> multivector_t
1436 {
1437 if (frm.is_contiguous())
1438 return *this;
1439 else
1440 {
1441 auto result = multivector_t();
1442 for (auto& this_term : *this)
1443 result.insert(term_t(this_term.first.fold(frm), this_term.second));
1444 return result;
1445 }
1446 }
1447
1449 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1450 auto
1452 unfold(const index_set_t frm) const -> multivector_t
1453 {
1454 if (frm.is_contiguous())
1455 return *this;
1456 else
1457 {
1458 auto result = multivector_t();
1459 for (auto& this_term : *this)
1460 result.insert(term_t(this_term.first.unfold(frm), this_term.second));
1461 return result;
1462 }
1463 }
1464
1466 // Reference: [L] 16.4 Periodicity of 8, p216
1467 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1468 auto
1471 {
1472 // We add 4 to q by subtracting 4 from p
1473 if (q+4 > -LO)
1474 throw error_t("centre_pm4_qp4(p,q): LO is too high to represent this value");
1475 if (this->frame().max() > p-4)
1476 {
1477 using index_pair_t = typename index_set_t::index_pair_t;
1478 const auto pm3210 = index_set_t(index_pair_t(p-3,p), true);
1479 const auto qm4321 = index_set_t(index_pair_t(-q-4,-q-1), true);
1480 const auto& tqm4321 = term_t(qm4321, Scalar_T(1));
1481 auto result = multivector_t();
1482 for (auto& this_term : *this)
1483 {
1484 const auto ist = this_term.first;
1485 if (ist.max() > p-4)
1486 {
1487 auto var_term = var_term_t();
1488 for (auto
1489 n = index_t(0);
1490 n != index_t(4);
1491 ++n)
1492 if (ist[n+p-3])
1493 var_term *= term_t(index_set_t(n-q-4), Scalar_T(1)) * tqm4321;
1494 // Mask out {p-3}..{p}
1495 result.insert(term_t(ist & ~pm3210, this_term.second) *
1496 term_t(var_term.first, var_term.second));
1497 }
1498 else
1499 result.insert(this_term);
1500 }
1501 *this = result;
1502 }
1503 p -=4; q += 4;
1504 return *this;
1505 }
1506
1508 // Reference: [L] 16.4 Periodicity of 8, p216
1509 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1510 auto
1513 {
1514 // We add 4 to p by subtracting 4 from q
1515 if (p+4 > HI)
1516 throw error_t("centre_pp4_qm4(p,q): HI is too low to represent this value");
1517 if (this->frame().min() < -q+4)
1518 {
1519 using index_pair_t = typename index_set_t::index_pair_t;
1520 const auto qp0123 = index_set_t(index_pair_t(-q,-q+3), true);
1521 const auto pp1234 = index_set_t(index_pair_t(p+1,p+4), true);
1522 const auto& tpp1234 = term_t(pp1234, Scalar_T(1));
1523 auto result = multivector_t();
1524 for (auto& this_term : *this)
1525 {
1526 index_set_t ist = this_term.first;
1527 if (ist.min() < -q+4)
1528 {
1529 auto var_term = var_term_t();
1530 for (auto
1531 n = index_t(0);
1532 n != index_t(4);
1533 ++n)
1534 if (ist[n-q])
1535 var_term *= term_t(index_set_t(n+p+1), Scalar_T(1)) * tpp1234;
1536 // Mask out {-q}..{-q+3}
1537 result.insert(term_t(var_term.first, var_term.second) *
1538 term_t(ist & ~qp0123, this_term.second));
1539 }
1540 else
1541 result.insert(this_term);
1542 }
1543 *this = result;
1544 }
1545 p +=4; q -= 4;
1546 return *this;
1547 }
1548
1550 // Reference: [P] Proposition 15.20, p 131
1551 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1552 auto
1555 {
1556 if (q+1 > HI)
1557 throw error_t("centre_qp1_pm1(p,q): HI is too low to represent this value");
1558 if (p-1 > -LO)
1559 throw error_t("centre_qp1_pm1(p,q): LO is too high to represent this value");
1560 const auto qp1 = index_set_t(q+1);
1561 const auto& tqp1 = term_t(qp1, Scalar_T(1));
1562 auto result = multivector_t();
1563 for (auto& this_term : *this)
1564 {
1565 const auto ist = this_term.first;
1566 auto var_term = var_term_t(index_set_t(), this_term.second);
1567 for (auto
1568 n = -q;
1569 n != p;
1570 ++n)
1571 if (n != 0 && ist[n])
1572 var_term *= term_t(index_set_t(-n) | qp1, Scalar_T(1));
1573 if (p != 0 && ist[p])
1574 var_term *= tqp1;
1575 result.insert(term_t(var_term.first, var_term.second));
1576 }
1577 index_t orig_p = p;
1578 p = q+1;
1579 q = orig_p-1;
1580 return *this = result;
1581 }
1582
1584 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1585 auto
1587 divide(const index_set_t ist) const -> const framed_pair_t
1588 {
1589 auto quo = multivector_t();
1590 auto rem = multivector_t();
1591 for (auto& this_term : *this)
1592 if ((this_term.first | ist) == this_term.first)
1593 quo.insert(term_t(this_term.first ^ ist, this_term.second));
1594 else
1595 rem.insert(this_term);
1596 return framed_pair_t(quo, rem);
1597 }
1598
1600 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1601 auto
1603 fast(const index_t level, const bool odd) const -> const matrix_t
1604 {
1605 // Assume val is already folded and centred
1606 if (this->empty())
1607 {
1608 using matrix_index_t = typename matrix_multi_t::matrix_index_t;
1609 const auto dim = matrix_index_t(1) << level;
1610 auto result =matrix_t(dim, dim);
1611 result.clear();
1612 return result;
1613 }
1614 if (level == 0)
1615 return matrix::unit<matrix_t>(1) * this->scalar();
1616
1617 using basis_matrix_t = typename matrix_multi_t::basis_matrix_t;
1618 using basis_scalar_t = typename basis_matrix_t::value_type;
1619
1620 const auto& I = matrix::unit<basis_matrix_t>(2);
1621 auto J = basis_matrix_t(2,2,2);
1622 J.clear();
1623 J(0,1) = basis_scalar_t(-1);
1624 J(1,0) = basis_scalar_t( 1);
1625 auto K = J;
1626 K(0,1) = basis_scalar_t( 1);
1627 auto JK = I;
1628 JK(0,0) = basis_scalar_t(-1);
1629
1630 const auto ist_mn = index_set_t(-level);
1631 const auto ist_pn = index_set_t(level);
1632 if (level == 1)
1633 {
1634 if (odd)
1635 return matrix_t(J) * (*this)[ist_mn] + matrix_t(K) * (*this)[ist_pn];
1636 else
1637 return matrix_t(I) * this->scalar() + matrix_t(JK) * (*this)[ist_mn ^ ist_pn];
1638 }
1639 else
1640 {
1641 const auto& pair_mn = this->divide(ist_mn);
1642 const auto& quo_mn = pair_mn.first;
1643 const auto& rem_mn = pair_mn.second;
1644 const auto& pair_quo_mnpn = quo_mn.divide(ist_pn);
1645 const auto& val_mnpn = pair_quo_mnpn.first;
1646 const auto& val_mn = pair_quo_mnpn.second;
1647 const auto& pair_rem_mnpn = rem_mn.divide(ist_pn);
1648 const auto& val_pn = pair_rem_mnpn.first;
1649 const auto& val_1 = pair_rem_mnpn.second;
1650 using matrix::kron;
1651 if (odd)
1652 return - kron(JK, val_1.fast (level-1, 1))
1653 + kron(I, val_mnpn.fast(level-1, 1))
1654 + kron(J, val_mn.fast (level-1, 0))
1655 + kron(K, val_pn.fast (level-1, 0));
1656 else
1657 return kron(I, val_1.fast (level-1, 0))
1658 + kron(JK, val_mnpn.fast(level-1, 0))
1659 + kron(K, val_mn.fast (level-1, 1))
1660 - kron(J, val_pn.fast (level-1, 1));
1661 }
1662 }
1663
1665 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1666 template< typename Other_Scalar_T >
1667 auto
1670 {
1671 // Fold val
1672 auto val = this->fold(frm);
1673 auto p = frm.count_pos();
1674 auto q = frm.count_neg();
1675 const auto bott_offset = gen::offset_to_super[pos_mod(p - q, 8)];
1676 p += std::max(bott_offset,index_t(0));
1677 q -= std::min(bott_offset,index_t(0));
1678 if (p > HI)
1679 throw error_t("fast_matrix_multi(frm): HI is too low to represent this value");
1680 if (q > -LO)
1681 throw error_t("fast_matrix_multi(frm): LO is too high to represent this value");
1682 // Centre val
1683 while (p - q > 4)
1684 val.centre_pm4_qp4(p, q);
1685 while (p - q < -3)
1686 val.centre_pp4_qm4(p, q);
1687 if (p - q > 1)
1688 val.centre_qp1_pm1(p, q);
1689 const index_t level = (p + q)/2;
1690
1691 // Do the fast transform
1692 const auto& ev_val = val.even();
1693 const auto& od_val = val.odd();
1694 return matrix_multi<Other_Scalar_T,LO,HI,Tune_P>(ev_val.fast(level, 0) + od_val.fast(level, 1), frm);
1695 }
1696
1697 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1698 inline
1699 auto
1703
1705 template< typename Scalar_T, const index_t LO, const index_t HI >
1706 inline
1707 static
1708 auto
1709 crd_of_mult(const std::pair<const index_set<LO,HI>, Scalar_T>& lhs,
1710 const std::pair<const index_set<LO,HI>, Scalar_T>& rhs) -> Scalar_T
1711 { return lhs.first.sign_of_mult(rhs.first) * lhs.second * rhs.second; }
1712
1714 template< typename Scalar_T, const index_t LO, const index_t HI >
1715 inline
1716 auto
1717 operator* (const std::pair<const index_set<LO,HI>, Scalar_T>& lhs,
1718 const std::pair<const index_set<LO,HI>, Scalar_T>& rhs) -> const std::pair<const index_set<LO,HI>, Scalar_T>
1719 {
1720 using term_t = std::pair<const index_set<LO,HI>, Scalar_T>;
1721 return term_t(lhs.first ^ rhs.first, crd_of_mult(lhs, rhs));
1722 }
1723
1725 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1726 auto
1728 {
1729 using traits_t = numeric_traits<Scalar_T>;
1730 if (val.isnan())
1731 return traits_t::NaN();
1732
1733 check_complex(val, i, prechecked);
1734
1735 const auto realval = val.scalar();
1736 if (val == realval)
1737 {
1738 if (realval < Scalar_T(0))
1739 return i * traits_t::sqrt(-realval);
1740 else
1741 return traits_t::sqrt(realval);
1742 }
1743 using matrix_multi_t = typename framed_multi<Scalar_T,LO,HI,Tune_P>::matrix_multi_t;
1744 return sqrt(matrix_multi_t(val), matrix_multi_t(i), prechecked);
1745 }
1746
1748 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1749 auto
1751 {
1752 using traits_t = numeric_traits<Scalar_T>;
1753 if (val.isnan())
1754 return traits_t::NaN();
1755
1756 const auto s = scalar(val);
1757 if (val == s)
1758 return traits_t::exp(s);
1759
1760 const double size = val.size();
1761 const auto frm_count = val.frame().count();
1762 const auto algebra_dim = set_value_t(1) << frm_count;
1763
1764 if( (size * size <= double(algebra_dim)) || (frm_count < Tune_P::mult_matrix_threshold))
1765 {
1766 switch (Tune_P::function_precision)
1767 {
1768 case precision_demoted:
1769 {
1770 using demoted_scalar_t = typename traits_t::demoted::type;
1771 using demoted_multivector_t = framed_multi<demoted_scalar_t,LO,HI,Tune_P>;
1772
1773 const auto& demoted_val = demoted_multivector_t(val);
1774 return clifford_exp(demoted_val);
1775 }
1776 break;
1777 case precision_promoted:
1778 {
1779 using promoted_scalar_t = typename traits_t::promoted::type;
1780 using promoted_multivector_t = framed_multi<promoted_scalar_t,LO,HI,Tune_P>;
1781
1782 const auto& promoted_val = promoted_multivector_t(val);
1783 return clifford_exp(promoted_val);
1784 }
1785 break;
1786 default:
1787 return clifford_exp(val);
1788 }
1789 }
1790 else
1791 {
1792 using matrix_multi_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
1793 return exp(matrix_multi_t(val));
1794 }
1795 }
1796
1798 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1799 auto
1801 {
1802 using traits_t = numeric_traits<Scalar_T>;
1803 if (val == Scalar_T(0) || val.isnan())
1804 return traits_t::NaN();
1805
1806 check_complex(val, i, prechecked);
1807
1808 const auto realval = val.scalar();
1809 if (val == realval)
1810 {
1811 if (realval < Scalar_T(0))
1812 return i * traits_t::pi() + traits_t::log(-realval);
1813 else
1814 return traits_t::log(realval);
1815 }
1816 using matrix_multi_t = typename framed_multi<Scalar_T,LO,HI,Tune_P>::matrix_multi_t;
1817 return log(matrix_multi_t(val), matrix_multi_t(i), prechecked);
1818 }
1819}
1820#endif // _GLUCAT_FRAMED_MULTI_IMP_H
virtual auto frame() const -> const index_set_t=0
Subalgebra generated by all generators of terms of given multivector.
virtual auto truncated(const Scalar_T &limit=default_truncation) const -> const multivector_t=0
virtual auto norm() const -> Scalar_T=0
Scalar_T norm == sum of norm of coordinates.
Specific exception class.
Definition errors.h:57
A framed_multi<Scalar_T,LO,HI,Tune_P> is a framed approximation to a multivector.
auto divide(const index_set_t ist) const -> const framed_pair_t
Divide multivector into part divisible by index_set and remainder.
auto fast(const index_t level, const bool odd) const -> const matrix_t
Generalized FFT from multivector_t to matrix_t.
auto operator+=(const term_t &term) -> multivector_t &
Add a term, if non-zero.
friend class framed_multi
framed_multi multivector_t
std::vector< Scalar_T > vector_t
auto fast_matrix_multi(const index_set_t frm) const -> const matrix_multi< Other_Scalar_T, LO, HI, Tune_P >
Use generalized FFT to construct a matrix_multi_t.
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}.
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.
std::pair< index_t, index_t > index_pair_t
Definition index_set.h:85
auto max() const -> index_t
Maximum member.
A matrix_multi<Scalar_T,LO,HI,Tune_P> is a matrix approximation to a multivector.
typename matrix_t::size_type matrix_index_t
auto basis_element(const index_set< LO, HI > &ist) const -> const basis_matrix_t
Create a basis element matrix within the current frame.
matrix_t m_matrix
Matrix value representing the multivector within the folded frame.
ublas::compressed_matrix< int, orientation_t > basis_matrix_t
Extra traits which extend numeric limits.
Definition scalar.h:48
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:373
auto isinf(const Matrix_T &m) -> bool
Infinite.
Definition matrix_imp.h:275
auto nnz(const Matrix_T &m) -> typename Matrix_T::size_type
Number of non-zeros.
Definition matrix_imp.h:258
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:83
auto unit(const typename Matrix_T::size_type n) -> const Matrix_T
Unit matrix - as per Matlab eye.
Definition matrix_imp.h:310
auto isnan(const Matrix_T &m) -> bool
Not a Number.
Definition matrix_imp.h:292
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 odd(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Odd part.
auto vector_part(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const std::vector< Scalar_T >
Vector part of multivector, as a vector_t with respect to frame()