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 
65 namespace 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 >
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
927  random(const index_set<LO,HI> frm, Scalar_T fill) -> const multivector_t
928  {
930  }
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
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 >
1162  class basis_table :
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
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 
1192  using basis_table_t = basis_table<Scalar_T, LO, HI, basis_matrix_t>;
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
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
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 
1380 namespace 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< >
1416  struct pade_sqrt_numer<float>
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< >
1428  struct pade_sqrt_denom<float>
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< >
1441  struct pade_sqrt_numer<long double>
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< >
1455  struct pade_sqrt_denom<long double>
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< >
1471  struct pade_sqrt_numer<dd_real>
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< >
1491  struct pade_sqrt_denom<dd_real>
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< >
1512  struct pade_sqrt_numer<qd_real>
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< >
1538  struct pade_sqrt_denom<qd_real>
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"),
1553 qd_real("4923689695575")/qd_real("4294967296"), qd_real("6499270398159")/qd_real("17179869184"),
1554  qd_real("456864812569")/qd_real("4294967296"), qd_real("3486599885395")/qd_real("137438953472"),
1555 qd_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 
1566 namespace 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 
1710 namespace 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< >
1746  struct pade_log_numer<float>
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< >
1758  struct pade_log_denom<float>
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< >
1771  struct pade_log_numer<long double>
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< >
1785  struct pade_log_denom<long double>
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< >
1800  struct pade_log_numer<dd_real>
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< >
1820  struct pade_log_denom<dd_real>
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< >
1841  struct pade_log_numer<qd_real>
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< >
1867  struct pade_log_denom<qd_real>
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 
1895 namespace 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
std::array< float, 10 > array
error< multivector_t > error_t
Definition: matrix_multi.h:148
framed_multi< Scalar_T, LO, HI, Tune_P > framed_multi_t
Definition: matrix_multi.h:149
Coefficients of denominator polynomials of Pade approximations produced by Pade1(log(1+x),x,n,n)
virtual auto operator*=(const Scalar_T &scr) -> multivector_t &=0
Product of multivector and scalar.
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.
std::array< Scalar_T, 14 > array
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.
virtual auto quad() const -> Scalar_T=0
Scalar_T quadratic form == (rev(x)*x)(0)
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 unfold(const index_set_t frm) const -> multivector_t
Subalgebra isomorphism: unfold each term within the given frame.
auto trace(const Matrix_T &val) -> typename Matrix_T::value_type
Matrix trace.
Definition: matrix_imp.h:416
static auto fast(const Matrix_T &X, index_t level) -> Multivector_T
Inverse generalized Fast Fourier Transform.
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 exp(const framed_multi< Scalar_T, LO, HI, Tune_P > &val) -> const framed_multi< Scalar_T, LO, HI, Tune_P >
Exponential of multivector.
static auto pade_log(const matrix_multi< Scalar_T, LO, HI, Tune_P > &val) -> const matrix_multi< Scalar_T, LO, HI, Tune_P >
Pade&#39; approximation of log.
static const array denom
virtual auto outer_pow(int m) const -> const multivector_t=0
Outer product power.
std::array< dd_real, 22 > array
A framed_multi<Scalar_T,LO,HI,Tune_P> is a framed approximation to a multivector. ...
Definition: framed_multi.h:56
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.
virtual auto operator%=(const multivector_t &rhs) -> multivector_t &=0
Contraction.
static auto NaN() -> Scalar_T
Smart NaN.
Definition: scalar.h:115
float pi
Definition: PyClical.pyx:1907
auto operator=(const basis_table &) -> basis_table &=delete
double scalar_t
Definition: PyClical.h:147
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
std::array< long double, 18 > array
virtual auto norm() const -> Scalar_T=0
Scalar_T norm == sum of norm of coordinates.
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.
index_set_t m_frame
Index set representing the frame for the subalgebra which contains the multivector.
Definition: matrix_multi.h:278
auto centre_pm4_qp4(index_t &p, index_t &q) -> multivector_t &
Subalgebra isomorphism: R_{p,q} to R_{p-4,q+4}.
auto norm(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> Scalar_T
Scalar_T norm == sum of norm of coordinates.
auto scalar(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> Scalar_T
Scalar part.
virtual auto operator[](const index_set_t ist) const -> Scalar_T=0
Subscripting: map from index set to scalar coordinate.
Extra traits which extend numeric limits.
Definition: scalar.h:47
Table of basis elements used as a cache by basis_element()
std::array< float, 10 > array
virtual auto inv() const -> const multivector_t=0
Geometric multiplicative inverse.
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.
virtual auto conj() const -> const multivector_t=0
Conjugation, reverse o involute == involute o reverse.
virtual auto even() const -> const multivector_t=0
Even part of multivector, sum of even grade terms.
std::pair< const index_set_t, Scalar_T > term_t
Definition: matrix_multi.h:146
virtual auto operator-=(const multivector_t &rhs) -> multivector_t &=0
Geometric difference.
std::array< long double, 18 > array
auto pos_mod(LHS_T lhs, RHS_T rhs) -> LHS_T
Modulo function which works reliably for lhs < 0.
Definition: global.h:117
static const array numer
index_set< LO, HI > index_set_t
Definition: matrix_multi.h:145
std::array< dd_real, 22 > array
auto classify_eigenvalues(const Matrix_T &val) -> eig_genus< Matrix_T >
Classify the eigenvalues of a matrix.
Definition: matrix_imp.h:548
auto clifford_exp(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Exponential of multivector.
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 isnan(const Matrix_T &m) -> bool
Not a Number.
Definition: matrix_imp.h:292
const scalar_t epsilon
Definition: PyClical.h:150
std::array< float, 10 > array
Coefficients of numerator polynomials of Pade approximations produced by Pade1(sqrt(1+x),x,n,n)
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.
A matrix_multi<Scalar_T,LO,HI,Tune_P> is a matrix approximation to a multivector. ...
Definition: framed_multi.h:59
static auto generator() -> generator_table< Matrix_T > &
Single instance of generator table.
virtual void write(const std::string &msg="") const=0
Write formatted multivector to output.
virtual auto operator^=(const multivector_t &rhs) -> multivector_t &=0
Outer product.
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&#39; approximation.
ublas::matrix< Scalar_T, orientation_t > matrix_t
Definition: matrix_multi.h:158
auto operator>>(std::istream &s, framed_multi< Scalar_T, LO, HI, Tune_P > &val) -> std::istream &
Read multivector from input.
auto basis_element(const index_set< LO, HI > &ist) const -> const basis_matrix_t
Create a basis element matrix within the current frame.
virtual auto reverse() const -> const multivector_t=0
Reversion, eg. {1}*{2} -> {2}*{1}.
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&#39; approximation of log.
virtual auto operator()(index_t grade) const -> const multivector_t=0
Pure grade-vector part.
virtual auto frame() const -> const index_set_t=0
Subalgebra generated by all generators of terms of given multivector.
static auto folded_dim(const index_set< LO, HI > &sub) -> Matrix_Index_T
Determine the matrix dimension of the fold of a subalegbra.
std::array< Scalar_T, 14 > array
virtual auto vector_part() const -> const vector_t=0
Vector part of multivector, as a vector_t with respect to frame()
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.
~basis_table()=default
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.
Abstract exception class.
Definition: errors.h:41
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.
Coefficients of denominator polynomials of Pade approximations produced by Pade1(sqrt(1+x),x,n,n)
virtual auto truncated(const Scalar_T &limit=default_truncation) const -> const multivector_t=0
Remove all terms with relative size smaller than limit.
auto centre_qp1_pm1(index_t &p, index_t &q) -> multivector_t &
Subalgebra isomorphism: R_{p,q} to R_{q+1,p-1}.
Coefficients of numerator polynomials of Pade approximations produced by Pade1(log(1+x),x,n,n)
static const array denom
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.
std::array< qd_real, 34 > array
virtual auto odd() const -> const multivector_t=0
Odd part of multivector, sum of odd grade terms.
auto fast_matrix_multi(const index_set_t frm) const -> const matrix_multi_t
Use generalized FFT to construct a matrix_multi_t.
std::array< long double, 18 > array
std::vector< Scalar_T > vector_t
Definition: matrix_multi.h:147
auto centre_pp4_qm4(index_t &p, index_t &q) -> multivector_t &
Subalgebra isomorphism: R_{p,q} to R_{p+4,q-4}.
typename matrix_t::size_type matrix_index_t
Definition: matrix_multi.h:159
static auto classname() -> const std::string
Class name used in messages.
friend class friend_for_private_destructor
std::array< dd_real, 22 > array
std::array< long double, 18 > array
std::array< Scalar_T, 14 > array
int index_t
Size of index_t should be enough to represent LO, HI.
Definition: global.h:77
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 max() const -> index_t
Maximum member.
std::array< qd_real, 34 > array
virtual auto grade() const -> index_t=0
Maximum of the grades of each term.
auto sqrt(const Multivector< Scalar_T, LO, HI, Tune_P > &val, const Multivector< Scalar_T, LO, HI, Tune_P > &i, const bool prechecked=false) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Square root of multivector with specified complexifier.
auto norm_frob2(const Matrix_T &val) -> typename Matrix_T::value_type
Square of Frobenius norm.
Definition: matrix_imp.h:395
Index set class based on std::bitset<> in Gnu standard C++ library.
Definition: index_set.h:45
virtual auto isinf() const -> bool=0
Check if a multivector contains any infinite values.
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.
virtual auto operator==(const multivector_t &val) const -> bool=0
Test for equality of multivectors.
auto operator+=(const term_t &rhs) -> multivector_t &
Add a term, if non-zero.
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.
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()
virtual auto operator/=(const Scalar_T &scr) -> multivector_t &=0
Quotient of multivector and scalar.
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.
static auto random(const index_set_t frm, Scalar_T fill=Scalar_T(1)) -> const matrix_multi_t
Random multivector within a frame.
ublas::compressed_matrix< int, orientation_t > basis_matrix_t
Definition: matrix_multi.h:157
matrix_multi()
Default constructor.
virtual auto max_abs() const -> Scalar_T=0
Maximum of absolute values of components of multivector: multivector infinity norm.
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.
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.
std::array< dd_real, 22 > array
static auto basis() -> basis_table &
Single instance of basis table.
matrix_t m_matrix
Matrix value representing the multivector within the folded frame.
Definition: matrix_multi.h:280
virtual auto operator|=(const multivector_t &rhs) -> multivector_t &=0
Transformation via twisted adjoint action.
def e(obj)
Definition: PyClical.pyx:1936
virtual auto involute() const -> const multivector_t=0
Main involution, each {i} is replaced by -{i} in each term, eg. {1} -> -{1}.
static auto random(const index_set_t frm, Scalar_T fill=Scalar_T(1)) -> const multivector_t
Random multivector within a frame.
static const array numer
virtual auto isnan() const -> bool=0
Check if a multivector contains any IEEE NaN values.
auto min() const -> index_t
Minimum member.
std::array< qd_real, 34 > array
std::array< qd_real, 34 > array
auto offset_level(const index_t p, const index_t q) -> index_t
Determine the log2 dim corresponding to signature p, q.
_GLUCAT_CLIFFORD_ALGEBRA_OPERATIONS auto operator=(const multivector_t &rhs) -> multivector_t &
Assignment operator.
auto isinf(const Matrix_T &m) -> bool
Infinite.
Definition: matrix_imp.h:275
matrix_multi multivector_t
Definition: matrix_multi.h:141
std::array< float, 10 > array
virtual auto operator&=(const multivector_t &rhs) -> multivector_t &=0
Inner product.
auto count() const -> index_t
Cardinality: Number of indices included in set.
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.
std::array< Scalar_T, 14 > array