glucat  0.8.4
framed_multi_imp.h
Go to the documentation of this file.
1 #ifndef _GLUCAT_FRAMED_MULTI_IMP_H
2 #define _GLUCAT_FRAMED_MULTI_IMP_H
3 /***************************************************************************
4  GluCat : Generic library of universal Clifford algebra templates
5  framed_multi_imp.h : Implement the coordinate map representation of a
6  Clifford algebra element
7  -------------------
8  begin : Sun 2001-12-09
9  copyright : (C) 2001-2016 by Paul C. Leopardi
10  ***************************************************************************
11 
12  This library is free software: you can redistribute it and/or modify
13  it under the terms of the GNU Lesser General Public License as published
14  by the Free Software Foundation, either version 3 of the License, or
15  (at your option) any later version.
16 
17  This library is distributed in the hope that it will be useful,
18  but WITHOUT ANY WARRANTY; without even the implied warranty of
19  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  GNU Lesser General Public License for more details.
21 
22  You should have received a copy of the GNU Lesser General Public License
23  along with this library. If not, see <http://www.gnu.org/licenses/>.
24 
25  ***************************************************************************
26  This library is based on a prototype written by Arvind Raja and was
27  licensed under the LGPL with permission of the author. See Arvind Raja,
28  "Object-oriented implementations of Clifford algebras in C++: a prototype",
29  in Ablamowicz, Lounesto and Parra (eds.)
30  "Clifford algebras with numeric and symbolic computations", Birkhauser, 1996.
31  ***************************************************************************
32  See also Arvind Raja's original header comments in glucat.h
33  ***************************************************************************/
34 
35 #include "glucat/framed_multi.h"
36 
37 #include "glucat/random.h"
38 
39 #if defined(_GLUCAT_USE_BOOST_POOL_ALLOC)
40 // Use the Boost pool allocator
41 #include <boost/pool/pool_alloc.hpp>
42 #endif
43 
44 #include <sstream>
45 #include <fstream>
46 
47 namespace glucat
48 {
50  template< typename Scalar_T, const index_t LO, const index_t HI >
51  const std::string
53  classname()
54  { return "framed_multi"; }
55 
56 #if defined(_GLUCAT_MAP_IS_HASH)
57 #define _GLUCAT_HASH_N(x) (x)
58 #define _GLUCAT_HASH_SIZE_T(x) (typename multivector_t::hash_size_t)(x)
59 #else
60 #define _GLUCAT_HASH_N(x)
61 #define _GLUCAT_HASH_SIZE_T(x)
62 #endif
63 
65  template< typename Scalar_T, const index_t LO, const index_t HI >
68  : map_t(_GLUCAT_HASH_N(0))
69  { }
70 
72  template< typename Scalar_T, const index_t LO, const index_t HI >
73  framed_multi<Scalar_T,LO,HI>::
74  framed_multi(const hash_size_t& hash_size)
75  : map_t(_GLUCAT_HASH_N(hash_size()))
76  { }
77 
79  template< typename Scalar_T, const index_t LO, const index_t HI >
80  template< typename Other_Scalar_T >
81  framed_multi<Scalar_T,LO,HI>::
82  framed_multi(const framed_multi<Other_Scalar_T,LO,HI>& val)
83  : map_t(_GLUCAT_HASH_N(val.size()))
84  {
85  typedef framed_multi<Other_Scalar_T,LO,HI> other_framed_multi_t;
86  typedef typename other_framed_multi_t::const_iterator other_const_iterator;
87  const other_const_iterator val_begin = val.begin();
88  const other_const_iterator val_end = val.end();
89  for (other_const_iterator val_it = val_begin; val_it != val_end; ++val_it)
90  this->insert(term_t(val_it->first, numeric_traits<Scalar_T>::to_scalar_t(val_it->second)));
91  }
92 
94  template< typename Scalar_T, const index_t LO, const index_t HI >
95  template< typename Other_Scalar_T >
98  const index_set_t frm, const bool prechecked)
99  : map_t(_GLUCAT_HASH_N(val.size()))
100  {
101  typedef framed_multi<Other_Scalar_T,LO,HI> other_framed_multi_t;
102  typedef typename other_framed_multi_t::const_iterator other_const_iterator;
103  const other_const_iterator val_begin = val.begin();
104  const other_const_iterator val_end = val.end();
105  for (other_const_iterator val_it = val_begin; val_it != val_end; ++val_it)
106  this->insert(term_t(val_it->first, static_cast<Scalar_T>(val_it->second)));
107  }
108 
110  template< typename Scalar_T, const index_t LO, const index_t HI >
112  framed_multi(const framed_multi_t& val,
113  const index_set_t frm, const bool prechecked)
114  : map_t(val)
115  { }
116 
118  template< typename Scalar_T, const index_t LO, const index_t HI >
120  framed_multi(const index_set_t ist, const Scalar_T& crd)
121  : map_t(_GLUCAT_HASH_N(1))
122  {
123  if (crd != Scalar_T(0))
124  this->insert(term_t(ist, crd));
125  }
126 
128  template< typename Scalar_T, const index_t LO, const index_t HI >
130  framed_multi(const index_set_t ist, const Scalar_T& crd,
131  const index_set_t frm, const bool prechecked)
132  : map_t(_GLUCAT_HASH_N(1))
133  {
134  if (!prechecked && (ist | frm) != frm)
135  throw error_t("multivector_t(ist,crd,frm): cannot initialize with value outside of frame");
136  if (crd != Scalar_T(0))
137  this->insert(term_t(ist, crd));
138  }
139 
141  template< typename Scalar_T, const index_t LO, const index_t HI >
143  framed_multi(const Scalar_T& scr, const index_set_t frm)
144  : map_t(_GLUCAT_HASH_N(1))
145  {
146  if (scr != Scalar_T(0))
147  this->insert(term_t(index_set_t(), scr));
148  }
149 
151  template< typename Scalar_T, const index_t LO, const index_t HI >
153  framed_multi(const int scr, const index_set_t frm)
154  : map_t(_GLUCAT_HASH_N(1))
155  {
156  if (scr != Scalar_T(0))
157  this->insert(term_t(index_set_t(), Scalar_T(scr)));
158  }
159 
161  template< typename Scalar_T, const index_t LO, const index_t HI >
163  framed_multi(const vector_t& vec,
164  const index_set_t frm, const bool prechecked)
165  : map_t(_GLUCAT_HASH_N(vec.size()))
166  {
167  if (!prechecked && index_t(vec.size()) != frm.count())
168  throw error_t("multivector_t(vec,frm): cannot initialize with vector not matching frame");
169  typename vector_t::const_iterator vec_it = vec.begin();
170  const index_t begin_index = frm.min();
171  const index_t end_index = frm.max()+1;
172  for (index_t
173  idx = begin_index;
174  idx != end_index;
175  ++idx)
176  if (frm[idx])
177  {
178  *this += term_t(index_set_t(idx), *vec_it);
179  ++vec_it;
180  }
181  }
182 
184  template< typename Scalar_T, const index_t LO, const index_t HI >
186  framed_multi(const std::string& str)
187  : map_t(_GLUCAT_HASH_N(0))
188  {
189  std::istringstream ss(str);
190  ss >> *this;
191  if (!ss)
192  throw error_t("multivector_t(str): could not parse string");
193  // Peek to see if the end of the string has been reached.
194  ss.peek();
195  if (!ss.eof())
196  throw error_t("multivector_t(str): could not parse entire string");
197  }
198 
200  template< typename Scalar_T, const index_t LO, const index_t HI >
202  framed_multi(const std::string& str, const index_set_t frm, const bool prechecked)
203  : map_t(_GLUCAT_HASH_N(0))
204  {
205  if (prechecked)
206  *this = multivector_t(str);
207  else
208  *this = multivector_t(multivector_t(str), frm, false);
209  }
210 
212  template< typename Scalar_T, const index_t LO, const index_t HI >
213  template< typename Other_Scalar_T >
216  : map_t(_GLUCAT_HASH_N(1))
217  {
218  if (val == Other_Scalar_T(0))
219  return;
220 
221  typedef typename matrix_multi_t::matrix_index_t matrix_index_t;
222  const matrix_index_t dim = val.m_matrix.size1();
223  typedef numeric_traits<Scalar_T> traits_t;
224  if (dim == 1)
225  {
226  this->insert(term_t(index_set_t(), traits_t::to_scalar_t(val.m_matrix(0, 0))));
227  return;
228  }
230  try
231  {
232  *this = val.template fast_framed_multi<Scalar_T>();
233  return;
234  }
235  catch (const glucat_error& e)
236  { }
237 
238  const index_set_t frm = val.frame();
239  const set_value_t algebra_dim = 1 << frm.count();
240  const Scalar_T val_norm = traits_t::to_scalar_t(val.norm());
241  if (traits_t::isNaN_or_isInf(val_norm))
242  {
243  *this = traits_t::NaN();
244  return;
245  }
246  const Scalar_T eps = std::numeric_limits<Scalar_T>::epsilon();
247  const Scalar_T tol = traits_t::abs(val_norm * eps * eps);
248 #if defined(_GLUCAT_MAP_IS_HASH)
249  const size_t max_size = std::min<size_t>(algebra_dim, matrix::nnz(val.m_matrix));
250  *this = multivector_t(_GLUCAT_HASH_SIZE_T(max_size));
251 #endif
252  for (set_value_t
253  stv = 0;
254  stv != algebra_dim;
255  stv++)
256  {
257  const index_set_t ist = index_set_t(stv, frm, true);
258  const Scalar_T crd =
259  traits_t::to_scalar_t(matrix::inner<Other_Scalar_T>(val.basis_element(ist), val.m_matrix));
260  const Scalar_T abs_crd = traits_t::abs(crd);
261  if ((abs_crd * abs_crd) > tol)
262  this->insert(term_t(ist, crd));
263  }
264  }
265 
267  template< typename Scalar_T, const index_t LO, const index_t HI >
268  bool
269  framed_multi<Scalar_T,LO,HI>::
270  operator== (const multivector_t& rhs) const
271  {
272  if (this->size() != rhs.size())
273  return false;
274  const const_iterator this_begin = this->begin();
275  const const_iterator this_end = this->end();
276  const const_iterator rhs_end = rhs.end();
277 #if defined(_GLUCAT_MAP_IS_ORDERED)
278  const_iterator this_it = this_begin;
279  const_iterator rhs_it = rhs.begin();
280  for (;
281  (this_it != this_end) && (rhs_it != rhs_end);
282  this_it++, rhs_it++)
283  if (*this_it != *rhs_it)
284  return false;
285  return (this_it == this_end) && (rhs_it == rhs_end);
286 #else
287  for (const_iterator
288  this_it = this_begin;
289  this_it != this_end;
290  this_it++)
291  {
292  const const_iterator& rhs_it = rhs.find(this_it->first);
293  if (rhs_it == rhs_end)
294  return false;
295  else if (rhs_it->second != this_it->second)
296  return false;
297  }
298  return true;
299 #endif
300  }
301 
303  template< typename Scalar_T, const index_t LO, const index_t HI >
304  inline
305  bool
306  framed_multi<Scalar_T,LO,HI>::
307  operator== (const Scalar_T& scr) const
308  {
309  switch (this->size())
310  {
311  case 0:
312  return scr == Scalar_T(0);
313  case 1:
314  {
315  const const_iterator& this_it = this->begin();
316  return this_it->first == index_set_t() && this_it->second == scr;
317  }
318  default:
319  return false;
320  }
321  }
322 
323 
325  template< typename Scalar_T, const index_t LO, const index_t HI >
326  inline
327  framed_multi<Scalar_T,LO,HI> &
329  operator+= (const Scalar_T& scr)
330  {
331  *this += term_t(index_set_t(),scr);
332  return *this;
333  }
334 
336  template< typename Scalar_T, const index_t LO, const index_t HI >
337  inline
338  framed_multi<Scalar_T,LO,HI> &
340  operator+= (const multivector_t& rhs)
341  { // simply add terms
342  for (const_iterator
343  rhs_it = rhs.begin();
344  rhs_it != rhs.end();
345  ++rhs_it)
346  *this += *rhs_it;
347  return *this;
348  }
349 
351  template< typename Scalar_T, const index_t LO, const index_t HI >
352  inline
353  framed_multi<Scalar_T,LO,HI> &
354  framed_multi<Scalar_T,LO,HI>::
355  operator-= (const multivector_t& rhs)
356  {
357  for (const_iterator
358  rhs_it = rhs.begin();
359  rhs_it != rhs.end();
360  ++rhs_it)
361  *this += term_t(rhs_it->first, -(rhs_it->second));
362  return *this;
363  }
364 
366  template< typename Scalar_T, const index_t LO, const index_t HI >
367  inline
370  operator- () const
371  { return *this * Scalar_T(-1); }
372 
374  template< typename Scalar_T, const index_t LO, const index_t HI >
375  framed_multi<Scalar_T,LO,HI> &
376  framed_multi<Scalar_T,LO,HI>::
377  operator*= (const Scalar_T& scr)
378  { // multiply coordinates of all terms by scalar
379  typedef numeric_traits<Scalar_T> traits_t;
380 
381  if (traits_t::isNaN_or_isInf(scr))
382  return *this = traits_t::NaN();
383  if (scr == Scalar_T(0))
384  if (this->isnan())
385  *this = traits_t::NaN();
386  else
387  this->clear();
388  else
389  for (iterator
390  this_it = this->begin();
391  this_it != this->end();
392  ++this_it)
393  this_it->second *= scr;
394  return *this;
395  }
396 
398  template< typename Scalar_T, const index_t LO, const index_t HI >
399  const framed_multi<Scalar_T,LO,HI>
400  operator* (const framed_multi<Scalar_T,LO,HI>& lhs, const framed_multi<Scalar_T,LO,HI>& rhs)
401  {
402  typedef framed_multi<Scalar_T,LO,HI> multivector_t;
403  typedef numeric_traits<Scalar_T> traits_t;
404  typedef typename multivector_t::index_set_t index_set_t;
405  typedef typename multivector_t::term_t term_t;
406  typedef typename multivector_t::map_t map_t;
407  typedef typename map_t::const_iterator const_iterator;
408 
409  if (lhs.isnan() || rhs.isnan())
410  return traits_t::NaN();
411 
412  const double lhs_size = lhs.size();
413  const double rhs_size = rhs.size();
414  const index_set_t our_frame = lhs.frame() | rhs.frame();
415  const index_t frm_count = our_frame.count();
416  const set_value_t algebra_dim = 1 << frm_count;
417  const bool direct_mult = lhs_size * rhs_size <= double(algebra_dim)
418 #if defined(_GLUCAT_MAP_IS_HASH)
419  || frm_count < Tune_P::mult_matrix_threshold
420 #endif
421  ;
422  if (direct_mult)
423  { // If we have a sparse multiply, store the result directly
424  multivector_t result =
425  multivector_t(_GLUCAT_HASH_SIZE_T(size_t(std::min(lhs_size * rhs_size, double(algebra_dim)))));
426  const const_iterator lhs_begin = lhs.begin();
427  const const_iterator lhs_end = lhs.end();
428  const const_iterator rhs_begin = rhs.begin();
429  const const_iterator rhs_end = rhs.end();
430 
431  for (const_iterator
432  lhs_it = lhs_begin;
433  lhs_it != lhs_end;
434  ++lhs_it)
435  {
436  const term_t& lhs_term = *lhs_it;
437  for (const_iterator
438  rhs_it = rhs_begin;
439  rhs_it != rhs_end;
440  ++rhs_it)
441  result += lhs_term * *rhs_it;
442  }
443  return result;
444  }
445 #if !defined(_GLUCAT_MAP_IS_HASH)
446  else if (frm_count < Tune_P::mult_matrix_threshold)
447  { // Fastest dense algorithm in low dimensions stores result in array
448  typedef std::vector<Scalar_T> array_t;
449  array_t result_array(algebra_dim, Scalar_T(0));
450 
451  const const_iterator lhs_begin = lhs.begin();
452  const const_iterator lhs_end = lhs.end();
453  const const_iterator rhs_begin = rhs.begin();
454  const const_iterator rhs_end = rhs.end();
455 
456  for (const_iterator
457  lhs_it = lhs_begin;
458  lhs_it != lhs_end;
459  ++lhs_it)
460  {
461  const term_t& lhs_term = *lhs_it;
462  for (const_iterator
463  rhs_it = rhs_begin;
464  rhs_it != rhs_end;
465  ++rhs_it)
466  {
467  const term_t& term = lhs_term * *rhs_it;
468  const set_value_t stv = term.first.value_of_fold(our_frame);
469  result_array[stv] += term.second;
470  }
471  }
472  multivector_t result;
473  for (set_value_t
474  stv = 0;
475  stv != algebra_dim;
476  ++stv)
477  if (result_array[stv] != Scalar_T(0))
478  result.insert(term_t(index_set_t(stv, our_frame, true), result_array[stv]));
479  return result;
480  }
481 #endif
482  else
483  { // Past a certain threshold, the matrix algorithm is fastest
484  typedef typename multivector_t::matrix_multi_t matrix_multi_t;
485  return matrix_multi_t(lhs, our_frame, true) *
486  matrix_multi_t(rhs, our_frame, true);
487  }
488  }
489 
491  template< typename Scalar_T, const index_t LO, const index_t HI >
492  inline
493  framed_multi<Scalar_T,LO,HI> &
494  framed_multi<Scalar_T,LO,HI>::
495  operator*= (const multivector_t& rhs)
496  { return *this = *this * rhs; }
497 
499  template< typename Scalar_T, const index_t LO, const index_t HI >
500  const framed_multi<Scalar_T,LO,HI>
501  operator^ (const framed_multi<Scalar_T,LO,HI>& lhs, const framed_multi<Scalar_T,LO,HI>& rhs)
502  { // Arvind Raja's original reference:
503  // "old clical, outerproduct(p,q:pterm):pterm in file compmod.pas"
504 
505  if (lhs.empty() || rhs.empty())
506  return Scalar_T(0);
507 
508  typedef framed_multi<Scalar_T,LO,HI> multivector_t;
509  typedef typename multivector_t::index_set_t index_set_t;
510  typedef typename multivector_t::term_t term_t;
511  typedef typename multivector_t::map_t map_t;
512  typedef typename map_t::const_iterator const_iterator;
513 
514  multivector_t result;
515  const index_set_t empty_set = index_set_t();
516 
517  const const_iterator lhs_begin = lhs.begin();
518  const const_iterator lhs_end = lhs.end();
519  const const_iterator rhs_begin = rhs.begin();
520  const const_iterator rhs_end = rhs.end();
521  const double lhs_size = lhs.size();
522  const double rhs_size = rhs.size();
523 
524  if (lhs_size * rhs_size > double(Tune_P::products_size_threshold))
525  {
526  const index_set_t lhs_frame = lhs.frame();
527  const index_set_t rhs_frame = rhs.frame();
528 
529  const index_set_t our_frame = lhs_frame | rhs_frame;
530  const set_value_t algebra_dim = 1 << our_frame.count();
531  multivector_t result =
532  multivector_t(_GLUCAT_HASH_SIZE_T(size_t(std::min(lhs_size * rhs_size, double(algebra_dim)))));
533  for (set_value_t
534  result_stv = 0;
535  result_stv != algebra_dim;
536  ++result_stv)
537  {
538  const index_set_t result_ist = index_set_t(result_stv, our_frame, true);
539  const index_set_t lhs_result_frame = lhs_frame & result_ist;
540  const set_value_t lhs_result_dim = 1 << lhs_result_frame.count();
541  Scalar_T result_crd = Scalar_T(0);
542  for (set_value_t
543  lhs_stv = 0;
544  lhs_stv != lhs_result_dim;
545  ++lhs_stv)
546  {
547  const index_set_t lhs_ist = index_set_t(lhs_stv, lhs_result_frame, true);
548  const index_set_t rhs_ist = result_ist ^ lhs_ist;
549  if ((rhs_ist | rhs_frame) == rhs_frame)
550  {
551  const const_iterator lhs_it = lhs.find(lhs_ist);
552  if (lhs_it != lhs_end)
553  {
554  const const_iterator rhs_it = rhs.find(rhs_ist);
555  if (rhs_it != rhs_end)
556  result_crd += crd_of_mult(*lhs_it, *rhs_it);
557  }
558  }
559  }
560  if (result_crd != Scalar_T(0))
561  result.insert(term_t(result_ist, result_crd));
562  }
563  return result;
564  }
565  else
566  {
567  multivector_t result;
568  for (const_iterator
569  rhs_it = rhs_begin;
570  rhs_it != rhs_end;
571  ++rhs_it)
572  {
573  const term_t& rhs_term = *rhs_it;
574  const index_set_t rhs_ist = rhs_term.first;
575  for (const_iterator
576  lhs_it = lhs_begin;
577  lhs_it != lhs_end;
578  ++lhs_it)
579  {
580  const term_t& lhs_term = *lhs_it;
581  const index_set_t lhs_ist = lhs_term.first;
582  if ((lhs_ist & rhs_ist) == empty_set)
583  result += lhs_term * rhs_term;
584  }
585  }
586  return result;
587  }
588  }
589 
591  template< typename Scalar_T, const index_t LO, const index_t HI >
592  inline
593  framed_multi<Scalar_T,LO,HI> &
594  framed_multi<Scalar_T,LO,HI>::
595  operator^= (const multivector_t& rhs)
596  { return *this = *this ^ rhs; }
597 
599  template< typename Scalar_T, const index_t LO, const index_t HI >
600  const framed_multi<Scalar_T,LO,HI>
601  operator& (const framed_multi<Scalar_T,LO,HI>& lhs, const framed_multi<Scalar_T,LO,HI>& rhs)
602  { // Arvind Raja's original reference:
603  // "old clical, innerproduct(p,q:pterm):pterm in file compmod.pas"
604 
605  if (lhs.empty() || rhs.empty())
606  return Scalar_T(0);
607 
608  typedef framed_multi<Scalar_T,LO,HI> multivector_t;
609  typedef typename multivector_t::index_set_t index_set_t;
610  typedef typename multivector_t::term_t term_t;
611  typedef typename multivector_t::map_t map_t;
612  typedef typename map_t::const_iterator const_iterator;
613 
614  const const_iterator lhs_end = lhs.end();
615  const const_iterator rhs_end = rhs.end();
616  const double lhs_size = lhs.size();
617  const double rhs_size = rhs.size();
618 
619  if (lhs_size * rhs_size > double(Tune_P::products_size_threshold))
620  {
621  const index_set_t lhs_frame = lhs.frame();
622  const index_set_t rhs_frame = rhs.frame();
623 
624  const index_set_t our_frame = lhs_frame | rhs_frame;
625  const set_value_t algebra_dim = 1 << our_frame.count();
626  multivector_t result =
627  multivector_t(_GLUCAT_HASH_SIZE_T(size_t(std::min(lhs_size * rhs_size, double(algebra_dim)))));
628  for (set_value_t
629  result_stv = 0;
630  result_stv != algebra_dim;
631  ++result_stv)
632  {
633  const index_set_t result_ist = index_set_t(result_stv, our_frame, true);
634  const index_set_t comp_frame = our_frame & ~result_ist;
635  const set_value_t comp_dim = 1 << comp_frame.count();
636  Scalar_T result_crd = Scalar_T(0);
637  for (set_value_t
638  comp_stv = 1;
639  comp_stv != comp_dim;
640  ++comp_stv)
641  {
642  const index_set_t comp_ist = index_set_t(comp_stv, comp_frame, true);
643  const index_set_t our_ist = result_ist ^ comp_ist;
644  if ((our_ist | lhs_frame) == lhs_frame)
645  {
646  const const_iterator lhs_it = lhs.find(our_ist);
647  if (lhs_it != lhs_end)
648  {
649  const const_iterator rhs_it = rhs.find(comp_ist);
650  if (rhs_it != rhs_end)
651  result_crd += crd_of_mult(*lhs_it, *rhs_it);
652  }
653  }
654  if (result_stv != 0)
655  {
656  if ((our_ist | rhs_frame) == rhs_frame)
657  {
658  const const_iterator rhs_it = rhs.find(our_ist);
659  if (rhs_it != rhs_end)
660  {
661  const const_iterator lhs_it = lhs.find(comp_ist);
662  if (lhs_it != lhs_end)
663  result_crd += crd_of_mult(*lhs_it, *rhs_it);
664  }
665  }
666  }
667  }
668  if (result_crd != Scalar_T(0))
669  result.insert(term_t(result_ist, result_crd));
670  }
671  return result;
672  }
673  else
674  {
675  const index_set_t empty_set = index_set_t();
676 
677  const const_iterator lhs_begin = lhs.begin();
678  const const_iterator rhs_begin = rhs.begin();
679 
680  multivector_t result;
681  for (const_iterator
682  lhs_it = lhs_begin;
683  lhs_it != lhs_end;
684  ++lhs_it)
685  {
686  const term_t& lhs_term = *lhs_it;
687  const index_set_t lhs_ist = lhs_term.first;
688  if (lhs_ist != empty_set)
689  for (const_iterator
690  rhs_it = rhs_begin;
691  rhs_it != rhs_end;
692  ++rhs_it)
693  {
694  const term_t& rhs_term = *rhs_it;
695  const index_set_t rhs_ist = rhs_term.first;
696  if (rhs_ist != empty_set)
697  {
698  const index_set_t our_ist = lhs_ist | rhs_ist;
699  if ((lhs_ist == our_ist) || (rhs_ist == our_ist))
700  result += lhs_term * rhs_term;
701  }
702  }
703  }
704  return result;
705  }
706  }
707 
709  template< typename Scalar_T, const index_t LO, const index_t HI >
710  inline
711  framed_multi<Scalar_T,LO,HI> &
712  framed_multi<Scalar_T,LO,HI>::
713  operator&= (const multivector_t& rhs)
714  { return *this = *this & rhs; }
715 
717  template< typename Scalar_T, const index_t LO, const index_t HI >
718  const framed_multi<Scalar_T,LO,HI>
719  operator% (const framed_multi<Scalar_T,LO,HI>& lhs, const framed_multi<Scalar_T,LO,HI>& rhs)
720  {
721  // Reference: Leo Dorst, "Honing geometric algebra for its use in the computer sciences",
722  // in Geometric Computing with Clifford Algebras, ed. G. Sommer,
723  // Springer 2001, Chapter 6, pp. 127-152.
724  // http://staff.science.uva.nl/~leo/clifford/index.html
725 
726  if (lhs.empty() || rhs.empty())
727  return Scalar_T(0);
728 
729  typedef framed_multi<Scalar_T,LO,HI> multivector_t;
730  typedef typename multivector_t::index_set_t index_set_t;
731  typedef typename multivector_t::term_t term_t;
732  typedef typename multivector_t::map_t map_t;
733  typedef typename map_t::const_iterator const_iterator;
734 
735 #if defined(_GLUCAT_MAP_IS_ORDERED)
736  // Both lhs and rhs are sorted by increasing grade, then lexicographically,
737  // and a "larger" index set cannot be a subset of a "smaller" one.
738 
739  const const_iterator lhs_begin = lhs.begin();
740 
741  typedef typename map_t::const_reverse_iterator const_reverse_iterator;
742  const const_reverse_iterator rhs_rbegin = rhs.rbegin();
743  const const_reverse_iterator rhs_rlower_bound =
744  static_cast<const_reverse_iterator>(rhs.lower_bound(lhs_begin->first));
745 
746  multivector_t result;
747 
748  for (const_reverse_iterator
749  rhs_it = rhs_rbegin;
750  rhs_it != rhs_rlower_bound;
751  ++rhs_it)
752  {
753  const term_t& rhs_term = *rhs_it;
754  const index_set_t rhs_ist = rhs_term.first;
755  const const_iterator lhs_upper_bound = lhs.upper_bound(rhs_ist);
756  for (const_iterator
757  lhs_it = lhs_begin;
758  lhs_it != lhs_upper_bound;
759  ++lhs_it)
760  {
761  const term_t& lhs_term = *lhs_it;
762  const index_set_t lhs_ist = lhs_term.first;
763  if ((lhs_ist | rhs_ist) == rhs_ist)
764  result += lhs_term * rhs_term;
765  }
766  }
767  return result;
768 #else
769  const const_iterator lhs_end = lhs.end();
770  const const_iterator rhs_end = rhs.end();
771  const double lhs_size = lhs.size();
772  const double rhs_size = rhs.size();
773 
774  if (lhs_size * rhs_size > double(Tune_P::products_size_threshold))
775  {
776  const index_set_t lhs_frame = lhs.frame();
777  const index_set_t rhs_frame = rhs.frame();
778 
779  const index_set_t our_frame = lhs_frame | rhs_frame;
780  const set_value_t algebra_dim = 1 << our_frame.count();
781  multivector_t result =
782  multivector_t(_GLUCAT_HASH_SIZE_T(size_t(std::min(lhs_size * rhs_size, double(algebra_dim)))));
783  for (set_value_t
784  result_stv = 0;
785  result_stv != algebra_dim;
786  ++result_stv)
787  {
788  const index_set_t result_ist = index_set_t(result_stv, our_frame, true);
789  const index_set_t comp_frame = lhs_frame & ~result_ist;
790  const set_value_t comp_dim = 1 << comp_frame.count();
791  Scalar_T result_crd = Scalar_T(0);
792  for (set_value_t
793  comp_stv = 0;
794  comp_stv != comp_dim;
795  ++comp_stv)
796  {
797  const index_set_t comp_ist = index_set_t(comp_stv, comp_frame, true);
798  const index_set_t rhs_ist = result_ist ^ comp_ist;
799  if ((rhs_ist | rhs_frame) == rhs_frame)
800  {
801  const const_iterator rhs_it = rhs.find(rhs_ist);
802  if (rhs_it != rhs_end)
803  {
804  const const_iterator lhs_it = lhs.find(comp_ist);
805  if (lhs_it != lhs_end)
806  result_crd += crd_of_mult(*lhs_it, *rhs_it);
807  }
808  }
809  }
810  if (result_crd != Scalar_T(0))
811  result.insert(term_t(result_ist, result_crd));
812  }
813  return result;
814  }
815  else
816  {
817  const const_iterator rhs_begin = rhs.begin();
818  const const_iterator lhs_begin = lhs.begin();
819 
820  multivector_t result;
821  for (const_iterator
822  rhs_it = rhs_begin;
823  rhs_it != rhs_end;
824  ++rhs_it)
825  {
826  const term_t& rhs_term = *rhs_it;
827  const index_set_t rhs_ist = rhs_term.first;
828  for (const_iterator
829  lhs_it = lhs_begin;
830  lhs_it != lhs_end;
831  ++lhs_it)
832  {
833  const term_t& lhs_term = *lhs_it;
834  const index_set_t lhs_ist = lhs_term.first;
835  if ((lhs_ist | rhs_ist) == rhs_ist)
836  result += lhs_term * rhs_term;
837  }
838  }
839  return result;
840  }
841 #endif
842  }
843 
845  template< typename Scalar_T, const index_t LO, const index_t HI >
846  inline
847  framed_multi<Scalar_T,LO,HI> &
848  framed_multi<Scalar_T,LO,HI>::
849  operator%= (const multivector_t& rhs)
850  { return *this = *this % rhs; }
851 
853  template< typename Scalar_T, const index_t LO, const index_t HI >
854  Scalar_T
855  star(const framed_multi<Scalar_T,LO,HI>& lhs, const framed_multi<Scalar_T,LO,HI>& rhs)
856  {
857  typedef framed_multi<Scalar_T,LO,HI> multivector_t;
858  typedef typename multivector_t::map_t map_t;
859  typedef typename map_t::const_iterator const_iterator;
860  typedef typename multivector_t::index_set_t index_set_t;
861 
862  Scalar_T result = Scalar_T(0);
863  const bool small_star_large = lhs.size() < rhs.size();
864  const multivector_t* smallp =
865  small_star_large
866  ? &lhs
867  : &rhs;
868  const multivector_t* largep =
869  small_star_large
870  ? &rhs
871  : &lhs;
872 
873  for (const_iterator
874  small_it = smallp->begin();
875  small_it != smallp->end();
876  ++small_it)
877  {
878  const index_set_t small_ist = small_it->first;
879  const Scalar_T large_crd = (*largep)[small_ist];
880  if (large_crd != Scalar_T(0))
881  result += small_ist.sign_of_square() * small_it->second * large_crd;
882  }
883  return result;
884  }
885 
887  template< typename Scalar_T, const index_t LO, const index_t HI >
890  operator/= (const Scalar_T& scr)
891  { // Divide coordinates of all terms by scr
892  typedef numeric_traits<Scalar_T> traits_t;
893 
894  if (traits_t::isNaN(scr))
895  return *this = traits_t::NaN();
896  if (traits_t::isInf(scr))
897  if (this->isnan())
898  *this = traits_t::NaN();
899  else
900  this->clear();
901  else
902  for (iterator
903  this_it = this->begin();
904  this_it != this->end();
905  ++this_it)
906  this_it->second /= scr;
907  return *this;
908  }
909 
911  template< typename Scalar_T, const index_t LO, const index_t HI >
912  inline
913  const framed_multi<Scalar_T,LO,HI>
914  operator/ (const framed_multi<Scalar_T,LO,HI>& lhs, const framed_multi<Scalar_T,LO,HI>& rhs)
915  {
916  typedef framed_multi<Scalar_T,LO,HI> multivector_t;
917  typedef numeric_traits<Scalar_T> traits_t;
918  typedef typename multivector_t::index_set_t index_set_t;
919  typedef typename multivector_t::matrix_multi_t matrix_multi_t;
920 
921  if (rhs == Scalar_T(0))
922  return traits_t::NaN();
923 
924  const index_set_t our_frame = lhs.frame() | rhs.frame();
925  return matrix_multi_t(lhs, our_frame, true) / matrix_multi_t(rhs, our_frame, true);
926  }
927 
929  template< typename Scalar_T, const index_t LO, const index_t HI >
930  inline
931  framed_multi<Scalar_T,LO,HI> &
932  framed_multi<Scalar_T,LO,HI>::
933  operator/= (const multivector_t& rhs)
934  { return *this = *this / rhs; }
935 
937  template< typename Scalar_T, const index_t LO, const index_t HI >
938  inline
939  const framed_multi<Scalar_T,LO,HI>
940  operator| (const framed_multi<Scalar_T,LO,HI>& lhs, const framed_multi<Scalar_T,LO,HI>& rhs)
941  {
942  typedef framed_multi<Scalar_T,LO,HI> multivector_t;
943  typedef typename multivector_t::matrix_multi_t matrix_multi_t;
944 
945  return matrix_multi_t(rhs) * matrix_multi_t(lhs) / matrix_multi_t(rhs.involute());
946  }
947 
949  template< typename Scalar_T, const index_t LO, const index_t HI >
950  inline
951  framed_multi<Scalar_T,LO,HI>&
952  framed_multi<Scalar_T,LO,HI>::
953  operator|= (const multivector_t& rhs)
954  { return *this = *this | rhs; }
955 
957  template< typename Scalar_T, const index_t LO, const index_t HI >
958  inline
959  const framed_multi<Scalar_T,LO,HI>
961  inv() const
962  {
963  matrix_multi_t result = matrix_multi_t(Scalar_T(1), this->frame());
964  return result /= matrix_multi_t(*this);
965  }
966 
968  template< typename Scalar_T, const index_t LO, const index_t HI >
969  const framed_multi<Scalar_T,LO,HI>
971  pow(int m) const
972  { return glucat::pow(*this, m); }
973 
975  template< typename Scalar_T, const index_t LO, const index_t HI >
976  const
977  framed_multi<Scalar_T,LO,HI>
979  outer_pow(int m) const
980  {
981  if (m < 0)
982  throw error_t("outer_pow(int): negative exponent");
983  multivector_t result = Scalar_T(1);
984  multivector_t a = *this;
985  for (;
986  m != 0;
987  m >>= 1, a ^= a)
988  if (m & 1)
989  result ^= a;
990  return result;
991  }
992 
994  template< typename Scalar_T, const index_t LO, const index_t HI >
995  inline
996  const index_set<LO,HI>
997  framed_multi<Scalar_T,LO,HI>::
998  frame() const
999  {
1000  index_set_t result;
1001  for (const_iterator
1002  this_it = this->begin();
1003  this_it != this->end();
1004  ++this_it)
1005  result |= this_it->first;
1006  return result;
1007  }
1008 
1010  template< typename Scalar_T, const index_t LO, const index_t HI >
1011  inline
1012  index_t
1013  framed_multi<Scalar_T,LO,HI>::
1014  grade() const
1015  {
1016  index_t result = 0;
1017  for (const_iterator
1018  this_it = this->begin();
1019  this_it != this->end();
1020  ++this_it)
1021  result = std::max( result, this_it->first.count() );
1022  return result;
1023  }
1024 
1026  template< typename Scalar_T, const index_t LO, const index_t HI >
1027  inline
1028  Scalar_T
1029  framed_multi<Scalar_T,LO,HI>::
1030  operator[] (const index_set_t ist) const
1031  {
1032  const const_iterator& this_it = this->find(ist);
1033  if (this_it == this->end())
1034  return Scalar_T(0);
1035  else
1036  return this_it->second;
1037  }
1038 
1040  template< typename Scalar_T, const index_t LO, const index_t HI >
1041  const framed_multi<Scalar_T,LO,HI>
1042  framed_multi<Scalar_T,LO,HI>::
1043  operator() (index_t grade) const
1044  {
1045  if ((grade < 0) || (grade > HI-LO))
1046  return Scalar_T(0);
1047  else
1048  {
1049  multivector_t result;
1050  for (const_iterator
1051  this_it = this->begin();
1052  this_it != this->end();
1053  ++this_it)
1054  if (this_it->first.count() == grade)
1055  result += *this_it;
1056  return result;
1057  }
1058  }
1059 
1061  template< typename Scalar_T, const index_t LO, const index_t HI >
1062  inline
1063  Scalar_T
1065  scalar() const
1066  { return (*this)[index_set_t()]; }
1067 
1069  template< typename Scalar_T, const index_t LO, const index_t HI >
1070  inline
1071  const framed_multi<Scalar_T,LO,HI>
1073  pure() const
1074  { return *this - this->scalar(); }
1075 
1077  template< typename Scalar_T, const index_t LO, const index_t HI >
1078  const framed_multi<Scalar_T,LO,HI>
1080  even() const
1081  { // even part of x, sum of the pure(count) with even count
1082  multivector_t result;
1083  for (const_iterator
1084  this_it = this->begin();
1085  this_it != this->end();
1086  ++this_it)
1087  if ((this_it->first.count() % 2) == 0)
1088  result.insert(*this_it);
1089  return result;
1090  }
1091 
1093  template< typename Scalar_T, const index_t LO, const index_t HI >
1094  const framed_multi<Scalar_T,LO,HI>
1096  odd() const
1097  { // even part of x, sum of the pure(count) with even count
1098  multivector_t result;
1099  for (const_iterator
1100  this_it = this->begin();
1101  this_it != this->end();
1102  ++this_it)
1103  if ((this_it->first.count() % 2) == 1)
1104  result.insert(*this_it);
1105  return result;
1106  }
1107 
1109  template< typename Scalar_T, const index_t LO, const index_t HI >
1112  vector_part() const
1113  { return this->vector_part(this->frame(), true); }
1114 
1116  template< typename Scalar_T, const index_t LO, const index_t HI >
1119  vector_part(const index_set_t frm, const bool prechecked) const
1120  {
1121  if (!prechecked && (this->frame() | frm) != frm)
1122  throw error_t("vector_part(frm): value is outside of requested frame");
1123  vector_t result;
1124  result.reserve(frm.count());
1125  const index_t frm_end = frm.max()+1;
1126  for (index_t
1127  idx = frm.min();
1128  idx != frm_end;
1129  ++idx)
1130  // Frame may contain indices which do not correspond to a grade 1 term but
1131  // frame cannot omit any index corresponding to a grade 1 term
1132  if (frm[idx])
1133  result.push_back((*this)[index_set_t(idx)]);
1134  return result;
1135  }
1136 
1138  template< typename Scalar_T, const index_t LO, const index_t HI >
1139  const framed_multi<Scalar_T,LO,HI>
1141  involute() const
1142  {
1143  multivector_t result = *this;
1144  for (iterator
1145  result_it = result.begin();
1146  result_it != result.end();
1147  ++result_it)
1148  { // for a k-vector u, involute(u) == (-1)^k * u
1149  if ((result_it->first.count() % 2) == 1)
1150  result_it->second *= Scalar_T(-1);
1151  }
1152  return result;
1153  }
1154 
1156  template< typename Scalar_T, const index_t LO, const index_t HI >
1157  const framed_multi<Scalar_T,LO,HI>
1159  reverse() const
1160  {
1161  multivector_t result = *this;
1162  for (iterator
1163  result_it = result.begin();
1164  result_it != result.end();
1165  ++result_it)
1166  // For a k-vector u, reverse(u) = { -u, k == 2,3 (mod 4)
1167  // { u, k == 0,1 (mod 4)
1168  switch (result_it->first.count() % 4)
1169  {
1170  case 2:
1171  case 3:
1172  result_it->second *= Scalar_T(-1);
1173  break;
1174  default:
1175  break;
1176  }
1177  return result;
1178  }
1179 
1181  template< typename Scalar_T, const index_t LO, const index_t HI >
1182  const framed_multi<Scalar_T,LO,HI>
1184  conj() const
1185  {
1186  multivector_t result = *this;
1187  for (iterator
1188  result_it = result.begin();
1189  result_it != result.end();
1190  ++result_it)
1191  // For a k-vector u, conj(u) = { -u, k == 1,2 (mod 4)
1192  // { u, k == 0,3 (mod 4)
1193  switch (result_it->first.count() % 4)
1194  {
1195  case 1:
1196  case 2:
1197  result_it->second *= Scalar_T(-1);
1198  break;
1199  default:
1200  break;
1201  }
1202  return result;
1203  }
1204 
1206  template< typename Scalar_T, const index_t LO, const index_t HI >
1207  Scalar_T
1209  quad() const
1210  {
1211  // scalar(conj(x)*x) = 2*quad(even(x)) - quad(x)
1212  // ref: old clical: quadfunction(p:pter):pterm in file compmod.pas
1213  Scalar_T result = Scalar_T(0);
1214  for (const_iterator
1215  this_it = this->begin();
1216  this_it != this->end();
1217  ++this_it)
1218  {
1219  const Scalar_T sign =
1220  (this_it->first.count_neg() % 2)
1221  ? -Scalar_T(1)
1222  : Scalar_T(1);
1223  result += sign * (this_it->second) * (this_it->second);
1224  }
1225  return result;
1226  }
1227 
1229  template< typename Scalar_T, const index_t LO, const index_t HI >
1230  Scalar_T
1232  norm() const
1233  {
1234  typedef numeric_traits<Scalar_T> traits_t;
1235 
1236  Scalar_T result = Scalar_T(0);
1237  for (const_iterator
1238  this_it = this->begin();
1239  this_it != this->end();
1240  ++this_it)
1241  {
1242  const Scalar_T abs_crd = traits_t::abs(this_it->second);
1243  result += abs_crd * abs_crd;
1244  }
1245  return result;
1246  }
1247 
1249  template< typename Scalar_T, const index_t LO, const index_t HI >
1250  Scalar_T
1252  max_abs() const
1253  {
1254  typedef numeric_traits<Scalar_T> traits_t;
1255 
1256  Scalar_T result = Scalar_T(0);
1257  for (const_iterator
1258  this_it = this->begin();
1259  this_it != this->end();
1260  ++this_it)
1261  {
1262  const Scalar_T abs_crd = traits_t::abs(this_it->second);
1263  if (abs_crd > result)
1264  result = abs_crd;
1265  }
1266  return result;
1267  }
1268 
1270  template< typename Scalar_T, const index_t LO, const index_t HI >
1271  const framed_multi<Scalar_T,LO,HI>
1273  random(const index_set<LO,HI> frm, Scalar_T fill)
1274  {
1275  typedef framed_multi<Scalar_T,LO,HI> multivector_t;
1276  typedef typename multivector_t::index_set_t index_set_t;
1277  typedef typename multivector_t::term_t term_t;
1278 
1279  typedef random_generator<Scalar_T> random_generator_t;
1280  random_generator_t& generator = random_generator_t::generator();
1281 
1282  fill =
1283  (fill < Scalar_T(0))
1284  ? Scalar_T(0)
1285  : (fill > Scalar_T(1))
1286  ? Scalar_T(1)
1287  : fill;
1288  const set_value_t algebra_dim = 1 << frm.count();
1289  typedef numeric_traits<Scalar_T> traits_t;
1290  const Scalar_T mean_abs = traits_t::sqrt(Scalar_T(double(algebra_dim)));
1291  multivector_t result;
1292  for (set_value_t
1293  stv = 0;
1294  stv != algebra_dim;
1295  ++stv)
1296  if (generator.uniform() < fill)
1297  {
1298  const Scalar_T& result_crd = generator.normal() / mean_abs;
1299  result.insert(term_t(index_set_t(stv, frm, true), result_crd));
1300  }
1301  return result;
1302  }
1303 
1305  template< typename Scalar_T, const index_t LO, const index_t HI >
1306  inline
1307  void
1309  write(const std::string& msg) const
1310  { std::cout << msg << std::endl << " " << (*this) << std::endl; }
1311 
1313  template< typename Scalar_T, const index_t LO, const index_t HI >
1314  inline
1315  void
1316  framed_multi<Scalar_T,LO,HI>::
1317  write(std::ofstream& ofile, const std::string& msg) const
1318  {
1319  if (!ofile)
1320  throw error_t("write(ofile,msg): cannot write to output file");
1321  ofile << msg << std::endl << " " << (*this) << std::endl;
1322  }
1323 
1325  template< typename Map_T,typename Sorted_Map_T >
1326  class sorted_range
1327  {
1328  public:
1329  typedef Map_T map_t;
1330  typedef Sorted_Map_T sorted_map_t;
1331  typedef typename Sorted_Map_T::const_iterator sorted_iterator;
1332 
1333  sorted_range (Sorted_Map_T &sorted_val, const Map_T& val)
1334  {
1335  for (typename map_t::const_iterator
1336  val_it = val.begin();
1337  val_it != val.end();
1338  ++val_it)
1339  sorted_val.insert(*val_it);
1340  sorted_begin = sorted_val.begin();
1341  sorted_end = sorted_val.end();
1342  }
1345  };
1346 
1347  template< typename Sorted_Map_T >
1348  class sorted_range< Sorted_Map_T, Sorted_Map_T >
1349  {
1350  public:
1351  typedef Sorted_Map_T map_t;
1352  typedef Sorted_Map_T sorted_map_t;
1353  typedef typename Sorted_Map_T::const_iterator sorted_iterator;
1354 
1355  sorted_range (Sorted_Map_T &sorted_val, const Sorted_Map_T& val)
1356  : sorted_begin( val.begin() ),
1357  sorted_end( val.end() )
1358  { }
1361  };
1362 
1364  template< typename Scalar_T, const index_t LO, const index_t HI >
1365  std::ostream&
1366  operator<< (std::ostream& os, const framed_multi<Scalar_T,LO,HI>& val)
1367  {
1368  if (val.empty())
1369  os << 0;
1370  else if (val.isnan())
1371  os << std::numeric_limits<Scalar_T>::quiet_NaN();
1372  else
1373  {
1374  typedef framed_multi<Scalar_T,LO,HI> multivector_t;
1375  typedef typename multivector_t::map_t map_t;
1376  typedef typename multivector_t::sorted_map_t sorted_map_t;
1377  typedef typename sorted_map_t::const_iterator sorted_iterator;
1378  sorted_map_t sorted_val;
1379  sorted_range< map_t, sorted_map_t > sorted_val_range(sorted_val, val);
1380  sorted_iterator sorted_it = sorted_val_range.sorted_begin;
1381  os << *sorted_it;
1382  for (++sorted_it;
1383  sorted_it != sorted_val_range.sorted_end;
1384  ++sorted_it)
1385  {
1386  const Scalar_T& scr = sorted_it->second;
1387  if (scr >= 0.0)
1388  os << '+';
1389  os << *sorted_it;
1390  }
1391  }
1392  return os;
1393  }
1394 
1396  template< typename Scalar_T, const index_t LO, const index_t HI >
1397  std::ostream&
1398  operator<< (std::ostream& os, const std::pair< const index_set<LO,HI>, Scalar_T >& term)
1399  {
1400  const double second_as_double = numeric_traits<Scalar_T>::to_double(term.second);
1401  const bool use_double =
1402  (os.precision() <= std::numeric_limits<double>::digits10) ||
1403  (term.second == Scalar_T(second_as_double));
1404  if (term.first.count() == 0)
1405  if (use_double)
1406  os << second_as_double;
1407  else
1408  os << term.second;
1409  else if (term.second == Scalar_T(-1))
1410  {
1411  os << '-';
1412  os << term.first;
1413  }
1414  else if (term.second != Scalar_T(1))
1415  {
1416  if (use_double)
1417  {
1418  double tol = std::pow(10.0,-os.precision());
1419  if ( std::fabs(second_as_double + 1.0) < tol )
1420  os << '-';
1421  else if ( std::fabs(second_as_double - 1.0) >= tol )
1422  os << second_as_double;
1423  }
1424  else
1425  os << term.second;
1426  os << term.first;
1427  }
1428  else
1429  os << term.first;
1430  return os;
1431  }
1432 
1434  template< typename Scalar_T, const index_t LO, const index_t HI >
1435  std::istream&
1436  operator>> (std::istream& s, framed_multi<Scalar_T,LO,HI> & val)
1437  { // Input looks like 1.0-2.0{1,2}+3.2{3,4}.
1438  typedef framed_multi<Scalar_T,LO,HI> multivector_t;
1439  // Parsing variables.
1440  multivector_t local_val;
1441  int c = 0;
1442  // Parsing control variables.
1443  bool negative = false;
1444  bool expect_term = true;
1445  // The multivector may begin with '+' or '-'. Check for this.
1446  c = s.peek();
1447  if (s.good() && (c == int('+') || c == int('-')))
1448  { // A '-' here negates the following term.
1449  negative = (c == int('-'));
1450  // Consume the '+' or '-'.
1451  s.get();
1452  }
1453  while (s.good())
1454  { // Parse a term.
1455  // A term consists of an optional scalar, followed by an optional index set.
1456  // At least one of the two must be present.
1457  // Default coordinate is Scalar_T(1).
1458  Scalar_T coordinate = Scalar_T(1);
1459  // Default index set is empty.
1460  index_set<LO,HI> ist;
1461  // First, check for an opening brace.
1462  c = s.peek();
1463  if (s.good())
1464  { // If the character is not an opening brace,
1465  // a coordinate value is expected here.
1466  if (c != int('{'))
1467  { // Try to read a coordinate value.
1468  double coordinate_as_double;
1469  s >> coordinate_as_double;
1470  // Reading the coordinate may have resulted in an end of file condition.
1471  // This is not a failure.
1472  if (s)
1473  coordinate = Scalar_T(coordinate_as_double);
1474  }
1475  }
1476  else
1477  { // End of file here ends parsing while a term may still be expected.
1478  break;
1479  }
1480  // Coordinate is now Scalar_T(1) or a Scalar_T value.
1481  // Parse an optional index set.
1482  if (s.good())
1483  {
1484  c = s.peek();
1485  if (s.good() && c == int('{'))
1486  { // Try to read index set.
1487  s >> ist;
1488  }
1489  }
1490  // Reading the term may have resulted in an end of file condition.
1491  // This is not a failure.
1492  if (s)
1493  {
1494  // Immediately after parsing a term, another term is not expected.
1495  expect_term = false;
1496  if (coordinate != Scalar_T(0))
1497  {
1498  // Add the term to the local multivector.
1499  coordinate =
1500  negative
1501  ? -coordinate
1502  : coordinate;
1503  typedef typename multivector_t::term_t term_t;
1504  local_val += term_t(ist, coordinate);
1505  }
1506  }
1507  // Check if anything follows the current term.
1508  if (s.good())
1509  {
1510  c = s.peek();
1511  if (s.good())
1512  { // Only '+' and '-' are valid here.
1513  if (c == int('+') || c == int('-'))
1514  { // A '-' here negates the following term.
1515  negative = (c == int('-'));
1516  // Consume the '+' or '-'.
1517  s.get();
1518  // Immediately after '+' or '-',
1519  // expect another term.
1520  expect_term = true;
1521  }
1522  else
1523  { // Any other character here is a not failure,
1524  // but still ends the parsing of the multivector.
1525  break;
1526  }
1527  }
1528  }
1529  }
1530  // If a term is still expected, this is a failure.
1531  if (expect_term)
1532  s.clear(std::istream::failbit);
1533  // End of file is not a failure.
1534  if (s)
1535  { // The multivector has been successfully parsed.
1536  val = local_val;
1537  }
1538  return s;
1539  }
1540 
1542  template< typename Scalar_T, const index_t LO, const index_t HI >
1543  unsigned long
1545  nbr_terms () const
1546  { return this->size(); }
1547 
1549  template< typename Scalar_T, const index_t LO, const index_t HI >
1550  inline
1551  framed_multi<Scalar_T,LO,HI> &
1553  operator+= (const term_t& term)
1554  { // Do not insert terms with 0 coordinate
1555  if (term.second != Scalar_T(0))
1556  {
1557  const iterator& this_it = this->find(term.first);
1558  if (this_it == this->end())
1559  this->insert(term);
1560  else if (this_it->second + term.second == Scalar_T(0))
1561  // Erase term if resulting coordinate is 0
1562  this->erase(this_it);
1563  else
1564  this_it->second += term.second;
1565  }
1566  return *this;
1567  }
1568 
1570  template< typename Scalar_T, const index_t LO, const index_t HI >
1571  bool
1573  isnan() const
1574  {
1575  typedef numeric_traits<Scalar_T> traits_t;
1576 
1577  if (std::numeric_limits<Scalar_T>::has_quiet_NaN)
1578  for (const_iterator
1579  this_it = this->begin();
1580  this_it != this->end();
1581  ++this_it)
1582  if (traits_t::isNaN(this_it->second))
1583  return true;
1584  return false;
1585  }
1586 
1588  template< typename Scalar_T, const index_t LO, const index_t HI >
1591  truncated(const Scalar_T& limit) const
1592  {
1593  typedef numeric_traits<Scalar_T> traits_t;
1594 
1595  const Scalar_T abs_limit = traits_t::abs(limit);
1596  if (this->isnan())
1597  return *this;
1598  Scalar_T top = max_abs();
1599  multivector_t result;
1600  if (top != Scalar_T(0))
1601  for (const_iterator
1602  this_it = this->begin();
1603  this_it != this->end();
1604  ++this_it)
1605  if (traits_t::abs(this_it->second / top) > abs_limit)
1606  result.insert(term_t(this_it->first, this_it->second));
1607  return result;
1608  }
1609 
1611  template< typename Scalar_T, const index_t LO, const index_t HI >
1612  framed_multi<Scalar_T,LO,HI>
1614  fold(const index_set_t frm) const
1615  {
1616  if (frm.is_contiguous())
1617  return *this;
1618  else
1619  {
1620  multivector_t result;
1621  for (const_iterator
1622  this_it = this->begin();
1623  this_it != this->end();
1624  ++this_it)
1625  result.insert(term_t(this_it->first.fold(frm), this_it->second));
1626  return result;
1627  }
1628  }
1629 
1631  template< typename Scalar_T, const index_t LO, const index_t HI >
1632  framed_multi<Scalar_T,LO,HI>
1634  unfold(const index_set_t frm) const
1635  {
1636  if (frm.is_contiguous())
1637  return *this;
1638  else
1639  {
1640  multivector_t result;
1641  for (const_iterator
1642  this_it = this->begin();
1643  this_it != this->end();
1644  ++this_it)
1645  result.insert(term_t(this_it->first.unfold(frm, true), this_it->second));
1646  return result;
1647  }
1648  }
1649 
1651  // Reference: [L] 16.4 Periodicity of 8, p216
1652  template< typename Scalar_T, const index_t LO, const index_t HI >
1653  framed_multi<Scalar_T,LO,HI>&
1656  {
1657  // We add 4 to q by subtracting 4 from p
1658  if (q+4 > -LO)
1659  throw error_t("centre_pm4_qp4(p,q): LO is too high to represent this value");
1660  if (this->frame().max() > p-4)
1661  {
1662  typedef typename index_set_t::index_pair_t index_pair_t;
1663  const index_set_t pm3210(index_pair_t(p-3,p), true);
1664  const index_set_t qm4321(index_pair_t(-q-4,-q-1), true);
1665  const term_t& tqm4321 = term_t(qm4321, Scalar_T(1));
1666  multivector_t result;
1667  for (const_iterator
1668  this_it = this->begin();
1669  this_it != this->end();
1670  ++this_it)
1671  {
1672  index_set_t ist = this_it->first;
1673  if (ist.max() > p-4)
1674  {
1675  var_term_t term;
1676  for (index_t
1677  n = 0;
1678  n != 4;
1679  ++n)
1680  if (ist[n+p-3])
1681  term *= term_t(index_set_t(n-q-4), Scalar_T(1)) * tqm4321;
1682  // Mask out {p-3}..{p}
1683  result.insert(term_t(ist & ~pm3210, this_it->second) *
1684  term_t(term.first, term.second));
1685  }
1686  else
1687  result.insert(*this_it);
1688  }
1689  *this = result;
1690  }
1691  p -=4; q += 4;
1692  return *this;
1693  }
1694 
1696  // Reference: [L] 16.4 Periodicity of 8, p216
1697  template< typename Scalar_T, const index_t LO, const index_t HI >
1698  framed_multi<Scalar_T,LO,HI>&
1701  {
1702  // We add 4 to p by subtracting 4 from q
1703  if (p+4 > HI)
1704  throw error_t("centre_pp4_qm4(p,q): HI is too low to represent this value");
1705  if (this->frame().min() < -q+4)
1706  {
1707  typedef typename index_set_t::index_pair_t index_pair_t;
1708  const index_set_t qp0123(index_pair_t(-q,-q+3), true);
1709  const index_set_t pp1234(index_pair_t(p+1,p+4), true);
1710  const term_t& tpp1234 = term_t(pp1234, Scalar_T(1));
1711  multivector_t result;
1712  for (const_iterator
1713  this_it = this->begin();
1714  this_it != this->end();
1715  ++this_it)
1716  {
1717  index_set_t ist = this_it->first;
1718  if (ist.min() < -q+4)
1719  {
1720  var_term_t term;
1721  for (index_t
1722  n = 0;
1723  n != 4;
1724  ++n)
1725  if (ist[n-q])
1726  term *= term_t(index_set_t(n+p+1), Scalar_T(1)) * tpp1234;
1727  // Mask out {-q}..{-q+3}
1728  result.insert(term_t(term.first, term.second) *
1729  term_t(ist & ~qp0123, this_it->second));
1730  }
1731  else
1732  result.insert(*this_it);
1733  }
1734  *this = result;
1735  }
1736  p +=4; q -= 4;
1737  return *this;
1738  }
1739 
1741  // Reference: [P] Proposition 15.20, p 131
1742  template< typename Scalar_T, const index_t LO, const index_t HI >
1743  framed_multi<Scalar_T,LO,HI>&
1746  {
1747  if (q+1 > HI)
1748  throw error_t("centre_qp1_pm1(p,q): HI is too low to represent this value");
1749  if (p-1 > -LO)
1750  throw error_t("centre_qp1_pm1(p,q): LO is too high to represent this value");
1751  const index_set_t qp1 = index_set_t(q+1);
1752  const term_t& tqp1 = term_t(qp1, Scalar_T(1));
1753  multivector_t result;
1754  for (const_iterator
1755  this_it = this->begin();
1756  this_it != this->end();
1757  ++this_it)
1758  {
1759  const index_set_t ist = this_it->first;
1760  var_term_t term = var_term_t(index_set_t(), this_it->second);
1761  for (index_t
1762  n = -q;
1763  n != p;
1764  ++n)
1765  if (n != 0 && ist[n])
1766  term *= term_t(index_set_t(-n) | qp1, Scalar_T(1));
1767  if (p != 0 && ist[p])
1768  term *= tqp1;
1769  result.insert(term_t(term.first, term.second));
1770  }
1771  index_t orig_p = p;
1772  p = q+1;
1773  q = orig_p-1;
1774  return *this = result;
1775  }
1776 
1778  template< typename Scalar_T, const index_t LO, const index_t HI >
1779  const std::pair< const framed_multi<Scalar_T,LO,HI>, const framed_multi<Scalar_T,LO,HI> >
1781  divide(const index_set_t ist) const
1782  {
1783  multivector_t quo;
1784  multivector_t rem;
1785  for (const_iterator
1786  this_it = this->begin();
1787  this_it != this->end();
1788  ++this_it)
1789  if ((this_it->first | ist) == this_it->first)
1790  quo.insert(term_t(this_it->first ^ ist, this_it->second));
1791  else
1792  rem.insert(*this_it);
1793  return framed_pair_t(quo, rem);
1794  }
1795 
1797  template< typename Scalar_T, const index_t LO, const index_t HI >
1800  fast(const index_t level, const bool odd) const
1801  {
1802  // Assume val is already folded and centred
1803  if (this->empty())
1804  {
1805  typedef typename matrix_multi_t::matrix_index_t matrix_index_t;
1806  const matrix_index_t dim = 1 << level;
1807  matrix_t result(dim, dim);
1808  result.clear();
1809  return result;
1810  }
1811  if (level == 0)
1812  return matrix::unit<matrix_t>(1) * this->scalar();
1813 
1814  typedef typename matrix_multi_t::basis_matrix_t basis_matrix_t;
1815  typedef typename basis_matrix_t::value_type basis_scalar_t;
1816 
1817  const basis_matrix_t& I = matrix::unit<basis_matrix_t>(2);
1818  basis_matrix_t J(2,2,2);
1819  J.clear();
1820  J(0,1) = basis_scalar_t(-1);
1821  J(1,0) = basis_scalar_t( 1);
1822  basis_matrix_t K = J;
1823  K(0,1) = basis_scalar_t( 1);
1824  basis_matrix_t JK = I;
1825  JK(0,0) = basis_scalar_t(-1);
1826 
1827  const index_set_t ist_mn = index_set_t(-level);
1828  const index_set_t ist_pn = index_set_t(level);
1829  if (level == 1)
1830  {
1831  if (odd)
1832  return matrix_t(J) * (*this)[ist_mn] + matrix_t(K) * (*this)[ist_pn];
1833  else
1834  return matrix_t(I) * this->scalar() + matrix_t(JK) * (*this)[ist_mn ^ ist_pn];
1835  }
1836  else
1837  {
1838  const framed_pair_t& pair_mn = this->divide(ist_mn);
1839  const framed_multi_t& quo_mn = pair_mn.first;
1840  const framed_multi_t& rem_mn = pair_mn.second;
1841  const framed_pair_t& pair_quo_mnpn = quo_mn.divide(ist_pn);
1842  const framed_multi_t& val_mnpn = pair_quo_mnpn.first;
1843  const framed_multi_t& val_mn = pair_quo_mnpn.second;
1844  const framed_pair_t& pair_rem_mnpn = rem_mn.divide(ist_pn);
1845  const framed_multi_t& val_pn = pair_rem_mnpn.first;
1846  const framed_multi_t& val_1 = pair_rem_mnpn.second;
1847  using matrix::kron;
1848  if (odd)
1849  return - kron(JK, val_1.fast (level-1, 1))
1850  + kron(I, val_mnpn.fast(level-1, 1))
1851  + kron(J, val_mn.fast (level-1, 0))
1852  + kron(K, val_pn.fast (level-1, 0));
1853  else
1854  return kron(I, val_1.fast (level-1, 0))
1855  + kron(JK, val_mnpn.fast(level-1, 0))
1856  + kron(K, val_mn.fast (level-1, 1))
1857  - kron(J, val_pn.fast (level-1, 1));
1858  }
1859  }
1860 
1862  template< typename Scalar_T, const index_t LO, const index_t HI >
1863  template< typename Other_Scalar_T >
1864  const matrix_multi<Other_Scalar_T,LO,HI>
1866  fast_matrix_multi(const index_set_t frm) const
1867  {
1868  // Fold val
1869  multivector_t val = this->fold(frm);
1870  index_t p = frm.count_pos();
1871  index_t q = frm.count_neg();
1872  const index_t bott_offset = gen::offset_to_super[pos_mod(p - q, 8)];
1873  p += std::max(bott_offset,index_t(0));
1874  q -= std::min(bott_offset,index_t(0));
1875  if (p > HI)
1876  throw error_t("fast_matrix_multi(frm): HI is too low to represent this value");
1877  if (q > -LO)
1878  throw error_t("fast_matrix_multi(frm): LO is too high to represent this value");
1879  // Centre val
1880  while (p - q > 4)
1881  val.centre_pm4_qp4(p, q);
1882  while (p - q < -3)
1883  val.centre_pp4_qm4(p, q);
1884  if (p - q > 1)
1885  val.centre_qp1_pm1(p, q);
1886  const index_t level = (p + q)/2;
1887 
1888  // Do the fast transform
1889  const multivector_t& ev_val = val.even();
1890  const multivector_t& od_val = val.odd();
1891  return matrix_multi<Other_Scalar_T,LO,HI>(ev_val.fast(level, 0) + od_val.fast(level, 1), frm);
1892  }
1893 
1894  template< typename Scalar_T, const index_t LO, const index_t HI >
1895  inline
1898  fast_framed_multi() const
1899  { return *this; }
1900 
1902  template< typename Scalar_T, const index_t LO, const index_t HI >
1903  inline
1904  static
1905  Scalar_T
1906  crd_of_mult(const std::pair<const index_set<LO,HI>, Scalar_T>& lhs,
1907  const std::pair<const index_set<LO,HI>, Scalar_T>& rhs)
1908  { return lhs.first.sign_of_mult(rhs.first) * lhs.second * rhs.second; }
1909 
1911  template< typename Scalar_T, const index_t LO, const index_t HI >
1912  inline
1913  const std::pair<const index_set<LO,HI>, Scalar_T>
1914  operator* (const std::pair<const index_set<LO,HI>, Scalar_T>& lhs,
1915  const std::pair<const index_set<LO,HI>, Scalar_T>& rhs)
1916  {
1917  typedef std::pair<const index_set<LO,HI>, Scalar_T> term_t;
1918  return term_t(lhs.first ^ rhs.first, crd_of_mult(lhs, rhs));
1919  }
1920 
1922  template< typename Scalar_T, const index_t LO, const index_t HI >
1923  const framed_multi<Scalar_T,LO,HI>
1924  sqrt(const framed_multi<Scalar_T,LO,HI>& val, const framed_multi<Scalar_T,LO,HI>& i, bool prechecked)
1925  {
1926  typedef numeric_traits<Scalar_T> traits_t;
1927  if (val.isnan())
1928  return traits_t::NaN();
1929 
1930  check_complex(val, i, prechecked);
1931 
1932  const Scalar_T realval = val.scalar();
1933  if (val == realval)
1934  {
1935  if (realval < Scalar_T(0))
1936  return i * traits_t::sqrt(-realval);
1937  else
1938  return traits_t::sqrt(realval);
1939  }
1940  typedef typename framed_multi<Scalar_T,LO,HI>::matrix_multi_t matrix_multi_t;
1941  return sqrt(matrix_multi_t(val), matrix_multi_t(i), prechecked);
1942  }
1943 
1945  template< typename Scalar_T, const index_t LO, const index_t HI >
1946  const framed_multi<Scalar_T,LO,HI>
1947  exp(const framed_multi<Scalar_T,LO,HI>& val)
1948  {
1949  typedef numeric_traits<Scalar_T> traits_t;
1950  if (val.isnan())
1951  return traits_t::NaN();
1952 
1953  const Scalar_T s = scalar(val);
1954  if (val == s)
1955  return traits_t::exp(s);
1956 
1957  const double size = val.size();
1958  const index_t frm_count = val.frame().count();
1959  const set_value_t algebra_dim = 1 << frm_count;
1960 
1961  if( (size * size <= double(algebra_dim)) || (frm_count < Tune_P::mult_matrix_threshold))
1962  {
1964  {
1965  case precision_demoted:
1966  {
1967  typedef typename traits_t::demoted::type demoted_scalar_t;
1968  typedef framed_multi<demoted_scalar_t,LO,HI> demoted_multivector_t;
1969 
1970  const demoted_multivector_t& demoted_val = demoted_multivector_t(val);
1971  return clifford_exp(demoted_val);
1972  }
1973  break;
1974  case precision_promoted:
1975  {
1976  typedef typename traits_t::promoted::type promoted_scalar_t;
1977  typedef framed_multi<promoted_scalar_t,LO,HI> promoted_multivector_t;
1978 
1979  const promoted_multivector_t& promoted_val = promoted_multivector_t(val);
1980  return clifford_exp(promoted_val);
1981  }
1982  break;
1983  default:
1984  return clifford_exp(val);
1985  }
1986  }
1987  else
1988  {
1989  typedef matrix_multi<Scalar_T,LO,HI> matrix_multi_t;
1990  return exp(matrix_multi_t(val));
1991  }
1992  }
1993 
1995  template< typename Scalar_T, const index_t LO, const index_t HI >
1996  const framed_multi<Scalar_T,LO,HI>
1997  log(const framed_multi<Scalar_T,LO,HI>& val, const framed_multi<Scalar_T,LO,HI>& i, bool prechecked)
1998  {
1999  typedef numeric_traits<Scalar_T> traits_t;
2000  if (val == Scalar_T(0) || val.isnan())
2001  return traits_t::NaN();
2002 
2003  check_complex(val, i, prechecked);
2004 
2005  const Scalar_T realval = val.scalar();
2006  if (val == realval)
2007  {
2008  if (realval < Scalar_T(0))
2009  return i * traits_t::pi() + traits_t::log(-realval);
2010  else
2011  return traits_t::log(realval);
2012  }
2013  typedef typename framed_multi<Scalar_T,LO,HI>::matrix_multi_t matrix_multi_t;
2014  return log(matrix_multi_t(val), matrix_multi_t(i), prechecked);
2015  }
2016 }
2017 #endif // _GLUCAT_FRAMED_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::framed_multi::var_term_t
class var_term var_term_t
Definition: framed_multi.h:164
glucat::precision_demoted
@ precision_demoted
Definition: global.h:148
glucat::framed_multi::index_set_t
index_set< LO, HI > index_set_t
Definition: framed_multi.h:152
PyClical.e
def e(obj)
Definition: PyClical.pyx:1886
glucat::scalar
Scalar_T scalar(const Multivector< Scalar_T, LO, HI > &val)
Scalar part.
Definition: clifford_algebra_imp.h:421
glucat::index_set::min
index_t min() const
Minimum member.
Definition: index_set_imp.h:490
glucat::framed_multi::centre_qp1_pm1
multivector_t & centre_qp1_pm1(index_t &p, index_t &q)
Subalgebra isomorphism: R_{p,q} to R_{q+1,p-1}.
Definition: framed_multi_imp.h:1775
glucat::framed_multi::multivector_t
framed_multi multivector_t
Definition: framed_multi.h:149
glucat::sorted_range::map_t
Map_T map_t
Definition: framed_multi_imp.h:1359
glucat::sorted_range::sorted_end
sorted_iterator sorted_end
Definition: framed_multi_imp.h:1374
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::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
glucat::framed_multi::matrix_t
matrix_multi_t::matrix_t matrix_t
Definition: framed_multi.h:165
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::tuning::products_size_threshold
@ products_size_threshold
Definition: global.h:206
glucat::framed_multi::classname
static const std::string classname()
Class name used in messages.
Definition: framed_multi_imp.h:83
glucat::index_t
int index_t
Size of index_t should be enough to represent LO, HI.
Definition: global.h:106
_GLUCAT_HASH_N
#define _GLUCAT_HASH_N(x)
Definition: framed_multi_imp.h:90
glucat::set_value_t
unsigned long set_value_t
Size of set_value_t should be enough to contain index_set<LO,HI>
Definition: global.h:108
glucat::sorted_range::sorted_map_t
Sorted_Map_T sorted_map_t
Definition: framed_multi_imp.h:1360
glucat::sorted_range< Sorted_Map_T, Sorted_Map_T >::sorted_iterator
Sorted_Map_T::const_iterator sorted_iterator
Definition: framed_multi_imp.h:1383
glucat::framed_multi::vector_t
std::vector< Scalar_T > vector_t
Definition: framed_multi.h:154
glucat::precision_promoted
@ precision_promoted
Definition: global.h:150
glucat::index_set
Index set class based on std::bitset<> in Gnu standard C++ library.
Definition: index_set.h:74
random.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::even
const Multivector< Scalar_T, LO, HI > even(const Multivector< Scalar_T, LO, HI > &val)
Even part.
Definition: clifford_algebra_imp.h:456
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::tuning::inv_fast_dim_threshold
@ inv_fast_dim_threshold
Definition: global.h:203
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::framed_multi::map_t
std::unordered_map< index_set_t, Scalar_T, index_set_hash< LO, HI > > map_t
Definition: framed_multi.h:175
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::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::framed_multi::centre_pp4_qm4
multivector_t & centre_pp4_qm4(index_t &p, index_t &q)
Subalgebra isomorphism: R_{p,q} to R_{p+4,q-4}.
Definition: framed_multi_imp.h:1730
glucat::crd_of_mult
static Scalar_T crd_of_mult(const std::pair< const index_set< LO, HI >, Scalar_T > &lhs, const std::pair< const index_set< LO, HI >, Scalar_T > &rhs)
Coordinate of product of terms.
glucat::sorted_range::sorted_iterator
Sorted_Map_T::const_iterator sorted_iterator
Definition: framed_multi_imp.h:1361
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::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
framed_multi.h
PyClical.fill
fill
Definition: PyClical.pyx:1814
glucat::sorted_range
Sorted range for use with output.
Definition: framed_multi_imp.h:1357
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::framed_multi::term_t
std::pair< const index_set_t, Scalar_T > term_t
Definition: framed_multi.h:153
glucat::framed_multi::unfold
multivector_t unfold(const index_set_t frm) const
Subalgebra isomorphism: unfold each term within the given frame.
Definition: framed_multi_imp.h:1664
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::framed_multi::fold
multivector_t fold(const index_set_t frm) const
Subalgebra isomorphism: fold each term within the given frame.
Definition: framed_multi_imp.h:1644
glucat::matrix::kron
const RHS_T kron(const LHS_T &lhs, const RHS_T &rhs)
Kronecker tensor product of matrices - as per Matlab kron.
Definition: matrix_imp.h:135
glucat::matrix_multi
A matrix_multi<Scalar_T,LO,HI> is a matrix approximation to a multivector.
Definition: framed_multi.h:68
glucat::matrix::nnz
Matrix_T::size_type nnz(const Matrix_T &m)
Number of non-zeros.
Definition: matrix_imp.h:331
glucat::framed_multi
A framed_multi<Scalar_T,LO,HI> is a framed approximation to a multivector.
Definition: framed_multi.h:65
glucat::framed_multi::framed_multi
friend class framed_multi
Definition: framed_multi.h:160
PyClical.pi
float pi
Definition: PyClical.pyx:1857
glucat::sorted_range::sorted_range
sorted_range(Sorted_Map_T &sorted_val, const Map_T &val)
Definition: framed_multi_imp.h:1363
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_HASH_SIZE_T
#define _GLUCAT_HASH_SIZE_T(x)
Definition: framed_multi_imp.h:91
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::tuning::mult_matrix_threshold
@ mult_matrix_threshold
Definition: global.h:184
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::sorted_range::sorted_begin
sorted_iterator sorted_begin
Definition: framed_multi_imp.h:1373
glucat::numeric_traits
Extra traits which extend numeric limits.
Definition: scalar.h:76
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::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::numeric_traits::to_double
static double to_double(const Scalar_T &val)
Cast to double.
Definition: scalar.h:190
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
Definition: clifford_algebra.h:39
glucat::framed_multi::error_t
error< multivector_t > error_t
Definition: framed_multi.h:155
glucat::framed_multi::matrix_multi_t
matrix_multi< Scalar_T, LO, HI > matrix_multi_t
Definition: framed_multi.h:156
glucat::matrix_multi::matrix_index_t
matrix_t::size_type matrix_index_t
Definition: matrix_multi.h:188
glucat::framed_multi::fast
const matrix_t fast(const index_t level, const bool odd) const
Generalized FFT from framed_multi_t to matrix_t.
Definition: framed_multi_imp.h:1830
glucat::framed_multi::fast_framed_multi
const framed_multi_t fast_framed_multi() const
Use inverse generalized FFT to construct a framed_multi_t.
Definition: framed_multi_imp.h:1928
glucat::framed_multi::fast_matrix_multi
const matrix_multi< Other_Scalar_T, LO, HI > fast_matrix_multi(const index_set_t frm) const
Use generalized FFT to construct a matrix_multi_t.
Definition: framed_multi_imp.h:1896
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::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::sorted_range< Sorted_Map_T, Sorted_Map_T >::sorted_map_t
Sorted_Map_T sorted_map_t
Definition: framed_multi_imp.h:1382
glucat::framed_multi::nbr_terms
_GLUCAT_CLIFFORD_ALGEBRA_OPERATIONS unsigned long nbr_terms() const
Number of terms.
Definition: framed_multi_imp.h:1575
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::framed_multi::operator+=
multivector_t & operator+=(const term_t &term)
Add a term, if non-zero.
Definition: framed_multi_imp.h:359
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
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::framed_multi::centre_pm4_qp4
multivector_t & centre_pm4_qp4(index_t &p, index_t &q)
Subalgebra isomorphism: R_{p,q} to R_{p-4,q+4}.
Definition: framed_multi_imp.h:1685
glucat::framed_multi::divide
const framed_pair_t divide(const index_set_t ist) const
Divide multivector into part divisible by index_set and remainder.
Definition: framed_multi_imp.h:1811
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