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