glucat 0.13.0
matrix_multi_imp.h
Go to the documentation of this file.
1#ifndef _GLUCAT_MATRIX_MULTI_IMP_H
2#define _GLUCAT_MATRIX_MULTI_IMP_H
3/***************************************************************************
4 GluCat : Generic library of universal Clifford algebra templates
5 matrix_multi_imp.h : Implement the matrix representation of a multivector
6 -------------------
7 begin : Sun 2001-12-09
8 copyright : (C) 2001-2021 by Paul C. Leopardi
9 ***************************************************************************
10
11 This library is free software: you can redistribute it and/or modify
12 it under the terms of the GNU Lesser General Public License as published
13 by the Free Software Foundation, either version 3 of the License, or
14 (at your option) any later version.
15
16 This library is distributed in the hope that it will be useful,
17 but WITHOUT ANY WARRANTY; without even the implied warranty of
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 GNU Lesser General Public License for more details.
20
21 You should have received a copy of the GNU Lesser General Public License
22 along with this library. If not, see <http://www.gnu.org/licenses/>.
23
24 ***************************************************************************
25 This library is based on a prototype written by Arvind Raja and was
26 licensed under the LGPL with permission of the author. See Arvind Raja,
27 "Object-oriented implementations of Clifford algebras in C++: a prototype",
28 in Ablamowicz, Lounesto and Parra (eds.)
29 "Clifford algebras with numeric and symbolic computations", Birkhauser, 1996.
30 ***************************************************************************
31 See also Arvind Raja's original header comments in glucat.h
32 ***************************************************************************/
33
34#include "glucat/matrix_multi.h"
35
36#include "glucat/scalar.h"
37#include "glucat/generation.h"
38#include "glucat/matrix.h"
39
40# if defined(_GLUCAT_GCC_IGNORE_UNUSED_LOCAL_TYPEDEFS)
41# pragma GCC diagnostic push
42# pragma GCC diagnostic ignored "-Wunused-local-typedefs"
43# endif
44# if defined(_GLUCAT_HAVE_BOOST_SERIALIZATION_ARRAY_WRAPPER_H)
45# include <boost/serialization/array_wrapper.hpp>
46# endif
47#include <boost/numeric/ublas/matrix.hpp>
48#include <boost/numeric/ublas/matrix_expression.hpp>
49#include <boost/numeric/ublas/matrix_proxy.hpp>
50#include <boost/numeric/ublas/matrix_sparse.hpp>
51#include <boost/numeric/ublas/operation.hpp>
52#include <boost/numeric/ublas/operation_sparse.hpp>
53#include <boost/numeric/ublas/triangular.hpp>
54#include <boost/numeric/ublas/lu.hpp>
55#include <boost/numeric/ublas/io.hpp>
56# if defined(_GLUCAT_GCC_IGNORE_UNUSED_LOCAL_TYPEDEFS)
57# pragma GCC diagnostic pop
58# endif
59
60#include <fstream>
61#include <iomanip>
62#include <array>
63#include <iostream>
64
65namespace glucat
66{
67 // References for algorithms:
68 // [CHKL]:
69 // [L]: Pertti Lounesto, "Clifford algebras and spinors", Cambridge UP, 1997.
70 // [MB]: Beatrice Meini, "The Matrix Square Root From a New Functional Perspective:
71 // Theoretical Results and Computational Issues", SIAM Journal on
72 // Matrix Analysis and Applications 26(2):362-376, 2004.
73 // [P]: Ian R. Porteous, "Clifford algebras and the classical groups", Cambridge UP, 1995.
74
76 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
77 auto
79 classname() -> const std::string
80 { return "matrix_multi"; }
81
83 // Reference: [P] Table 15.27, p 133
84 inline
85 auto
86 offset_level(const index_t p, const index_t q) -> index_t
87 {
88 // Offsets between the log2 of the matrix dimension for the current signature
89 // and that of the real superalgebra
90 static const std::array<int, 8> offset_log2_dim = {0, 1, 0, 1, 1, 2, 1, 1};
91 const index_t bott = pos_mod(p-q, 8);
92 return (p+q)/2 + offset_log2_dim[bott];
93 }
94
96 // Reference: [P] Table 15.27, p 133
97 template< typename Matrix_Index_T, const index_t LO, const index_t HI >
98 inline
99 static
100 auto
101 folded_dim( const index_set<LO,HI>& sub ) -> Matrix_Index_T
102 { return 1 << offset_level(sub.count_pos(), sub.count_neg()); }
103
105 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
111
113 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
114 template< typename Other_Scalar_T >
117 : m_frame( val.m_frame ), m_matrix( val.m_matrix.size1(), val.m_matrix.size2() )
118 {
119 this->m_matrix.clear();
120 for (auto
121 val_it1 = val.m_matrix.begin1();
122 val_it1 != val.m_matrix.end1();
123 ++val_it1)
124 for (auto
125 val_it2 = val_it1.begin();
126 val_it2 != val_it1.end();
127 ++val_it2)
128 this->m_matrix(val_it2.index1(), val_it2.index2()) = numeric_traits<Scalar_T>::to_scalar_t(*val_it2);
129 }
130
132 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
133 template< typename Other_Scalar_T >
135 matrix_multi(const matrix_multi<Other_Scalar_T,LO,HI,Tune_P>& val, const index_set_t frm, const bool prechecked)
136 : m_frame( frm )
137 {
138 if (frm != val.m_frame)
139 *this = multivector_t(framed_multi_t(val), frm);
140 else
141 {
143 this->m_matrix.resize(dim, dim, false);
144 this->m_matrix.clear();
145 for (auto
146 val_it1 = val.m_matrix.begin1();
147 val_it1 != val.m_matrix.end1();
148 ++val_it1)
149 for (auto
150 val_it2 = val_it1.begin();
151 val_it2 != val_it1.end();
152 ++val_it2)
153 this->m_matrix(val_it2.index1(), val_it2.index2()) = numeric_traits<Scalar_T>::to_scalar_t(*val_it2);
154 }
155 }
156
158 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
160 matrix_multi(const multivector_t& val, const index_set_t frm, const bool prechecked)
161 : m_frame( frm )
162 {
163 if (frm != val.m_frame)
164 *this = multivector_t(framed_multi_t(val), frm);
165 else
166 this->m_matrix = val.m_matrix;
168
170 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
172 matrix_multi(const index_set_t ist, const Scalar_T& crd)
173 : m_frame( ist )
174 {
175 const auto dim = folded_dim<matrix_index_t>(this->m_frame);
176 this->m_matrix.resize(dim, dim, false);
177 this->m_matrix.clear();
178 *this += term_t(ist, crd);
180
181 /// Construct a multivector, within a given frame, from an index set and a scalar coordinate
182 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
184 matrix_multi(const index_set_t ist, const Scalar_T& crd, const index_set_t frm, const bool prechecked)
185 : m_frame( frm )
187 if (!prechecked && (ist | frm) != frm)
188 throw error_t("multivector_t(ist,crd,frm): cannot initialize with value outside of frame");
190 this->m_matrix.resize(dim, dim, false);
191 this->m_matrix.clear();
192 *this += term_t(ist, crd);
194
196 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
198 matrix_multi(const Scalar_T& scr, const index_set_t frm)
199 : m_frame( frm )
200 {
201 const auto dim = folded_dim<matrix_index_t>(frm);
202 this->m_matrix.resize(dim, dim, false);
203 this->m_matrix.clear();
204 *this += term_t(index_set_t(), scr);
205 }
206
207 /// Construct a multivector from an int (within a frame, if given)
208 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
210 matrix_multi(const int scr, const index_set_t frm)
211 { *this = multivector_t(Scalar_T(scr), frm); }
212
213 /// Construct a multivector, within a given frame, from a given vector
214 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
216 matrix_multi(const vector_t& vec,
217 const index_set_t frm, const bool prechecked)
218 : m_frame( frm )
219 {
220 if (!prechecked && index_t(vec.size()) != frm.count())
221 throw error_t("multivector_t(vec,frm): cannot initialize with vector not matching frame");
222 const auto dim = folded_dim<matrix_index_t>(frm);
223 this->m_matrix.resize(dim, dim, false);
224 this->m_matrix.clear();
225 auto idx = frm.min();
226 const auto frm_end = frm.max()+1;
227 for (auto& crd : vec)
229 *this += term_t(index_set_t(idx), crd);
230 for (
231 ++idx;
232 idx != frm_end && !frm[idx];
233 ++idx)
234 ;
235 }
236 }
237
239 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
241 matrix_multi(const std::string& str)
242 { *this = framed_multi_t(str); }
243
245 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
247 matrix_multi(const std::string& str, const index_set_t frm, const bool prechecked)
248 { *this = multivector_t(framed_multi_t(str), frm, prechecked); }
249
251 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
252 template< typename Other_Scalar_T >
255 : m_frame( val.frame() )
256 {
257 if (val.size() >= Tune_P::fast_size_threshold)
258 try
259 {
260 *this = val.template fast_matrix_multi<Scalar_T>(this->m_frame);
261 return;
262 }
263 catch (const glucat_error& e)
264 { }
265 const auto dim = folded_dim<matrix_index_t>(this->m_frame);
266 this->m_matrix.resize(dim, dim, false);
267 this->m_matrix.clear();
268
269 for (auto& val_term : val)
270 *this += val_term;
271 }
274 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
275 template< typename Other_Scalar_T >
277 matrix_multi(const framed_multi<Other_Scalar_T,LO,HI,Tune_P>& framed_val, const index_set_t frm, const bool prechecked)
278 {
279 const auto val = framed_val.truncated();
280 const auto our_frame = val.frame() | frm;
281 if (val.size() >= Tune_P::fast_size_threshold)
282 try
283 {
284 *this = val.template fast_matrix_multi<Scalar_T>(our_frame);
285 return;
286 }
287 catch (const glucat_error& e)
288 { }
289 this->m_frame = our_frame;
290 const auto dim = folded_dim<matrix_index_t>(our_frame);
291 this->m_matrix.resize(dim, dim, false);
292 this->m_matrix.clear();
293
294 for (auto& val_term : val)
295 *this += val_term;
296 }
297
299 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
300 template< typename Matrix_T >
302 matrix_multi(const Matrix_T& mtx, const index_set_t frm)
303 : m_frame( frm ), m_matrix( mtx.size1(), mtx.size2() )
304 {
305 this->m_matrix.clear();
306
307 for (auto
308 mtx_it1 = mtx.begin1();
309 mtx_it1 != mtx.end1();
310 ++mtx_it1)
311 for (auto
312 mtx_it2 = mtx_it1.begin();
313 mtx_it2 != mtx_it1.end();
314 ++mtx_it2)
315 this->m_matrix(mtx_it2.index1(), mtx_it2.index2()) = numeric_traits<Scalar_T>::to_scalar_t(*mtx_it2);
316 }
317
319 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
321 matrix_multi(const matrix_t& mtx, const index_set_t frm)
322 : m_frame( frm ), m_matrix( mtx )
323 { }
324
326 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
327 auto
330 {
331 // Check for assignment to self
332 if (this == &rhs)
333 return *this;
334 this->m_frame = rhs.m_frame;
335 this->m_matrix = rhs.m_matrix;
336 return *this;
337 }
338
340 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
341 inline
342 auto
345 {
346 using index_set_t = index_set<LO, HI>;
347 using multivector_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
348 using framed_multi_t = typename multivector_t::framed_multi_t;
349 // Determine the initial common frame
350 index_set_t our_frame = lhs.m_frame | rhs.m_frame;
351 framed_multi_t framed_lhs;
352 framed_multi_t framed_rhs;
353 if ((lhs.m_frame != our_frame) || (rhs.m_frame != our_frame))
354 {
355 // The common frame may expand as a result of the transform to framed_multi_t
356 framed_lhs = framed_multi_t(lhs);
357 framed_rhs = framed_multi_t(rhs);
358 our_frame |= framed_lhs.frame() | framed_rhs.frame();
359 }
360 // Do the reframing only where necessary
361 if (lhs.m_frame != our_frame)
362 lhs_reframed = multivector_t(framed_lhs, our_frame, true);
363 if (rhs.m_frame != our_frame)
364 rhs_reframed = multivector_t(framed_rhs, our_frame, true);
365 return our_frame;
366 }
367
369 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
370 auto
372 operator== (const multivector_t& rhs) const -> bool
373 {
374 // Ensure that there is no aliasing
375 if (this == &rhs)
376 return true;
377
378 // Operate only within a common frame
379 multivector_t lhs_reframed;
380 multivector_t rhs_reframed;
381 const index_set_t our_frame = reframe(*this, rhs, lhs_reframed, rhs_reframed);
382 const multivector_t& lhs_ref = (this->m_frame == our_frame)
383 ? *this
384 : lhs_reframed;
385 const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
386 ? rhs
387 : rhs_reframed;
388
389 return ublas::norm_inf(lhs_ref.m_matrix - rhs_ref.m_matrix) == 0;
390 }
391
392 // Test for equality of multivector and scalar
393 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
394 inline
395 auto
397 operator== (const Scalar_T& scr) const -> bool
398 {
399 if (scr != Scalar_T(0))
400 return *this == multivector_t(framed_multi_t(scr), this->m_frame, true);
401 else if (ublas::norm_inf(this->m_matrix) != 0)
402 return false;
403 else
404 {
405 const matrix_index_t dim = this->m_matrix.size1();
406 return !(dim == 1 && this->isnan());
407 }
408 }
409
411 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
412 inline
413 auto
415 operator+= (const Scalar_T& scr) -> multivector_t&
416 { return *this += term_t(index_set_t(), scr); }
417
419 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
420 inline
421 auto
423 operator+= (const multivector_t& rhs) -> multivector_t&
424 {
425 // Ensure that there is no aliasing
426 if (this == &rhs)
427 return *this *= Scalar_T(2);
428
429 // Operate only within a common frame
430 multivector_t rhs_reframed;
431 const index_set_t our_frame = reframe(*this, rhs, *this, rhs_reframed);
432 const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
433 ? rhs
434 : rhs_reframed;
435
436 noalias(this->m_matrix) += rhs_ref.m_matrix;
437 return *this;
438 }
439
441 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
442 inline
443 auto
445 operator-= (const Scalar_T& scr) -> multivector_t&
446 { return *this += term_t(index_set_t(), -scr); }
447
449 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
450 inline
451 auto
453 operator-= (const multivector_t& rhs) -> multivector_t&
454 {
455 // Ensure that there is no aliasing
456 if (this == &rhs)
457 return *this = Scalar_T(0);
458
459 // Operate only within a common frame
460 multivector_t rhs_reframed;
461 const index_set_t our_frame = reframe(*this, rhs, *this, rhs_reframed);
462 const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
463 ? rhs
464 : rhs_reframed;
465
466 noalias(this->m_matrix) -= rhs_ref.m_matrix;
467 return *this;
468 }
469
471 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
472 inline
473 auto
475 operator- () const -> const multivector_t
476 { return multivector_t(-(this->m_matrix), this->m_frame); }
477
479 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
480 inline
481 auto
483 operator*= (const Scalar_T& scr) -> multivector_t&
484 { // multiply coordinates of all terms by scalar
485
486 using traits_t = numeric_traits<Scalar_T>;
487 if (traits_t::isNaN_or_isInf(scr) || this->isnan())
488 return *this = traits_t::NaN();
489 if (scr == Scalar_T(0))
490 *this = Scalar_T(0);
491 else
492 this->m_matrix *= scr;
493 return *this;
494 }
495
497 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
498 inline
499 auto
501 {
502 using multivector_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
503 using index_set_t = typename multivector_t::index_set_t;
504
505 if (lhs.isnan() || rhs.isnan())
507
508 // Operate only within a common frame
509 multivector_t lhs_reframed;
510 multivector_t rhs_reframed;
511 const index_set_t our_frame = reframe(lhs, rhs, lhs_reframed, rhs_reframed);
512 const multivector_t& lhs_ref = (lhs.m_frame == our_frame)
513 ? lhs
514 : lhs_reframed;
515 const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
516 ? rhs
517 : rhs_reframed;
518
519 using matrix_t = typename multivector_t::matrix_t;
520 using matrix_index_t = typename matrix_t::size_type;
521
522 const matrix_index_t dim = lhs_ref.m_matrix.size1();
523 multivector_t result = multivector_t(matrix_t(dim, dim), our_frame);
524 result.m_matrix.clear();
525 ublas::axpy_prod(lhs_ref.m_matrix, rhs_ref.m_matrix, result.m_matrix, true);
526 return result;
527 }
528
530 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
531 inline
532 auto
534 operator*= (const multivector_t& rhs) -> multivector_t&
535 { return *this = *this * rhs; }
536
538 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
539 inline
540 auto
542 {
543 using multivector_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
544 using framed_multi_t = typename multivector_t::framed_multi_t;
545 return framed_multi_t(lhs) ^ framed_multi_t(rhs);
546 }
547
549 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
550 inline
551 auto
553 operator^= (const multivector_t& rhs) -> multivector_t&
554 { return *this = *this ^ rhs; }
555
557 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
558 inline
559 auto
561 {
562 using multivector_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
563 using framed_multi_t = typename multivector_t::framed_multi_t;
564 return framed_multi_t(lhs) & framed_multi_t(rhs);
565 }
566
568 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
569 inline
570 auto
572 operator&= (const multivector_t& rhs) -> multivector_t&
573 { return *this = *this & rhs; }
574
576 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
577 inline
578 auto
580 {
581 using multivector_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
582 using framed_multi_t = typename multivector_t::framed_multi_t;
583 return framed_multi_t(lhs) % framed_multi_t(rhs);
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 inline
597 auto
599 { return (lhs * rhs).scalar(); }
600
602 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
603 inline
604 auto
606 operator/= (const Scalar_T& scr) -> multivector_t&
607 { return *this *= Scalar_T(1)/scr; }
608
610 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
611 auto
613 {
614 using traits_t = numeric_traits<Scalar_T>;
615
616 if (lhs.isnan() || rhs.isnan())
617 return traits_t::NaN();
618
619 if (rhs == Scalar_T(0))
620 return traits_t::NaN();
621
622 using multivector_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
623
624 // Operate only within a common frame
625 multivector_t lhs_reframed;
626 multivector_t rhs_reframed;
627 const auto our_frame = reframe(lhs, rhs, lhs_reframed, rhs_reframed);
628 const auto& lhs_ref = (lhs.m_frame == our_frame)
629 ? lhs
630 : lhs_reframed;
631 const auto& rhs_ref = (rhs.m_frame == our_frame)
632 ? rhs
633 : rhs_reframed;
634
635 // Solve result == lhs_ref/rhs_ref <=> result*rhs_ref == lhs_ref
636 // We now solve X == B/A
637 // (where X == result, B == lhs_ref.m_matrix and A == rhs_ref.m_matrix)
638 // X == B/A <=> X*A == B <=> AT*XT == BT
639 // So, we solve AT*XT == BT
640
641 using matrix_t = typename multivector_t::matrix_t;
642 using matrix_index_t = typename matrix_t::size_type;
643
644 const auto& AT = matrix_t(ublas::trans(rhs_ref.m_matrix));
645 auto LU = AT;
646
647 using permutation_t = ublas::permutation_matrix<matrix_index_t>;
648
649 auto pvector = permutation_t(AT.size1());
650 if (! ublas::lu_factorize(LU, pvector))
651 {
652 const auto& BT = matrix_t(ublas::trans(lhs_ref.m_matrix));
653 auto XT = BT;
654 ublas::lu_substitute(LU, pvector, XT);
655 if (matrix::isnan(XT))
656 return traits_t::NaN();
657
658 // Iterative refinement.
659 // Reference: Nicholas J. Higham, "Accuracy and Stability of Numerical Algorithms",
660 // SIAM, 1996, ISBN 0-89871-355-2, Chapter 11
661 if (Tune_P::div_max_steps > 0)
662 {
663 // matrix_t R = ublas::prod(AT, XT) - BT;
664 auto R = matrix_t(-BT);
665 ublas::axpy_prod(AT, XT, R, false);
666 if (matrix::isnan(R))
667 return traits_t::NaN();
668
669 auto nr = Scalar_T(ublas::norm_inf(R));
670 if ( nr != Scalar_T(0) && !traits_t::isNaN_or_isInf(nr) )
671 {
672 auto XTnew = XT;
673 auto nrold = nr + Scalar_T(1);
674 for (auto
675 step = 0;
676 step != Tune_P::div_max_steps &&
677 nr < nrold &&
678 nr != Scalar_T(0) &&
679 nr == nr;
680 ++step)
681 {
682 nrold = nr;
683 if (step != 0)
684 XT = XTnew;
685 auto& D = R;
686 ublas::lu_substitute(LU, pvector, D);
687 XTnew -= D;
688 // noalias(R) = ublas::prod(AT, XTnew) - BT;
689 R = -BT;
690 ublas::axpy_prod(AT, XTnew, R, false);
691 nr = ublas::norm_inf(R);
692 }
693 }
694 }
695 return multivector_t(ublas::trans(XT), our_frame);
696 }
697 else
698 // AT is singular. Return NaN
699 return traits_t::NaN();
700 }
701
703 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
704 inline
705 auto
707 operator/= (const multivector_t& rhs) -> multivector_t&
708 { return *this = *this / rhs; }
709
711 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
712 inline
713 auto
715 { return rhs * lhs / rhs.involute(); }
716
718 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
719 inline
720 auto
722 operator|= (const multivector_t& rhs) -> multivector_t&
723 { return *this = rhs * *this / rhs.involute(); }
724
726 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
727 inline
728 auto
730 inv() const -> const multivector_t
731 { return multivector_t(Scalar_T(1), this->m_frame) / *this; }
732
734 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
735 inline
736 auto
738 pow(int m) const -> const multivector_t
739 { return glucat::pow(*this, m); }
740
742 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
743 auto
745 outer_pow(int m) const -> const multivector_t
746 {
747 if (m < 0)
748 throw error_t("outer_pow(m): negative exponent");
749 framed_multi_t a = *this;
750 return a.outer_pow(m);
751 }
752
754 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
755 inline
756 auto
758 grade() const -> index_t
759 { return framed_multi_t(*this).grade(); }
760
762 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
763 inline
764 auto
766 frame() const -> const index_set_t
767 { return this->m_frame; }
768
770 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
771 inline
772 auto
774 operator[] (const index_set_t ist) const -> Scalar_T
775 {
776 // Use matrix inner product only if ist is in frame
777 if ( (ist | this->m_frame) == this->m_frame)
778 return matrix::inner<Scalar_T>(this->basis_element(ist), this->m_matrix);
779 else
780 return Scalar_T(0);
781 }
782
784 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
785 inline
786 auto
788 operator() (index_t grade) const -> const multivector_t
789 {
790 if ((grade < 0) || (grade > HI-LO))
791 return 0;
792 else
793 return (framed_multi_t(*this))(grade);
794 }
795
797 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
798 inline
799 auto
801 scalar() const -> Scalar_T
802 {
803 const matrix_index_t dim = this->m_matrix.size1();
804 return matrix::trace(this->m_matrix) / Scalar_T( double(dim) );
805 }
806
808 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
809 inline
810 auto
812 pure() const -> const multivector_t
813 { return *this - this->scalar(); }
814
816 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
817 inline
818 auto
820 even() const -> const multivector_t
821 { return framed_multi_t(*this).even(); }
822
824 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
825 inline
826 auto
828 odd() const -> const multivector_t
829 { return framed_multi_t(*this).odd(); }
830
832 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
833 auto
835 vector_part() const -> const vector_t
836 { return this->vector_part(this->frame(), true); }
837
839 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
840 auto
842 vector_part(const index_set_t frm, const bool prechecked) const -> const vector_t
843 {
844 if (!prechecked && (this->frame() | frm) != frm)
845 throw error_t("vector_part(frm): value is outside of requested frame");
846 vector_t result;
847 // If we need to enlarge the frame we may as well use a framed_multi_t
848 if (this->frame() != frm)
849 return framed_multi_t(*this).vector_part(frm, true);
850
851 const auto begin_index = frm.min();
852 const auto end_index = frm.max()+1;
853 for (auto
854 idx = begin_index;
855 idx != end_index;
856 ++idx)
857 if (frm[idx])
858 // Frame may contain indices which do not correspond to a grade 1 term but
859 // frame cannot omit any index corresponding to a grade 1 term
860 result.push_back(
861 matrix::inner<Scalar_T>(this->basis_element(index_set_t(idx)),
862 this->m_matrix));
863 return result;
864 }
865
867 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
868 inline
869 auto
871 involute() const -> const multivector_t
872 { return framed_multi_t(*this).involute(); }
873
875 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
876 inline
877 auto
879 reverse() const -> const multivector_t
880 { return framed_multi_t(*this).reverse(); }
881
883 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
884 inline
885 auto
887 conj() const -> const multivector_t
888 { return framed_multi_t(*this).conj(); }
889
891 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
892 inline
893 auto
895 quad() const -> Scalar_T
896 { // scalar(conj(x)*x) = 2*quad(even(x)) - quad(x)
897 // Arvind Raja ref: "old clical: quadfunction(p:pter):pterm in file compmod.pas"
898 return framed_multi_t(*this).quad();
899 }
900
902 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
903 inline
904 auto
906 norm() const -> Scalar_T
907 {
908 const matrix_index_t dim = this->m_matrix.size1();
909 return matrix::norm_frob2(this->m_matrix) / Scalar_T( double(dim) );
910 }
911
913 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
914 inline
915 auto
917 max_abs() const -> Scalar_T
918 { return framed_multi_t(*this).max_abs(); }
919
921 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
922 auto
928
930 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
931 inline
932 void
934 write(const std::string& msg) const
935 { framed_multi_t(*this).write(msg); }
936
938 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
939 inline
940 void
942 write(std::ofstream& ofile, const std::string& msg) const
943 {
944 if (!ofile)
945 throw error_t("write(ofile,msg): cannot write to output file");
946 framed_multi_t(*this).write(ofile, msg);
947 }
948
950 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
951 inline
952 auto
953 operator<< (std::ostream& os, const matrix_multi<Scalar_T,LO,HI,Tune_P>& val) -> std::ostream&
954 {
955 os << typename matrix_multi<Scalar_T,LO,HI,Tune_P>::framed_multi_t(val);
956 return os;
957 }
958
960 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
961 inline
962 auto
963 operator>> (std::istream& s, matrix_multi<Scalar_T,LO,HI,Tune_P>& val) -> std::istream&
964 { // Input looks like 1.0-2.0{1,2}+3.2{3,4}
966 s >> local;
967 // If s.bad() then we have a corrupt input
968 // otherwise we are fine and can copy the resulting matrix_multi
969 if (!s.bad())
970 val = local;
971 return s;
972 }
973
975 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
976 inline
977 auto
979 isinf() const -> bool
980 {
981 if (std::numeric_limits<Scalar_T>::has_infinity)
982 return matrix::isinf(this->m_matrix);
983 else
984 return false;
985 }
986
988 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
989 inline
990 auto
992 isnan() const -> bool
993 {
994 if (std::numeric_limits<Scalar_T>::has_quiet_NaN)
995 return matrix::isnan(this->m_matrix);
996 else
997 return false;
998 }
999
1001 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1002 inline
1003 auto
1005 truncated(const Scalar_T& limit) const -> const multivector_t
1006 { return framed_multi_t(*this).truncated(limit); }
1007
1009 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1010 inline
1011 auto
1013 operator+= (const term_t& term) -> multivector_t&
1014 {
1015 if (term.second != Scalar_T(0))
1016 this->m_matrix.plus_assign(matrix_t(this->basis_element(term.first)) * term.second);
1017 return *this;
1018 }
1019
1021 template< typename Multivector_T, typename Matrix_T, typename Basis_Matrix_T >
1022 static
1023 auto
1024 fast(const Matrix_T& X, index_t level) -> Multivector_T
1025 {
1026 using framed_multi_t = Multivector_T;
1027
1028 using index_set_t = typename framed_multi_t::index_set_t;
1029 using Scalar_T = typename framed_multi_t::scalar_t;
1030 using matrix_t = Matrix_T;
1031 using basis_matrix_t = Basis_Matrix_T;
1032 using basis_scalar_t = typename basis_matrix_t::value_type;
1033 using traits_t = numeric_traits<Scalar_T>;
1034
1035 if (level == 0)
1036 return framed_multi_t(traits_t::to_scalar_t(X(0,0)));
1037
1038 if (ublas::norm_inf(X) == 0)
1039 return Scalar_T(0);
1040
1041 const basis_matrix_t& I = matrix::unit<basis_matrix_t>(2);
1042 basis_matrix_t J(2,2,2);
1043 J.clear();
1044 J(0,1) = basis_scalar_t(-1);
1045 J(1,0) = basis_scalar_t( 1);
1046 basis_matrix_t K = J;
1047 K(0,1) = basis_scalar_t( 1);
1048 basis_matrix_t JK = I;
1049 JK(0,0) = basis_scalar_t(-1);
1050
1052 const index_set_t ist_mn = index_set_t(-level);
1053 const index_set_t ist_pn = index_set_t(level);
1054 const index_set_t ist_mnpn = ist_mn | ist_pn;
1055 if (level == 1)
1056 {
1057 using term_t = typename framed_multi_t::term_t;
1058 const Scalar_T i_x = traits_t::to_scalar_t(signed_perm_nork(I, X)(0, 0));
1059 const Scalar_T j_x = traits_t::to_scalar_t(signed_perm_nork(J, X)(0, 0));
1060 const Scalar_T k_x = traits_t::to_scalar_t(signed_perm_nork(K, X)(0, 0));
1061 const Scalar_T jk_x = traits_t::to_scalar_t(signed_perm_nork(JK,X)(0, 0));
1062 framed_multi_t
1063 result = i_x;
1064 result += term_t(ist_mn, j_x); // j_x * mn;
1065 result += term_t(ist_pn, k_x); // k_x * pn;
1066 return result += term_t(ist_mnpn, jk_x); // jk_x * mnpn;
1067 }
1068 else
1069 {
1070 const framed_multi_t& mn = framed_multi_t(ist_mn);
1071 const framed_multi_t& pn = framed_multi_t(ist_pn);
1072 const framed_multi_t& mnpn = framed_multi_t(ist_mnpn);
1073 const framed_multi_t& i_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1074 (signed_perm_nork(I, X), level-1);
1075 const framed_multi_t& j_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1076 (signed_perm_nork(J, X), level-1);
1077 const framed_multi_t& k_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1078 (signed_perm_nork(K, X), level-1);
1079 const framed_multi_t& jk_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1080 (signed_perm_nork(JK,X), level-1);
1081 framed_multi_t
1082 result = i_x.even() - jk_x.odd();
1083 result += (j_x.even() - k_x.odd()) * mn;
1084 result += (k_x.even() - j_x.odd()) * pn;
1085 return result += (jk_x.even() - i_x.odd()) * mnpn;
1086 }
1087 }
1088
1090 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1091 inline
1092 auto
1094 fast_matrix_multi(const index_set_t frm) const -> const multivector_t
1095 {
1096 if (this->m_frame == frm)
1097 return *this;
1098 else
1099 return (this->template fast_framed_multi<Scalar_T>()).template fast_matrix_multi<Scalar_T>(frm);
1100 }
1101
1103 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1104 template <typename Other_Scalar_T>
1105 auto
1107 fast_framed_multi() const -> const framed_multi<Other_Scalar_T,LO,HI,Tune_P>
1108 {
1109 // Determine the amount of off-centering needed
1110 index_t p = this->m_frame.count_pos();
1111 index_t q = this->m_frame.count_neg();
1112
1113 const index_t bott = pos_mod(p-q, 8);
1114 p += std::max(gen::offset_to_super[bott],index_t(0));
1115 q -= std::min(gen::offset_to_super[bott],index_t(0));
1116
1117 const index_t orig_p = p;
1118 const index_t orig_q = q;
1119 while (p-q > 4)
1120 { p -= 4; q += 4; }
1121 while (p-q < -3)
1122 { p += 4; q -= 4; }
1123 if (p-q > 1)
1124 {
1125 index_t old_p = p;
1126 p = q+1;
1127 q = old_p-1;
1128 }
1129 const index_t level = (p+q)/2;
1130
1131 // Do the inverse fast transform
1134
1135 // Off-centre val
1136 switch (pos_mod(orig_p-orig_q, 8))
1137 {
1138 case 2:
1139 case 3:
1140 case 4:
1141 val.centre_qp1_pm1(p, q);
1142 break;
1143 default:
1144 break;
1145 }
1146 if (orig_p-orig_q > 4)
1147 while (p != orig_p)
1148 val.centre_pp4_qm4(p, q);
1149 if (orig_p-orig_q < -3)
1150 while (p != orig_p)
1151 val.centre_pm4_qp4(p, q);
1152
1153 // Return unfolded val
1154 return val.unfold(this->m_frame);
1155 }
1156
1158 template< typename Scalar_T, const index_t LO, const index_t HI, typename Matrix_T >
1160 public std::map< std::pair< const index_set<LO,HI>, const index_set<LO,HI> >,
1161 Matrix_T* >
1162 {
1163 public:
1165 static auto basis() -> basis_table& { static basis_table b; return b;}
1166 private:
1171 // Enforce singleton
1172 // Reference: A. Alexandrescu, "Modern C++ Design", Chapter 6
1173 basis_table() = default;
1174 ~basis_table() = default;
1175 public:
1176 basis_table(const basis_table&) = delete;
1177 auto operator= (const basis_table&) -> basis_table& = delete;
1178 };
1179
1181 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1182 auto
1184 basis_element(const index_set_t& ist) const -> const basis_matrix_t
1185 {
1186 using index_set_pair_t = std::pair<const index_set_t, const index_set_t>;
1187 const auto& unfolded_pair = index_set_pair_t(ist, this->m_frame);
1188
1190 auto& basis_cache = basis_table_t::basis();
1191
1192 const auto frame_count = this->m_frame.count();
1193 const auto use_cache = frame_count <= index_t(Tune_P::basis_max_count);
1194
1195 if (use_cache)
1196 {
1197 const auto basis_it = basis_cache.find(unfolded_pair);
1198 if (basis_it != basis_cache.end())
1199 return *(basis_it->second);
1200 }
1201 const auto folded_set = ist.fold(this->m_frame);
1202 const auto folded_frame = this->m_frame.fold();
1203 const auto& folded_pair = index_set_pair_t(folded_set, folded_frame);
1204 using basis_pair_t = std::pair<const index_set_pair_t, basis_matrix_t *>;
1205 if (use_cache)
1206 {
1207 const auto basis_it = basis_cache.find(folded_pair);
1208 if (basis_it != basis_cache.end())
1209 {
1210 auto* result_ptr = basis_it->second;
1211 basis_cache.insert(basis_pair_t(unfolded_pair, result_ptr));
1212 return *result_ptr;
1213 }
1214 }
1215 const auto folded_max = folded_frame.max();
1216 const auto folded_min = folded_frame.min();
1217 const auto p = std::max(folded_max, index_t(0));
1218 const auto q = std::max(index_t(-folded_min), index_t(0));
1219 const auto* e = (gen::generator_table<basis_matrix_t>::generator())(p, q);
1220 const auto dim = matrix_index_t(1) << offset_level(p, q);
1221 auto result = matrix::unit<basis_matrix_t>(dim);
1222 for (auto
1223 k = folded_min;
1224 k <= folded_max;
1225 ++k)
1226 if (folded_set[k])
1227 result = matrix::mono_prod(result, e[k]);
1228 if (use_cache)
1229 {
1230 auto* result_ptr = new basis_matrix_t(result);
1231 basis_cache.insert(basis_pair_t(folded_pair, result_ptr));
1232 basis_cache.insert(basis_pair_t(unfolded_pair, result_ptr));
1233 }
1234 return result;
1235 }
1236
1238 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P, const size_t Size >
1239 inline
1240 static
1241 auto
1243 const std::array<Scalar_T, Size>& numer,
1244 const std::array<Scalar_T, Size>& denom,
1246 {
1247 // Pade' approximation
1248 // Reference: [GW], Section 4.3, pp318-322
1249 // Reference: [GL], Section 11.3, p572-576.
1250
1251 using multivector_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
1252 using traits_t = numeric_traits<Scalar_T>;
1253
1254 if (X.isnan())
1255 return traits_t::NaN();
1256
1257 // Array size is assumed to be even
1258 const auto nbr_even_powers = Size/2 - 1;
1259
1260 // Create an array of even powers
1261 auto XX = std::vector<multivector_t>(nbr_even_powers);
1262 XX[0] = X * X;
1263 XX[1] = XX[0] * XX[0];
1264 for (auto
1265 k = size_t(2);
1266 k != nbr_even_powers;
1267 ++k)
1268 XX[k] = XX[k-2] * XX[1];
1269
1270 // Calculate numerator N and denominator D
1271 auto N = multivector_t(numer[1]);
1272 for (auto
1273 k = size_t(0);
1274 k != nbr_even_powers;
1275 ++k)
1276 N += XX[k] * numer[2*k + 3];
1277 N *= X;
1278 N += numer[0];
1279 for (auto
1280 k = size_t(0);
1281 k != nbr_even_powers;
1282 ++k)
1283 N += XX[k] * numer[2*k + 2];
1284 auto D = multivector_t(denom[1]);
1285 for (auto
1286 k = size_t(0);
1287 k != nbr_even_powers;
1288 ++k)
1289 D += XX[k] * denom[2*k + 3];
1290 D *= X;
1291 D += denom[0];
1292 for (auto
1293 k = size_t(0);
1294 k != nbr_even_powers;
1295 ++k)
1296 D += XX[k] * denom[2*k + 2];
1297 return N / D;
1298 }
1299
1301 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1302 inline
1303 static
1304 void
1306 {
1307 // Reference: [CHKL]
1308 const auto& invM = inv(M);
1309 M = ((M + invM)/Scalar_T(2) + Scalar_T(1)) / Scalar_T(2);
1310 Y *= (invM + Scalar_T(1)) / Scalar_T(2);
1311 }
1312
1314 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1315 static
1316 auto
1318 Scalar_T norm_tol=std::pow(std::numeric_limits<Scalar_T>::epsilon(), 4)) -> const matrix_multi<Scalar_T,LO,HI,Tune_P>
1319 {
1320 // Reference: [CHKL]
1321 if (val == Scalar_T(0))
1322 return val;
1323
1324 static const auto sqrt_max_steps = Tune_P::db_sqrt_max_steps;
1325 auto M = val;
1326 auto Y = val;
1327
1328 for (auto
1329 step = 0;
1330 step != sqrt_max_steps && norm(M - Scalar_T(1)) > norm_tol;
1331 ++step)
1332 {
1333 if (Y.isnan())
1335 db_step(M, Y);
1336 }
1337 return Y;
1338 }
1339
1341 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1342 static
1343 auto
1345 Scalar_T norm_Y_tol=std::pow(std::numeric_limits<Scalar_T>::epsilon(), 1)) -> const matrix_multi<Scalar_T,LO,HI,Tune_P>
1346 {
1347 // Reference: [MB]
1348 if (val == Scalar_T(0))
1349 return val;
1350
1351 static const auto sqrt_max_steps = Tune_P::cr_sqrt_max_steps;
1352 auto Z = Scalar_T(2) * (Scalar_T(1) + val);
1353 auto Y = Scalar_T(1) - val;
1354 auto norm_Y = norm(Y);
1355 for (auto
1356 step = 0;
1357 step != sqrt_max_steps && norm_Y > norm_Y_tol;
1358 ++step)
1359 {
1360 const auto old_norm_Y = norm_Y;
1361 Y = (-Y / Z) * Y;
1362 norm_Y = norm(Y);
1363 if (Y.isnan() || (norm_Y > old_norm_Y * Scalar_T(2)))
1365
1366 Z += Y * Scalar_T(2);
1367 }
1368 return Z / Scalar_T(4);
1369 }
1370}
1371
1372namespace pade {
1374 // Reference: [Z], Pade1
1375 template< typename Scalar_T >
1377 {
1378 using array = std::array<Scalar_T, 14>;
1379 static const array numer;
1380 };
1381 template< typename Scalar_T >
1383 {
1384 1.0, 27.0/4.0, 81.0/4.0, 2277.0/64.0,
1385 10395.0/256.0, 32319.0/1024.0, 8721.0/512.0, 26163.0/4096.0,
1386 53703.0/32768.0, 36465.0/131072.0, 3861.0/131072.0, 7371.0/4194304.0,
1387 819.0/16777216.0, 27.0/67108864.0
1388 };
1389
1391 // Reference: [Z], Pade1
1392 template< typename Scalar_T >
1394 {
1395 using array = std::array<Scalar_T, 14>;
1396 static const array denom;
1397 };
1398 template< typename Scalar_T >
1400 {
1401 1.0, 25.0/4.0, 69.0/4.0, 1771.0/64.0,
1402 7315.0/256.0, 20349.0/1024.0, 4845.0/512.0, 12597.0/4096.0,
1403 21879.0/32768.0, 12155.0/131072.0, 1001.0/131072.0, 1365.0/4194304.0,
1404 91.0/16777216.0, 1.0/67108864.0
1405 };
1406
1407 template< >
1408 struct pade_sqrt_numer<float>
1409 {
1410 using array = std::array<float, 10>;
1411 static const array numer;
1412 };
1414 {
1415 1.0, 19.0/4.0, 19.0/2.0, 665.0/64.0,
1416 1729.0/256.0, 2717.0/1024.0, 627.0/1024.0, 627.0/8192.0,
1417 285.0/65536.0, 19.0/262144.0
1418 };
1419 template< >
1420 struct pade_sqrt_denom<float>
1421 {
1422 using array = std::array<float, 10>;
1423 static const array denom;
1424 };
1426 {
1427 1.0, 17.0/4.0, 15.0/2.0, 455.0/64.0,
1428 1001.0/256.0, 1287.0/1024.0, 231.0/1024.0, 165.0/8192.0,
1429 45.0/65536, 1.0/262144.0
1430 };
1431
1432 template< >
1433 struct pade_sqrt_numer<long double>
1434 {
1435 using array = std::array<long double, 18>;
1436 static const array numer;
1437 };
1439 {
1440 1.0L, 35.0L/4.0L, 35.0L, 5425.0L/64.0L,
1441 35525.0L/256.0L, 166257.0L/1024.0L, 143325.0L/1024.0L, 740025.0L/8192.0L,
1442 2877875.0L/65536.0L, 4206125.0L/262144.0L, 572033.0L/131072.0L, 1820105.0L/2097152.0L,
1443 1028755.0L/8388608.0L, 395675.0L/33554432.0L, 24225.0L/33554432.0L, 6783.0L/268435456.0L,
1444 1785.0L/4294967296.0L, 35.0L/17179869184.0L
1445 };
1446 template< >
1447 struct pade_sqrt_denom<long double>
1448 {
1449 using array = std::array<long double, 18>;
1450 static const array denom;
1451 };
1453 {
1454 1.0L, 33.0L/4.0L, 31.0L, 4495.0L/64.0L,
1455 27405.0L/256.0L, 118755.0L/1024.0L, 94185.0L/1024.0L, 444015.0L/8192.0L,
1456 1562275.0L/65536.0L, 2042975.0L/262144.0L, 245157.0L/131072.0L, 676039.0L/2097152.0L,
1457 323323.0L/8388608.0L, 101745.0L/33554432.0L, 4845.0L/33554432.0L, 969.0L/268435456.0L,
1458 153.0L/4294967296.0L, 1.0L/17179869184.0L
1459 };
1460
1461#if defined(_GLUCAT_USE_QD)
1462 template< >
1463 struct pade_sqrt_numer<dd_real>
1464 {
1465 using array = std::array<dd_real, 22>;
1466 static const array numer;
1467 };
1469 {
1470 dd_real("1"), dd_real("43")/dd_real("4"),
1471 dd_real("215")/dd_real("4"), dd_real("10621")/dd_real("64"),
1472 dd_real("90687")/dd_real("256"), dd_real("567987")/dd_real("1024"),
1473 dd_real("168861")/dd_real("256"), dd_real("1246355")/dd_real("2048"),
1474 dd_real("7228859")/dd_real("16384"), dd_real("16583853")/dd_real("65536"),
1475 dd_real("7538115")/dd_real("65536"), dd_real("173376645")/dd_real("4194304"),
1476 dd_real("195747825")/dd_real("16777216"), dd_real("171655785")/dd_real("67108864"),
1477 dd_real("14375115")/dd_real("33554432"), dd_real("14375115")/dd_real("268435456"),
1478 dd_real("20764055")/dd_real("4294967296"), dd_real("5167525")/dd_real("17179869184"),
1479 dd_real("206701")/dd_real("17179869184"), dd_real("76153")/dd_real("274877906944"),
1480 dd_real("3311")/dd_real("1099511627776") , dd_real("43")/dd_real("4398046511104")
1481 };
1482 template< >
1483 struct pade_sqrt_denom<dd_real>
1484 {
1485 using array = std::array<dd_real, 22>;
1486 static const array denom;
1487 };
1489 {
1490 dd_real("1"), dd_real("41")/dd_real("4"),
1491 dd_real("195")/dd_real("4"), dd_real("9139")/dd_real("64"),
1492 dd_real("73815")/dd_real("256"), dd_real("435897")/dd_real("1024"),
1493 dd_real("121737")/dd_real("256"), dd_real("840565")/dd_real("2048"),
1494 dd_real("4539051")/dd_real("16384"), dd_real("9641775")/dd_real("65536"),
1495 dd_real("4032015")/dd_real("65536"), dd_real("84672315")/dd_real("4194304"),
1496 dd_real("86493225")/dd_real("16777216"), dd_real("67863915")/dd_real("67108864"),
1497 dd_real("5014575")/dd_real("33554432"), dd_real("4345965")/dd_real("268435456"),
1498 dd_real("5311735")/dd_real("4294967296"), dd_real("1081575")/dd_real("17179869184"),
1499 dd_real("33649")/dd_real("17179869184"), dd_real("8855")/dd_real("274877906944"),
1500 dd_real("231")/dd_real("1099511627776"), dd_real("1")/dd_real("4398046511104")
1501 };
1502
1503 template< >
1504 struct pade_sqrt_numer<qd_real>
1505 {
1506 using array = std::array<qd_real, 34>;
1507 static const array numer;
1508 };
1510 {
1511 qd_real("1"), qd_real("67")/qd_real("4"),
1512 qd_real("134"), qd_real("43617")/qd_real("64"),
1513 qd_real("633485")/qd_real("256"), qd_real("6992857")/qd_real("1024"),
1514 qd_real("15246721")/qd_real("1024"), qd_real("215632197")/qd_real("8192"),
1515 qd_real("2518145487")/qd_real("65536"), qd_real("12301285425")/qd_real("262144"),
1516 qd_real("6344873535")/qd_real("131072"), qd_real("89075432355")/qd_real("2097152"),
1517 qd_real("267226297065")/qd_real("8388608"), qd_real("687479618945")/qd_real("33554432"),
1518 qd_real("379874182975")/qd_real("33554432"), qd_real("1443521895305")/qd_real("268435456"),
1519 qd_real("9425348845815")/qd_real("4294967296"), qd_real("13195488384141")/qd_real("17179869184"),
1520 qd_real("987417498133")/qd_real("4294967296"), qd_real("8055248011085")/qd_real("137438953472"),
1521 qd_real("6958363175533")/qd_real("549755813888"), qd_real("5056698705201")/qd_real("2199023255552"),
1522 qd_real("766166470485")/qd_real("2199023255552"), qd_real("766166470485")/qd_real("17592186044416"),
1523 qd_real("623623871325")/qd_real("140737488355328"), qd_real("203123203803")/qd_real("562949953421312"),
1524 qd_real("6478601247")/qd_real("281474976710656"), qd_real("5038912081")/qd_real("4503599627370496"),
1525 qd_real("719844583")/qd_real("18014398509481984"), qd_real("71853815")/qd_real("72057594037927936"),
1526 qd_real("1165197")/qd_real("72057594037927936"), qd_real("87703")/qd_real("576460752303423488"),
1527 qd_real("12529")/qd_real("18446744073709551616"), qd_real("67")/qd_real("73786976294838206464")
1528 };
1529 template< >
1530 struct pade_sqrt_denom<qd_real>
1531 {
1532 using array = std::array<qd_real, 34>;
1533 static const array denom;
1534 };
1536 {
1537 qd_real("1"), qd_real("65")/qd_real("4"),
1538 qd_real("126"), qd_real("39711")/qd_real("64"),
1539 qd_real("557845")/qd_real("256"), qd_real("5949147")/qd_real("1024"),
1540 qd_real("12515965")/qd_real("1024"), qd_real("170574723")/qd_real("8192"),
1541 qd_real("1916797311")/qd_real("65536"), qd_real("8996462475")/qd_real("262144"),
1542 qd_real("4450881435")/qd_real("131072"), qd_real("59826782925")/qd_real("2097152"),
1543 qd_real("171503444385")/qd_real("8388608"), qd_real("420696483235")/qd_real("33554432"),
1544 qd_real("221120793075")/qd_real("33554432"), qd_real("797168807855")/qd_real("268435456"),
1545qd_real("4923689695575")/qd_real("4294967296"), qd_real("6499270398159")/qd_real("17179869184"),
1546 qd_real("456864812569")/qd_real("4294967296"), qd_real("3486599885395")/qd_real("137438953472"),
1547qd_real("2804116503573")/qd_real("549755813888"), qd_real("1886827875075")/qd_real("2199023255552"),
1548 qd_real("263012370465")/qd_real("2199023255552"), qd_real("240141729555")/qd_real("17592186044416"),
1549 qd_real("176848560525")/qd_real("140737488355328"), qd_real("51538723353")/qd_real("562949953421312"),
1550 qd_real("1450433115")/qd_real("281474976710656"), qd_real("977699359")/qd_real("4503599627370496"),
1551 qd_real("118183439")/qd_real("18014398509481984"), qd_real("9652005")/qd_real("72057594037927936"),
1552 qd_real("121737")/qd_real("72057594037927936"), qd_real("6545")/qd_real("576460752303423488"),
1553 qd_real("561")/qd_real("18446744073709551616"), qd_real("1")/qd_real("73786976294838206464")
1554 };
1555#endif
1556}
1557
1558namespace glucat
1559{
1561 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1562 auto
1565 const index_t level) -> const matrix_multi<Scalar_T,LO,HI,Tune_P>
1566 {
1567 // Reference: [GW], Section 4.3, pp318-322
1568 // Reference: [GL], Section 11.3, p572-576
1569 // Reference: [Z], Pade1
1570
1571 using traits_t = numeric_traits<Scalar_T>;
1572
1573 if (val.isnan())
1574 return traits_t::NaN();
1575
1576 const auto scr_val = val.scalar();
1577 if (val == scr_val)
1578 {
1579 if (scr_val < Scalar_T(0))
1580 return i * traits_t::sqrt(-scr_val);
1581 else
1582 return traits_t::sqrt(scr_val);
1583 }
1584
1585 // Scale val towards abs(A) == 1 or towards A == 1 as appropriate
1586 const auto scale =
1587 (scr_val != Scalar_T(0) && norm(val/scr_val - Scalar_T(1)) < Scalar_T(1))
1588 ? scr_val
1589 : (scr_val < Scalar_T(0))
1590 ? -abs(val)
1591 : abs(val);
1592 const auto sqrt_scale = traits_t::sqrt(traits_t::abs(scale));
1593 if (traits_t::isNaN_or_isInf(sqrt_scale))
1594 return traits_t::NaN();
1595
1596 using multivector_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
1597 auto rescale = multivector_t(sqrt_scale);
1598 if (scale < Scalar_T(0))
1599 rescale = i * sqrt_scale;
1600
1601 const auto& unitval = val / scale;
1602 static const auto max_norm = Scalar_T(1.0/4.0);
1603 auto use_approx_sqrt = true;
1604 auto use_cr_sqrt = false;
1605 auto scaled_result = multivector_t();
1606#if defined(_GLUCAT_USE_EIGENVALUES)
1607 static const auto sqrt_2 = traits_t::sqrt(Scalar_T(2));
1608 if (level == 0)
1609 {
1610 // What kind of eigenvalues does the matrix contain?
1611 const auto genus = matrix::classify_eigenvalues(unitval.m_matrix);
1612 const index_t next_level =
1613 (genus.m_is_singular)
1614 ? level
1615 : level + 1;
1616 switch (genus.m_eig_case)
1617 {
1618 case matrix::neg_real_eigs:
1619 scaled_result = matrix_sqrt(-i * unitval, i, next_level) * (i + Scalar_T(1)) / sqrt_2;
1620 use_approx_sqrt = false;
1621 break;
1622 case matrix::both_eigs:
1623 {
1624 const auto safe_arg = genus.m_safe_arg;
1625 scaled_result = matrix_sqrt(exp(i*safe_arg) * unitval, i, next_level) * exp(-i*safe_arg / Scalar_T(2));
1626 }
1627 use_approx_sqrt = false;
1628 break;
1629 default:
1630 break;
1631 }
1632 use_cr_sqrt = genus.m_is_singular;
1633 }
1634#endif
1635 if (use_approx_sqrt)
1636 {
1637 scaled_result =
1638 (norm(unitval - Scalar_T(1)) < max_norm)
1639 // Pade' approximation of square root
1642 unitval - Scalar_T(1))
1643 // Product form of Denman-Beavers square root iteration
1644 : (use_cr_sqrt)
1645 ? cr_sqrt(unitval)
1646 : db_sqrt(unitval);
1647 }
1648 return (scaled_result.isnan() ||
1649 !approx_equal(pow(scaled_result, 2), unitval))
1650 ? traits_t::NaN()
1651 : scaled_result * rescale;
1652 }
1653
1655 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1656 auto
1658 {
1659 // Reference: [GW], Section 4.3, pp318-322
1660 // Reference: [GL], Section 11.3, p572-576
1661 // Reference: [Z], Pade1
1662
1663 using traits_t = numeric_traits<Scalar_T>;
1664
1665 if (val.isnan())
1666 return traits_t::NaN();
1667
1668 check_complex(val, i, prechecked);
1669
1670 switch (Tune_P::function_precision)
1671 {
1672 case precision_demoted:
1673 {
1674 using demoted_scalar_t = typename traits_t::demoted::type;
1675 using demoted_multivector_t = matrix_multi<demoted_scalar_t,LO,HI,Tune_P>;
1676
1677 const auto& demoted_val = demoted_multivector_t(val);
1678 const auto& demoted_i = demoted_multivector_t(i);
1679
1680 return matrix_sqrt(demoted_val, demoted_i, 0);
1681 }
1682 break;
1683 case precision_promoted:
1684 {
1685 using promoted_scalar_t = typename traits_t::promoted::type;
1686 using promoted_multivector_t = matrix_multi<promoted_scalar_t,LO,HI,Tune_P>;
1687
1688 const auto& promoted_val = promoted_multivector_t(val);
1689 const auto& promoted_i = promoted_multivector_t(i);
1690
1691 return matrix_sqrt(promoted_val, promoted_i, 0);
1692 }
1693 break;
1694 default:
1695 return matrix_sqrt(val, i, 0);
1696 }
1697 }
1698}
1699
1700namespace pade {
1702 // Reference: [Z], Pade1
1703 template< typename Scalar_T >
1705 {
1706 using array = std::array<Scalar_T, 14>;
1707 static const array numer;
1708 };
1709 template< typename Scalar_T >
1711 {
1712 0.0, 1.0, 6.0, 4741.0/300.0,
1713 1441.0/60.0, 107091.0/4600.0, 8638.0/575.0, 263111.0/40250.0,
1714 153081.0/80500.0, 395243.0/1101240.0, 28549.0/688275.0, 605453.0/228813200.0,
1715 785633.0/10296594000.0, 1145993.0/1873980108000.0
1716 };
1717
1719 // Reference: [Z], Pade1
1720 template< typename Scalar_T >
1722 {
1723 using array = std::array<Scalar_T, 14>;
1724 static const array denom;
1725 };
1726 template< typename Scalar_T >
1728 {
1729 1.0, 13.0/2.0, 468.0/25.0, 1573.0/50.0,
1730 1573.0/46.0, 11583.0/460.0, 10296.0/805.0, 2574.0/575.0,
1731 11583.0/10925.0, 143.0/874.0, 572.0/37145.0, 117.0/148580.0,
1732 13.0/742900.0, 1.0/10400600.0
1733 };
1734
1735 template< >
1736 struct pade_log_numer<float>
1737 {
1738 using array = std::array<float, 10>;
1739 static const array numer;
1740 };
1742 {
1743 0.0, 1.0, 4.0, 1337.0/204.0,
1744 385.0/68.0, 1879.0/680.0, 193.0/255.0, 197.0/1820.0,
1745 419.0/61880.0, 7129.0/61261200.0
1746 };
1747 template< >
1748 struct pade_log_denom<float>
1749 {
1750 using array = std::array<float, 10>;
1751 static const array denom;
1752 };
1754 {
1755 1.0, 9.0/2.0, 144.0/17.0, 147.0/17.0,
1756 441.0/85.0, 63.0/34.0, 84.0/221.0, 9.0/221.0,
1757 9.0/4862.0, 1.0/48620.0
1758 };
1759
1760 template< >
1761 struct pade_log_numer<long double>
1762 {
1763 using array = std::array<long double, 18>;
1764 static const array numer;
1765 };
1767 {
1768 0.0L, 1.0L, 8.0L, 3835.0L/132.0L,
1769 8365.0L/132.0L, 11363807.0L/122760.0L, 162981.0L/1705.0L, 9036157.0L/125860.0L,
1770 18009875.0L/453096.0L, 44211925.0L/2718576.0L, 4149566.0L/849555.0L, 16973929.0L/16020180.0L,
1771 172459.0L/1068012.0L, 116317061.0L/7025382936.0L, 19679783.0L/18441630207.0L, 23763863.0L/614721006900.0L,
1772 50747.0L/79318839600.0L, 42142223.0L/14295951736466400.0L
1773 };
1774 template< >
1775 struct pade_log_denom<long double>
1776 {
1777 using array = std::array<long double, 18>;
1778 static const array denom;
1779 };
1781 {
1782 1.0L, 17.0L/2.0L, 1088.0L/33.0L, 850.0L/11.0L,
1783 41650.0L/341.0L, 140777.0L/1023.0L, 1126216.0L/9889.0L, 63206.0L/899.0L,
1784 790075.0L/24273.0L, 60775.0L/5394.0L, 38896.0L/13485.0L, 21658.0L/40455.0L,
1785 21658.0L/310155.0L, 4165.0L/682341.0L, 680.0L/2047023.0L, 34.0L/3411705.0L,
1786 17.0L/129644790.0L, 1.0L/2333606220
1787 };
1788#if defined(_GLUCAT_USE_QD)
1789 template< >
1790 struct pade_log_numer<dd_real>
1791 {
1792 using array = std::array<dd_real, 22>;
1793 static const array numer;
1794 };
1796 {
1797 dd_real("0"), dd_real("1"),
1798 dd_real("10"), dd_real("22781")/dd_real("492"),
1799 dd_real("21603")/dd_real("164"), dd_real("5492649")/dd_real("21320"),
1800 dd_real("978724")/dd_real("2665"), dd_real("4191605")/dd_real("10619"),
1801 dd_real("12874933")/dd_real("39442"), dd_real("11473457")/dd_real("54612"),
1802 dd_real("2406734")/dd_real("22755"), dd_real("166770367")/dd_real("4004880"),
1803 dd_real("30653165")/dd_real("2402928"), dd_real("647746389")/dd_real("215195552"),
1804 dd_real("25346331")/dd_real("47074027"), dd_real("278270613")/dd_real("3900419380"),
1805 dd_real("105689791")/dd_real("15601677520"), dd_real("606046475")/dd_real("1379188292768"),
1806 dd_real("969715")/dd_real("53502994116"), dd_real("11098301")/dd_real("26204577562592"),
1807 dd_real("118999")/dd_real("26204577562592"), dd_real("18858053")/dd_real("1392249205900512960")
1808 };
1809 template< >
1810 struct pade_log_denom<dd_real>
1811 {
1812 using array = std::array<dd_real, 22>;
1813 static const array denom;
1814 };
1816 {
1817 dd_real("1"), dd_real("21")/dd_real("2"),
1818 dd_real("2100")/dd_real("41"), dd_real("12635")/dd_real("82"),
1819 dd_real("341145")/dd_real("1066"), dd_real("1037799")/dd_real("2132"),
1820 dd_real("11069856")/dd_real("19721"), dd_real("9883800")/dd_real("19721"),
1821 dd_real("6918660")/dd_real("19721"), dd_real("293930")/dd_real("1517"),
1822 dd_real("1410864")/dd_real("16687"), dd_real("88179")/dd_real("3034"),
1823 dd_real("734825")/dd_real("94054"), dd_real("305235")/dd_real("188108"),
1824 dd_real("348840")/dd_real("1363783"), dd_real("40698")/dd_real("1363783"),
1825 dd_real("6783")/dd_real("2727566"), dd_real("9975")/dd_real("70916716"),
1826 dd_real("266")/dd_real("53187537"), dd_real("7")/dd_real("70916716"),
1827 dd_real("7")/dd_real("8155422340"), dd_real("1")/dd_real("538257874440")
1828 };
1829
1830 template< >
1831 struct pade_log_numer<qd_real>
1832 {
1833 using array = std::array<qd_real, 34>;
1834 static const array numer;
1835 };
1837 {
1838 qd_real("0"), qd_real("1"),
1839 qd_real("16"), qd_real("95201")/qd_real("780"),
1840 qd_real("30721")/qd_real("52"), qd_real("7416257")/qd_real("3640"),
1841 qd_real("1039099")/qd_real("195"), qd_real("6097772319")/qd_real("555100"),
1842 qd_real("1564058073")/qd_real("85400"), qd_real("30404640205")/qd_real("1209264"),
1843 qd_real("725351278")/qd_real("25193"), qd_real("4092322670789")/qd_real("147429436"),
1844 qd_real("4559713849589")/qd_real("201040140"), qd_real("5049361751189")/qd_real("320023080"),
1845 qd_real("74979677195")/qd_real("8000577"), qd_real("16569850691873")/qd_real("3481514244"),
1846 qd_real("1065906022369")/qd_real("515779888"), qd_real("335956770855841")/qd_real("438412904800"),
1847 qd_real("1462444287585964")/qd_real("6041877844275"), qd_real("397242326339851")/qd_real("6122436215532"),
1848 qd_real("64211291334131")/qd_real("4373168725380"), qd_real("142322343550859")/qd_real("51080680851480"),
1849 qd_real("154355972958659")/qd_real("351179680853925"), qd_real("167483568676259")/qd_real("2937139148960100"),
1850 qd_real("4230788929433")/qd_real("704913395750424"), qd_real("197968763176019")/qd_real("392923948371995600"),
1851 qd_real("10537522306718")/qd_real("319250708052246425"), qd_real("236648286272519")/qd_real("144249197475035425500"),
1852 qd_real("260715545088119")/qd_real("4375558990076074573500"), qd_real("289596255666839")/qd_real("192874640282553367199880"),
1853 qd_real("8802625510547")/qd_real("361639950529787563499775"), qd_real("373831661521439")/qd_real("1659204093030665341336967700"),
1854 qd_real("446033437968239")/qd_real("464577146048586295574350956000"), qd_real("53676090078349")/qd_real("47386868896955802148583797512000")
1855 };
1856 template< >
1857 struct pade_log_denom<qd_real>
1858 {
1859 using array = std::array<qd_real, 34>;
1860 static const array denom;
1861 };
1863 {
1864 qd_real("1"), qd_real("33")/qd_real("2"),
1865 qd_real("8448")/qd_real("65"), qd_real("42284")/qd_real("65"),
1866 qd_real("211420")/qd_real("91"), qd_real("573562")/qd_real("91"),
1867 qd_real("32119472")/qd_real("2379"), qd_real("92917044")/qd_real("3965"),
1868 qd_real("603960786")/qd_real("17995"), qd_real("144626625")/qd_real("3599"),
1869 qd_real("2776831200")/qd_real("68381"), qd_real("16692542100")/qd_real("478667"),
1870 qd_real("12241197540")/qd_real("478667"), qd_real("1098569010")/qd_real("68381"),
1871 qd_real("31387686000")/qd_real("3624193"), qd_real("9939433900")/qd_real("2479711"),
1872 qd_real("67091178825")/qd_real("42155087"), qd_real("2683647153")/qd_real("4959422"),
1873 qd_real("19083713088")/qd_real("121505839"), qd_real("4708152900")/qd_real("121505839"),
1874 qd_real("941630580")/qd_real("116546417"), qd_real("88704330")/qd_real("62755763"),
1875 qd_real("12902448")/qd_real("62755763"), qd_real("1542684")/qd_real("62755763"),
1876 qd_real("6427850")/qd_real("2698497809"), qd_real("3471039")/qd_real("18889484663"),
1877 qd_real("8544096")/qd_real("774468871183"), qd_real("39556")/qd_real("79027435835"),
1878 qd_real("118668")/qd_real("7191496660985"), qd_real("10230")/qd_real("27327687311743"),
1879 qd_real("5456")/qd_real("1011124430534491"), qd_real("44")/qd_real("1011124430534491"),
1880 qd_real("11")/qd_real("70778710137414370"), qd_real("1")/qd_real("7219428434016265740")
1881 };
1882#endif
1883}
1884
1885namespace glucat{
1887 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1888 static
1889 auto
1891 {
1892 // Reference: [GW], Section 4.3, pp318-322
1893 // Reference: [CHKL]
1894 // Reference: [GL], Section 11.3, p572-576
1895 // Reference: [Z], Pade1
1896
1897 using traits_t = numeric_traits<Scalar_T>;
1898 if (val == Scalar_T(0) || val.isnan())
1899 return traits_t::NaN();
1900 else
1903 val - Scalar_T(1));
1904 }
1905
1907 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1908 static
1909 auto
1911 {
1912 // Reference: [CHKL]
1913 using multivector_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
1914 using traits_t = numeric_traits<Scalar_T>;
1915 if (val == Scalar_T(0) || val.isnan())
1916 return traits_t::NaN();
1917
1918 using limits_t = std::numeric_limits<Scalar_T>;
1919 static const auto epsilon = limits_t::epsilon();
1920 static const auto max_inner_norm = traits_t::pow(epsilon, 2);
1921 static const auto max_outer_norm = Scalar_T(6.0/limits_t::digits);
1922 auto Y = val;
1923 auto E = multivector_t(Scalar_T(0));
1924 Scalar_T norm_Y_1;
1925 auto pow_2_outer_step = Scalar_T(1);
1926 auto pow_4_outer_step = Scalar_T(1);
1927 int outer_step;
1928 for (outer_step = 0, norm_Y_1 = norm(Y - Scalar_T(1));
1929 outer_step != Tune_P::log_max_outer_steps && norm_Y_1 * pow_2_outer_step > max_outer_norm;
1930 ++outer_step, norm_Y_1 = norm(Y - Scalar_T(1)))
1931 {
1932 if (Y == Scalar_T(0) || Y.isnan())
1933 return traits_t::NaN();
1934
1935 // Incomplete product form of Denman-Beavers square root iteration
1936 auto M = Y;
1937 for (auto
1938 inner_step = 0;
1939 inner_step != Tune_P::log_max_inner_steps &&
1940 norm(M - Scalar_T(1)) * pow_4_outer_step > max_inner_norm;
1941 ++inner_step)
1942 db_step(M, Y);
1943
1944 E += (M - Scalar_T(1)) * pow_2_outer_step;
1945 pow_2_outer_step *= Scalar_T(2);
1946 pow_4_outer_step *= Scalar_T(4);
1947 }
1948 if (outer_step == Tune_P::log_max_outer_steps && norm_Y_1 * pow_2_outer_step > max_outer_norm)
1949 return traits_t::NaN();
1950 else
1951 return pade_log(Y) * pow_2_outer_step - E;
1952 }
1953
1955 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1956 auto
1959 const index_t level) -> const matrix_multi<Scalar_T,LO,HI,Tune_P>
1960 {
1961 // Scaled incomplete square root cascade and scaled Pade' approximation of log
1962 // Reference: [CHKL]
1963
1964 using traits_t = numeric_traits<Scalar_T>;
1965 if (val == Scalar_T(0) || val.isnan())
1966 return traits_t::NaN();
1967
1968 static const auto pi = traits_t::pi();
1969 const auto scr_val = val.scalar();
1970 if (val == scr_val)
1971 {
1972 if (scr_val < Scalar_T(0))
1973 return i * pi + traits_t::log(-scr_val);
1974 else
1975 return traits_t::log(scr_val);
1976 }
1977
1978 // Scale val towards abs(A) == 1 or towards A == 1 as appropriate
1979 const auto max_norm = Scalar_T(1.0/9.0);
1980 const auto scale =
1981 (scr_val != Scalar_T(0) && norm(val/scr_val - Scalar_T(1)) < max_norm)
1982 ? scr_val
1983 : (scr_val < Scalar_T(0))
1984 ? -abs(val)
1985 : abs(val);
1986 if (scale == Scalar_T(0))
1987 return traits_t::NaN();
1988
1989 using multivector_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
1990 const auto log_scale = traits_t::log(traits_t::abs(scale));
1991 auto rescale = multivector_t(log_scale);
1992 if (scale < Scalar_T(0))
1993 rescale = i * pi + log_scale;
1994 const auto unitval = val/scale;
1995 if (inv(unitval).isnan())
1996 return traits_t::NaN();
1997
1998#if defined(_GLUCAT_USE_EIGENVALUES)
1999 auto scaled_result = multivector_t();
2000 if (level == 0)
2001 {
2002 // What kind of eigenvalues does the matrix contain?
2003 auto genus = matrix::classify_eigenvalues(unitval.m_matrix);
2004 switch (genus.m_eig_case)
2005 {
2006 case matrix::neg_real_eigs:
2007 scaled_result = matrix_log(-i * unitval, i, level + 1) + i * pi/Scalar_T(2);
2008 break;
2009 case matrix::both_eigs:
2010 {
2011 const Scalar_T safe_arg = genus.m_safe_arg;
2012 scaled_result = matrix_log(exp(i*safe_arg) * unitval, i, level + 1) - i * safe_arg;
2013 }
2014 break;
2015 default:
2016 scaled_result = cascade_log(unitval);
2017 break;
2018 }
2019 }
2020 else
2021 scaled_result = cascade_log(unitval);
2022#else
2023 auto scaled_result = cascade_log(unitval);
2024#endif
2025 return (scaled_result.isnan())
2026 ? traits_t::NaN()
2027 : scaled_result + rescale;
2028 }
2029
2031 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
2032 auto
2034 {
2035 using traits_t = numeric_traits<Scalar_T>;
2036
2037 if (val == Scalar_T(0) || val.isnan())
2038 return traits_t::NaN();
2039
2040 check_complex(val, i, prechecked);
2041
2042 switch (Tune_P::function_precision)
2043 {
2044 case precision_demoted:
2045 {
2046 using demoted_scalar_t = typename traits_t::demoted::type;
2047 using demoted_multivector_t = matrix_multi<demoted_scalar_t,LO,HI,Tune_P>;
2048
2049 const auto& demoted_val = demoted_multivector_t(val);
2050 const auto& demoted_i = demoted_multivector_t(i);
2051
2052 return matrix_log(demoted_val, demoted_i, 0);
2053 }
2054 break;
2055 case precision_promoted:
2056 {
2057 using promoted_scalar_t = typename traits_t::promoted::type;
2058 using promoted_multivector_t = matrix_multi<promoted_scalar_t,LO,HI,Tune_P>;
2059
2060 const auto& promoted_val = promoted_multivector_t(val);
2061 const auto& promoted_i = promoted_multivector_t(i);
2062
2063 return matrix_log(promoted_val, promoted_i, 0);
2064 }
2065 break;
2066 default:
2067 return matrix_log(val, i, 0);
2068 }
2069 }
2070
2072 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
2073 auto
2075 {
2076 using traits_t = numeric_traits<Scalar_T>;
2077 if (val.isnan())
2078 return traits_t::NaN();
2079
2080 const auto scr_val = val.scalar();
2081 if (val == scr_val)
2082 return traits_t::exp(scr_val);
2083
2084 switch (Tune_P::function_precision)
2085 {
2086 case precision_demoted:
2087 {
2088 using demoted_scalar_t = typename traits_t::demoted::type;
2089 using demoted_multivector_t = matrix_multi<demoted_scalar_t,LO,HI,Tune_P>;
2090
2091 const auto& demoted_val = demoted_multivector_t(val);
2092 return clifford_exp(demoted_val);
2093 }
2094 break;
2095 case precision_promoted:
2096 {
2097 using promoted_scalar_t = typename traits_t::promoted::type;
2098 using promoted_multivector_t = matrix_multi<promoted_scalar_t,LO,HI,Tune_P>;
2099
2100 const auto& promoted_val = promoted_multivector_t(val);
2101 return clifford_exp(promoted_val);
2102 }
2103 break;
2104 default:
2105 return clifford_exp(val);
2106 }
2107 }
2108}
2109#endif // _GLUCAT_MATRIX_MULTI_IMP_H
const scalar_t epsilon
Definition PyClical.h:150
Table of basis elements used as a cache by basis_element()
auto operator=(const basis_table &) -> basis_table &=delete
static auto basis() -> basis_table &
Single instance of basis table.
~basis_table()=default
friend class friend_for_private_destructor
basis_table(const basis_table &)=delete
virtual auto truncated(const Scalar_T &limit=default_truncation) const -> const multivector_t=0
Remove all terms with relative size smaller than limit.
A framed_multi<Scalar_T,LO,HI,Tune_P> is a framed approximation to a multivector.
auto centre_pm4_qp4(index_t &p, index_t &q) -> multivector_t &
Subalgebra isomorphism: R_{p,q} to R_{p-4,q+4}.
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 unfold(const index_set_t frm) const -> multivector_t
Subalgebra isomorphism: unfold each term within the given frame.
auto centre_pp4_qm4(index_t &p, index_t &q) -> multivector_t &
Subalgebra isomorphism: R_{p,q} to R_{p+4,q-4}.
static auto generator() -> generator_table< Matrix_T > &
Single instance of generator table.
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 max() const -> index_t
Maximum member.
A matrix_multi<Scalar_T,LO,HI,Tune_P> is a matrix approximation to a multivector.
_GLUCAT_CLIFFORD_ALGEBRA_OPERATIONS auto operator=(const multivector_t &rhs) -> multivector_t &
Assignment operator.
typename matrix_t::size_type matrix_index_t
friend class framed_multi
static auto classname() -> const std::string
Class name used in messages.
index_set< LO, HI > index_set_t
auto basis_element(const index_set< LO, HI > &ist) const -> const basis_matrix_t
Create a basis element matrix within the current frame.
std::pair< const index_set_t, Scalar_T > term_t
auto operator+=(const term_t &rhs) -> multivector_t &
Add a term, if non-zero.
matrix_t m_matrix
Matrix value representing the multivector within the folded frame.
auto fast_framed_multi() const -> const framed_multi< Other_Scalar_T, LO, HI, Tune_P >
Use inverse generalized FFT to construct a framed_multi_t.
auto fast_matrix_multi(const index_set_t frm) const -> const matrix_multi_t
Use generalized FFT to construct a matrix_multi_t.
ublas::compressed_matrix< int, orientation_t > basis_matrix_t
index_set_t m_frame
Index set representing the frame for the subalgebra which contains the multivector.
friend class matrix_multi
ublas::matrix< Scalar_T, orientation_t > matrix_t
matrix_multi multivector_t
std::vector< Scalar_T > vector_t
static auto random(const index_set_t frm, Scalar_T fill=Scalar_T(1)) -> const matrix_multi_t
Random multivector within a frame.
framed_multi< Scalar_T, LO, HI, Tune_P > framed_multi_t
Extra traits which extend numeric limits.
Definition scalar.h:48
static auto NaN() -> Scalar_T
Smart NaN.
Definition scalar.h:115
static auto to_scalar_t(const Other_Scalar_T &val) -> Scalar_T
Cast to Scalar_T.
Definition scalar.h:141
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 signed_perm_nork(const LHS_T &lhs, const RHS_T &rhs) -> const RHS_T
Left inverse of Kronecker product where lhs is a signed permutation matrix.
Definition matrix_imp.h:223
auto trace(const Matrix_T &val) -> typename Matrix_T::value_type
Matrix trace.
Definition matrix_imp.h:411
auto isinf(const Matrix_T &m) -> bool
Infinite.
Definition matrix_imp.h:270
auto norm_frob2(const Matrix_T &val) -> typename Matrix_T::value_type
Square of Frobenius norm.
Definition matrix_imp.h:390
auto classify_eigenvalues(const Matrix_T &val) -> eig_genus< Matrix_T >
Classify the eigenvalues of a matrix.
Definition matrix_imp.h:492
auto mono_prod(const ublas::matrix_expression< LHS_T > &lhs, const ublas::matrix_expression< RHS_T > &rhs) -> const typename RHS_T::expression_type
Product of monomial matrices.
Definition matrix_imp.h:315
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.
auto approx_equal(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const RHS< Scalar_T, LO, HI, Tune_P > &rhs, const Scalar_T threshold, const Scalar_T tolerance) -> bool
Test for approximate equality of multivectors.
static auto pade_approx(const std::array< Scalar_T, Size > &numer, const std::array< Scalar_T, Size > &denom, const matrix_multi< Scalar_T, LO, HI, Tune_P > &X) -> const matrix_multi< Scalar_T, LO, HI, Tune_P >
Pade' approximation.
static auto fast(const Matrix_T &X, index_t level) -> Multivector_T
Inverse generalized Fast Fourier Transform.
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 cr_sqrt(const matrix_multi< Scalar_T, LO, HI, Tune_P > &val, Scalar_T norm_Y_tol=std::pow(std::numeric_limits< Scalar_T >::epsilon(), 1)) -> const matrix_multi< Scalar_T, LO, HI, Tune_P >
Cyclic reduction square root iteration.
static auto db_sqrt(const matrix_multi< Scalar_T, LO, HI, Tune_P > &val, Scalar_T norm_tol=std::pow(std::numeric_limits< Scalar_T >::epsilon(), 4)) -> const matrix_multi< Scalar_T, LO, HI, Tune_P >
Product form of Denman-Beavers square root iteration.
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.
static void db_step(matrix_multi< Scalar_T, LO, HI, Tune_P > &M, matrix_multi< Scalar_T, LO, HI, Tune_P > &Y)
Single step of product form of Denman-Beavers square root iteration.
auto reframe(const matrix_multi< Scalar_T, LO, HI, Tune_P > &lhs, const matrix_multi< Scalar_T, LO, HI, Tune_P > &rhs, matrix_multi< Scalar_T, LO, HI, Tune_P > &lhs_reframed, matrix_multi< Scalar_T, LO, HI, Tune_P > &rhs_reframed) -> const index_set< LO, HI >
Find a common frame for operands of a binary operator.
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.
auto abs(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> Scalar_T
Absolute value == sqrt(norm)
auto inv(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Geometric multiplicative inverse.
static auto pade_log(const matrix_multi< Scalar_T, LO, HI, Tune_P > &val) -> const matrix_multi< Scalar_T, LO, HI, Tune_P >
Pade' approximation of log.
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.
static auto cascade_log(const matrix_multi< Scalar_T, LO, HI, Tune_P > &val) -> const matrix_multi< Scalar_T, LO, HI, Tune_P >
Incomplete square root cascade and Pade' approximation of log.
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 norm(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> Scalar_T
Scalar_T norm == sum of norm of coordinates.
auto matrix_sqrt(const matrix_multi< Scalar_T, LO, HI, Tune_P > &val, const matrix_multi< Scalar_T, LO, HI, Tune_P > &i, const index_t level) -> const matrix_multi< Scalar_T, LO, HI, Tune_P >
Square root 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 matrix_log(const matrix_multi< Scalar_T, LO, HI, Tune_P > &val, const matrix_multi< Scalar_T, LO, HI, Tune_P > &i, const index_t level) -> const matrix_multi< Scalar_T, LO, HI, Tune_P >
Natural logarithm of multivector with specified complexifier.
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.
static auto folded_dim(const index_set< LO, HI > &sub) -> Matrix_Index_T
Determine the matrix dimension of the fold of a subalegbra.
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()
auto offset_level(const index_t p, const index_t q) -> index_t
Determine the log2 dim corresponding to signature p, q.
std::array< dd_real, 22 > array
std::array< float, 10 > array
std::array< long double, 18 > array
std::array< qd_real, 34 > array
Coefficients of denominator polynomials of Pade approximations produced by Pade1(log(1+x),...
static const array denom
std::array< Scalar_T, 14 > array
std::array< dd_real, 22 > array
std::array< float, 10 > array
std::array< long double, 18 > array
std::array< qd_real, 34 > array
Coefficients of numerator polynomials of Pade approximations produced by Pade1(log(1+x),...
static const array numer
std::array< Scalar_T, 14 > array
std::array< dd_real, 22 > array
std::array< float, 10 > array
std::array< long double, 18 > array
std::array< qd_real, 34 > array
Coefficients of denominator polynomials of Pade approximations produced by Pade1(sqrt(1+x),...
static const array denom
std::array< Scalar_T, 14 > array
std::array< dd_real, 22 > array
std::array< float, 10 > array
std::array< long double, 18 > array
std::array< qd_real, 34 > array
Coefficients of numerator polynomials of Pade approximations produced by Pade1(sqrt(1+x),...
std::array< Scalar_T, 14 > array
static const array numer