glucat  0.8.2
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-2016 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/matrix.h"
37 #include "glucat/generation.h"
38 
39 # if defined(_GLUCAT_GCC_IGNORE_UNUSED_LOCAL_TYPEDEFS)
40 # pragma GCC diagnostic push
41 # pragma GCC diagnostic ignored "-Wunused-local-typedefs"
42 # endif
43 #include <boost/numeric/ublas/matrix.hpp>
44 #include <boost/numeric/ublas/matrix_expression.hpp>
45 #include <boost/numeric/ublas/matrix_proxy.hpp>
46 #include <boost/numeric/ublas/matrix_sparse.hpp>
47 #include <boost/numeric/ublas/operation.hpp>
48 #include <boost/numeric/ublas/operation_sparse.hpp>
49 #include <boost/numeric/ublas/triangular.hpp>
50 #include <boost/numeric/ublas/lu.hpp>
51 # if defined(_GLUCAT_GCC_IGNORE_UNUSED_LOCAL_TYPEDEFS)
52 # pragma GCC diagnostic pop
53 # endif
54 
55 #include <fstream>
56 #include <iomanip>
57 
58 namespace glucat
59 {
60  // References for algorithms:
61  // [M]: Scott Meyers, "Effective C++" Second Edition, Addison-Wesley, 1998.
62  // [P]: Ian R. Porteous, "Clifford algebras and the classical groups", Cambridge UP, 1995.
63  // [L]: Pertti Lounesto, "Clifford algebras and spinors", Cambridge UP, 1997.
64 
66  template< typename Scalar_T, const index_t LO, const index_t HI >
67  const std::string
70  { return "matrix_multi"; }
71 
73  // Reference: [P] Table 15.27, p 133
74  inline
75  index_t
76  offset_level(const index_t p, const index_t q)
77  {
78  // Offsets between the log2 of the matrix dimension for the current signature
79  // and that of the real superalgebra
80  static const int offset_log2_dim[] = {0, 1, 0, 1, 1, 2, 1, 1};
81  const index_t bott = pos_mod(p-q, 8);
82  return (p+q)/2 + offset_log2_dim[bott];
83  }
84 
86  // Reference: [P] Table 15.27, p 133
87  template< typename Matrix_Index_T, const index_t LO, const index_t HI >
88  inline
89  static
90  Matrix_Index_T
92  { return 1 << offset_level(sub.count_pos(), sub.count_neg()); }
93 
95  template< typename Scalar_T, const index_t LO, const index_t HI >
98  : m_frame( index_set_t() ),
99  m_matrix( matrix_t( 1, 1 ) )
100  { this->m_matrix.clear(); }
101 
103  template< typename Scalar_T, const index_t LO, const index_t HI >
104  template< typename Other_Scalar_T >
107  : m_frame( val.m_frame ), m_matrix( val.m_matrix.size1(), val.m_matrix.size2() )
108  {
109  this->m_matrix.clear();
110  typedef typename matrix_multi<Other_Scalar_T,LO,HI>::matrix_t other_matrix_t;
111  typedef typename other_matrix_t::const_iterator1 other_const_iterator1;
112  typedef typename other_matrix_t::const_iterator2 other_const_iterator2;
113  for (other_const_iterator1
114  val_it1 = val.m_matrix.begin1();
115  val_it1 != val.m_matrix.end1();
116  ++val_it1)
117  for (other_const_iterator2
118  val_it2 = val_it1.begin();
119  val_it2 != val_it1.end();
120  ++val_it2)
121  this->m_matrix(val_it2.index1(), val_it2.index2()) = numeric_traits<Scalar_T>::to_scalar_t(*val_it2);
122  }
123 
125  template< typename Scalar_T, const index_t LO, const index_t HI >
126  template< typename Other_Scalar_T >
128  matrix_multi(const matrix_multi<Other_Scalar_T,LO,HI>& val, const index_set_t frm, const bool prechecked)
129  : m_frame( frm )
130  {
131  if (frm != val.m_frame)
132  *this = multivector_t(framed_multi_t(val), frm);
133  else
134  {
135  const matrix_index_t dim = folded_dim<matrix_index_t>(frm);
136  this->m_matrix.resize(dim, dim, false);
137  this->m_matrix.clear();
138  typedef typename matrix_multi<Other_Scalar_T,LO,HI>::matrix_t other_matrix_t;
139  typedef typename other_matrix_t::const_iterator1 other_const_iterator1;
140  typedef typename other_matrix_t::const_iterator2 other_const_iterator2;
141  for (other_const_iterator1
142  val_it1 = val.m_matrix.begin1();
143  val_it1 != val.m_matrix.end1();
144  ++val_it1)
145  for (other_const_iterator2
146  val_it2 = val_it1.begin();
147  val_it2 != val_it1.end();
148  ++val_it2)
149  this->m_matrix(val_it2.index1(), val_it2.index2()) = numeric_traits<Scalar_T>::to_scalar_t(*val_it2);
150  }
151  }
152 
154  template< typename Scalar_T, const index_t LO, const index_t HI >
156  matrix_multi(const multivector_t& val, const index_set_t frm, const bool prechecked)
157  : m_frame( frm )
158  {
159  if (frm != val.m_frame)
160  *this = multivector_t(framed_multi_t(val), frm);
161  else
162  this->m_matrix = val.m_matrix;
163  }
164 
166  template< typename Scalar_T, const index_t LO, const index_t HI >
168  matrix_multi(const index_set_t ist, const Scalar_T& crd)
169  : m_frame( ist )
170  {
171  const matrix_index_t dim = folded_dim<matrix_index_t>(this->m_frame);
172  this->m_matrix.resize(dim, dim, false);
173  this->m_matrix.clear();
174  *this += term_t(ist, crd);
175  }
176 
178  template< typename Scalar_T, const index_t LO, const index_t HI >
180  matrix_multi(const index_set_t ist, const Scalar_T& crd, const index_set_t frm, const bool prechecked)
181  : m_frame( frm )
182  {
183  if (!prechecked && (ist | frm) != frm)
184  throw error_t("multivector_t(ist,crd,frm): cannot initialize with value outside of frame");
185  const matrix_index_t dim = folded_dim<matrix_index_t>(frm);
186  this->m_matrix.resize(dim, dim, false);
187  this->m_matrix.clear();
188  *this += term_t(ist, crd);
189  }
190 
192  template< typename Scalar_T, const index_t LO, const index_t HI >
194  matrix_multi(const Scalar_T& scr, const index_set_t frm)
195  : m_frame( frm )
196  {
197  const matrix_index_t dim = folded_dim<matrix_index_t>(frm);
198  this->m_matrix.resize(dim, dim, false);
199  this->m_matrix.clear();
200  *this += term_t(index_set_t(), scr);
201  }
202 
204  template< typename Scalar_T, const index_t LO, const index_t HI >
206  matrix_multi(const int scr, const index_set_t frm)
207  { *this = multivector_t(Scalar_T(scr), frm); }
208 
210  template< typename Scalar_T, const index_t LO, const index_t HI >
213  const index_set_t frm, const bool prechecked)
214  : m_frame( frm )
215  {
216  if (!prechecked && index_t(vec.size()) != frm.count())
217  throw error_t("multivector_t(vec,frm): cannot initialize with vector not matching frame");
218  const matrix_index_t dim = folded_dim<matrix_index_t>(frm);
219  this->m_matrix.resize(dim, dim, false);
220  this->m_matrix.clear();
221  typename vector_t::const_iterator vec_it = vec.begin();
222  const index_t begin_index = frm.min();
223  const index_t end_index = frm.max()+1;
224 
225  for (index_t
226  idx = begin_index;
227  idx != end_index;
228  ++idx)
229  if (frm[idx])
230  {
231  *this += term_t(index_set_t(idx), *vec_it);
232  ++vec_it;
233  }
234  }
235 
237  template< typename Scalar_T, const index_t LO, const index_t HI >
239  matrix_multi(const std::string& str)
240  { *this = framed_multi_t(str); }
241 
243  template< typename Scalar_T, const index_t LO, const index_t HI >
245  matrix_multi(const std::string& str, const index_set_t frm, const bool prechecked)
246  { *this = multivector_t(framed_multi_t(str), frm, prechecked); }
247 
249  template< typename Scalar_T, const index_t LO, const index_t HI >
250  template< typename Other_Scalar_T >
253  : m_frame( val.frame() )
254  {
255  if (val.size() >= Tune_P::fast_size_threshold)
256  try
257  {
258  *this = val.template fast_matrix_multi<Scalar_T>(this->m_frame);
259  return;
260  }
261  catch (const glucat_error& e)
262  { }
263  const matrix_index_t dim = folded_dim<matrix_index_t>(this->m_frame);
264  this->m_matrix.resize(dim, dim, false);
265  this->m_matrix.clear();
266 
268  for (typename framed_multi_t::const_iterator
269  val_it = val.begin();
270  val_it != val.end();
271  ++val_it)
272  *this += *val_it;
273  }
274 
276  template< typename Scalar_T, const index_t LO, const index_t HI >
277  template< typename Other_Scalar_T >
279  matrix_multi(const framed_multi<Other_Scalar_T,LO,HI>& val, const index_set_t frm, const bool prechecked)
280  {
281  const index_set_t 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 matrix_index_t dim = folded_dim<matrix_index_t>(our_frame);
292  this->m_matrix.resize(dim, dim, false);
293  this->m_matrix.clear();
294 
296  for (typename framed_multi_t::const_iterator
297  val_it = val.begin();
298  val_it != val.end();
299  ++val_it)
300  *this += *val_it;
301  }
302 
304  template< typename Scalar_T, const index_t LO, const index_t HI >
305  template< typename Matrix_T >
307  matrix_multi(const Matrix_T& mtx, const index_set_t frm)
308  : m_frame( frm ), m_matrix( mtx.size1(), mtx.size2() )
309  {
310  this->m_matrix.clear();
311 
312  typedef typename Matrix_T::const_iterator1 const_iterator1;
313  typedef typename Matrix_T::const_iterator2 const_iterator2;
314  for (const_iterator1
315  mtx_it1 = mtx.begin1();
316  mtx_it1 != mtx.end1();
317  ++mtx_it1)
318  for (const_iterator2
319  mtx_it2 = mtx_it1.begin();
320  mtx_it2 != mtx_it1.end();
321  ++mtx_it2)
322  this->m_matrix(mtx_it2.index1(), mtx_it2.index2()) = numeric_traits<Scalar_T>::to_scalar_t(*mtx_it2);
323  }
324 
326  template< typename Scalar_T, const index_t LO, const index_t HI >
328  matrix_multi(const matrix_t& mtx, const index_set_t frm)
329  : m_frame( frm ), m_matrix( mtx )
330  { }
331 
333  template< typename Scalar_T, const index_t LO, const index_t HI >
337  {
338  // Check for assignment to self
339  if (this == &rhs)
340  return *this;
341  this->m_frame = rhs.m_frame;
342  this->m_matrix = rhs.m_matrix;
343  return *this;
344  }
345 
347  template< typename Scalar_T, const index_t LO, const index_t HI >
348  inline
349  const index_set<LO,HI>
351  matrix_multi<Scalar_T,LO,HI>& lhs_reframed, matrix_multi<Scalar_T,LO,HI>& rhs_reframed)
352  {
356  // Determine the initial common frame
357  index_set_t our_frame = lhs.m_frame | rhs.m_frame;
358  framed_multi_t framed_lhs;
359  framed_multi_t framed_rhs;
360  if ((lhs.m_frame != our_frame) || (rhs.m_frame != our_frame))
361  {
362  // The common frame may expand as a result of the transform to framed_multi_t
363  framed_lhs = framed_multi_t(lhs);
364  framed_rhs = framed_multi_t(rhs);
365  our_frame |= framed_lhs.frame() | framed_rhs.frame();
366  }
367  // Do the reframing only where necessary
368  if (lhs.m_frame != our_frame)
369  lhs_reframed = multivector_t(framed_lhs, our_frame, true);
370  if (rhs.m_frame != our_frame)
371  rhs_reframed = multivector_t(framed_rhs, our_frame, true);
372  return our_frame;
373  }
374 
376  template< typename Scalar_T, const index_t LO, const index_t HI >
377  bool
379  operator== (const multivector_t& rhs) const
380  {
381  // Ensure that there is no aliasing
382  if (this == &rhs)
383  return true;
384 
385  // Operate only within a common frame
386  multivector_t lhs_reframed;
387  multivector_t rhs_reframed;
388  const index_set_t our_frame = reframe(*this, rhs, lhs_reframed, rhs_reframed);
389  const multivector_t& lhs_ref = (this->m_frame == our_frame)
390  ? *this
391  : lhs_reframed;
392  const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
393  ? rhs
394  : rhs_reframed;
395 
396 #if defined(_GLUCAT_USE_DENSE_MATRICES)
397  return ublas::norm_inf(lhs_ref.m_matrix - rhs_ref.m_matrix) == 0;
398 #else
399  typedef typename matrix_t::const_iterator1 const_iterator1;
400  typedef typename matrix_t::const_iterator2 const_iterator2;
401  // If either matrix contains zero entries,
402  // compare using subtraction and ublas::norm_inf
403  for (const_iterator1
404  it1 = lhs_ref.m_matrix.begin1();
405  it1 != lhs_ref.m_matrix.end1();
406  ++it1)
407  for (const_iterator2
408  it2 = it1.begin();
409  it2 != it1.end();
410  ++it2)
411  if (*it2 == 0)
412  return ublas::norm_inf(lhs_ref.m_matrix - rhs_ref.m_matrix) == 0;
413  for (const_iterator1
414  it1 = rhs_ref.m_matrix.begin1();
415  it1 != rhs_ref.m_matrix.end1();
416  ++it1)
417  for (const_iterator2
418  it2 = it1.begin();
419  it2 != it1.end();
420  ++it2)
421  if (*it2 == 0)
422  return ublas::norm_inf(lhs_ref.m_matrix - rhs_ref.m_matrix) == 0;
423  // Neither matrix contains zero entries.
424  // Compare by iterating over both matrices in lock step.
425  const_iterator1 this_it1 = lhs_ref.m_matrix.begin1();
426  const_iterator1 it1 = rhs_ref.m_matrix.begin1();
427  for (;
428  (this_it1 != lhs_ref.m_matrix.end1()) && (it1 != rhs_ref.m_matrix.end1());
429  ++this_it1, ++it1)
430  {
431  if ( this_it1.index1() != it1.index1() )
432  return false;
433  const_iterator2 this_it2 = this_it1.begin();
434  const_iterator2 it2 = it1.begin();
435  for (;
436  (this_it2 != this_it1.end()) && (it2 != it1.end());
437  ++this_it2, ++it2)
438  if ( (this_it2.index2() != it2.index2()) || (*this_it2 != *it2) )
439  return false;
440  if ( (this_it2 != this_it1.end()) || (it2 != it1.end()) )
441  return false;
442  }
443  return (this_it1 == lhs_ref.m_matrix.end1()) && (it1 == rhs_ref.m_matrix.end1());
444 #endif
445  }
446 
447  // Test for equality of multivector and scalar
448  template< typename Scalar_T, const index_t LO, const index_t HI >
449  inline
450  bool
452  operator== (const Scalar_T& scr) const
453  {
454  if (scr != Scalar_T(0))
455  return *this == multivector_t(framed_multi_t(scr), this->m_frame, true);
456  else if (ublas::norm_inf(this->m_matrix) != 0)
457  return false;
458  else
459  {
460  const matrix_index_t dim = this->m_matrix.size1();
461  return !(dim == 1 && this->isnan());
462  }
463  }
464 
466  template< typename Scalar_T, const index_t LO, const index_t HI >
467  inline
470  operator+= (const Scalar_T& scr)
471  { return *this += term_t(index_set_t(), scr); }
472 
474  template< typename Scalar_T, const index_t LO, const index_t HI >
475  inline
478  operator+= (const multivector_t& rhs)
479  {
480  // Ensure that there is no aliasing
481  if (this == &rhs)
482  return *this *= Scalar_T(2);
483 
484  // Operate only within a common frame
485  multivector_t rhs_reframed;
486  const index_set_t our_frame = reframe(*this, rhs, *this, rhs_reframed);
487  const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
488  ? rhs
489  : rhs_reframed;
490 
491  noalias(this->m_matrix) += rhs_ref.m_matrix;
492  return *this;
493  }
494 
496  template< typename Scalar_T, const index_t LO, const index_t HI >
497  inline
500  operator-= (const multivector_t& rhs)
501  {
502  // Ensure that there is no aliasing
503  if (this == &rhs)
504  return *this = Scalar_T(0);
505 
506  // Operate only within a common frame
507  multivector_t rhs_reframed;
508  const index_set_t our_frame = reframe(*this, rhs, *this, rhs_reframed);
509  const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
510  ? rhs
511  : rhs_reframed;
512 
513  noalias(this->m_matrix) -= rhs_ref.m_matrix;
514  return *this;
515  }
516 
518  template< typename Scalar_T, const index_t LO, const index_t HI >
519  inline
522  operator- () const
523  { return multivector_t(-(this->m_matrix), this->m_frame); }
524 
526  template< typename Scalar_T, const index_t LO, const index_t HI >
527  inline
530  operator*= (const Scalar_T& scr)
531  { // multiply coordinates of all terms by scalar
532 
533  typedef numeric_traits<Scalar_T> traits_t;
534  if (traits_t::isNaN_or_isInf(scr) || this->isnan())
535  return *this = traits_t::NaN();
536  if (scr == Scalar_T(0))
537  *this = Scalar_T(0);
538  else
539  this->m_matrix *= scr;
540  return *this;
541  }
542 
544  template< typename Scalar_T, const index_t LO, const index_t HI >
545  inline
548  {
550  typedef typename multivector_t::index_set_t index_set_t;
551 
552 #if defined(_GLUCAT_CHECK_ISNAN)
553  if (lhs.isnan() || rhs.isnan())
555 #endif
556 
557  // Operate only within a common frame
558  multivector_t lhs_reframed;
559  multivector_t rhs_reframed;
560  const index_set_t our_frame = reframe(lhs, rhs, lhs_reframed, rhs_reframed);
561  const multivector_t& lhs_ref = (lhs.m_frame == our_frame)
562  ? lhs
563  : lhs_reframed;
564  const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
565  ? rhs
566  : rhs_reframed;
567 
568  typedef typename multivector_t::matrix_t matrix_t;
569 #if defined(_GLUCAT_USE_DENSE_MATRICES)
570  typedef typename matrix_t::size_type matrix_index_t;
571 
572  const matrix_index_t dim = lhs_ref.m_matrix.size1();
573  multivector_t result = multivector_t(matrix_t(dim, dim), our_frame);
574  result.m_matrix.clear();
575  ublas::axpy_prod(lhs_ref.m_matrix, rhs_ref.m_matrix, result.m_matrix, true);
576  return result;
577 #else
578  typedef typename matrix_t::expression_type expression_t;
579 
580  return
581  multivector_t(ublas::sparse_prod<expression_t>(lhs_ref.m_matrix, rhs_ref.m_matrix), our_frame);
582 #endif
583  }
584 
586  template< typename Scalar_T, const index_t LO, const index_t HI >
587  inline
590  operator*= (const multivector_t& rhs)
591  { return *this = *this * rhs; }
592 
594  template< typename Scalar_T, const index_t LO, const index_t HI >
595  inline
598  {
601  return framed_multi_t(lhs) ^ framed_multi_t(rhs);
602  }
603 
605  template< typename Scalar_T, const index_t LO, const index_t HI >
606  inline
609  operator^= (const multivector_t& rhs)
610  { return *this = *this ^ rhs; }
611 
613  template< typename Scalar_T, const index_t LO, const index_t HI >
614  inline
617  {
620  return framed_multi_t(lhs) & framed_multi_t(rhs);
621  }
622 
624  template< typename Scalar_T, const index_t LO, const index_t HI >
625  inline
628  operator&= (const multivector_t& rhs)
629  { return *this = *this & rhs; }
630 
632  template< typename Scalar_T, const index_t LO, const index_t HI >
633  inline
636  {
639  return framed_multi_t(lhs) % framed_multi_t(rhs);
640  }
641 
643  template< typename Scalar_T, const index_t LO, const index_t HI >
644  inline
647  operator%= (const multivector_t& rhs)
648  { return *this = *this % rhs; }
649 
651  template< typename Scalar_T, const index_t LO, const index_t HI >
652  inline
653  Scalar_T
655  { return (lhs * rhs).scalar(); }
656 
658  template< typename Scalar_T, const index_t LO, const index_t HI >
659  inline
662  operator/= (const Scalar_T& scr)
663  { return *this *= Scalar_T(1)/scr; }
664 
666  template< typename Scalar_T, const index_t LO, const index_t HI >
669  {
670  typedef numeric_traits<Scalar_T> traits_t;
671 
672 #if defined(_GLUCAT_CHECK_ISNAN)
673  if (lhs.isnan() || rhs.isnan())
674  return traits_t::NaN();
675 #endif
676 
677  if (rhs == Scalar_T(0))
678  return traits_t::NaN();
679 
681  typedef typename multivector_t::index_set_t index_set_t;
682 
683  // Operate only within a common frame
684  multivector_t lhs_reframed;
685  multivector_t rhs_reframed;
686  const index_set_t our_frame = reframe(lhs, rhs, lhs_reframed, rhs_reframed);
687  const multivector_t& lhs_ref = (lhs.m_frame == our_frame)
688  ? lhs
689  : lhs_reframed;
690  const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
691  ? rhs
692  : rhs_reframed;
693 
694  // Solve result == lhs_ref/rhs_ref <=> result*rhs_ref == lhs_ref
695  // We now solve X == B/A
696  // (where X == result, B == lhs_ref.m_matrix and A == rhs_ref.m_matrix)
697  // X == B/A <=> X*A == B <=> AT*XT == BT
698  // So, we solve AT*XT == BT
699 
700  typedef typename multivector_t::matrix_t matrix_t;
701  typedef typename matrix_t::size_type matrix_index_t;
702 
703  const matrix_t& AT = ublas::trans(rhs_ref.m_matrix);
704  matrix_t LU = AT;
705 
706  typedef ublas::permutation_matrix<matrix_index_t> permutation_t;
707 
708  permutation_t pvector(AT.size1());
709  if (! ublas::lu_factorize(LU, pvector))
710  {
711  const matrix_t& BT = ublas::trans(lhs_ref.m_matrix);
712  matrix_t XT = BT;
713  ublas::lu_substitute(LU, pvector, XT);
714 #if defined(_GLUCAT_CHECK_ISNAN)
715  if (matrix::isnan(XT))
716  return traits_t::NaN();
717 #endif
718 
719  // Iterative refinement.
720  // Reference: Nicholas J. Higham, "Accuracy and Stability of Numerical Algorithms",
721  // SIAM, 1996, ISBN 0-89871-355-2, Chapter 11
722  if (Tune_P::div_max_steps > 0)
723  {
724  // matrix_t R = ublas::prod(AT, XT) - BT;
725  matrix_t R = -BT;
726  ublas::axpy_prod(AT, XT, R, false);
727 #if defined(_GLUCAT_CHECK_ISNAN)
728  if (matrix::isnan(R))
729  return traits_t::NaN();
730 #endif
731 
732  Scalar_T nr = ublas::norm_inf(R);
733  if ( nr != Scalar_T(0) && !traits_t::isNaN_or_isInf(nr) )
734  {
735  matrix_t XTnew = XT;
736  Scalar_T nrold = nr + Scalar_T(1);
737  for (int
738  step = 0;
739  step != Tune_P::div_max_steps &&
740  nr < nrold &&
741  nr != Scalar_T(0) &&
742  nr == nr;
743  ++step)
744  {
745  nrold = nr;
746  if (step != 0)
747  XT = XTnew;
748  matrix_t& D = R;
749  ublas::lu_substitute(LU, pvector, D);
750  XTnew -= D;
751  // noalias(R) = ublas::prod(AT, XTnew) - BT;
752  R = -BT;
753  ublas::axpy_prod(AT, XTnew, R, false);
754  nr = ublas::norm_inf(R);
755  }
756  }
757  }
758  return multivector_t(ublas::trans(XT), our_frame);
759  }
760  else
761  // AT is singular. Return NaN
762  return traits_t::NaN();
763  }
764 
766  template< typename Scalar_T, const index_t LO, const index_t HI >
767  inline
770  operator/= (const multivector_t& rhs)
771  { return *this = *this / rhs; }
772 
774  template< typename Scalar_T, const index_t LO, const index_t HI >
775  inline
778  { return rhs * lhs / rhs.involute(); }
779 
781  template< typename Scalar_T, const index_t LO, const index_t HI >
782  inline
785  operator|= (const multivector_t& rhs)
786  { return *this = rhs * *this / rhs.involute(); }
787 
789  template< typename Scalar_T, const index_t LO, const index_t HI >
790  inline
793  inv() const
794  { return multivector_t(Scalar_T(1), this->m_frame) / *this; }
795 
797  template< typename Scalar_T, const index_t LO, const index_t HI >
798  inline
801  pow(int m) const
802  { return glucat::pow(*this, m); }
803 
805  template< typename Scalar_T, const index_t LO, const index_t HI >
808  outer_pow(int m) const
809  {
810  if (m < 0)
811  throw error_t("outer_pow(m): negative exponent");
812  framed_multi_t result = Scalar_T(1);
813  framed_multi_t a = *this;
814  for (;
815  m != 0;
816  m >>= 1)
817  {
818  if (m & 1)
819  result ^= a;
820  a ^= a;
821  }
822  return result;
823  }
824 
826  template< typename Scalar_T, const index_t LO, const index_t HI >
827  inline
828  index_t
830  grade() const
831  { return framed_multi_t(*this).grade(); }
832 
834  template< typename Scalar_T, const index_t LO, const index_t HI >
835  inline
836  const index_set<LO,HI>
838  frame() const
839  { return this->m_frame; }
840 
842  template< typename Scalar_T, const index_t LO, const index_t HI >
843  inline
844  Scalar_T
846  operator[] (const index_set_t ist) const
847  {
848  // Use matrix inner product only if ist is in frame
849  if ( (ist | this->m_frame) == this->m_frame)
850  return matrix::inner<Scalar_T>(this->basis_element(ist), this->m_matrix);
851  else
852  return Scalar_T(0);
853  }
854 
856  template< typename Scalar_T, const index_t LO, const index_t HI >
857  inline
861  {
862  if ((grade < 0) || (grade > HI-LO))
863  return 0;
864  else
865  return (framed_multi_t(*this))(grade);
866  }
867 
869  template< typename Scalar_T, const index_t LO, const index_t HI >
870  inline
871  Scalar_T
873  scalar() const
874  {
875  const matrix_index_t dim = this->m_matrix.size1();
876  return matrix::trace(this->m_matrix) / Scalar_T( double(dim) );
877  }
878 
880  template< typename Scalar_T, const index_t LO, const index_t HI >
881  inline
884  pure() const
885  { return *this - this->scalar(); }
886 
888  template< typename Scalar_T, const index_t LO, const index_t HI >
889  inline
892  even() const
893  { return framed_multi_t(*this).even(); }
894 
896  template< typename Scalar_T, const index_t LO, const index_t HI >
897  inline
900  odd() const
901  { return framed_multi_t(*this).odd(); }
902 
904  template< typename Scalar_T, const index_t LO, const index_t HI >
907  vector_part() const
908  { return this->vector_part(this->frame(), true); }
909 
911  template< typename Scalar_T, const index_t LO, const index_t HI >
914  vector_part(const index_set_t frm, const bool prechecked) const
915  {
916  if (!prechecked && (this->frame() | frm) != frm)
917  throw error_t("vector_part(frm): value is outside of requested frame");
918  vector_t result;
919  // If we need to enlarge the frame we may as well use a framed_multi_t
920  if (this->frame() != frm)
921  return framed_multi_t(*this).vector_part(frm, true);
922 
923  const index_t begin_index = frm.min();
924  const index_t end_index = frm.max()+1;
925  for (index_t
926  idx = begin_index;
927  idx != end_index;
928  ++idx)
929  if (frm[idx])
930  // Frame may contain indices which do not correspond to a grade 1 term but
931  // frame cannot omit any index corresponding to a grade 1 term
932  result.push_back(
933  matrix::inner<Scalar_T>(this->basis_element(index_set_t(idx)),
934  this->m_matrix));
935  return result;
936  }
937 
939  template< typename Scalar_T, const index_t LO, const index_t HI >
940  inline
943  involute() const
944  { return framed_multi_t(*this).involute(); }
945 
947  template< typename Scalar_T, const index_t LO, const index_t HI >
948  inline
951  reverse() const
952  { return framed_multi_t(*this).reverse(); }
953 
955  template< typename Scalar_T, const index_t LO, const index_t HI >
956  inline
959  conj() const
960  { return framed_multi_t(*this).conj(); }
961 
963  template< typename Scalar_T, const index_t LO, const index_t HI >
964  inline
965  Scalar_T
967  quad() const
968  { // scalar(conj(x)*x) = 2*quad(even(x)) - quad(x)
969  // Arvind Raja ref: "old clical: quadfunction(p:pter):pterm in file compmod.pas"
970  return framed_multi_t(*this).quad();
971  }
972 
974  template< typename Scalar_T, const index_t LO, const index_t HI >
975  inline
976  Scalar_T
978  norm() const
979  {
980  const matrix_index_t dim = this->m_matrix.size1();
981  return matrix::norm_frob2(this->m_matrix) / Scalar_T( double(dim) );
982  }
983 
985  template< typename Scalar_T, const index_t LO, const index_t HI >
986  inline
987  Scalar_T
989  max_abs() const
990  { return framed_multi_t(*this).max_abs(); }
991 
993  template< typename Scalar_T, const index_t LO, const index_t HI >
996  random(const index_set<LO,HI> frm, Scalar_T fill)
997  {
998  return framed_multi<Scalar_T,LO,HI>::random(frm, fill);
999  }
1000 
1002  template< typename Scalar_T, const index_t LO, const index_t HI >
1003  inline
1004  void
1006  write(const std::string& msg) const
1007  { framed_multi_t(*this).write(msg); }
1008 
1010  template< typename Scalar_T, const index_t LO, const index_t HI >
1011  inline
1012  void
1014  write(std::ofstream& ofile, const std::string& msg) const
1015  {
1016  if (!ofile)
1017  throw error_t("write(ofile,msg): cannot write to output file");
1018  framed_multi_t(*this).write(ofile, msg);
1019  }
1020 
1022  template< typename Scalar_T, const index_t LO, const index_t HI >
1023  inline
1024  std::ostream&
1025  operator<< (std::ostream& os, const matrix_multi<Scalar_T,LO,HI>& val)
1026  {
1027  os << typename matrix_multi<Scalar_T,LO,HI>::framed_multi_t(val);
1028  return os;
1029  }
1030 
1032  template< typename Scalar_T, const index_t LO, const index_t HI >
1033  inline
1034  std::istream&
1036  { // Input looks like 1.0-2.0{1,2}+3.2{3,4}
1038  s >> local;
1039  // If s.bad() then we have a corrupt input
1040  // otherwise we are fine and can copy the resulting matrix_multi
1041  if (!s.bad())
1042  val = local;
1043  return s;
1044  }
1045 
1047  template< typename Scalar_T, const index_t LO, const index_t HI >
1048  inline
1049  bool
1051  isnan() const
1052  {
1053  if (std::numeric_limits<Scalar_T>::has_quiet_NaN)
1054  return matrix::isnan(this->m_matrix);
1055  else
1056  return false;
1057  }
1058 
1060  template< typename Scalar_T, const index_t LO, const index_t HI >
1061  inline
1064  truncated(const Scalar_T& limit) const
1065  { return framed_multi_t(*this).truncated(limit); }
1066 
1068  template< typename Scalar_T, const index_t LO, const index_t HI >
1069  inline
1072  operator+= (const term_t& term)
1073  {
1074  if (term.second != Scalar_T(0))
1075  this->m_matrix.plus_assign(matrix_t(this->basis_element(term.first)) * term.second);
1076  return *this;
1077  }
1078 
1080  template< typename Multivector_T, typename Matrix_T, typename Basis_Matrix_T >
1081  static
1082  Multivector_T
1083  fast(const Matrix_T& X, index_t level)
1084  {
1085  typedef Multivector_T framed_multi_t;
1086 
1087  typedef typename framed_multi_t::index_set_t index_set_t;
1088  typedef typename framed_multi_t::scalar_t Scalar_T;
1089  typedef Matrix_T matrix_t;
1090  typedef Basis_Matrix_T basis_matrix_t;
1091  typedef typename basis_matrix_t::value_type basis_scalar_t;
1092  typedef numeric_traits<Scalar_T> traits_t;
1093 
1094  if (level == 0)
1095  return framed_multi_t(traits_t::to_scalar_t(X(0,0)));
1096 
1097  if (ublas::norm_inf(X) == 0)
1098  return Scalar_T(0);
1099 
1100  const basis_matrix_t& I = matrix::unit<basis_matrix_t>(2);
1101  basis_matrix_t J(2,2,2);
1102  J.clear();
1103  J(0,1) = basis_scalar_t(-1);
1104  J(1,0) = basis_scalar_t( 1);
1105  basis_matrix_t K = J;
1106  K(0,1) = basis_scalar_t( 1);
1107  basis_matrix_t JK = I;
1108  JK(0,0) = basis_scalar_t(-1);
1109 
1111  const index_set_t ist_mn = index_set_t(-level);
1112  const index_set_t ist_pn = index_set_t(level);
1113  const index_set_t ist_mnpn = ist_mn | ist_pn;
1114  if (level == 1)
1115  {
1116  typedef typename framed_multi_t::term_t term_t;
1117  const Scalar_T i_x = traits_t::to_scalar_t(signed_perm_nork(I, X)(0, 0));
1118  const Scalar_T j_x = traits_t::to_scalar_t(signed_perm_nork(J, X)(0, 0));
1119  const Scalar_T k_x = traits_t::to_scalar_t(signed_perm_nork(K, X)(0, 0));
1120  const Scalar_T jk_x = traits_t::to_scalar_t(signed_perm_nork(JK,X)(0, 0));
1121  framed_multi_t
1122  result = i_x;
1123  result += term_t(ist_mn, j_x); // j_x * mn;
1124  result += term_t(ist_pn, k_x); // k_x * pn;
1125  return result += term_t(ist_mnpn, jk_x); // jk_x * mnpn;
1126  }
1127  else
1128  {
1129  const framed_multi_t& mn = framed_multi_t(ist_mn);
1130  const framed_multi_t& pn = framed_multi_t(ist_pn);
1131  const framed_multi_t& mnpn = framed_multi_t(ist_mnpn);
1132  const framed_multi_t& i_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1133  (signed_perm_nork(I, X), level-1);
1134  const framed_multi_t& j_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1135  (signed_perm_nork(J, X), level-1);
1136  const framed_multi_t& k_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1137  (signed_perm_nork(K, X), level-1);
1138  const framed_multi_t& jk_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1139  (signed_perm_nork(JK,X), level-1);
1140  framed_multi_t
1141  result = i_x.even() - jk_x.odd();
1142  result += (j_x.even() - k_x.odd()) * mn;
1143  result += (k_x.even() - j_x.odd()) * pn;
1144  return result += (jk_x.even() - i_x.odd()) * mnpn;
1145  }
1146  }
1147 
1149  template< typename Scalar_T, const index_t LO, const index_t HI >
1150  inline
1154  {
1155  if (this->m_frame == frm)
1156  return *this;
1157  else
1158  return (this->template fast_framed_multi<Scalar_T>()).template fast_matrix_multi<Scalar_T>(frm);
1159  }
1160 
1162  template< typename Scalar_T, const index_t LO, const index_t HI >
1163  template <typename Other_Scalar_T>
1167  {
1168  // Determine the amount of off-centering needed
1169  index_t p = this->m_frame.count_pos();
1170  index_t q = this->m_frame.count_neg();
1171 
1172  const index_t bott = pos_mod(p-q, 8);
1173  p += std::max(gen::offset_to_super[bott],index_t(0));
1174  q -= std::min(gen::offset_to_super[bott],index_t(0));
1175 
1176  const index_t orig_p = p;
1177  const index_t orig_q = q;
1178  while (p-q > 4)
1179  { p -= 4; q += 4; }
1180  while (p-q < -3)
1181  { p += 4; q -= 4; }
1182  if (p-q > 1)
1183  {
1184  index_t old_p = p;
1185  p = q+1;
1186  q = old_p-1;
1187  }
1188  const index_t level = (p+q)/2;
1189 
1190  // Do the inverse fast transform
1192  framed_multi_t val = fast<framed_multi_t, matrix_t, basis_matrix_t>(this->m_matrix, level);
1193 
1194  // Off-centre val
1195  switch (pos_mod(orig_p-orig_q, 8))
1196  {
1197  case 2:
1198  case 3:
1199  case 4:
1200  val.centre_qp1_pm1(p, q);
1201  break;
1202  default:
1203  break;
1204  }
1205  if (orig_p-orig_q > 4)
1206  while (p != orig_p)
1207  val.centre_pp4_qm4(p, q);
1208  if (orig_p-orig_q < -3)
1209  while (p != orig_p)
1210  val.centre_pm4_qp4(p, q);
1211 
1212  // Return unfolded val
1213  return val.unfold(this->m_frame);
1214  }
1215 
1217  template< typename Scalar_T, const index_t LO, const index_t HI, typename Matrix_T >
1218  class basis_table :
1219  public std::map< std::pair< const index_set<LO,HI>, const index_set<LO,HI> >,
1220  Matrix_T* >
1221  {
1222  public:
1224  static basis_table& basis() { static basis_table b; return b;}
1225  private:
1226  // Enforce singleton
1227  // Reference: A. Alexandrescu, "Modern C++ Design", Chapter 6
1230  basis_table(const basis_table&);
1232 
1236  friend class friend_for_private_destructor;
1237  };
1238 
1240  template< typename Scalar_T, const index_t LO, const index_t HI >
1243  basis_element(const index_set_t& ist) const
1244  {
1245  typedef std::pair<const index_set_t, const index_set_t> index_set_pair_t;
1246  const index_set_pair_t& unfolded_pair = index_set_pair_t(ist, this->m_frame);
1247 
1248  typedef basis_table<Scalar_T,LO,HI,basis_matrix_t> basis_table_t;
1249  typedef typename basis_table_t::const_iterator basis_const_iterator_t;
1250  basis_table_t& basis_cache = basis_table_t::basis();
1251 
1252  const index_t frame_count = this->m_frame.count();
1253  const bool use_cache = frame_count <= index_t(Tune_P::basis_max_count);
1254 
1255  if (use_cache)
1256  {
1257  const basis_const_iterator_t basis_it = basis_cache.find(unfolded_pair);
1258  if (basis_it != basis_cache.end())
1259  return *(basis_it->second);
1260  }
1261  const index_set_t folded_set = ist.fold(this->m_frame);
1262  const index_set_t folded_frame = this->m_frame.fold();
1263  const index_set_pair_t& folded_pair = index_set_pair_t(folded_set, folded_frame);
1264  typedef std::pair<const index_set_pair_t, basis_matrix_t*> basis_pair_t;
1265  if (use_cache)
1266  {
1267  const basis_const_iterator_t basis_it = basis_cache.find(folded_pair);
1268  if (basis_it != basis_cache.end())
1269  {
1270  basis_matrix_t* result_ptr = basis_it->second;
1271  basis_cache.insert(basis_pair_t(unfolded_pair, result_ptr));
1272  return *result_ptr;
1273  }
1274  }
1275  const index_t folded_max = folded_frame.max();
1276  const index_t folded_min = folded_frame.min();
1277  const index_t p = std::max(folded_max, index_t(0));
1278  const index_t q = std::max(index_t(-folded_min), index_t(0));
1280  const matrix_index_t dim = 1 << offset_level(p, q);
1281  basis_matrix_t result = matrix::unit<basis_matrix_t>(dim);
1282  for (index_t
1283  k = folded_min;
1284  k <= folded_max;
1285  ++k)
1286  if (folded_set[k])
1287  result = matrix::mono_prod(result, e[k]);
1288  if (use_cache)
1289  {
1290  basis_matrix_t* result_ptr = new basis_matrix_t(result);
1291  basis_cache.insert(basis_pair_t(folded_pair, result_ptr));
1292  basis_cache.insert(basis_pair_t(unfolded_pair, result_ptr));
1293  }
1294  return result;
1295  }
1296 
1298  template< typename Scalar_T, const index_t LO, const index_t HI >
1299  inline
1300  static
1302  pade_approx(const int array_size, const Scalar_T a[], const Scalar_T b[], const matrix_multi<Scalar_T,LO,HI>& X)
1303  {
1304  // Pade' approximation
1305  // Reference: [GW], Section 4.3, pp318-322
1306  // Reference: [GL], Section 11.3, p572-576.
1307 
1309  typedef numeric_traits<Scalar_T> traits_t;
1310 
1311  if (X.isnan())
1312  return traits_t::NaN();
1313 
1314  // Array size is assumed to be even
1315  const int nbr_even_powers = array_size/2 - 1;
1316 
1317  // Create an array of even powers
1318  std::vector<multivector_t> XX(nbr_even_powers);
1319  XX[0] = X * X;
1320  XX[1] = XX[0] * XX[0];
1321  for (int
1322  k = 2;
1323  k != nbr_even_powers;
1324  ++k)
1325  XX[k] = XX[k-2] * XX[1];
1326 
1327  // Calculate numerator N and denominator D
1328  multivector_t N = a[1];
1329  for (int
1330  k = 0;
1331  k != nbr_even_powers;
1332  ++k)
1333  N += XX[k] * a[2*k + 3];
1334  N *= X;
1335  N += a[0];
1336  for (int
1337  k = 0;
1338  k != nbr_even_powers;
1339  ++k)
1340  N += XX[k] * a[2*k + 2];
1341  multivector_t D = b[1];
1342  for (int
1343  k = 0;
1344  k != nbr_even_powers;
1345  ++k)
1346  D += XX[k] * b[2*k + 3];
1347  D *= X;
1348  D += b[0];
1349  for (int
1350  k = 0;
1351  k != nbr_even_powers;
1352  ++k)
1353  D += XX[k] * b[2*k + 2];
1354  return N / D;
1355  }
1356 
1358  template< typename Scalar_T, const index_t LO, const index_t HI >
1359  inline
1360  static
1361  void
1363  {
1364  // Reference: [CHKL]
1366  const multivector_t& invM = inv(M);
1367  M = ((M + invM)/Scalar_T(2) + Scalar_T(1)) / Scalar_T(2);
1368  Y *= (invM + Scalar_T(1)) / Scalar_T(2);
1369  }
1370 
1372  template< typename Scalar_T, const index_t LO, const index_t HI >
1373  static
1376  {
1377  // Reference: [CHKL]
1379 
1380  if (val == Scalar_T(0))
1381  return val;
1382 
1383  typedef std::numeric_limits<Scalar_T> limits_t;
1384  static const Scalar_T tol = std::pow(limits_t::epsilon(), 2);
1385  static const Scalar_T tol2 = tol * tol;
1386  static const int sqrt_max_steps = Tune_P::sqrt_max_steps;
1387  multivector_t M = val;
1388  multivector_t Y = val;
1389  Scalar_T norm_M_1 = norm(M - Scalar_T(1));
1390  typedef numeric_traits<Scalar_T> traits_t;
1391 
1392  for (int step = 0;
1393  step != sqrt_max_steps && norm_M_1 > tol2;
1394  ++step)
1395  {
1396  if (Y.isnan())
1397  return traits_t::NaN();
1398  db_step(M, Y);
1399  norm_M_1 = norm(M - Scalar_T(1));
1400  }
1401  if (norm_M_1 > tol2)
1402  return traits_t::NaN();
1403  else
1404  return Y;
1405  }
1406 }
1407 
1408 namespace {
1410  // Reference: [Z], Pade1
1411  template< typename Scalar_T >
1412  struct pade_sqrt_a
1413  {
1414  static const int array_size = 14;
1415  static const Scalar_T array[array_size];
1416  };
1417 
1419  // Reference: [Z], Pade1
1420  template< typename Scalar_T >
1421  struct pade_sqrt_b
1422  {
1423  static const int array_size = 14;
1424  static const Scalar_T array[array_size];
1425  };
1426 
1427  template< typename Scalar_T >
1428  const Scalar_T pade_sqrt_a<Scalar_T>::array[pade_sqrt_a<Scalar_T>::array_size] =
1429  {
1430  1.0, 27.0/4.0, 81.0/4.0, 2277.0/64.0,
1431  10395.0/256.0, 32319.0/1024.0, 8721.0/512.0, 26163.0/4096.0,
1432  53703.0/32768.0, 36465.0/131072.0, 3861.0/131072.0, 7371.0/4194304.0,
1433  819.0/16777216.0, 27.0/67108864.0
1434  };
1435  template< typename Scalar_T >
1436  const Scalar_T pade_sqrt_b<Scalar_T>::array[pade_sqrt_b<Scalar_T>::array_size] =
1437  {
1438  1.0, 25.0/4.0, 69.0/4.0, 1771.0/64.0,
1439  7315.0/256.0, 20349.0/1024.0, 4845.0/512.0, 12597.0/4096.0,
1440  21879.0/32768.0, 12155.0/131072.0, 1001.0/131072.0, 1365.0/4194304.0,
1441  91.0/16777216.0, 1.0/67108864.0
1442  };
1443 
1444  template< >
1445  struct pade_sqrt_a<float>
1446  {
1447  static const int array_size = 10;
1448  static const float array[array_size];
1449  };
1450  template< >
1451  struct pade_sqrt_b<float>
1452  {
1453  static const int array_size = 10;
1454  static const float array[array_size];
1455  };
1456  const float pade_sqrt_a<float>::array[pade_sqrt_a<float>::array_size] =
1457  {
1458  1.0, 19.0/4.0, 19.0/2.0, 665.0/64.0,
1459  1729.0/256.0, 2717.0/1024.0, 627.0/1024.0, 627.0/8192.0,
1460  285.0/65536.0, 19.0/262144.0
1461  };
1462  const float pade_sqrt_b<float>::array[pade_sqrt_a<float>::array_size] =
1463  {
1464  1.0, 17.0/4.0, 15.0/2.0, 455.0/64.0,
1465  1001.0/256.0, 1287.0/1024.0, 231.0/1024.0, 165.0/8192.0,
1466  45.0/65536, 1.0/262144.0
1467  };
1468 
1469  template< >
1470  struct pade_sqrt_a<long double>
1471  {
1472  static const int array_size = 18;
1473  static const long double array[array_size];
1474  };
1475  template< >
1476  struct pade_sqrt_b<long double>
1477  {
1478  static const int array_size = 18;
1479  static const long double array[array_size];
1480  };
1481  const long double pade_sqrt_a<long double>::array[pade_sqrt_a<long double>::array_size] =
1482  {
1483  1.0L, 35.0L/4.0L, 35.0L, 5425.0L/64.0L,
1484  35525.0L/256.0L, 166257.0L/1024.0L, 143325.0L/1024.0L, 740025.0L/8192.0L,
1485  2877875.0L/65536.0L, 4206125.0L/262144.0L, 572033.0L/131072.0L, 1820105.0L/2097152.0L,
1486  1028755.0L/8388608.0L, 395675.0L/33554432.0L, 24225.0L/33554432.0L, 6783.0L/268435456.0L,
1487  1785.0L/4294967296.0L, 35.0L/17179869184.0L
1488  };
1489  const long double pade_sqrt_b<long double>::array[pade_sqrt_a<long double>::array_size] =
1490  {
1491  1.0L, 33.0L/4.0L, 31.0L, 4495.0L/64.0L,
1492  27405.0L/256.0L, 118755.0L/1024.0L, 94185.0L/1024.0L, 444015.0L/8192.0L,
1493  1562275.0L/65536.0L, 2042975.0L/262144.0L, 245157.0L/131072.0L, 676039.0L/2097152.0L,
1494  323323.0L/8388608.0L, 101745.0L/33554432.0L, 4845.0L/33554432.0L, 969.0L/268435456.0L,
1495  153.0L/4294967296.0L, 1.0L/17179869184.0L
1496  };
1497 
1498 #if defined(_GLUCAT_USE_QD)
1499  template< >
1500  struct pade_sqrt_a<dd_real>
1501  {
1502  static const int array_size = 22;
1503  static const dd_real array[array_size];
1504  };
1505  template< >
1506  struct pade_sqrt_b<dd_real>
1507  {
1508  static const int array_size = 22;
1509  static const dd_real array[array_size];
1510  };
1511  const dd_real pade_sqrt_a<dd_real>::array[pade_sqrt_a<dd_real>::array_size] =
1512  {
1513  dd_real("1"), dd_real("43")/dd_real("4"),
1514  dd_real("215")/dd_real("4"), dd_real("10621")/dd_real("64"),
1515  dd_real("90687")/dd_real("256"), dd_real("567987")/dd_real("1024"),
1516  dd_real("168861")/dd_real("256"), dd_real("1246355")/dd_real("2048"),
1517  dd_real("7228859")/dd_real("16384"), dd_real("16583853")/dd_real("65536"),
1518  dd_real("7538115")/dd_real("65536"), dd_real("173376645")/dd_real("4194304"),
1519  dd_real("195747825")/dd_real("16777216"), dd_real("171655785")/dd_real("67108864"),
1520  dd_real("14375115")/dd_real("33554432"), dd_real("14375115")/dd_real("268435456"),
1521  dd_real("20764055")/dd_real("4294967296"), dd_real("5167525")/dd_real("17179869184"),
1522  dd_real("206701")/dd_real("17179869184"), dd_real("76153")/dd_real("274877906944"),
1523  dd_real("3311")/dd_real("1099511627776") , dd_real("43")/dd_real("4398046511104")
1524  };
1525  const dd_real pade_sqrt_b<dd_real>::array[pade_sqrt_a<dd_real>::array_size] =
1526 {
1527  dd_real("1"), dd_real("41")/dd_real("4"),
1528  dd_real("195")/dd_real("4"), dd_real("9139")/dd_real("64"),
1529  dd_real("73815")/dd_real("256"), dd_real("435897")/dd_real("1024"),
1530  dd_real("121737")/dd_real("256"), dd_real("840565")/dd_real("2048"),
1531  dd_real("4539051")/dd_real("16384"), dd_real("9641775")/dd_real("65536"),
1532  dd_real("4032015")/dd_real("65536"), dd_real("84672315")/dd_real("4194304"),
1533  dd_real("86493225")/dd_real("16777216"), dd_real("67863915")/dd_real("67108864"),
1534  dd_real("5014575")/dd_real("33554432"), dd_real("4345965")/dd_real("268435456"),
1535  dd_real("5311735")/dd_real("4294967296"), dd_real("1081575")/dd_real("17179869184"),
1536  dd_real("33649")/dd_real("17179869184"), dd_real("8855")/dd_real("274877906944"),
1537  dd_real("231")/dd_real("1099511627776"), dd_real("1")/dd_real("4398046511104")
1538  };
1539 
1540  template< >
1541  struct pade_sqrt_a<qd_real>
1542  {
1543  static const int array_size = 34;
1544  static const qd_real array[array_size];
1545  };
1546  template< >
1547  struct pade_sqrt_b<qd_real>
1548  {
1549  static const int array_size = 34;
1550  static const qd_real array[array_size];
1551  };
1552  const qd_real pade_sqrt_a<qd_real>::array[pade_sqrt_a<qd_real>::array_size] =
1553  {
1554  qd_real("1"), qd_real("67")/qd_real("4"),
1555  qd_real("134"), qd_real("43617")/qd_real("64"),
1556  qd_real("633485")/qd_real("256"), qd_real("6992857")/qd_real("1024"),
1557  qd_real("15246721")/qd_real("1024"), qd_real("215632197")/qd_real("8192"),
1558  qd_real("2518145487")/qd_real("65536"), qd_real("12301285425")/qd_real("262144"),
1559  qd_real("6344873535")/qd_real("131072"), qd_real("89075432355")/qd_real("2097152"),
1560  qd_real("267226297065")/qd_real("8388608"), qd_real("687479618945")/qd_real("33554432"),
1561  qd_real("379874182975")/qd_real("33554432"), qd_real("1443521895305")/qd_real("268435456"),
1562  qd_real("9425348845815")/qd_real("4294967296"), qd_real("13195488384141")/qd_real("17179869184"),
1563  qd_real("987417498133")/qd_real("4294967296"), qd_real("8055248011085")/qd_real("137438953472"),
1564  qd_real("6958363175533")/qd_real("549755813888"), qd_real("5056698705201")/qd_real("2199023255552"),
1565  qd_real("766166470485")/qd_real("2199023255552"), qd_real("766166470485")/qd_real("17592186044416"),
1566  qd_real("623623871325")/qd_real("140737488355328"), qd_real("203123203803")/qd_real("562949953421312"),
1567  qd_real("6478601247")/qd_real("281474976710656"), qd_real("5038912081")/qd_real("4503599627370496"),
1568  qd_real("719844583")/qd_real("18014398509481984"), qd_real("71853815")/qd_real("72057594037927936"),
1569  qd_real("1165197")/qd_real("72057594037927936"), qd_real("87703")/qd_real("576460752303423488"),
1570  qd_real("12529")/qd_real("18446744073709551616"), qd_real("67")/qd_real("73786976294838206464")
1571  };
1572  const qd_real pade_sqrt_b<qd_real>::array[pade_sqrt_a<qd_real>::array_size] =
1573  {
1574  qd_real("1"), qd_real("65")/qd_real("4"),
1575  qd_real("126"), qd_real("39711")/qd_real("64"),
1576  qd_real("557845")/qd_real("256"), qd_real("5949147")/qd_real("1024"),
1577  qd_real("12515965")/qd_real("1024"), qd_real("170574723")/qd_real("8192"),
1578  qd_real("1916797311")/qd_real("65536"), qd_real("8996462475")/qd_real("262144"),
1579  qd_real("4450881435")/qd_real("131072"), qd_real("59826782925")/qd_real("2097152"),
1580  qd_real("171503444385")/qd_real("8388608"), qd_real("420696483235")/qd_real("33554432"),
1581  qd_real("221120793075")/qd_real("33554432"), qd_real("797168807855")/qd_real("268435456"),
1582 qd_real("4923689695575")/qd_real("4294967296"), qd_real("6499270398159")/qd_real("17179869184"),
1583  qd_real("456864812569")/qd_real("4294967296"), qd_real("3486599885395")/qd_real("137438953472"),
1584 qd_real("2804116503573")/qd_real("549755813888"), qd_real("1886827875075")/qd_real("2199023255552"),
1585  qd_real("263012370465")/qd_real("2199023255552"), qd_real("240141729555")/qd_real("17592186044416"),
1586  qd_real("176848560525")/qd_real("140737488355328"), qd_real("51538723353")/qd_real("562949953421312"),
1587  qd_real("1450433115")/qd_real("281474976710656"), qd_real("977699359")/qd_real("4503599627370496"),
1588  qd_real("118183439")/qd_real("18014398509481984"), qd_real("9652005")/qd_real("72057594037927936"),
1589  qd_real("121737")/qd_real("72057594037927936"), qd_real("6545")/qd_real("576460752303423488"),
1590  qd_real("561")/qd_real("18446744073709551616"), qd_real("1")/qd_real("73786976294838206464")
1591  };
1592 #endif
1593 }
1594 
1595 namespace glucat
1596 {
1598  template< typename Scalar_T, const index_t LO, const index_t HI >
1601  {
1602  // Reference: [GW], Section 4.3, pp318-322
1603  // Reference: [GL], Section 11.3, p572-576
1604  // Reference: [Z], Pade1
1605 
1606  typedef numeric_traits<Scalar_T> traits_t;
1607 
1608  if (val.isnan())
1609  return traits_t::NaN();
1610 
1611  check_complex(val, i, prechecked);
1612 
1614  {
1615  case precision_demoted:
1616  {
1617  typedef typename traits_t::demoted::type demoted_scalar_t;
1618  typedef matrix_multi<demoted_scalar_t,LO,HI> demoted_multivector_t;
1619 
1620  const demoted_multivector_t& demoted_val = demoted_multivector_t(val);
1621  const demoted_multivector_t& demoted_i = demoted_multivector_t(i);
1622 
1623  return matrix_sqrt(demoted_val, demoted_i);
1624  }
1625  break;
1626  case precision_promoted:
1627  {
1628  typedef typename traits_t::promoted::type promoted_scalar_t;
1629  typedef matrix_multi<promoted_scalar_t,LO,HI> promoted_multivector_t;
1630 
1631  const promoted_multivector_t& promoted_val = promoted_multivector_t(val);
1632  const promoted_multivector_t& promoted_i = promoted_multivector_t(i);
1633 
1634  return matrix_sqrt(promoted_val, promoted_i);
1635  }
1636  break;
1637  default:
1638  return matrix_sqrt(val, i);
1639  }
1640  }
1641 
1643  template< typename Scalar_T, const index_t LO, const index_t HI >
1646  {
1647  // Reference: [GW], Section 4.3, pp318-322
1648  // Reference: [GL], Section 11.3, p572-576
1649  // Reference: [Z], Pade1
1650 
1651  typedef numeric_traits<Scalar_T> traits_t;
1652 
1653  if (val.isnan())
1654  return traits_t::NaN();
1655 
1657 
1658  const Scalar_T realval = val.scalar();
1659  if (val == realval)
1660  {
1661  if (realval < Scalar_T(0))
1662  return i * traits_t::sqrt(-realval);
1663  else
1664  return traits_t::sqrt(realval);
1665  }
1666 
1667  static const Scalar_T sqrt_2 = traits_t::sqrt(Scalar_T(2));
1668 
1669 #if !defined(_GLUCAT_USE_EIGENVALUES)
1670  const multivector_t val2 = val*val;
1671  const Scalar_T real_val2 = val2.scalar();
1672  if (val2 == real_val2 && real_val2 > Scalar_T(0))
1673  return matrix_sqrt(-i * val, i) * (i + Scalar_T(1)) / sqrt_2;
1674 #endif
1675 
1676  // Scale val towards abs(A) == 1 or towards A == 1 as appropriate
1677  const Scalar_T scale =
1678  (realval != Scalar_T(0) && norm(val/realval - Scalar_T(1)) < Scalar_T(1))
1679  ? realval
1680  : (realval < Scalar_T(0))
1681  ? -abs(val)
1682  : abs(val);
1683  const Scalar_T sqrt_scale = traits_t::sqrt(traits_t::abs(scale));
1684  if (traits_t::isNaN_or_isInf(sqrt_scale))
1685  return traits_t::NaN();
1686 
1688  multivector_t rescale = sqrt_scale;
1689  if (scale < Scalar_T(0))
1690  rescale = i * sqrt_scale;
1691 
1692  const multivector_t& unitval = val / scale;
1693  const Scalar_T max_norm = Scalar_T(1.0/4.0);
1694 
1695 #if defined(_GLUCAT_USE_EIGENVALUES)
1696  multivector_t scaled_result;
1697  typedef typename multivector_t::matrix_t matrix_t;
1698 
1699  // What kind of eigenvalues does the matrix contain?
1701  switch (genus.m_eig_case)
1702  {
1704  scaled_result = matrix_sqrt(-i * unitval, i) * (i + Scalar_T(1)) / sqrt_2;
1705  break;
1706  case matrix::both_eig_case:
1707  {
1708  const Scalar_T safe_arg = genus.m_safe_arg;
1709  scaled_result = matrix_sqrt(exp(i*safe_arg) * unitval, i) * exp(-i*safe_arg/Scalar_T(2));
1710  }
1711  break;
1712  default:
1713  scaled_result =
1714  (norm(unitval - Scalar_T(1)) < max_norm)
1715  // Pade' approximation of square root
1716  ? pade_approx(pade_sqrt_a<Scalar_T>::array_size,
1717  pade_sqrt_a<Scalar_T>::array,
1718  pade_sqrt_b<Scalar_T>::array,
1719  unitval - Scalar_T(1))
1720  // Product form of Denman-Beavers square root iteration
1721  : db_sqrt(unitval);
1722  break;
1723  }
1724  if (scaled_result.isnan())
1725  return traits_t::NaN();
1726  else
1727  return scaled_result * rescale;
1728 #else
1729  const multivector_t& scaled_result =
1730  (norm(unitval - Scalar_T(1)) < max_norm)
1731  // Pade' approximation of square root
1732  ? pade_approx(pade_sqrt_a<Scalar_T>::array_size,
1733  pade_sqrt_a<Scalar_T>::array,
1734  pade_sqrt_b<Scalar_T>::array,
1735  unitval - Scalar_T(1))
1736  // Product form of Denman-Beavers square root iteration
1737  : db_sqrt(unitval);
1738  if (scaled_result.isnan())
1739  {
1740  if (inv(unitval).isnan())
1741  return traits_t::NaN();
1742 
1743  const multivector_t& mi_unitval = -i * unitval;
1744 
1745  const multivector_t& scaled_mi_result =
1746  (norm(mi_unitval - Scalar_T(1)) < max_norm)
1747  // Pade' approximation of square root
1748  ? pade_approx(pade_sqrt_a<Scalar_T>::array_size,
1749  pade_sqrt_a<Scalar_T>::array,
1750  pade_sqrt_b<Scalar_T>::array,
1751  mi_unitval - Scalar_T(1))
1752  // Product form of Denman-Beavers square root iteration
1753  : db_sqrt(mi_unitval);
1754  if (scaled_mi_result.isnan())
1755  return traits_t::NaN();
1756  else
1757  return scaled_mi_result * rescale * (i + Scalar_T(1)) / sqrt_2;
1758  }
1759  else
1760  return scaled_result * rescale;
1761 #endif
1762  }
1763 }
1764 
1765 namespace {
1767  // Reference: [Z], Pade1
1768  template< typename Scalar_T >
1769  struct pade_log_a
1770  {
1771  static const int array_size = 14;
1772  static const Scalar_T array[array_size];
1773  };
1774 
1776  // Reference: [Z], Pade1
1777  template< typename Scalar_T >
1778  struct pade_log_b
1779  {
1780  static const int array_size = 14;
1781  static const Scalar_T array[array_size];
1782  };
1783  template< typename Scalar_T >
1784  const Scalar_T pade_log_a<Scalar_T>::array[pade_log_a<Scalar_T>::array_size] =
1785  {
1786  0.0, 1.0, 6.0, 4741.0/300.0,
1787  1441.0/60.0, 107091.0/4600.0, 8638.0/575.0, 263111.0/40250.0,
1788  153081.0/80500.0, 395243.0/1101240.0, 28549.0/688275.0, 605453.0/228813200.0,
1789  785633.0/10296594000.0, 1145993.0/1873980108000.0
1790  };
1791  template< typename Scalar_T >
1792  const Scalar_T pade_log_b<Scalar_T>::array[pade_log_b<Scalar_T>::array_size] =
1793  {
1794  1.0, 13.0/2.0, 468.0/25.0, 1573.0/50.0,
1795  1573.0/46.0, 11583.0/460.0, 10296.0/805.0, 2574.0/575.0,
1796  11583.0/10925.0, 143.0/874.0, 572.0/37145.0, 117.0/148580.0,
1797  13.0/742900.0, 1.0/10400600.0
1798  };
1799 
1800  template< >
1801  struct pade_log_a<float>
1802  {
1803  static const int array_size = 10;
1804  static const float array[array_size];
1805  };
1806  template< >
1807  struct pade_log_b<float>
1808  {
1809  static const int array_size = 10;
1810  static const float array[array_size];
1811  };
1812  const float pade_log_a<float>::array[pade_log_a<float>::array_size] =
1813  {
1814  0.0, 1.0, 4.0, 1337.0/204.0,
1815  385.0/68.0, 1879.0/680.0, 193.0/255.0, 197.0/1820.0,
1816  419.0/61880.0, 7129.0/61261200.0
1817  };
1818  const float pade_log_b<float>::array[pade_log_a<float>::array_size] =
1819  {
1820  1.0, 9.0/2.0, 144.0/17.0, 147.0/17.0,
1821  441.0/85.0, 63.0/34.0, 84.0/221.0, 9.0/221.0,
1822  9.0/4862.0, 1.0/48620.0
1823  };
1824 
1825  template< >
1826  struct pade_log_a<long double>
1827  {
1828  static const int array_size = 18;
1829  static const long double array[array_size];
1830  };
1831  template< >
1832  struct pade_log_b<long double>
1833  {
1834  static const int array_size = 18;
1835  static const long double array[array_size];
1836  };
1837  const long double pade_log_a<long double>::array[pade_log_a<long double>::array_size] =
1838  {
1839  0.0L, 1.0L, 8.0L, 3835.0L/132.0L,
1840  8365.0L/132.0L, 11363807.0L/122760.0L, 162981.0L/1705.0L, 9036157.0L/125860.0L,
1841  18009875.0L/453096.0L, 44211925.0L/2718576.0L, 4149566.0L/849555.0L, 16973929.0L/16020180.0L,
1842  172459.0L/1068012.0L, 116317061.0L/7025382936.0L, 19679783.0L/18441630207.0L, 23763863.0L/614721006900.0L,
1843  50747.0L/79318839600.0L, 42142223.0L/14295951736466400.0L
1844  };
1845  const long double pade_log_b<long double>::array[pade_log_a<long double>::array_size] =
1846  {
1847  1.0L, 17.0L/2.0L, 1088.0L/33.0L, 850.0L/11.0L,
1848  41650.0L/341.0L, 140777.0L/1023.0L, 1126216.0L/9889.0L, 63206.0L/899.0L,
1849  790075.0L/24273.0L, 60775.0L/5394.0L, 38896.0L/13485.0L, 21658.0L/40455.0L,
1850  21658.0L/310155.0L, 4165.0L/682341.0L, 680.0L/2047023.0L, 34.0L/3411705.0L,
1851  17.0L/129644790.0L, 1.0L/2333606220
1852  };
1853 #if defined(_GLUCAT_USE_QD)
1854  template< >
1855  struct pade_log_a<dd_real>
1856  {
1857  static const int array_size = 22;
1858  static const dd_real array[array_size];
1859  };
1860  template< >
1861  struct pade_log_b<dd_real>
1862  {
1863  static const int array_size = 22;
1864  static const dd_real array[array_size];
1865  };
1866  const dd_real pade_log_a<dd_real>::array[pade_log_a<dd_real>::array_size] =
1867  {
1868  dd_real("0"), dd_real("1"),
1869  dd_real("10"), dd_real("22781")/dd_real("492"),
1870  dd_real("21603")/dd_real("164"), dd_real("5492649")/dd_real("21320"),
1871  dd_real("978724")/dd_real("2665"), dd_real("4191605")/dd_real("10619"),
1872  dd_real("12874933")/dd_real("39442"), dd_real("11473457")/dd_real("54612"),
1873  dd_real("2406734")/dd_real("22755"), dd_real("166770367")/dd_real("4004880"),
1874  dd_real("30653165")/dd_real("2402928"), dd_real("647746389")/dd_real("215195552"),
1875  dd_real("25346331")/dd_real("47074027"), dd_real("278270613")/dd_real("3900419380"),
1876  dd_real("105689791")/dd_real("15601677520"), dd_real("606046475")/dd_real("1379188292768"),
1877  dd_real("969715")/dd_real("53502994116"), dd_real("11098301")/dd_real("26204577562592"),
1878  dd_real("118999")/dd_real("26204577562592"), dd_real("18858053")/dd_real("1392249205900512960")
1879  };
1880  const dd_real pade_log_b<dd_real>::array[pade_log_a<dd_real>::array_size] =
1881  {
1882  dd_real("1"), dd_real("21")/dd_real("2"),
1883  dd_real("2100")/dd_real("41"), dd_real("12635")/dd_real("82"),
1884  dd_real("341145")/dd_real("1066"), dd_real("1037799")/dd_real("2132"),
1885  dd_real("11069856")/dd_real("19721"), dd_real("9883800")/dd_real("19721"),
1886  dd_real("6918660")/dd_real("19721"), dd_real("293930")/dd_real("1517"),
1887  dd_real("1410864")/dd_real("16687"), dd_real("88179")/dd_real("3034"),
1888  dd_real("734825")/dd_real("94054"), dd_real("305235")/dd_real("188108"),
1889  dd_real("348840")/dd_real("1363783"), dd_real("40698")/dd_real("1363783"),
1890  dd_real("6783")/dd_real("2727566"), dd_real("9975")/dd_real("70916716"),
1891  dd_real("266")/dd_real("53187537"), dd_real("7")/dd_real("70916716"),
1892  dd_real("7")/dd_real("8155422340"), dd_real("1")/dd_real("538257874440")
1893  };
1894 
1895  template< >
1896  struct pade_log_a<qd_real>
1897  {
1898  static const int array_size = 34;
1899  static const qd_real array[array_size];
1900  };
1901  template< >
1902  struct pade_log_b<qd_real>
1903  {
1904  static const int array_size = 34;
1905  static const qd_real array[array_size];
1906  };
1907  const qd_real pade_log_a<qd_real>::array[pade_log_a<qd_real>::array_size] =
1908  {
1909  qd_real("0"), qd_real("1"),
1910  qd_real("16"), qd_real("95201")/qd_real("780"),
1911  qd_real("30721")/qd_real("52"), qd_real("7416257")/qd_real("3640"),
1912  qd_real("1039099")/qd_real("195"), qd_real("6097772319")/qd_real("555100"),
1913  qd_real("1564058073")/qd_real("85400"), qd_real("30404640205")/qd_real("1209264"),
1914  qd_real("725351278")/qd_real("25193"), qd_real("4092322670789")/qd_real("147429436"),
1915  qd_real("4559713849589")/qd_real("201040140"), qd_real("5049361751189")/qd_real("320023080"),
1916  qd_real("74979677195")/qd_real("8000577"), qd_real("16569850691873")/qd_real("3481514244"),
1917  qd_real("1065906022369")/qd_real("515779888"), qd_real("335956770855841")/qd_real("438412904800"),
1918  qd_real("1462444287585964")/qd_real("6041877844275"), qd_real("397242326339851")/qd_real("6122436215532"),
1919  qd_real("64211291334131")/qd_real("4373168725380"), qd_real("142322343550859")/qd_real("51080680851480"),
1920  qd_real("154355972958659")/qd_real("351179680853925"), qd_real("167483568676259")/qd_real("2937139148960100"),
1921  qd_real("4230788929433")/qd_real("704913395750424"), qd_real("197968763176019")/qd_real("392923948371995600"),
1922  qd_real("10537522306718")/qd_real("319250708052246425"), qd_real("236648286272519")/qd_real("144249197475035425500"),
1923  qd_real("260715545088119")/qd_real("4375558990076074573500"), qd_real("289596255666839")/qd_real("192874640282553367199880"),
1924  qd_real("8802625510547")/qd_real("361639950529787563499775"), qd_real("373831661521439")/qd_real("1659204093030665341336967700"),
1925  qd_real("446033437968239")/qd_real("464577146048586295574350956000"), qd_real("53676090078349")/qd_real("47386868896955802148583797512000")
1926  };
1927  const qd_real pade_log_b<qd_real>::array[pade_log_a<qd_real>::array_size] =
1928  {
1929  qd_real("1"), qd_real("33")/qd_real("2"),
1930  qd_real("8448")/qd_real("65"), qd_real("42284")/qd_real("65"),
1931  qd_real("211420")/qd_real("91"), qd_real("573562")/qd_real("91"),
1932  qd_real("32119472")/qd_real("2379"), qd_real("92917044")/qd_real("3965"),
1933  qd_real("603960786")/qd_real("17995"), qd_real("144626625")/qd_real("3599"),
1934  qd_real("2776831200")/qd_real("68381"), qd_real("16692542100")/qd_real("478667"),
1935  qd_real("12241197540")/qd_real("478667"), qd_real("1098569010")/qd_real("68381"),
1936  qd_real("31387686000")/qd_real("3624193"), qd_real("9939433900")/qd_real("2479711"),
1937  qd_real("67091178825")/qd_real("42155087"), qd_real("2683647153")/qd_real("4959422"),
1938  qd_real("19083713088")/qd_real("121505839"), qd_real("4708152900")/qd_real("121505839"),
1939  qd_real("941630580")/qd_real("116546417"), qd_real("88704330")/qd_real("62755763"),
1940  qd_real("12902448")/qd_real("62755763"), qd_real("1542684")/qd_real("62755763"),
1941  qd_real("6427850")/qd_real("2698497809"), qd_real("3471039")/qd_real("18889484663"),
1942  qd_real("8544096")/qd_real("774468871183"), qd_real("39556")/qd_real("79027435835"),
1943  qd_real("118668")/qd_real("7191496660985"), qd_real("10230")/qd_real("27327687311743"),
1944  qd_real("5456")/qd_real("1011124430534491"), qd_real("44")/qd_real("1011124430534491"),
1945  qd_real("11")/qd_real("70778710137414370"), qd_real("1")/qd_real("7219428434016265740")
1946  };
1947 #endif
1948 }
1949 
1950 namespace glucat{
1952  template< typename Scalar_T, const index_t LO, const index_t HI >
1953  static
1956  {
1957  // Reference: [GW], Section 4.3, pp318-322
1958  // Reference: [CHKL]
1959  // Reference: [GL], Section 11.3, p572-576
1960  // Reference: [Z], Pade1
1961 
1962  typedef numeric_traits<Scalar_T> traits_t;
1963  if (val == Scalar_T(0) || val.isnan())
1964  return traits_t::NaN();
1965  else
1966  return pade_approx(pade_log_a<Scalar_T>::array_size,
1967  pade_log_a<Scalar_T>::array,
1968  pade_log_b<Scalar_T>::array,
1969  val - Scalar_T(1));
1970  }
1971 
1973  template< typename Scalar_T, const index_t LO, const index_t HI >
1974  static
1977  {
1978  // Reference: [CHKL]
1980  typedef numeric_traits<Scalar_T> traits_t;
1981  if (val == Scalar_T(0) || val.isnan())
1982  return traits_t::NaN();
1983 
1984  typedef std::numeric_limits<Scalar_T> limits_t;
1985  static const Scalar_T epsilon = limits_t::epsilon();
1986  static const Scalar_T max_inner_norm = traits_t::pow(epsilon, 2);
1987  static const Scalar_T max_outer_norm = Scalar_T(6.0/limits_t::digits);
1988  multivector_t Y = val;
1989  multivector_t E = Scalar_T(0);
1990  Scalar_T norm_Y_1;
1991  int outer_step;
1992  Scalar_T pow_2_outer_step = Scalar_T(1);
1993  Scalar_T pow_4_outer_step = Scalar_T(1);
1994  for (outer_step = 0, norm_Y_1 = norm(Y - Scalar_T(1));
1995  outer_step != Tune_P::log_max_outer_steps && norm_Y_1 * pow_2_outer_step > max_outer_norm;
1996  ++outer_step, norm_Y_1 = norm(Y - Scalar_T(1)))
1997  {
1998  if (Y == Scalar_T(0) || Y.isnan())
1999  return traits_t::NaN();
2000 
2001  // Incomplete product form of Denman-Beavers square root iteration
2002  multivector_t M = Y;
2003  for (int
2004  inner_step = 0;
2005  inner_step != Tune_P::log_max_inner_steps &&
2006  norm(M - Scalar_T(1)) * pow_4_outer_step > max_inner_norm;
2007  ++inner_step)
2008  db_step(M, Y);
2009 
2010  E += (M - Scalar_T(1)) * pow_2_outer_step;
2011  pow_2_outer_step *= Scalar_T(2);
2012  pow_4_outer_step *= Scalar_T(4);
2013  }
2014  if (outer_step == Tune_P::log_max_outer_steps && norm_Y_1 * pow_2_outer_step > max_outer_norm)
2015  return traits_t::NaN();
2016  else
2017  return pade_log(Y) * pow_2_outer_step - E;
2018  }
2019 
2021  template< typename Scalar_T, const index_t LO, const index_t HI >
2024  {
2025  typedef numeric_traits<Scalar_T> traits_t;
2026 
2027  if (val == Scalar_T(0) || val.isnan())
2028  return traits_t::NaN();
2029 
2030  check_complex(val, i, prechecked);
2031 
2033  {
2034  case precision_demoted:
2035  {
2036  typedef typename traits_t::demoted::type demoted_scalar_t;
2037  typedef matrix_multi<demoted_scalar_t,LO,HI> demoted_multivector_t;
2038 
2039  const demoted_multivector_t& demoted_val = demoted_multivector_t(val);
2040  const demoted_multivector_t& demoted_i = demoted_multivector_t(i);
2041 
2042  return matrix_log(demoted_val, demoted_i);
2043  }
2044  break;
2045  case precision_promoted:
2046  {
2047  typedef typename traits_t::promoted::type promoted_scalar_t;
2048  typedef matrix_multi<promoted_scalar_t,LO,HI> promoted_multivector_t;
2049 
2050  const promoted_multivector_t& promoted_val = promoted_multivector_t(val);
2051  const promoted_multivector_t& promoted_i = promoted_multivector_t(i);
2052 
2053  return matrix_log(promoted_val, promoted_i);
2054  }
2055  break;
2056  default:
2057  return matrix_log(val, i);
2058  }
2059  }
2060 
2062  template< typename Scalar_T, const index_t LO, const index_t HI >
2065  {
2066  // Scaled incomplete square root cascade and scaled Pade' approximation of log
2067  // Reference: [CHKL]
2068 
2069  typedef numeric_traits<Scalar_T> traits_t;
2070  if (val == Scalar_T(0) || val.isnan())
2071  return traits_t::NaN();
2072 
2073  static const Scalar_T pi = traits_t::pi();
2074  const Scalar_T realval = val.scalar();
2075  if (val == realval)
2076  {
2077  if (realval < Scalar_T(0))
2078  return i * pi + traits_t::log(-realval);
2079  else
2080  return traits_t::log(realval);
2081  }
2083 #if !defined(_GLUCAT_USE_EIGENVALUES)
2084  const multivector_t val2 = val*val;
2085  const Scalar_T real_val2 = val2.scalar();
2086  if (val2 == real_val2 && real_val2 > 0)
2087  return matrix_log(-i * val, i) + i * pi/Scalar_T(2);
2088 #endif
2089  // Scale val towards abs(A) == 1 or towards A == 1 as appropriate
2090  const Scalar_T max_norm = Scalar_T(1.0/9.0);
2091  const Scalar_T scale =
2092  (realval != Scalar_T(0) && norm(val/realval - Scalar_T(1)) < max_norm)
2093  ? realval
2094  : (realval < Scalar_T(0))
2095  ? -abs(val)
2096  : abs(val);
2097  if (scale == Scalar_T(0))
2098  return traits_t::NaN();
2099 
2100  const Scalar_T log_scale = traits_t::log(traits_t::abs(scale));
2101  multivector_t rescale = log_scale;
2102  if (scale < Scalar_T(0))
2103  rescale = i * pi + log_scale;
2104  const multivector_t unitval = val/scale;
2105  if (inv(unitval).isnan())
2106  return traits_t::NaN();
2107 #if defined(_GLUCAT_USE_EIGENVALUES)
2108  multivector_t scaled_result;
2109  typedef typename multivector_t::matrix_t matrix_t;
2110 
2111  // What kind of eigenvalues does the matrix contain?
2113  switch (genus.m_eig_case)
2114  {
2116  scaled_result = matrix_log(-i * unitval, i) + i * pi/Scalar_T(2);
2117  break;
2118  case matrix::both_eig_case:
2119  {
2120  const Scalar_T safe_arg = genus.m_safe_arg;
2121  scaled_result = matrix_log(exp(i*safe_arg) * unitval, i) - i * safe_arg;
2122  }
2123  break;
2124  default:
2125  scaled_result = cascade_log(unitval);
2126  break;
2127  }
2128 #else
2129  multivector_t scaled_result = cascade_log(unitval);
2130 #endif
2131  if (scaled_result.isnan())
2132  return traits_t::NaN();
2133  else
2134  return scaled_result + rescale;
2135  }
2136 
2138  template< typename Scalar_T, const index_t LO, const index_t HI >
2141  {
2142  typedef numeric_traits<Scalar_T> traits_t;
2143  if (val.isnan())
2144  return traits_t::NaN();
2145 
2146  const Scalar_T s = scalar(val);
2147  if (val == s)
2148  return traits_t::exp(s);
2149 
2151  {
2152  case precision_demoted:
2153  {
2154  typedef typename traits_t::demoted::type demoted_scalar_t;
2155  typedef matrix_multi<demoted_scalar_t,LO,HI> demoted_multivector_t;
2156 
2157  const demoted_multivector_t& demoted_val = demoted_multivector_t(val);
2158  return clifford_exp(demoted_val);
2159  }
2160  break;
2161  case precision_promoted:
2162  {
2163  typedef typename traits_t::promoted::type promoted_scalar_t;
2164  typedef matrix_multi<promoted_scalar_t,LO,HI> promoted_multivector_t;
2165 
2166  const promoted_multivector_t& promoted_val = promoted_multivector_t(val);
2167  return clifford_exp(promoted_val);
2168  }
2169  break;
2170  default:
2171  return clifford_exp(val);
2172  }
2173  }
2174 }
2175 #endif // _GLUCAT_MATRIX_MULTI_IMP_H
multivector_t & centre_pm4_qp4(index_t &p, index_t &q)
Subalgebra isomorphism: R_{p,q} to R_{p-4,q+4}.
index_set< LO, HI > index_set_t
Definition: framed_multi.h:152
friend const matrix_multi_t operator|(const matrix_multi_t &lhs, const matrix_multi_t &rhs)
static Multivector_T fast(const Matrix_T &X, index_t level)
Inverse generalized Fast Fourier Transform.
std::vector< Scalar_T > vector_t
Definition: matrix_multi.h:141
multivector_t & centre_qp1_pm1(index_t &p, index_t &q)
Subalgebra isomorphism: R_{p,q} to R_{q+1,p-1}.
virtual index_t grade() const =0
Maximum of the grades of each term.
static const matrix_multi_t random(const index_set_t frm, Scalar_T fill=Scalar_T(1))
Random multivector within a frame.
const Multivector< Scalar_T, LO, HI > sqrt(const Multivector< Scalar_T, LO, HI > &val, const Multivector< Scalar_T, LO, HI > &i, const bool prechecked=false)
Square root of multivector with specified complexifier.
index_t count_pos() const
Number of positive indices included in set.
const index_set_t fold() const
Fold this index set within itself as a frame.
static const matrix_multi< Scalar_T, LO, HI > pade_log(const matrix_multi< Scalar_T, LO, HI > &val)
Pade&#39; approximation of log.
Scalar_T abs(const Multivector< Scalar_T, LO, HI > &val)
Absolute value == sqrt(norm)
A framed_multi<Scalar_T,LO,HI> is a framed approximation to a multivector.
Definition: framed_multi.h:65
const Multivector< Scalar_T, LO, HI > log(const Multivector< Scalar_T, LO, HI > &val, const Multivector< Scalar_T, LO, HI > &i, const bool prechecked=false)
Natural logarithm of multivector with specified complexifier.
static generator_table< Matrix_T > & generator()
Single instance of generator table.
map_t::const_iterator const_iterator
Definition: framed_multi.h:196
virtual multivector_t & operator|=(const multivector_t &rhs)=0
Transformation via twisted adjoint action.
friend const matrix_multi_t operator&(const matrix_multi_t &lhs, const matrix_multi_t &rhs)
index_t offset_level(const index_t p, const index_t q)
Determine the log2 dim corresponding to signature p, q.
virtual multivector_t & operator*=(const Scalar_T &scr)=0
Product of multivector and scalar.
static const matrix_multi< Scalar_T, LO, HI > pade_approx(const int array_size, const Scalar_T a[], const Scalar_T b[], const matrix_multi< Scalar_T, LO, HI > &X)
Pade&#39; approximation.
virtual const vector_t vector_part() const =0
Vector part of multivector, as a vector_t with respect to frame()
static void db_step(matrix_multi< Scalar_T, LO, HI > &M, matrix_multi< Scalar_T, LO, HI > &Y)
Single step of product form of Denman-Beavers square root iteration.
float pi
Definition: PyClical.pyx:1857
static const index_t offset_to_super[]
Offsets between the current signature and that of the real superalgebra.
Definition: generation.h:81
virtual multivector_t & operator/=(const Scalar_T &scr)=0
Quotient of multivector and scalar.
static const matrix_multi< Scalar_T, LO, HI > cascade_log(const matrix_multi< Scalar_T, LO, HI > &val)
Incomplete square root cascade and Pade&#39; approximation of log.
friend const matrix_multi_t operator^(const matrix_multi_t &lhs, const matrix_multi_t &rhs)
index_t max() const
Maximum member.
friend std::istream & operator>>(std::istream &s, multivector_t &val)
virtual const multivector_t inv() const =0
Geometric multiplicative inverse.
_GLUCAT_CLIFFORD_ALGEBRA_OPERATIONS multivector_t & operator=(const multivector_t &rhs)
Assignment operator.
virtual Scalar_T operator[](const index_set_t ist) const =0
Subscripting: map from index set to scalar coordinate.
friend const matrix_multi_t operator%(const matrix_multi_t &lhs, const matrix_multi_t &rhs)
Extra traits which extend numeric limits.
Definition: scalar.h:46
Table of basis elements used as a cache by basis_element()
virtual const multivector_t reverse() const =0
Reversion, eg. {1}*{2} -> {2}*{1}.
friend Scalar_T star(const matrix_multi_t &lhs, const matrix_multi_t &rhs)
const Multivector< Scalar_T, LO, HI > clifford_exp(const Multivector< Scalar_T, LO, HI > &val)
Exponential of multivector.
Matrix_T::value_type norm_frob2(const Matrix_T &val)
Square of Frobenius norm.
Definition: matrix_imp.h:413
const Multivector< Scalar_T, LO, HI > pow(const Multivector< Scalar_T, LO, HI > &lhs, int rhs)
Integer power of multivector.
virtual const multivector_t even() const =0
Even part of multivector, sum of even grade terms.
bool isnan(const Matrix_T &m)
Not a Number.
Definition: matrix_imp.h:292
index_t count() const
Cardinality: Number of indices included in set.
static void check_complex(const Multivector< Scalar_T, LO, HI > &val, const Multivector< Scalar_T, LO, HI > &i, const bool prechecked=false)
Check that i is a valid complexifier for val.
matrix_multi multivector_t
Definition: matrix_multi.h:136
index_set_t m_frame
Index set representing the frame for the subalgebra which contains the multivector.
Definition: matrix_multi.h:275
const framed_multi< Other_Scalar_T, LO, HI > fast_framed_multi() const
Use inverse generalized FFT to construct a framed_multi_t.
A matrix_multi<Scalar_T,LO,HI> is a matrix approximation to a multivector.
Definition: framed_multi.h:68
virtual Scalar_T norm() const =0
Scalar_T norm == sum of norm of coordinates.
Structure containing classification of eigenvalues.
Definition: matrix.h:131
index_t min() const
Minimum member.
index_set< LO, HI > index_set_t
Definition: matrix_multi.h:139
std::pair< const index_set_t, Scalar_T > term_t
Definition: matrix_multi.h:140
virtual Scalar_T max_abs() const =0
Maximum of absolute values of components of multivector: multivector infinity norm.
const framed_multi< Scalar_T, LO, HI > exp(const framed_multi< Scalar_T, LO, HI > &val)
Exponential of multivector.
virtual const multivector_t operator()(index_t grade) const =0
Pure grade-vector part.
error< multivector_t > error_t
Definition: matrix_multi.h:142
const matrix_multi_t fast_matrix_multi(const index_set_t frm) const
Use generalized FFT to construct a matrix_multi_t.
eig_case_t m_eig_case
What kind of eigenvalues does the matrix contain?
Definition: matrix.h:135
friend const matrix_multi< Other_Scalar_T, Other_LO, Other_HI > matrix_log(const matrix_multi< Other_Scalar_T, Other_LO, Other_HI > &val, const matrix_multi< Other_Scalar_T, Other_LO, Other_HI > &i)
Abstract exception class.
Definition: errors.h:41
std::pair< const index_set_t, Scalar_T > term_t
Definition: framed_multi.h:153
virtual const multivector_t involute() const =0
Main involution, each {i} is replaced by -{i} in each term, eg. {1} -> -{1}.
multivector_t & centre_pp4_qm4(index_t &p, index_t &q)
Subalgebra isomorphism: R_{p,q} to R_{p+4,q-4}.
static const precision_t function_precision
Precision used for exp, log and sqrt functions.
Definition: global.h:180
virtual const multivector_t outer_pow(int m) const =0
Outer product power.
virtual void write(const std::string &msg="") const =0
Write formatted multivector to output.
friend const matrix_multi_t operator/(const matrix_multi_t &lhs, const matrix_multi_t &rhs)
ublas::compressed_matrix< Scalar_T, orientation_t > matrix_t
Definition: matrix_multi.h:157
multivector_t unfold(const index_set_t frm) const
Subalgebra isomorphism: unfold each term within the given frame.
const basis_matrix_t basis_element(const index_set< LO, HI > &ist) const
Create a basis element matrix within the current frame.
friend const index_set< Other_LO, Other_HI > reframe(const matrix_multi< Other_Scalar_T, Other_LO, Other_HI > &lhs, const matrix_multi< Other_Scalar_T, Other_LO, Other_HI > &rhs, matrix_multi< Other_Scalar_T, Other_LO, Other_HI > &lhs_reframed, matrix_multi< Other_Scalar_T, Other_LO, Other_HI > &rhs_reframed)
virtual const multivector_t odd() const =0
Odd part of multivector, sum of odd grade terms.
static basis_table & basis()
Single instance of basis table.
framed_multi< Scalar_T, LO, HI > framed_multi_t
Definition: matrix_multi.h:143
matrix_t::size_type matrix_index_t
Definition: matrix_multi.h:159
virtual multivector_t & operator-=(const multivector_t &rhs)=0
Geometric difference.
friend const matrix_multi< Other_Scalar_T, Other_LO, Other_HI > matrix_sqrt(const matrix_multi< Other_Scalar_T, Other_LO, Other_HI > &val, const matrix_multi< Other_Scalar_T, Other_LO, Other_HI > &i)
static const std::string classname()
Class name used in messages.
matrix_t m_matrix
Matrix value representing the multivector within the folded frame.
Definition: matrix_multi.h:277
static const matrix_multi< Scalar_T, LO, HI > db_sqrt(const matrix_multi< Scalar_T, LO, HI > &val)
Product form of Denman-Beavers square root iteration.
Index set class based on std::bitset<> in Gnu standard C++ library.
Definition: index_set.h:45
virtual const multivector_t conj() const =0
Conjugation, reverse o involute == involute o reverse.
virtual multivector_t & operator^=(const multivector_t &rhs)=0
Outer product.
const RHS_T::expression_type mono_prod(const ublas::matrix_expression< LHS_T > &lhs, const ublas::matrix_expression< RHS_T > &rhs)
Product of monomial matrices.
Definition: matrix_imp.h:326
multivector_t & operator+=(const term_t &rhs)
Add a term, if non-zero.
virtual const multivector_t truncated(const Scalar_T &limit=Scalar_T(DEFAULT_TRUNCATION)) const =0
Remove all terms with relative size smaller than limit.
virtual bool operator==(const multivector_t &val) const =0
Test for equality of multivectors.
Matrix_T::value_type trace(const Matrix_T &val)
Matrix trace.
Definition: matrix_imp.h:437
Scalar_T m_safe_arg
Argument such that exp(pi-m_safe_arg) lies between arguments of eigenvalues.
Definition: matrix.h:137
ublas::compressed_matrix< int, orientation_t > basis_matrix_t
Definition: matrix_multi.h:152
int index_t
Size of index_t should be enough to represent LO, HI.
Definition: global.h:77
const RHS_T signed_perm_nork(const LHS_T &lhs, const RHS_T &rhs)
Left inverse of Kronecker product where lhs is a signed permutation matrix.
Definition: matrix_imp.h:237
def e(obj)
Definition: PyClical.pyx:1887
eig_genus< Matrix_T > classify_eigenvalues(const Matrix_T &val)
Classify the eigenvalues of a matrix.
Definition: matrix_imp.h:526
static Matrix_Index_T folded_dim(const index_set< LO, HI > &sub)
Determine the matrix dimension of the fold of a subalegbra.
LHS_T pos_mod(LHS_T lhs, RHS_T rhs)
Modulo function which works reliably for lhs < 0.
Definition: global.h:187
virtual multivector_t & operator&=(const multivector_t &rhs)=0
Inner product.
static const framed_multi_t random(const index_set_t frm, Scalar_T fill=Scalar_T(1))
Random multivector within a frame.
virtual Scalar_T quad() const =0
Scalar_T quadratic form == (rev(x)*x)(0)
index_t count_neg() const
Number of negative indices included in set.
virtual bool isnan() const =0
Check if a multivector contains any IEEE NaN values.
friend const matrix_multi_t operator*(const matrix_multi_t &lhs, const matrix_multi_t &rhs)
virtual const index_set_t frame() const =0
Subalgebra generated by all generators of terms of given multivector.
matrix_multi()
Default constructor.
virtual multivector_t & operator%=(const multivector_t &rhs)=0
Contraction.