glucat 0.12.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 >
108 : m_frame( index_set_t() ),
109 m_matrix( matrix_t( 1, 1 ) )
110 { this->m_matrix.clear(); }
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 {
142 const matrix_index_t dim = folded_dim<matrix_index_t>(frm);
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;
167 }
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);
179 }
180
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 )
186 {
187 if (!prechecked && (ist | frm) != frm)
188 throw error_t("multivector_t(ist,crd,frm): cannot initialize with value outside of frame");
189 const matrix_index_t dim = folded_dim<matrix_index_t>(frm);
190 this->m_matrix.resize(dim, dim, false);
191 this->m_matrix.clear();
192 *this += term_t(ist, crd);
193 }
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
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
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)
228 {
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
270 for (auto& val_term : val)
271 *this += val_term;
272 }
273
275 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
276 template< typename Other_Scalar_T >
278 matrix_multi(const framed_multi<Other_Scalar_T,LO,HI,Tune_P>& framed_val, const index_set_t frm, const bool prechecked)
279 {
280 const auto val = framed_val.truncated();
281 const auto our_frame = val.frame() | frm;
282 if (val.size() >= Tune_P::fast_size_threshold)
283 try
284 {
285 *this = val.template fast_matrix_multi<Scalar_T>(our_frame);
286 return;
287 }
288 catch (const glucat_error& e)
289 { }
290 this->m_frame = our_frame;
291 const auto dim = folded_dim<matrix_index_t>(our_frame);
292 this->m_matrix.resize(dim, dim, false);
293 this->m_matrix.clear();
294
296 for (auto& val_term : val)
297 *this += val_term;
298 }
299
301 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
302 template< typename Matrix_T >
304 matrix_multi(const Matrix_T& mtx, const index_set_t frm)
305 : m_frame( frm ), m_matrix( mtx.size1(), mtx.size2() )
306 {
307 this->m_matrix.clear();
308
309 for (auto
310 mtx_it1 = mtx.begin1();
311 mtx_it1 != mtx.end1();
312 ++mtx_it1)
313 for (auto
314 mtx_it2 = mtx_it1.begin();
315 mtx_it2 != mtx_it1.end();
316 ++mtx_it2)
317 this->m_matrix(mtx_it2.index1(), mtx_it2.index2()) = numeric_traits<Scalar_T>::to_scalar_t(*mtx_it2);
318 }
319
321 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
323 matrix_multi(const matrix_t& mtx, const index_set_t frm)
324 : m_frame( frm ), m_matrix( mtx )
325 { }
326
328 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
329 auto
332 {
333 // Check for assignment to self
334 if (this == &rhs)
335 return *this;
336 this->m_frame = rhs.m_frame;
337 this->m_matrix = rhs.m_matrix;
338 return *this;
339 }
340
342 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
343 inline
344 auto
347 {
348 using index_set_t = index_set<LO, HI>;
349 using multivector_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
350 using framed_multi_t = typename multivector_t::framed_multi_t;
351 // Determine the initial common frame
352 index_set_t our_frame = lhs.m_frame | rhs.m_frame;
353 framed_multi_t framed_lhs;
354 framed_multi_t framed_rhs;
355 if ((lhs.m_frame != our_frame) || (rhs.m_frame != our_frame))
356 {
357 // The common frame may expand as a result of the transform to framed_multi_t
358 framed_lhs = framed_multi_t(lhs);
359 framed_rhs = framed_multi_t(rhs);
360 our_frame |= framed_lhs.frame() | framed_rhs.frame();
361 }
362 // Do the reframing only where necessary
363 if (lhs.m_frame != our_frame)
364 lhs_reframed = multivector_t(framed_lhs, our_frame, true);
365 if (rhs.m_frame != our_frame)
366 rhs_reframed = multivector_t(framed_rhs, our_frame, true);
367 return our_frame;
368 }
369
371 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
372 auto
374 operator== (const multivector_t& rhs) const -> bool
375 {
376 // Ensure that there is no aliasing
377 if (this == &rhs)
378 return true;
379
380 // Operate only within a common frame
381 multivector_t lhs_reframed;
382 multivector_t rhs_reframed;
383 const index_set_t our_frame = reframe(*this, rhs, lhs_reframed, rhs_reframed);
384 const multivector_t& lhs_ref = (this->m_frame == our_frame)
385 ? *this
386 : lhs_reframed;
387 const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
388 ? rhs
389 : rhs_reframed;
390
391 return ublas::norm_inf(lhs_ref.m_matrix - rhs_ref.m_matrix) == 0;
392 }
393
394 // Test for equality of multivector and scalar
395 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
396 inline
397 auto
399 operator== (const Scalar_T& scr) const -> bool
400 {
401 if (scr != Scalar_T(0))
402 return *this == multivector_t(framed_multi_t(scr), this->m_frame, true);
403 else if (ublas::norm_inf(this->m_matrix) != 0)
404 return false;
405 else
406 {
407 const matrix_index_t dim = this->m_matrix.size1();
408 return !(dim == 1 && this->isnan());
409 }
410 }
411
413 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
414 inline
415 auto
417 operator+= (const Scalar_T& scr) -> multivector_t&
418 { return *this += term_t(index_set_t(), scr); }
419
421 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
422 inline
423 auto
425 operator+= (const multivector_t& rhs) -> multivector_t&
426 {
427 // Ensure that there is no aliasing
428 if (this == &rhs)
429 return *this *= Scalar_T(2);
430
431 // Operate only within a common frame
432 multivector_t rhs_reframed;
433 const index_set_t our_frame = reframe(*this, rhs, *this, rhs_reframed);
434 const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
435 ? rhs
436 : rhs_reframed;
437
438 noalias(this->m_matrix) += rhs_ref.m_matrix;
439 return *this;
440 }
441
443 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
444 inline
445 auto
447 operator-= (const Scalar_T& scr) -> multivector_t&
448 { return *this += term_t(index_set_t(), -scr); }
449
451 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
452 inline
453 auto
455 operator-= (const multivector_t& rhs) -> multivector_t&
456 {
457 // Ensure that there is no aliasing
458 if (this == &rhs)
459 return *this = Scalar_T(0);
460
461 // Operate only within a common frame
462 multivector_t rhs_reframed;
463 const index_set_t our_frame = reframe(*this, rhs, *this, rhs_reframed);
464 const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
465 ? rhs
466 : rhs_reframed;
467
468 noalias(this->m_matrix) -= rhs_ref.m_matrix;
469 return *this;
470 }
471
473 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
474 inline
475 auto
477 operator- () const -> const multivector_t
478 { return multivector_t(-(this->m_matrix), this->m_frame); }
479
481 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
482 inline
483 auto
485 operator*= (const Scalar_T& scr) -> multivector_t&
486 { // multiply coordinates of all terms by scalar
487
488 using traits_t = numeric_traits<Scalar_T>;
489 if (traits_t::isNaN_or_isInf(scr) || this->isnan())
490 return *this = traits_t::NaN();
491 if (scr == Scalar_T(0))
492 *this = Scalar_T(0);
493 else
494 this->m_matrix *= scr;
495 return *this;
496 }
497
499 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
500 inline
501 auto
503 {
504 using multivector_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
505 using index_set_t = typename multivector_t::index_set_t;
506
507 if (lhs.isnan() || rhs.isnan())
509
510 // Operate only within a common frame
511 multivector_t lhs_reframed;
512 multivector_t rhs_reframed;
513 const index_set_t our_frame = reframe(lhs, rhs, lhs_reframed, rhs_reframed);
514 const multivector_t& lhs_ref = (lhs.m_frame == our_frame)
515 ? lhs
516 : lhs_reframed;
517 const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
518 ? rhs
519 : rhs_reframed;
520
521 using matrix_t = typename multivector_t::matrix_t;
522 using matrix_index_t = typename matrix_t::size_type;
523
524 const matrix_index_t dim = lhs_ref.m_matrix.size1();
525 multivector_t result = multivector_t(matrix_t(dim, dim), our_frame);
526 result.m_matrix.clear();
527 ublas::axpy_prod(lhs_ref.m_matrix, rhs_ref.m_matrix, result.m_matrix, true);
528 return result;
529 }
530
532 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
533 inline
534 auto
536 operator*= (const multivector_t& rhs) -> multivector_t&
537 { return *this = *this * rhs; }
538
540 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
541 inline
542 auto
544 {
545 using multivector_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
546 using framed_multi_t = typename multivector_t::framed_multi_t;
547 return framed_multi_t(lhs) ^ framed_multi_t(rhs);
548 }
549
551 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
552 inline
553 auto
555 operator^= (const multivector_t& rhs) -> multivector_t&
556 { return *this = *this ^ rhs; }
557
559 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
560 inline
561 auto
563 {
564 using multivector_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
565 using framed_multi_t = typename multivector_t::framed_multi_t;
566 return framed_multi_t(lhs) & framed_multi_t(rhs);
567 }
568
570 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
571 inline
572 auto
574 operator&= (const multivector_t& rhs) -> multivector_t&
575 { return *this = *this & rhs; }
576
578 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
579 inline
580 auto
582 {
583 using multivector_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
584 using framed_multi_t = typename multivector_t::framed_multi_t;
585 return framed_multi_t(lhs) % framed_multi_t(rhs);
586 }
587
589 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
590 inline
591 auto
593 operator%= (const multivector_t& rhs) -> multivector_t&
594 { return *this = *this % rhs; }
595
597 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
598 inline
599 auto
601 { return (lhs * rhs).scalar(); }
602
604 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
605 inline
606 auto
608 operator/= (const Scalar_T& scr) -> multivector_t&
609 { return *this *= Scalar_T(1)/scr; }
610
612 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
613 auto
615 {
616 using traits_t = numeric_traits<Scalar_T>;
617
618 if (lhs.isnan() || rhs.isnan())
619 return traits_t::NaN();
620
621 if (rhs == Scalar_T(0))
622 return traits_t::NaN();
623
624 using multivector_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
625 using index_set_t = typename multivector_t::index_set_t;
626
627 // Operate only within a common frame
628 multivector_t lhs_reframed;
629 multivector_t rhs_reframed;
630 const auto our_frame = reframe(lhs, rhs, lhs_reframed, rhs_reframed);
631 const auto& lhs_ref = (lhs.m_frame == our_frame)
632 ? lhs
633 : lhs_reframed;
634 const auto& rhs_ref = (rhs.m_frame == our_frame)
635 ? rhs
636 : rhs_reframed;
637
638 // Solve result == lhs_ref/rhs_ref <=> result*rhs_ref == lhs_ref
639 // We now solve X == B/A
640 // (where X == result, B == lhs_ref.m_matrix and A == rhs_ref.m_matrix)
641 // X == B/A <=> X*A == B <=> AT*XT == BT
642 // So, we solve AT*XT == BT
643
644 using matrix_t = typename multivector_t::matrix_t;
645 using matrix_index_t = typename matrix_t::size_type;
646
647 const auto& AT = matrix_t(ublas::trans(rhs_ref.m_matrix));
648 auto LU = AT;
649
650 using permutation_t = ublas::permutation_matrix<matrix_index_t>;
651
652 auto pvector = permutation_t(AT.size1());
653 if (! ublas::lu_factorize(LU, pvector))
654 {
655 const auto& BT = matrix_t(ublas::trans(lhs_ref.m_matrix));
656 auto XT = BT;
657 ublas::lu_substitute(LU, pvector, XT);
658 if (matrix::isnan(XT))
659 return traits_t::NaN();
660
661 // Iterative refinement.
662 // Reference: Nicholas J. Higham, "Accuracy and Stability of Numerical Algorithms",
663 // SIAM, 1996, ISBN 0-89871-355-2, Chapter 11
664 if (Tune_P::div_max_steps > 0)
665 {
666 // matrix_t R = ublas::prod(AT, XT) - BT;
667 auto R = matrix_t(-BT);
668 ublas::axpy_prod(AT, XT, R, false);
669 if (matrix::isnan(R))
670 return traits_t::NaN();
671
672 auto nr = Scalar_T(ublas::norm_inf(R));
673 if ( nr != Scalar_T(0) && !traits_t::isNaN_or_isInf(nr) )
674 {
675 auto XTnew = XT;
676 auto nrold = nr + Scalar_T(1);
677 for (auto
678 step = 0;
679 step != Tune_P::div_max_steps &&
680 nr < nrold &&
681 nr != Scalar_T(0) &&
682 nr == nr;
683 ++step)
684 {
685 nrold = nr;
686 if (step != 0)
687 XT = XTnew;
688 auto& D = R;
689 ublas::lu_substitute(LU, pvector, D);
690 XTnew -= D;
691 // noalias(R) = ublas::prod(AT, XTnew) - BT;
692 R = -BT;
693 ublas::axpy_prod(AT, XTnew, R, false);
694 nr = ublas::norm_inf(R);
695 }
696 }
697 }
698 return multivector_t(ublas::trans(XT), our_frame);
699 }
700 else
701 // AT is singular. Return NaN
702 return traits_t::NaN();
703 }
704
706 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
707 inline
708 auto
710 operator/= (const multivector_t& rhs) -> multivector_t&
711 { return *this = *this / rhs; }
712
714 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
715 inline
716 auto
718 { return rhs * lhs / rhs.involute(); }
719
721 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
722 inline
723 auto
725 operator|= (const multivector_t& rhs) -> multivector_t&
726 { return *this = rhs * *this / rhs.involute(); }
727
729 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
730 inline
731 auto
733 inv() const -> const multivector_t
734 { return multivector_t(Scalar_T(1), this->m_frame) / *this; }
735
737 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
738 inline
739 auto
741 pow(int m) const -> const multivector_t
742 { return glucat::pow(*this, m); }
743
745 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
746 auto
748 outer_pow(int m) const -> const multivector_t
749 {
750 if (m < 0)
751 throw error_t("outer_pow(m): negative exponent");
752 framed_multi_t a = *this;
753 return a.outer_pow(m);
754 }
755
757 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
758 inline
759 auto
761 grade() const -> index_t
762 { return framed_multi_t(*this).grade(); }
763
765 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
766 inline
767 auto
769 frame() const -> const index_set_t
770 { return this->m_frame; }
771
773 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
774 inline
775 auto
777 operator[] (const index_set_t ist) const -> Scalar_T
778 {
779 // Use matrix inner product only if ist is in frame
780 if ( (ist | this->m_frame) == this->m_frame)
781 return matrix::inner<Scalar_T>(this->basis_element(ist), this->m_matrix);
782 else
783 return Scalar_T(0);
784 }
785
787 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
788 inline
789 auto
791 operator() (index_t grade) const -> const multivector_t
792 {
793 if ((grade < 0) || (grade > HI-LO))
794 return 0;
795 else
796 return (framed_multi_t(*this))(grade);
797 }
798
800 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
801 inline
802 auto
804 scalar() const -> Scalar_T
805 {
806 const matrix_index_t dim = this->m_matrix.size1();
807 return matrix::trace(this->m_matrix) / Scalar_T( double(dim) );
808 }
809
811 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
812 inline
813 auto
815 pure() const -> const multivector_t
816 { return *this - this->scalar(); }
817
819 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
820 inline
821 auto
823 even() const -> const multivector_t
824 { return framed_multi_t(*this).even(); }
825
827 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
828 inline
829 auto
831 odd() const -> const multivector_t
832 { return framed_multi_t(*this).odd(); }
833
835 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
836 auto
838 vector_part() const -> const vector_t
839 { return this->vector_part(this->frame(), true); }
840
842 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
843 auto
845 vector_part(const index_set_t frm, const bool prechecked) const -> const vector_t
846 {
847 if (!prechecked && (this->frame() | frm) != frm)
848 throw error_t("vector_part(frm): value is outside of requested frame");
849 vector_t result;
850 // If we need to enlarge the frame we may as well use a framed_multi_t
851 if (this->frame() != frm)
852 return framed_multi_t(*this).vector_part(frm, true);
853
854 const auto begin_index = frm.min();
855 const auto end_index = frm.max()+1;
856 for (auto
857 idx = begin_index;
858 idx != end_index;
859 ++idx)
860 if (frm[idx])
861 // Frame may contain indices which do not correspond to a grade 1 term but
862 // frame cannot omit any index corresponding to a grade 1 term
863 result.push_back(
864 matrix::inner<Scalar_T>(this->basis_element(index_set_t(idx)),
865 this->m_matrix));
866 return result;
867 }
868
870 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
871 inline
872 auto
874 involute() const -> const multivector_t
875 { return framed_multi_t(*this).involute(); }
876
878 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
879 inline
880 auto
882 reverse() const -> const multivector_t
883 { return framed_multi_t(*this).reverse(); }
884
886 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
887 inline
888 auto
890 conj() const -> const multivector_t
891 { return framed_multi_t(*this).conj(); }
892
894 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
895 inline
896 auto
898 quad() const -> Scalar_T
899 { // scalar(conj(x)*x) = 2*quad(even(x)) - quad(x)
900 // Arvind Raja ref: "old clical: quadfunction(p:pter):pterm in file compmod.pas"
901 return framed_multi_t(*this).quad();
902 }
903
905 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
906 inline
907 auto
909 norm() const -> Scalar_T
910 {
911 const matrix_index_t dim = this->m_matrix.size1();
912 return matrix::norm_frob2(this->m_matrix) / Scalar_T( double(dim) );
913 }
914
916 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
917 inline
918 auto
920 max_abs() const -> Scalar_T
921 { return framed_multi_t(*this).max_abs(); }
922
924 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
925 auto
931
933 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
934 inline
935 void
937 write(const std::string& msg) const
938 { framed_multi_t(*this).write(msg); }
939
941 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
942 inline
943 void
945 write(std::ofstream& ofile, const std::string& msg) const
946 {
947 if (!ofile)
948 throw error_t("write(ofile,msg): cannot write to output file");
949 framed_multi_t(*this).write(ofile, msg);
950 }
951
953 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
954 inline
955 auto
956 operator<< (std::ostream& os, const matrix_multi<Scalar_T,LO,HI,Tune_P>& val) -> std::ostream&
957 {
958 os << typename matrix_multi<Scalar_T,LO,HI,Tune_P>::framed_multi_t(val);
959 return os;
960 }
961
963 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
964 inline
965 auto
966 operator>> (std::istream& s, matrix_multi<Scalar_T,LO,HI,Tune_P>& val) -> std::istream&
967 { // Input looks like 1.0-2.0{1,2}+3.2{3,4}
969 s >> local;
970 // If s.bad() then we have a corrupt input
971 // otherwise we are fine and can copy the resulting matrix_multi
972 if (!s.bad())
973 val = local;
974 return s;
975 }
976
978 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
979 inline
980 auto
982 isinf() const -> bool
983 {
984 if (std::numeric_limits<Scalar_T>::has_infinity)
985 return matrix::isinf(this->m_matrix);
986 else
987 return false;
988 }
989
991 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
992 inline
993 auto
995 isnan() const -> bool
996 {
997 if (std::numeric_limits<Scalar_T>::has_quiet_NaN)
998 return matrix::isnan(this->m_matrix);
999 else
1000 return false;
1001 }
1002
1004 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1005 inline
1006 auto
1008 truncated(const Scalar_T& limit) const -> const multivector_t
1009 { return framed_multi_t(*this).truncated(limit); }
1010
1012 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1013 inline
1014 auto
1016 operator+= (const term_t& term) -> multivector_t&
1017 {
1018 if (term.second != Scalar_T(0))
1019 this->m_matrix.plus_assign(matrix_t(this->basis_element(term.first)) * term.second);
1020 return *this;
1021 }
1022
1024 template< typename Multivector_T, typename Matrix_T, typename Basis_Matrix_T >
1025 static
1026 auto
1027 fast(const Matrix_T& X, index_t level) -> Multivector_T
1028 {
1029 using framed_multi_t = Multivector_T;
1030
1031 using index_set_t = typename framed_multi_t::index_set_t;
1032 using Scalar_T = typename framed_multi_t::scalar_t;
1033 using matrix_t = Matrix_T;
1034 using basis_matrix_t = Basis_Matrix_T;
1035 using basis_scalar_t = typename basis_matrix_t::value_type;
1036 using traits_t = numeric_traits<Scalar_T>;
1037
1038 if (level == 0)
1039 return framed_multi_t(traits_t::to_scalar_t(X(0,0)));
1040
1041 if (ublas::norm_inf(X) == 0)
1042 return Scalar_T(0);
1043
1044 const basis_matrix_t& I = matrix::unit<basis_matrix_t>(2);
1045 basis_matrix_t J(2,2,2);
1046 J.clear();
1047 J(0,1) = basis_scalar_t(-1);
1048 J(1,0) = basis_scalar_t( 1);
1049 basis_matrix_t K = J;
1050 K(0,1) = basis_scalar_t( 1);
1051 basis_matrix_t JK = I;
1052 JK(0,0) = basis_scalar_t(-1);
1053
1055 const index_set_t ist_mn = index_set_t(-level);
1056 const index_set_t ist_pn = index_set_t(level);
1057 const index_set_t ist_mnpn = ist_mn | ist_pn;
1058 if (level == 1)
1059 {
1060 using term_t = typename framed_multi_t::term_t;
1061 const Scalar_T i_x = traits_t::to_scalar_t(signed_perm_nork(I, X)(0, 0));
1062 const Scalar_T j_x = traits_t::to_scalar_t(signed_perm_nork(J, X)(0, 0));
1063 const Scalar_T k_x = traits_t::to_scalar_t(signed_perm_nork(K, X)(0, 0));
1064 const Scalar_T jk_x = traits_t::to_scalar_t(signed_perm_nork(JK,X)(0, 0));
1065 framed_multi_t
1066 result = i_x;
1067 result += term_t(ist_mn, j_x); // j_x * mn;
1068 result += term_t(ist_pn, k_x); // k_x * pn;
1069 return result += term_t(ist_mnpn, jk_x); // jk_x * mnpn;
1070 }
1071 else
1072 {
1073 const framed_multi_t& mn = framed_multi_t(ist_mn);
1074 const framed_multi_t& pn = framed_multi_t(ist_pn);
1075 const framed_multi_t& mnpn = framed_multi_t(ist_mnpn);
1076 const framed_multi_t& i_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1077 (signed_perm_nork(I, X), level-1);
1078 const framed_multi_t& j_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1079 (signed_perm_nork(J, X), level-1);
1080 const framed_multi_t& k_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1081 (signed_perm_nork(K, X), level-1);
1082 const framed_multi_t& jk_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1083 (signed_perm_nork(JK,X), level-1);
1084 framed_multi_t
1085 result = i_x.even() - jk_x.odd();
1086 result += (j_x.even() - k_x.odd()) * mn;
1087 result += (k_x.even() - j_x.odd()) * pn;
1088 return result += (jk_x.even() - i_x.odd()) * mnpn;
1089 }
1090 }
1091
1093 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1094 inline
1095 auto
1097 fast_matrix_multi(const index_set_t frm) const -> const multivector_t
1098 {
1099 if (this->m_frame == frm)
1100 return *this;
1101 else
1102 return (this->template fast_framed_multi<Scalar_T>()).template fast_matrix_multi<Scalar_T>(frm);
1103 }
1104
1106 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1107 template <typename Other_Scalar_T>
1108 auto
1110 fast_framed_multi() const -> const framed_multi<Other_Scalar_T,LO,HI,Tune_P>
1111 {
1112 // Determine the amount of off-centering needed
1113 index_t p = this->m_frame.count_pos();
1114 index_t q = this->m_frame.count_neg();
1115
1116 const index_t bott = pos_mod(p-q, 8);
1117 p += std::max(gen::offset_to_super[bott],index_t(0));
1118 q -= std::min(gen::offset_to_super[bott],index_t(0));
1119
1120 const index_t orig_p = p;
1121 const index_t orig_q = q;
1122 while (p-q > 4)
1123 { p -= 4; q += 4; }
1124 while (p-q < -3)
1125 { p += 4; q -= 4; }
1126 if (p-q > 1)
1127 {
1128 index_t old_p = p;
1129 p = q+1;
1130 q = old_p-1;
1131 }
1132 const index_t level = (p+q)/2;
1133
1134 // Do the inverse fast transform
1136 framed_multi_t val = fast<framed_multi_t, matrix_t, basis_matrix_t>(this->m_matrix, level);
1137
1138 // Off-centre val
1139 switch (pos_mod(orig_p-orig_q, 8))
1140 {
1141 case 2:
1142 case 3:
1143 case 4:
1144 val.centre_qp1_pm1(p, q);
1145 break;
1146 default:
1147 break;
1148 }
1149 if (orig_p-orig_q > 4)
1150 while (p != orig_p)
1151 val.centre_pp4_qm4(p, q);
1152 if (orig_p-orig_q < -3)
1153 while (p != orig_p)
1154 val.centre_pm4_qp4(p, q);
1155
1156 // Return unfolded val
1157 return val.unfold(this->m_frame);
1158 }
1159
1161 template< typename Scalar_T, const index_t LO, const index_t HI, typename Matrix_T >
1163 public std::map< std::pair< const index_set<LO,HI>, const index_set<LO,HI> >,
1164 Matrix_T* >
1165 {
1166 public:
1168 static auto basis() -> basis_table& { static basis_table b; return b;}
1169 private:
1174 // Enforce singleton
1175 // Reference: A. Alexandrescu, "Modern C++ Design", Chapter 6
1176 basis_table() = default;
1177 ~basis_table() = default;
1178 public:
1179 basis_table(const basis_table&) = delete;
1180 auto operator= (const basis_table&) -> basis_table& = delete;
1181 };
1182
1184 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1185 auto
1187 basis_element(const index_set_t& ist) const -> const basis_matrix_t
1188 {
1189 using index_set_pair_t = std::pair<const index_set_t, const index_set_t>;
1190 const auto& unfolded_pair = index_set_pair_t(ist, this->m_frame);
1191
1193 auto& basis_cache = basis_table_t::basis();
1194
1195 const auto frame_count = this->m_frame.count();
1196 const auto use_cache = frame_count <= index_t(Tune_P::basis_max_count);
1197
1198 if (use_cache)
1199 {
1200 const auto basis_it = basis_cache.find(unfolded_pair);
1201 if (basis_it != basis_cache.end())
1202 return *(basis_it->second);
1203 }
1204 const auto folded_set = ist.fold(this->m_frame);
1205 const auto folded_frame = this->m_frame.fold();
1206 const auto& folded_pair = index_set_pair_t(folded_set, folded_frame);
1207 using basis_pair_t = std::pair<const index_set_pair_t, basis_matrix_t *>;
1208 if (use_cache)
1209 {
1210 const auto basis_it = basis_cache.find(folded_pair);
1211 if (basis_it != basis_cache.end())
1212 {
1213 auto* result_ptr = basis_it->second;
1214 basis_cache.insert(basis_pair_t(unfolded_pair, result_ptr));
1215 return *result_ptr;
1216 }
1217 }
1218 const auto folded_max = folded_frame.max();
1219 const auto folded_min = folded_frame.min();
1220 const auto p = std::max(folded_max, index_t(0));
1221 const auto q = std::max(index_t(-folded_min), index_t(0));
1222 const auto* e = (gen::generator_table<basis_matrix_t>::generator())(p, q);
1223 const auto dim = matrix_index_t(1) << offset_level(p, q);
1224 auto result = matrix::unit<basis_matrix_t>(dim);
1225 for (auto
1226 k = folded_min;
1227 k <= folded_max;
1228 ++k)
1229 if (folded_set[k])
1230 result = matrix::mono_prod(result, e[k]);
1231 if (use_cache)
1232 {
1233 auto* result_ptr = new basis_matrix_t(result);
1234 basis_cache.insert(basis_pair_t(folded_pair, result_ptr));
1235 basis_cache.insert(basis_pair_t(unfolded_pair, result_ptr));
1236 }
1237 return result;
1238 }
1239
1241 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P, const size_t Size >
1242 inline
1243 static
1244 auto
1246 const std::array<Scalar_T, Size>& numer,
1247 const std::array<Scalar_T, Size>& denom,
1249 {
1250 // Pade' approximation
1251 // Reference: [GW], Section 4.3, pp318-322
1252 // Reference: [GL], Section 11.3, p572-576.
1253
1254 using multivector_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
1255 using traits_t = numeric_traits<Scalar_T>;
1256
1257 if (X.isnan())
1258 return traits_t::NaN();
1259
1260 // Array size is assumed to be even
1261 const auto nbr_even_powers = Size/2 - 1;
1262
1263 // Create an array of even powers
1264 auto XX = std::vector<multivector_t>(nbr_even_powers);
1265 XX[0] = X * X;
1266 XX[1] = XX[0] * XX[0];
1267 for (auto
1268 k = size_t(2);
1269 k != nbr_even_powers;
1270 ++k)
1271 XX[k] = XX[k-2] * XX[1];
1272
1273 // Calculate numerator N and denominator D
1274 auto N = multivector_t(numer[1]);
1275 for (auto
1276 k = size_t(0);
1277 k != nbr_even_powers;
1278 ++k)
1279 N += XX[k] * numer[2*k + 3];
1280 N *= X;
1281 N += numer[0];
1282 for (auto
1283 k = size_t(0);
1284 k != nbr_even_powers;
1285 ++k)
1286 N += XX[k] * numer[2*k + 2];
1287 auto D = multivector_t(denom[1]);
1288 for (auto
1289 k = size_t(0);
1290 k != nbr_even_powers;
1291 ++k)
1292 D += XX[k] * denom[2*k + 3];
1293 D *= X;
1294 D += denom[0];
1295 for (auto
1296 k = size_t(0);
1297 k != nbr_even_powers;
1298 ++k)
1299 D += XX[k] * denom[2*k + 2];
1300 return N / D;
1301 }
1302
1304 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1305 inline
1306 static
1307 void
1309 {
1310 // Reference: [CHKL]
1311 const auto& invM = inv(M);
1312 M = ((M + invM)/Scalar_T(2) + Scalar_T(1)) / Scalar_T(2);
1313 Y *= (invM + Scalar_T(1)) / Scalar_T(2);
1314 }
1315
1317 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1318 static
1319 auto
1321 Scalar_T norm_tol=std::pow(std::numeric_limits<Scalar_T>::epsilon(), 4)) -> const matrix_multi<Scalar_T,LO,HI,Tune_P>
1322 {
1323 // Reference: [CHKL]
1324 using multivector_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
1325
1326 if (val == Scalar_T(0))
1327 return val;
1328
1329 static const auto sqrt_max_steps = Tune_P::db_sqrt_max_steps;
1330 auto M = val;
1331 auto Y = val;
1332
1333 for (auto
1334 step = 0;
1335 step != sqrt_max_steps && norm(M - Scalar_T(1)) > norm_tol;
1336 ++step)
1337 {
1338 if (Y.isnan())
1340 db_step(M, Y);
1341 }
1342 return Y;
1343 }
1344
1346 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1347 static
1348 auto
1350 Scalar_T norm_Y_tol=std::pow(std::numeric_limits<Scalar_T>::epsilon(), 1)) -> const matrix_multi<Scalar_T,LO,HI,Tune_P>
1351 {
1352 // Reference: [MB]
1353 using multivector_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
1354
1355 if (val == Scalar_T(0))
1356 return val;
1357
1358 static const auto sqrt_max_steps = Tune_P::cr_sqrt_max_steps;
1359 auto Z = Scalar_T(2) * (Scalar_T(1) + val);
1360 auto Y = Scalar_T(1) - val;
1361 using traits_t = numeric_traits<Scalar_T>;
1362 auto norm_Y = norm(Y);
1363 for (auto
1364 step = 0;
1365 step != sqrt_max_steps && norm_Y > norm_Y_tol;
1366 ++step)
1367 {
1368 const auto old_norm_Y = norm_Y;
1369 Y = (-Y / Z) * Y;
1370 norm_Y = norm(Y);
1371 if (Y.isnan() || (norm_Y > old_norm_Y * Scalar_T(2)))
1373
1374 Z += Y * Scalar_T(2);
1375 }
1376 return Z / Scalar_T(4);
1377 }
1378}
1379
1380namespace pade {
1382 // Reference: [Z], Pade1
1383 template< typename Scalar_T >
1385 {
1386 using array = std::array<Scalar_T, 14>;
1387 static const array numer;
1388 };
1389 template< typename Scalar_T >
1391 {
1392 1.0, 27.0/4.0, 81.0/4.0, 2277.0/64.0,
1393 10395.0/256.0, 32319.0/1024.0, 8721.0/512.0, 26163.0/4096.0,
1394 53703.0/32768.0, 36465.0/131072.0, 3861.0/131072.0, 7371.0/4194304.0,
1395 819.0/16777216.0, 27.0/67108864.0
1396 };
1397
1399 // Reference: [Z], Pade1
1400 template< typename Scalar_T >
1402 {
1403 using array = std::array<Scalar_T, 14>;
1404 static const array denom;
1405 };
1406 template< typename Scalar_T >
1408 {
1409 1.0, 25.0/4.0, 69.0/4.0, 1771.0/64.0,
1410 7315.0/256.0, 20349.0/1024.0, 4845.0/512.0, 12597.0/4096.0,
1411 21879.0/32768.0, 12155.0/131072.0, 1001.0/131072.0, 1365.0/4194304.0,
1412 91.0/16777216.0, 1.0/67108864.0
1413 };
1414
1415 template< >
1417 {
1418 using array = std::array<float, 10>;
1419 static const array numer;
1420 };
1422 {
1423 1.0, 19.0/4.0, 19.0/2.0, 665.0/64.0,
1424 1729.0/256.0, 2717.0/1024.0, 627.0/1024.0, 627.0/8192.0,
1425 285.0/65536.0, 19.0/262144.0
1426 };
1427 template< >
1429 {
1430 using array = std::array<float, 10>;
1431 static const array denom;
1432 };
1434 {
1435 1.0, 17.0/4.0, 15.0/2.0, 455.0/64.0,
1436 1001.0/256.0, 1287.0/1024.0, 231.0/1024.0, 165.0/8192.0,
1437 45.0/65536, 1.0/262144.0
1438 };
1439
1440 template< >
1442 {
1443 using array = std::array<long double, 18>;
1444 static const array numer;
1445 };
1447 {
1448 1.0L, 35.0L/4.0L, 35.0L, 5425.0L/64.0L,
1449 35525.0L/256.0L, 166257.0L/1024.0L, 143325.0L/1024.0L, 740025.0L/8192.0L,
1450 2877875.0L/65536.0L, 4206125.0L/262144.0L, 572033.0L/131072.0L, 1820105.0L/2097152.0L,
1451 1028755.0L/8388608.0L, 395675.0L/33554432.0L, 24225.0L/33554432.0L, 6783.0L/268435456.0L,
1452 1785.0L/4294967296.0L, 35.0L/17179869184.0L
1453 };
1454 template< >
1456 {
1457 using array = std::array<long double, 18>;
1458 static const array denom;
1459 };
1461 {
1462 1.0L, 33.0L/4.0L, 31.0L, 4495.0L/64.0L,
1463 27405.0L/256.0L, 118755.0L/1024.0L, 94185.0L/1024.0L, 444015.0L/8192.0L,
1464 1562275.0L/65536.0L, 2042975.0L/262144.0L, 245157.0L/131072.0L, 676039.0L/2097152.0L,
1465 323323.0L/8388608.0L, 101745.0L/33554432.0L, 4845.0L/33554432.0L, 969.0L/268435456.0L,
1466 153.0L/4294967296.0L, 1.0L/17179869184.0L
1467 };
1468
1469#if defined(_GLUCAT_USE_QD)
1470 template< >
1472 {
1473 using array = std::array<dd_real, 22>;
1474 static const array numer;
1475 };
1477 {
1478 dd_real("1"), dd_real("43")/dd_real("4"),
1479 dd_real("215")/dd_real("4"), dd_real("10621")/dd_real("64"),
1480 dd_real("90687")/dd_real("256"), dd_real("567987")/dd_real("1024"),
1481 dd_real("168861")/dd_real("256"), dd_real("1246355")/dd_real("2048"),
1482 dd_real("7228859")/dd_real("16384"), dd_real("16583853")/dd_real("65536"),
1483 dd_real("7538115")/dd_real("65536"), dd_real("173376645")/dd_real("4194304"),
1484 dd_real("195747825")/dd_real("16777216"), dd_real("171655785")/dd_real("67108864"),
1485 dd_real("14375115")/dd_real("33554432"), dd_real("14375115")/dd_real("268435456"),
1486 dd_real("20764055")/dd_real("4294967296"), dd_real("5167525")/dd_real("17179869184"),
1487 dd_real("206701")/dd_real("17179869184"), dd_real("76153")/dd_real("274877906944"),
1488 dd_real("3311")/dd_real("1099511627776") , dd_real("43")/dd_real("4398046511104")
1489 };
1490 template< >
1492 {
1493 using array = std::array<dd_real, 22>;
1494 static const array denom;
1495 };
1497 {
1498 dd_real("1"), dd_real("41")/dd_real("4"),
1499 dd_real("195")/dd_real("4"), dd_real("9139")/dd_real("64"),
1500 dd_real("73815")/dd_real("256"), dd_real("435897")/dd_real("1024"),
1501 dd_real("121737")/dd_real("256"), dd_real("840565")/dd_real("2048"),
1502 dd_real("4539051")/dd_real("16384"), dd_real("9641775")/dd_real("65536"),
1503 dd_real("4032015")/dd_real("65536"), dd_real("84672315")/dd_real("4194304"),
1504 dd_real("86493225")/dd_real("16777216"), dd_real("67863915")/dd_real("67108864"),
1505 dd_real("5014575")/dd_real("33554432"), dd_real("4345965")/dd_real("268435456"),
1506 dd_real("5311735")/dd_real("4294967296"), dd_real("1081575")/dd_real("17179869184"),
1507 dd_real("33649")/dd_real("17179869184"), dd_real("8855")/dd_real("274877906944"),
1508 dd_real("231")/dd_real("1099511627776"), dd_real("1")/dd_real("4398046511104")
1509 };
1510
1511 template< >
1513 {
1514 using array = std::array<qd_real, 34>;
1515 static const array numer;
1516 };
1518 {
1519 qd_real("1"), qd_real("67")/qd_real("4"),
1520 qd_real("134"), qd_real("43617")/qd_real("64"),
1521 qd_real("633485")/qd_real("256"), qd_real("6992857")/qd_real("1024"),
1522 qd_real("15246721")/qd_real("1024"), qd_real("215632197")/qd_real("8192"),
1523 qd_real("2518145487")/qd_real("65536"), qd_real("12301285425")/qd_real("262144"),
1524 qd_real("6344873535")/qd_real("131072"), qd_real("89075432355")/qd_real("2097152"),
1525 qd_real("267226297065")/qd_real("8388608"), qd_real("687479618945")/qd_real("33554432"),
1526 qd_real("379874182975")/qd_real("33554432"), qd_real("1443521895305")/qd_real("268435456"),
1527 qd_real("9425348845815")/qd_real("4294967296"), qd_real("13195488384141")/qd_real("17179869184"),
1528 qd_real("987417498133")/qd_real("4294967296"), qd_real("8055248011085")/qd_real("137438953472"),
1529 qd_real("6958363175533")/qd_real("549755813888"), qd_real("5056698705201")/qd_real("2199023255552"),
1530 qd_real("766166470485")/qd_real("2199023255552"), qd_real("766166470485")/qd_real("17592186044416"),
1531 qd_real("623623871325")/qd_real("140737488355328"), qd_real("203123203803")/qd_real("562949953421312"),
1532 qd_real("6478601247")/qd_real("281474976710656"), qd_real("5038912081")/qd_real("4503599627370496"),
1533 qd_real("719844583")/qd_real("18014398509481984"), qd_real("71853815")/qd_real("72057594037927936"),
1534 qd_real("1165197")/qd_real("72057594037927936"), qd_real("87703")/qd_real("576460752303423488"),
1535 qd_real("12529")/qd_real("18446744073709551616"), qd_real("67")/qd_real("73786976294838206464")
1536 };
1537 template< >
1539 {
1540 using array = std::array<qd_real, 34>;
1541 static const array denom;
1542 };
1544 {
1545 qd_real("1"), qd_real("65")/qd_real("4"),
1546 qd_real("126"), qd_real("39711")/qd_real("64"),
1547 qd_real("557845")/qd_real("256"), qd_real("5949147")/qd_real("1024"),
1548 qd_real("12515965")/qd_real("1024"), qd_real("170574723")/qd_real("8192"),
1549 qd_real("1916797311")/qd_real("65536"), qd_real("8996462475")/qd_real("262144"),
1550 qd_real("4450881435")/qd_real("131072"), qd_real("59826782925")/qd_real("2097152"),
1551 qd_real("171503444385")/qd_real("8388608"), qd_real("420696483235")/qd_real("33554432"),
1552 qd_real("221120793075")/qd_real("33554432"), qd_real("797168807855")/qd_real("268435456"),
1553qd_real("4923689695575")/qd_real("4294967296"), qd_real("6499270398159")/qd_real("17179869184"),
1554 qd_real("456864812569")/qd_real("4294967296"), qd_real("3486599885395")/qd_real("137438953472"),
1555qd_real("2804116503573")/qd_real("549755813888"), qd_real("1886827875075")/qd_real("2199023255552"),
1556 qd_real("263012370465")/qd_real("2199023255552"), qd_real("240141729555")/qd_real("17592186044416"),
1557 qd_real("176848560525")/qd_real("140737488355328"), qd_real("51538723353")/qd_real("562949953421312"),
1558 qd_real("1450433115")/qd_real("281474976710656"), qd_real("977699359")/qd_real("4503599627370496"),
1559 qd_real("118183439")/qd_real("18014398509481984"), qd_real("9652005")/qd_real("72057594037927936"),
1560 qd_real("121737")/qd_real("72057594037927936"), qd_real("6545")/qd_real("576460752303423488"),
1561 qd_real("561")/qd_real("18446744073709551616"), qd_real("1")/qd_real("73786976294838206464")
1562 };
1563#endif
1564}
1565
1566namespace glucat
1567{
1569 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1570 auto
1573 const index_t level) -> const matrix_multi<Scalar_T,LO,HI,Tune_P>
1574 {
1575 // Reference: [GW], Section 4.3, pp318-322
1576 // Reference: [GL], Section 11.3, p572-576
1577 // Reference: [Z], Pade1
1578
1579 using traits_t = numeric_traits<Scalar_T>;
1580
1581 if (val.isnan())
1582 return traits_t::NaN();
1583
1584 const auto scr_val = val.scalar();
1585 if (val == scr_val)
1586 {
1587 if (scr_val < Scalar_T(0))
1588 return i * traits_t::sqrt(-scr_val);
1589 else
1590 return traits_t::sqrt(scr_val);
1591 }
1592
1593 // Scale val towards abs(A) == 1 or towards A == 1 as appropriate
1594 const auto scale =
1595 (scr_val != Scalar_T(0) && norm(val/scr_val - Scalar_T(1)) < Scalar_T(1))
1596 ? scr_val
1597 : (scr_val < Scalar_T(0))
1598 ? -abs(val)
1599 : abs(val);
1600 const auto sqrt_scale = traits_t::sqrt(traits_t::abs(scale));
1601 if (traits_t::isNaN_or_isInf(sqrt_scale))
1602 return traits_t::NaN();
1603
1604 using multivector_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
1605 auto rescale = multivector_t(sqrt_scale);
1606 if (scale < Scalar_T(0))
1607 rescale = i * sqrt_scale;
1608
1609 const auto& unitval = val / scale;
1610 static const auto max_norm = Scalar_T(1.0/4.0);
1611 auto use_approx_sqrt = true;
1612 auto use_cr_sqrt = false;
1613 auto scaled_result = multivector_t();
1614#if defined(_GLUCAT_USE_EIGENVALUES)
1615 static const auto sqrt_2 = traits_t::sqrt(Scalar_T(2));
1616 if (level == 0)
1617 {
1618 using matrix_t = typename multivector_t::matrix_t;
1619
1620 // What kind of eigenvalues does the matrix contain?
1621 const auto genus = matrix::classify_eigenvalues(unitval.m_matrix);
1622 const index_t next_level =
1623 (genus.m_is_singular)
1624 ? level
1625 : level + 1;
1626 switch (genus.m_eig_case)
1627 {
1628 case matrix::neg_real_eigs:
1629 scaled_result = matrix_sqrt(-i * unitval, i, next_level) * (i + Scalar_T(1)) / sqrt_2;
1630 use_approx_sqrt = false;
1631 break;
1632 case matrix::both_eigs:
1633 {
1634 const auto safe_arg = genus.m_safe_arg;
1635 scaled_result = matrix_sqrt(exp(i*safe_arg) * unitval, i, next_level) * exp(-i*safe_arg / Scalar_T(2));
1636 }
1637 use_approx_sqrt = false;
1638 break;
1639 default:
1640 break;
1641 }
1642 use_cr_sqrt = genus.m_is_singular;
1643 }
1644#endif
1645 if (use_approx_sqrt)
1646 {
1647 scaled_result =
1648 (norm(unitval - Scalar_T(1)) < max_norm)
1649 // Pade' approximation of square root
1652 unitval - Scalar_T(1))
1653 // Product form of Denman-Beavers square root iteration
1654 : (use_cr_sqrt)
1655 ? cr_sqrt(unitval)
1656 : db_sqrt(unitval);
1657 }
1658 return (scaled_result.isnan() ||
1659 !approx_equal(pow(scaled_result, 2), unitval))
1660 ? traits_t::NaN()
1661 : scaled_result * rescale;
1662 }
1663
1665 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1666 auto
1668 {
1669 // Reference: [GW], Section 4.3, pp318-322
1670 // Reference: [GL], Section 11.3, p572-576
1671 // Reference: [Z], Pade1
1672
1673 using traits_t = numeric_traits<Scalar_T>;
1674
1675 if (val.isnan())
1676 return traits_t::NaN();
1677
1678 check_complex(val, i, prechecked);
1679
1680 switch (Tune_P::function_precision)
1681 {
1682 case precision_demoted:
1683 {
1684 using demoted_scalar_t = typename traits_t::demoted::type;
1685 using demoted_multivector_t = matrix_multi<demoted_scalar_t,LO,HI,Tune_P>;
1686
1687 const auto& demoted_val = demoted_multivector_t(val);
1688 const auto& demoted_i = demoted_multivector_t(i);
1689
1690 return matrix_sqrt(demoted_val, demoted_i, 0);
1691 }
1692 break;
1693 case precision_promoted:
1694 {
1695 using promoted_scalar_t = typename traits_t::promoted::type;
1696 using promoted_multivector_t = matrix_multi<promoted_scalar_t,LO,HI,Tune_P>;
1697
1698 const auto& promoted_val = promoted_multivector_t(val);
1699 const auto& promoted_i = promoted_multivector_t(i);
1700
1701 return matrix_sqrt(promoted_val, promoted_i, 0);
1702 }
1703 break;
1704 default:
1705 return matrix_sqrt(val, i, 0);
1706 }
1707 }
1708}
1709
1710namespace pade {
1712 // Reference: [Z], Pade1
1713 template< typename Scalar_T >
1715 {
1716 using array = std::array<Scalar_T, 14>;
1717 static const array numer;
1718 };
1719 template< typename Scalar_T >
1721 {
1722 0.0, 1.0, 6.0, 4741.0/300.0,
1723 1441.0/60.0, 107091.0/4600.0, 8638.0/575.0, 263111.0/40250.0,
1724 153081.0/80500.0, 395243.0/1101240.0, 28549.0/688275.0, 605453.0/228813200.0,
1725 785633.0/10296594000.0, 1145993.0/1873980108000.0
1726 };
1727
1729 // Reference: [Z], Pade1
1730 template< typename Scalar_T >
1732 {
1733 using array = std::array<Scalar_T, 14>;
1734 static const array denom;
1735 };
1736 template< typename Scalar_T >
1738 {
1739 1.0, 13.0/2.0, 468.0/25.0, 1573.0/50.0,
1740 1573.0/46.0, 11583.0/460.0, 10296.0/805.0, 2574.0/575.0,
1741 11583.0/10925.0, 143.0/874.0, 572.0/37145.0, 117.0/148580.0,
1742 13.0/742900.0, 1.0/10400600.0
1743 };
1744
1745 template< >
1747 {
1748 using array = std::array<float, 10>;
1749 static const array numer;
1750 };
1752 {
1753 0.0, 1.0, 4.0, 1337.0/204.0,
1754 385.0/68.0, 1879.0/680.0, 193.0/255.0, 197.0/1820.0,
1755 419.0/61880.0, 7129.0/61261200.0
1756 };
1757 template< >
1759 {
1760 using array = std::array<float, 10>;
1761 static const array denom;
1762 };
1764 {
1765 1.0, 9.0/2.0, 144.0/17.0, 147.0/17.0,
1766 441.0/85.0, 63.0/34.0, 84.0/221.0, 9.0/221.0,
1767 9.0/4862.0, 1.0/48620.0
1768 };
1769
1770 template< >
1772 {
1773 using array = std::array<long double, 18>;
1774 static const array numer;
1775 };
1777 {
1778 0.0L, 1.0L, 8.0L, 3835.0L/132.0L,
1779 8365.0L/132.0L, 11363807.0L/122760.0L, 162981.0L/1705.0L, 9036157.0L/125860.0L,
1780 18009875.0L/453096.0L, 44211925.0L/2718576.0L, 4149566.0L/849555.0L, 16973929.0L/16020180.0L,
1781 172459.0L/1068012.0L, 116317061.0L/7025382936.0L, 19679783.0L/18441630207.0L, 23763863.0L/614721006900.0L,
1782 50747.0L/79318839600.0L, 42142223.0L/14295951736466400.0L
1783 };
1784 template< >
1786 {
1787 using array = std::array<long double, 18>;
1788 static const array denom;
1789 };
1791 {
1792 1.0L, 17.0L/2.0L, 1088.0L/33.0L, 850.0L/11.0L,
1793 41650.0L/341.0L, 140777.0L/1023.0L, 1126216.0L/9889.0L, 63206.0L/899.0L,
1794 790075.0L/24273.0L, 60775.0L/5394.0L, 38896.0L/13485.0L, 21658.0L/40455.0L,
1795 21658.0L/310155.0L, 4165.0L/682341.0L, 680.0L/2047023.0L, 34.0L/3411705.0L,
1796 17.0L/129644790.0L, 1.0L/2333606220
1797 };
1798#if defined(_GLUCAT_USE_QD)
1799 template< >
1801 {
1802 using array = std::array<dd_real, 22>;
1803 static const array numer;
1804 };
1806 {
1807 dd_real("0"), dd_real("1"),
1808 dd_real("10"), dd_real("22781")/dd_real("492"),
1809 dd_real("21603")/dd_real("164"), dd_real("5492649")/dd_real("21320"),
1810 dd_real("978724")/dd_real("2665"), dd_real("4191605")/dd_real("10619"),
1811 dd_real("12874933")/dd_real("39442"), dd_real("11473457")/dd_real("54612"),
1812 dd_real("2406734")/dd_real("22755"), dd_real("166770367")/dd_real("4004880"),
1813 dd_real("30653165")/dd_real("2402928"), dd_real("647746389")/dd_real("215195552"),
1814 dd_real("25346331")/dd_real("47074027"), dd_real("278270613")/dd_real("3900419380"),
1815 dd_real("105689791")/dd_real("15601677520"), dd_real("606046475")/dd_real("1379188292768"),
1816 dd_real("969715")/dd_real("53502994116"), dd_real("11098301")/dd_real("26204577562592"),
1817 dd_real("118999")/dd_real("26204577562592"), dd_real("18858053")/dd_real("1392249205900512960")
1818 };
1819 template< >
1821 {
1822 using array = std::array<dd_real, 22>;
1823 static const array denom;
1824 };
1826 {
1827 dd_real("1"), dd_real("21")/dd_real("2"),
1828 dd_real("2100")/dd_real("41"), dd_real("12635")/dd_real("82"),
1829 dd_real("341145")/dd_real("1066"), dd_real("1037799")/dd_real("2132"),
1830 dd_real("11069856")/dd_real("19721"), dd_real("9883800")/dd_real("19721"),
1831 dd_real("6918660")/dd_real("19721"), dd_real("293930")/dd_real("1517"),
1832 dd_real("1410864")/dd_real("16687"), dd_real("88179")/dd_real("3034"),
1833 dd_real("734825")/dd_real("94054"), dd_real("305235")/dd_real("188108"),
1834 dd_real("348840")/dd_real("1363783"), dd_real("40698")/dd_real("1363783"),
1835 dd_real("6783")/dd_real("2727566"), dd_real("9975")/dd_real("70916716"),
1836 dd_real("266")/dd_real("53187537"), dd_real("7")/dd_real("70916716"),
1837 dd_real("7")/dd_real("8155422340"), dd_real("1")/dd_real("538257874440")
1838 };
1839
1840 template< >
1842 {
1843 using array = std::array<qd_real, 34>;
1844 static const array numer;
1845 };
1847 {
1848 qd_real("0"), qd_real("1"),
1849 qd_real("16"), qd_real("95201")/qd_real("780"),
1850 qd_real("30721")/qd_real("52"), qd_real("7416257")/qd_real("3640"),
1851 qd_real("1039099")/qd_real("195"), qd_real("6097772319")/qd_real("555100"),
1852 qd_real("1564058073")/qd_real("85400"), qd_real("30404640205")/qd_real("1209264"),
1853 qd_real("725351278")/qd_real("25193"), qd_real("4092322670789")/qd_real("147429436"),
1854 qd_real("4559713849589")/qd_real("201040140"), qd_real("5049361751189")/qd_real("320023080"),
1855 qd_real("74979677195")/qd_real("8000577"), qd_real("16569850691873")/qd_real("3481514244"),
1856 qd_real("1065906022369")/qd_real("515779888"), qd_real("335956770855841")/qd_real("438412904800"),
1857 qd_real("1462444287585964")/qd_real("6041877844275"), qd_real("397242326339851")/qd_real("6122436215532"),
1858 qd_real("64211291334131")/qd_real("4373168725380"), qd_real("142322343550859")/qd_real("51080680851480"),
1859 qd_real("154355972958659")/qd_real("351179680853925"), qd_real("167483568676259")/qd_real("2937139148960100"),
1860 qd_real("4230788929433")/qd_real("704913395750424"), qd_real("197968763176019")/qd_real("392923948371995600"),
1861 qd_real("10537522306718")/qd_real("319250708052246425"), qd_real("236648286272519")/qd_real("144249197475035425500"),
1862 qd_real("260715545088119")/qd_real("4375558990076074573500"), qd_real("289596255666839")/qd_real("192874640282553367199880"),
1863 qd_real("8802625510547")/qd_real("361639950529787563499775"), qd_real("373831661521439")/qd_real("1659204093030665341336967700"),
1864 qd_real("446033437968239")/qd_real("464577146048586295574350956000"), qd_real("53676090078349")/qd_real("47386868896955802148583797512000")
1865 };
1866 template< >
1868 {
1869 using array = std::array<qd_real, 34>;
1870 static const array denom;
1871 };
1873 {
1874 qd_real("1"), qd_real("33")/qd_real("2"),
1875 qd_real("8448")/qd_real("65"), qd_real("42284")/qd_real("65"),
1876 qd_real("211420")/qd_real("91"), qd_real("573562")/qd_real("91"),
1877 qd_real("32119472")/qd_real("2379"), qd_real("92917044")/qd_real("3965"),
1878 qd_real("603960786")/qd_real("17995"), qd_real("144626625")/qd_real("3599"),
1879 qd_real("2776831200")/qd_real("68381"), qd_real("16692542100")/qd_real("478667"),
1880 qd_real("12241197540")/qd_real("478667"), qd_real("1098569010")/qd_real("68381"),
1881 qd_real("31387686000")/qd_real("3624193"), qd_real("9939433900")/qd_real("2479711"),
1882 qd_real("67091178825")/qd_real("42155087"), qd_real("2683647153")/qd_real("4959422"),
1883 qd_real("19083713088")/qd_real("121505839"), qd_real("4708152900")/qd_real("121505839"),
1884 qd_real("941630580")/qd_real("116546417"), qd_real("88704330")/qd_real("62755763"),
1885 qd_real("12902448")/qd_real("62755763"), qd_real("1542684")/qd_real("62755763"),
1886 qd_real("6427850")/qd_real("2698497809"), qd_real("3471039")/qd_real("18889484663"),
1887 qd_real("8544096")/qd_real("774468871183"), qd_real("39556")/qd_real("79027435835"),
1888 qd_real("118668")/qd_real("7191496660985"), qd_real("10230")/qd_real("27327687311743"),
1889 qd_real("5456")/qd_real("1011124430534491"), qd_real("44")/qd_real("1011124430534491"),
1890 qd_real("11")/qd_real("70778710137414370"), qd_real("1")/qd_real("7219428434016265740")
1891 };
1892#endif
1893}
1894
1895namespace glucat{
1897 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1898 static
1899 auto
1901 {
1902 // Reference: [GW], Section 4.3, pp318-322
1903 // Reference: [CHKL]
1904 // Reference: [GL], Section 11.3, p572-576
1905 // Reference: [Z], Pade1
1906
1907 using traits_t = numeric_traits<Scalar_T>;
1908 if (val == Scalar_T(0) || val.isnan())
1909 return traits_t::NaN();
1910 else
1913 val - Scalar_T(1));
1914 }
1915
1917 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1918 static
1919 auto
1921 {
1922 // Reference: [CHKL]
1923 using multivector_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
1924 using traits_t = numeric_traits<Scalar_T>;
1925 if (val == Scalar_T(0) || val.isnan())
1926 return traits_t::NaN();
1927
1928 using limits_t = std::numeric_limits<Scalar_T>;
1929 static const auto epsilon = limits_t::epsilon();
1930 static const auto max_inner_norm = traits_t::pow(epsilon, 2);
1931 static const auto max_outer_norm = Scalar_T(6.0/limits_t::digits);
1932 auto Y = val;
1933 auto E = multivector_t(Scalar_T(0));
1934 Scalar_T norm_Y_1;
1935 auto pow_2_outer_step = Scalar_T(1);
1936 auto pow_4_outer_step = Scalar_T(1);
1937 int outer_step;
1938 for (outer_step = 0, norm_Y_1 = norm(Y - Scalar_T(1));
1939 outer_step != Tune_P::log_max_outer_steps && norm_Y_1 * pow_2_outer_step > max_outer_norm;
1940 ++outer_step, norm_Y_1 = norm(Y - Scalar_T(1)))
1941 {
1942 if (Y == Scalar_T(0) || Y.isnan())
1943 return traits_t::NaN();
1944
1945 // Incomplete product form of Denman-Beavers square root iteration
1946 auto M = Y;
1947 for (auto
1948 inner_step = 0;
1949 inner_step != Tune_P::log_max_inner_steps &&
1950 norm(M - Scalar_T(1)) * pow_4_outer_step > max_inner_norm;
1951 ++inner_step)
1952 db_step(M, Y);
1953
1954 E += (M - Scalar_T(1)) * pow_2_outer_step;
1955 pow_2_outer_step *= Scalar_T(2);
1956 pow_4_outer_step *= Scalar_T(4);
1957 }
1958 if (outer_step == Tune_P::log_max_outer_steps && norm_Y_1 * pow_2_outer_step > max_outer_norm)
1959 return traits_t::NaN();
1960 else
1961 return pade_log(Y) * pow_2_outer_step - E;
1962 }
1963
1965 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1966 auto
1969 const index_t level) -> const matrix_multi<Scalar_T,LO,HI,Tune_P>
1970 {
1971 // Scaled incomplete square root cascade and scaled Pade' approximation of log
1972 // Reference: [CHKL]
1973
1974 using traits_t = numeric_traits<Scalar_T>;
1975 if (val == Scalar_T(0) || val.isnan())
1976 return traits_t::NaN();
1977
1978 static const auto pi = traits_t::pi();
1979 const auto scr_val = val.scalar();
1980 if (val == scr_val)
1981 {
1982 if (scr_val < Scalar_T(0))
1983 return i * pi + traits_t::log(-scr_val);
1984 else
1985 return traits_t::log(scr_val);
1986 }
1987
1988 // Scale val towards abs(A) == 1 or towards A == 1 as appropriate
1989 const auto max_norm = Scalar_T(1.0/9.0);
1990 const auto scale =
1991 (scr_val != Scalar_T(0) && norm(val/scr_val - Scalar_T(1)) < max_norm)
1992 ? scr_val
1993 : (scr_val < Scalar_T(0))
1994 ? -abs(val)
1995 : abs(val);
1996 if (scale == Scalar_T(0))
1997 return traits_t::NaN();
1998
1999 using multivector_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
2000 const auto log_scale = traits_t::log(traits_t::abs(scale));
2001 auto rescale = multivector_t(log_scale);
2002 if (scale < Scalar_T(0))
2003 rescale = i * pi + log_scale;
2004 const auto unitval = val/scale;
2005 if (inv(unitval).isnan())
2006 return traits_t::NaN();
2007
2008#if defined(_GLUCAT_USE_EIGENVALUES)
2009 auto scaled_result = multivector_t();
2010 if (level == 0)
2011 {
2012 using matrix_t = typename multivector_t::matrix_t;
2013
2014 // What kind of eigenvalues does the matrix contain?
2015 auto genus = matrix::classify_eigenvalues(unitval.m_matrix);
2016 switch (genus.m_eig_case)
2017 {
2018 case matrix::neg_real_eigs:
2019 scaled_result = matrix_log(-i * unitval, i, level + 1) + i * pi/Scalar_T(2);
2020 break;
2021 case matrix::both_eigs:
2022 {
2023 const Scalar_T safe_arg = genus.m_safe_arg;
2024 scaled_result = matrix_log(exp(i*safe_arg) * unitval, i, level + 1) - i * safe_arg;
2025 }
2026 break;
2027 default:
2028 scaled_result = cascade_log(unitval);
2029 break;
2030 }
2031 }
2032 else
2033 scaled_result = cascade_log(unitval);
2034#else
2035 auto scaled_result = cascade_log(unitval);
2036#endif
2037 return (scaled_result.isnan())
2038 ? traits_t::NaN()
2039 : scaled_result + rescale;
2040 }
2041
2043 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
2044 auto
2046 {
2047 using traits_t = numeric_traits<Scalar_T>;
2048
2049 if (val == Scalar_T(0) || val.isnan())
2050 return traits_t::NaN();
2051
2052 check_complex(val, i, prechecked);
2053
2054 switch (Tune_P::function_precision)
2055 {
2056 case precision_demoted:
2057 {
2058 using demoted_scalar_t = typename traits_t::demoted::type;
2059 using demoted_multivector_t = matrix_multi<demoted_scalar_t,LO,HI,Tune_P>;
2060
2061 const auto& demoted_val = demoted_multivector_t(val);
2062 const auto& demoted_i = demoted_multivector_t(i);
2063
2064 return matrix_log(demoted_val, demoted_i, 0);
2065 }
2066 break;
2067 case precision_promoted:
2068 {
2069 using promoted_scalar_t = typename traits_t::promoted::type;
2070 using promoted_multivector_t = matrix_multi<promoted_scalar_t,LO,HI,Tune_P>;
2071
2072 const auto& promoted_val = promoted_multivector_t(val);
2073 const auto& promoted_i = promoted_multivector_t(i);
2074
2075 return matrix_log(promoted_val, promoted_i, 0);
2076 }
2077 break;
2078 default:
2079 return matrix_log(val, i, 0);
2080 }
2081 }
2082
2084 template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
2085 auto
2087 {
2088 using traits_t = numeric_traits<Scalar_T>;
2089 if (val.isnan())
2090 return traits_t::NaN();
2091
2092 const auto scr_val = val.scalar();
2093 if (val == scr_val)
2094 return traits_t::exp(scr_val);
2095
2096 switch (Tune_P::function_precision)
2097 {
2098 case precision_demoted:
2099 {
2100 using demoted_scalar_t = typename traits_t::demoted::type;
2101 using demoted_multivector_t = matrix_multi<demoted_scalar_t,LO,HI,Tune_P>;
2102
2103 const auto& demoted_val = demoted_multivector_t(val);
2104 return clifford_exp(demoted_val);
2105 }
2106 break;
2107 case precision_promoted:
2108 {
2109 using promoted_scalar_t = typename traits_t::promoted::type;
2110 using promoted_multivector_t = matrix_multi<promoted_scalar_t,LO,HI,Tune_P>;
2111
2112 const auto& promoted_val = promoted_multivector_t(val);
2113 return clifford_exp(promoted_val);
2114 }
2115 break;
2116 default:
2117 return clifford_exp(val);
2118 }
2119 }
2120}
2121#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 void write(const std::string &msg="") const =0
Write formatted multivector to output.
virtual auto involute() const -> const multivector_t=0
Main involution, each {i} is replaced by -{i} in each term, eg. {1} -> -{1}.
virtual auto outer_pow(int m) const -> const multivector_t=0
Outer product power.
virtual auto isinf() const -> bool=0
Check if a multivector contains any infinite values.
virtual auto quad() const -> Scalar_T=0
Scalar_T quadratic form == (rev(x)*x)(0)
virtual auto scalar() const -> Scalar_T=0
Scalar part.
virtual auto operator()(index_t grade) const -> const multivector_t=0
Pure grade-vector part.
virtual auto conj() const -> const multivector_t=0
Conjugation, reverse o involute == involute o reverse.
virtual auto operator-=(const multivector_t &rhs) -> multivector_t &=0
Geometric difference.
virtual auto grade() const -> index_t=0
Maximum of the grades of each term.
virtual auto even() const -> const multivector_t=0
Even part of multivector, sum of even grade terms.
virtual auto frame() const -> const index_set_t=0
Subalgebra generated by all generators of terms of given multivector.
virtual auto operator/=(const Scalar_T &scr) -> multivector_t &=0
Quotient of multivector and scalar.
virtual auto max_abs() const -> Scalar_T=0
Maximum of absolute values of components of multivector: multivector infinity norm.
virtual auto odd() const -> const multivector_t=0
Odd part of multivector, sum of odd grade terms.
virtual auto vector_part() const -> const vector_t=0
Vector part of multivector, as a vector_t with respect to frame()
virtual auto reverse() const -> const multivector_t=0
Reversion, eg. {1}*{2} -> {2}*{1}.
virtual auto operator%=(const multivector_t &rhs) -> multivector_t &=0
Contraction.
virtual auto operator-() const -> const multivector_t=0
Unary -.
virtual auto truncated(const Scalar_T &limit=default_truncation) const -> const multivector_t=0
Remove all terms with relative size smaller than limit.
virtual auto operator|=(const multivector_t &rhs) -> multivector_t &=0
Transformation via twisted adjoint action.
virtual auto pure() const -> const multivector_t=0
Pure part.
virtual auto operator==(const multivector_t &val) const -> bool=0
Test for equality of multivectors.
virtual auto operator^=(const multivector_t &rhs) -> multivector_t &=0
Outer product.
virtual auto isnan() const -> bool=0
Check if a multivector contains any IEEE NaN values.
virtual auto operator*=(const Scalar_T &scr) -> multivector_t &=0
Product of multivector and scalar.
virtual auto norm() const -> Scalar_T=0
Scalar_T norm == sum of norm of coordinates.
virtual auto pow(int m) const -> const multivector_t=0
*this to the m
virtual auto operator&=(const multivector_t &rhs) -> multivector_t &=0
Inner product.
virtual auto operator[](const index_set_t ist) const -> Scalar_T=0
Subscripting: map from index set to scalar coordinate.
virtual auto inv() const -> const multivector_t=0
Geometric multiplicative inverse.
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
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
error< multivector_t > error_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 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 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:228
auto trace(const Matrix_T &val) -> typename Matrix_T::value_type
Matrix trace.
Definition matrix_imp.h:416
auto isinf(const Matrix_T &m) -> bool
Infinite.
Definition matrix_imp.h:275
auto norm_frob2(const Matrix_T &val) -> typename Matrix_T::value_type
Square of Frobenius norm.
Definition matrix_imp.h:395
auto classify_eigenvalues(const Matrix_T &val) -> eig_genus< Matrix_T >
Classify the eigenvalues of a matrix.
Definition matrix_imp.h:548
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:320
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.
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),...
std::array< Scalar_T, 14 > array
static const array numer
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