glucat  0.12.0
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-2021 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/scalar.h"
38 #include "glucat/random.h"
39 #include "glucat/generation.h"
40 #include "glucat/matrix.h"
41 
42 #include <sstream>
43 #include <fstream>
44 
45 namespace glucat
46 {
48  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
49  auto
51  classname() -> const std::string
52  { return "framed_multi"; }
53 
54 #define _GLUCAT_HASH_N(x) (x)
55 #define _GLUCAT_HASH_SIZE_T(x) (typename multivector_t::hash_size_t)(x)
56 
58  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
61  : map_t(_GLUCAT_HASH_N(0))
62  { }
63 
65  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
67  framed_multi(const hash_size_t& hash_size)
68  : map_t(_GLUCAT_HASH_N(hash_size()))
69  { }
70 
72  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
73  template< typename Other_Scalar_T >
76  : map_t(_GLUCAT_HASH_N(val.size()))
77  {
78  for (auto& val_term : val)
79  this->insert(term_t(val_term.first, numeric_traits<Scalar_T>::to_scalar_t(val_term.second)));
80  }
81 
83  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
84  template< typename Other_Scalar_T >
87  const index_set_t frm, const bool prechecked)
88  : map_t(_GLUCAT_HASH_N(val.size()))
89  {
90  if (!prechecked && (val.frame() | frm) != frm)
91  throw error_t("multivector_t(val,frm): cannot initialize with value outside of frame");
92  for (auto& val_term : val)
93  this->insert(term_t(val_term.first, numeric_traits<Scalar_T>::to_scalar_t(val_term.second)));
94  }
95 
97  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
100  const index_set_t frm, const bool prechecked)
101  : map_t(_GLUCAT_HASH_N(val.size()))
102  {
103  if (!prechecked && (val.frame() | frm) != frm)
104  throw error_t("multivector_t(val,frm): cannot initialize with value outside of frame");
105  for (auto& val_term : val)
106  this->insert(val_term);
107  }
108 
110  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
112  framed_multi(const index_set_t ist, const Scalar_T& crd)
113  : map_t(_GLUCAT_HASH_N(1))
114  {
115  if (crd != Scalar_T(0))
116  this->insert(term_t(ist, crd));
117  }
118 
120  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
122  framed_multi(const index_set_t ist, const Scalar_T& crd,
123  const index_set_t frm, const bool prechecked)
124  : map_t(_GLUCAT_HASH_N(1))
125  {
126  if (!prechecked && (ist | frm) != frm)
127  throw error_t("multivector_t(ist,crd,frm): cannot initialize with value outside of frame");
128  if (crd != Scalar_T(0))
129  this->insert(term_t(ist, crd));
130  }
131 
133  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
135  framed_multi(const Scalar_T& scr, const index_set_t frm)
136  : map_t(_GLUCAT_HASH_N(1))
137  {
138  if (scr != Scalar_T(0))
139  this->insert(term_t(index_set_t(), scr));
140  }
141 
143  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
145  framed_multi(const int scr, const index_set_t frm)
146  : map_t(_GLUCAT_HASH_N(1))
147  {
148  if (scr != Scalar_T(0))
149  this->insert(term_t(index_set_t(), Scalar_T(scr)));
150  }
151 
153  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
156  const index_set_t frm, const bool prechecked)
157  : map_t(_GLUCAT_HASH_N(vec.size()))
158  {
159  if (!prechecked && index_t(vec.size()) != frm.count())
160  throw error_t("multivector_t(vec,frm): cannot initialize with vector not matching frame");
161  auto idx = frm.min();
162  const auto frm_end = frm.max()+1;
163  for (auto& crd : vec)
164  {
165  *this += term_t(index_set_t(idx), crd);
166  for (
167  ++idx;
168  idx != frm_end && !frm[idx];
169  ++idx)
170  ;
171  }
172  }
173 
175  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
177  framed_multi(const std::string& str)
178  : map_t(_GLUCAT_HASH_N(0))
179  {
180  std::istringstream ss(str);
181  ss >> *this;
182  if (!ss)
183  throw error_t("multivector_t(str): could not parse string");
184  // Peek to see if the end of the string has been reached.
185  ss.peek();
186  if (!ss.eof())
187  throw error_t("multivector_t(str): could not parse entire string");
188  }
189 
191  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
193  framed_multi(const std::string& str, const index_set_t frm, const bool prechecked)
194  : map_t(_GLUCAT_HASH_N(0))
195  {
196  if (prechecked)
197  *this = multivector_t(str);
198  else
199  *this = multivector_t(multivector_t(str), frm, false);
200  }
201 
203  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
204  template< typename Other_Scalar_T >
207  : map_t(_GLUCAT_HASH_N(1))
208  {
209  if (val == Other_Scalar_T(0))
210  return;
211 
212  const auto dim = val.m_matrix.size1();
213  using traits_t = numeric_traits<Scalar_T>;
214  if (dim == 1)
215  {
216  this->insert(term_t(index_set_t(), traits_t::to_scalar_t(val.m_matrix(0, 0))));
217  return;
218  }
219  if (dim >= Tune_P::inv_fast_dim_threshold)
220  try
221  {
222  *this = (val.template fast_framed_multi<Scalar_T>()).truncated();
223  return;
224  }
225  catch (const glucat_error& e)
226  { }
227 
228  const auto val_norm = traits_t::to_scalar_t(val.norm());
229  if (traits_t::isNaN_or_isInf(val_norm))
230  {
231  *this = traits_t::NaN();
232  return;
233  }
234  const auto frm = val.frame();
235  const auto algebra_dim = set_value_t(1) << frm.count();
236  auto result = multivector_t(
237  _GLUCAT_HASH_SIZE_T(std::min<size_t>(algebra_dim, matrix::nnz(val.m_matrix))));
238  for (auto
239  stv = set_value_t(0);
240  stv != algebra_dim;
241  stv++)
242  {
243  const auto ist = index_set_t(stv, frm, true);
244  const auto crd =
245  traits_t::to_scalar_t(matrix::inner<Other_Scalar_T>(val.basis_element(ist), val.m_matrix));
246  if (crd != Scalar_T(0))
247  result.insert(term_t(ist, crd));
248  }
249  *this = result.truncated();
250  }
251 
253  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
254  auto
256  operator== (const multivector_t& rhs) const -> bool
257  {
258  if (this->size() != rhs.size())
259  return false;
260  const auto rhs_end = rhs.end();
261  for (auto& this_term : *this)
262  {
263  const const_iterator& rhs_it = rhs.find(this_term.first);
264  if (rhs_it == rhs_end || rhs_it->second != this_term.second)
265  return false;
266  }
267  return true;
268  }
269 
271  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
272  inline
273  auto
275  operator== (const Scalar_T& scr) const -> bool
276  {
277  switch (this->size())
278  {
279  case 0:
280  return scr == Scalar_T(0);
281  case 1:
282  {
283  const auto& this_it = this->begin();
284  return this_it->first == index_set_t() && this_it->second == scr;
285  }
286  default:
287  return false;
288  }
289  }
290 
292  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
293  inline
294  auto
296  operator+= (const Scalar_T& scr) -> multivector_t&
297  {
298  *this += term_t(index_set_t(), scr);
299  return *this;
300  }
301 
303  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
304  inline
305  auto
307  operator+= (const multivector_t& rhs) -> multivector_t&
308  { // simply add terms
309  for (auto& rhs_term : rhs)
310  *this += rhs_term;
311  return *this;
312  }
313 
315  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
316  inline
317  auto
319  operator-= (const Scalar_T& scr) -> multivector_t&
320  {
321  *this += term_t(index_set_t(), -scr);
322  return *this;
323  }
324 
326  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
327  inline
328  auto
330  operator-= (const multivector_t& rhs) -> multivector_t&
331  {
332  for (auto& rhs_term : rhs)
333  *this += term_t(rhs_term.first, -(rhs_term.second));
334  return *this;
335  }
336 
338  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
339  inline
340  auto
342  operator- () const -> const multivector_t
343  { // multiply coordinates of all terms by -1
344  auto result = *this;
345  for (auto& result_term : result)
346  result_term.second *= Scalar_T(-1);
347  return result;
348  }
349 
351  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
352  auto
354  operator*= (const Scalar_T& scr) -> multivector_t&
355  { // multiply coordinates of all terms by scalar
356  using traits_t = numeric_traits<Scalar_T>;
357 
358  if (traits_t::isNaN_or_isInf(scr))
359  return *this = traits_t::NaN();
360  if (scr == Scalar_T(0))
361  if (this->isnan())
362  *this = traits_t::NaN();
363  else
364  this->clear();
365  else
366  for (auto& this_term : *this)
367  this_term.second *= scr;
368  return *this;
369  }
370 
372  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
373  auto
375  {
376  using multivector_t = framed_multi<Scalar_T,LO,HI,Tune_P>;
377  using traits_t = numeric_traits<Scalar_T>;
378 
379  if (lhs.isnan() || rhs.isnan())
380  return traits_t::NaN();
381 
382  const double lhs_size = lhs.size();
383  const double rhs_size = rhs.size();
384  const auto our_frame = lhs.frame() | rhs.frame();
385  const auto frm_count = our_frame.count();
386  const auto algebra_dim = set_value_t(1) << frm_count;
387  const auto direct_mult = lhs_size * rhs_size <= double(algebra_dim);
388  if (direct_mult)
389  { // If we have a sparse multiply, store the result directly
390  auto result = multivector_t(
391  _GLUCAT_HASH_SIZE_T(size_t(std::min(lhs_size * rhs_size, double(algebra_dim)))));
392  for (auto& lhs_term : lhs)
393  for (auto& rhs_term : rhs)
394  result += lhs_term * rhs_term;
395  return result;
396  }
397  else
398  { // Past a certain threshold, the matrix algorithm is fastest
399  using matrix_multi_t = typename multivector_t::matrix_multi_t;
400  return matrix_multi_t(lhs, our_frame, true) *
401  matrix_multi_t(rhs, our_frame, true);
402  }
403  }
404 
406  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
407  inline
408  auto
410  operator*= (const multivector_t& rhs) -> multivector_t&
411  { return *this = *this * rhs; }
412 
414  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
415  auto
417  { // Arvind Raja's original reference:
418  // "old clical, outerproduct(p,q:pterm):pterm in file compmod.pas"
419 
420  if (lhs.empty() || rhs.empty())
421  return Scalar_T(0);
422 
423  using multivector_t = framed_multi<Scalar_T,LO,HI,Tune_P>;
424  using index_set_t = typename multivector_t::index_set_t;
425  using term_t = typename multivector_t::term_t;
426 
427  const auto empty_set = index_set_t();
428 
429  const double lhs_size = lhs.size();
430  const double rhs_size = rhs.size();
431  const auto lhs_frame = lhs.frame();
432  const auto rhs_frame = rhs.frame();
433  const auto our_frame = lhs_frame | rhs_frame;
434  const auto algebra_dim = set_value_t(1) << our_frame.count();
435  auto result = multivector_t(
436  _GLUCAT_HASH_SIZE_T(size_t(std::min(lhs_size * rhs_size, double(algebra_dim)))));
437  const auto lhs_end = lhs.end();
438  const auto rhs_end = rhs.end();
439 
440  if (lhs_size * rhs_size > double(Tune_P::products_size_threshold))
441  {
442  for (auto
443  result_stv = set_value_t(0);
444  result_stv != algebra_dim;
445  ++result_stv)
446  {
447  const auto result_ist = index_set_t(result_stv, our_frame, true);
448  const auto lhs_result_frame = lhs_frame & result_ist;
449  const auto lhs_result_dim = set_value_t(1) << lhs_result_frame.count();
450  auto result_crd = Scalar_T(0);
451  for (auto
452  lhs_stv = set_value_t(0);
453  lhs_stv != lhs_result_dim;
454  ++lhs_stv)
455  {
456  const auto lhs_ist = index_set_t(lhs_stv, lhs_result_frame, true);
457  const auto rhs_ist = result_ist ^ lhs_ist;
458  if ((rhs_ist | rhs_frame) == rhs_frame)
459  {
460  const auto lhs_it = lhs.find(lhs_ist);
461  if (lhs_it != lhs_end)
462  {
463  const auto rhs_it = rhs.find(rhs_ist);
464  if (rhs_it != rhs_end)
465  result_crd += crd_of_mult(*lhs_it, *rhs_it);
466  }
467  }
468  }
469  if (result_crd != Scalar_T(0))
470  result.insert(term_t(result_ist, result_crd));
471  }
472  return result;
473  }
474  else
475  {
476  for (auto& lhs_term : lhs)
477  for (auto& rhs_term : rhs)
478  if ((lhs_term.first & rhs_term.first) == empty_set)
479  result += lhs_term * rhs_term;
480  return result;
481  }
482  }
483 
485  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
486  inline
487  auto
489  operator^= (const multivector_t& rhs) -> multivector_t&
490  { return *this = *this ^ rhs; }
491 
493  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
494  auto
496  { // Arvind Raja's original reference:
497  // "old clical, innerproduct(p,q:pterm):pterm in file compmod.pas"
498 
499  if (lhs.empty() || rhs.empty())
500  return Scalar_T(0);
501 
502  using multivector_t = framed_multi<Scalar_T,LO,HI,Tune_P>;
503  using index_set_t = typename multivector_t::index_set_t;
504  using term_t = typename multivector_t::term_t;
505 
506  const auto lhs_end = lhs.end();
507  const auto rhs_end = rhs.end();
508  const double lhs_size = lhs.size();
509  const double rhs_size = rhs.size();
510 
511  const auto lhs_frame = lhs.frame();
512  const auto rhs_frame = rhs.frame();
513 
514  const auto our_frame = lhs_frame | rhs_frame;
515  const auto algebra_dim = set_value_t(1) << our_frame.count();
516  auto result = multivector_t(
517  _GLUCAT_HASH_SIZE_T(size_t(std::min(lhs_size * rhs_size, double(algebra_dim)))));
518  if (lhs_size * rhs_size > double(Tune_P::products_size_threshold))
519  {
520  for (auto
521  result_stv = set_value_t(0);
522  result_stv != algebra_dim;
523  ++result_stv)
524  {
525  const auto result_ist = index_set_t(result_stv, our_frame, true);
526  const auto comp_frame = our_frame & ~result_ist;
527  const auto comp_dim = set_value_t(1) << comp_frame.count();
528  auto result_crd = Scalar_T(0);
529  for (auto
530  comp_stv = set_value_t(1);
531  comp_stv != comp_dim;
532  ++comp_stv)
533  {
534  const auto comp_ist = index_set_t(comp_stv, comp_frame, true);
535  const auto our_ist = result_ist ^ comp_ist;
536  if ((our_ist | lhs_frame) == lhs_frame)
537  {
538  const auto lhs_it = lhs.find(our_ist);
539  if (lhs_it != lhs_end)
540  {
541  const auto rhs_it = rhs.find(comp_ist);
542  if (rhs_it != rhs_end)
543  result_crd += crd_of_mult(*lhs_it, *rhs_it);
544  }
545  }
546  if (result_stv != 0)
547  {
548  if ((our_ist | rhs_frame) == rhs_frame)
549  {
550  const auto rhs_it = rhs.find(our_ist);
551  if (rhs_it != rhs_end)
552  {
553  const auto lhs_it = lhs.find(comp_ist);
554  if (lhs_it != lhs_end)
555  result_crd += crd_of_mult(*lhs_it, *rhs_it);
556  }
557  }
558  }
559  }
560  if (result_crd != Scalar_T(0))
561  result.insert(term_t(result_ist, result_crd));
562  }
563  }
564  else
565  {
566  const auto empty_set = index_set_t();
567  for (auto& lhs_term : lhs)
568  {
569  const auto lhs_ist = lhs_term.first;
570  if (lhs_ist != empty_set)
571  for (auto& rhs_term : rhs)
572  {
573  const auto rhs_ist = rhs_term.first;
574  if (rhs_ist != empty_set)
575  {
576  const auto our_ist = lhs_ist | rhs_ist;
577  if ((lhs_ist == our_ist) || (rhs_ist == our_ist))
578  result += lhs_term * rhs_term;
579  }
580  }
581  }
582  }
583  return result;
584  }
585 
587  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
588  inline
589  auto
591  operator&= (const multivector_t& rhs) -> multivector_t&
592  { return *this = *this & rhs; }
593 
595  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
596  auto
598  {
599  // Reference: Leo Dorst, "Honing geometric algebra for its use in the computer sciences",
600  // in Geometric Computing with Clifford Algebras, ed. G. Sommer,
601  // Springer 2001, Chapter 6, pp. 127-152.
602  // http://staff.science.uva.nl/~leo/clifford/index.html
603 
604  if (lhs.empty() || rhs.empty())
605  return Scalar_T(0);
606 
607  using multivector_t = framed_multi<Scalar_T,LO,HI,Tune_P>;
608  using index_set_t = typename multivector_t::index_set_t;
609  using term_t = typename multivector_t::term_t;
610  using map_t = typename multivector_t::map_t;
611 
612  const auto lhs_end = lhs.end();
613  const auto rhs_end = rhs.end();
614  const double lhs_size = lhs.size();
615  const double rhs_size = rhs.size();
616  const auto lhs_frame = lhs.frame();
617  const auto rhs_frame = rhs.frame();
618 
619  const auto our_frame = lhs_frame | rhs_frame;
620  const auto algebra_dim = set_value_t(1) << our_frame.count();
621  auto result = multivector_t(
622  _GLUCAT_HASH_SIZE_T(size_t(std::min(lhs_size * rhs_size, double(algebra_dim)))));
623 
624  if (lhs_size * rhs_size > double(Tune_P::products_size_threshold))
625  {
626  for (auto
627  result_stv = set_value_t(0);
628  result_stv != algebra_dim;
629  ++result_stv)
630  {
631  const auto result_ist = index_set_t(result_stv, our_frame, true);
632  const auto comp_frame = lhs_frame & ~result_ist;
633  const auto comp_dim = set_value_t(1) << comp_frame.count();
634  auto result_crd = Scalar_T(0);
635  for (auto
636  comp_stv = set_value_t(0);
637  comp_stv != comp_dim;
638  ++comp_stv)
639  {
640  const auto comp_ist = index_set_t(comp_stv, comp_frame, true);
641  const auto rhs_ist = result_ist ^ comp_ist;
642  if ((rhs_ist | rhs_frame) == rhs_frame)
643  {
644  const auto rhs_it = rhs.find(rhs_ist);
645  if (rhs_it != rhs_end)
646  {
647  const auto lhs_it = lhs.find(comp_ist);
648  if (lhs_it != lhs_end)
649  result_crd += crd_of_mult(*lhs_it, *rhs_it);
650  }
651  }
652  }
653  if (result_crd != Scalar_T(0))
654  result.insert(term_t(result_ist, result_crd));
655  }
656  }
657  else
658  {
659  for (auto& rhs_term : rhs)
660  {
661  const auto rhs_ist = rhs_term.first;
662  for (auto& lhs_term : lhs)
663  {
664  const index_set_t lhs_ist = lhs_term.first;
665  if ((lhs_ist | rhs_ist) == rhs_ist)
666  result += lhs_term * rhs_term;
667  }
668  }
669  }
670  return result;
671  }
672 
674  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
675  inline
676  auto
678  operator%= (const multivector_t& rhs) -> multivector_t&
679  { return *this = *this % rhs; }
680 
682  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
683  auto
685  {
686  using multivector_t = framed_multi<Scalar_T,LO,HI,Tune_P>;
687 
688  auto result = Scalar_T(0);
689  const auto small_star_large = lhs.size() < rhs.size();
690  const auto* smallp =
691  small_star_large
692  ? &lhs
693  : &rhs;
694  const auto* largep =
695  small_star_large
696  ? &rhs
697  : &lhs;
698 
699  for (auto& small_term : *smallp)
700  {
701  const auto small_ist = small_term.first;
702  const auto large_crd = (*largep)[small_ist];
703  if (large_crd != Scalar_T(0))
704  result += small_ist.sign_of_square() * small_term.second * large_crd;
705  }
706  return result;
707  }
708 
710  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
711  auto
713  operator/= (const Scalar_T& scr) -> multivector_t&
714  { // Divide coordinates of all terms by scr
715  using traits_t = numeric_traits<Scalar_T>;
716 
717  if (traits_t::isNaN(scr))
718  return *this = traits_t::NaN();
719  if (traits_t::isInf(scr))
720  if (this->isnan())
721  *this = traits_t::NaN();
722  else
723  this->clear();
724  else
725  for (auto& this_term : *this)
726  this_term.second /= scr;
727  return *this;
728  }
729 
731  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
732  inline
733  auto
735  {
736  using multivector_t = framed_multi<Scalar_T,LO,HI,Tune_P>;
737  using traits_t = numeric_traits<Scalar_T>;
738  using index_set_t = typename multivector_t::index_set_t;
739  using matrix_multi_t = typename multivector_t::matrix_multi_t;
740 
741  if (rhs == Scalar_T(0))
742  return traits_t::NaN();
743 
744  const auto our_frame = lhs.frame() | rhs.frame();
745  return matrix_multi_t(lhs, our_frame, true) / matrix_multi_t(rhs, our_frame, true);
746  }
747 
749  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
750  inline
751  auto
753  operator/= (const multivector_t& rhs) -> multivector_t&
754  { return *this = *this / rhs; }
755 
757  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
758  inline
759  auto
761  {
762  using multivector_t = framed_multi<Scalar_T,LO,HI,Tune_P>;
763  using matrix_multi_t = typename multivector_t::matrix_multi_t;
764 
765  return matrix_multi_t(rhs) * matrix_multi_t(lhs) / matrix_multi_t(rhs.involute());
766  }
767 
769  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
770  inline
771  auto
773  operator|= (const multivector_t& rhs) -> multivector_t&
774  { return *this = *this | rhs; }
775 
777  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
778  inline
779  auto
781  inv() const -> const multivector_t
782  {
783  auto result = matrix_multi_t(Scalar_T(1), this->frame());
784  return result /= matrix_multi_t(*this);
785  }
786 
788  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
789  auto
791  pow(int m) const -> const multivector_t
792  { return glucat::pow(*this, m); }
793 
795  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
796  auto
798  outer_pow(int m) const -> const multivector_t
799  {
800  if (m < 0)
801  throw error_t("outer_pow(int): negative exponent");
802  auto result = multivector_t(Scalar_T(1));
803  auto a = *this;
804  for (;
805  m != 0;
806  m >>= 1, a = a ^ a)
807  if (m & 1)
808  result ^= a;
809  return result;
810  }
811 
813  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
814  inline
815  auto
817  frame() const -> const index_set_t
818  {
819  auto result = index_set_t();
820  for (auto& this_term : *this)
821  result |= this_term.first;
822  return result;
823  }
824 
826  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
827  inline
828  auto
830  grade() const -> index_t
831  {
832  auto result = index_t(0);
833  for (auto& this_term : *this)
834  result = std::max(result, this_term.first.count());
835  return result;
836  }
837 
839  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
840  inline
841  auto
843  operator[] (const index_set_t ist) const -> Scalar_T
844  {
845  const auto& this_it = this->find(ist);
846  if (this_it == this->end())
847  return Scalar_T(0);
848  else
849  return this_it->second;
850  }
851 
853  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
854  auto
856  operator() (index_t grade) const -> const multivector_t
857  {
858  if ((grade < 0) || (grade > HI-LO))
859  return Scalar_T(0);
860  else
861  {
862  auto result = multivector_t();
863  for (auto& this_term : *this)
864  if (this_term.first.count() == grade)
865  result += this_term;
866  return result;
867  }
868  }
869 
871  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
872  inline
873  auto
875  scalar() const -> Scalar_T
876  { return (*this)[index_set_t()]; }
877 
879  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
880  inline
881  auto
883  pure() const -> const multivector_t
884  { return *this - this->scalar(); }
885 
887  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
888  auto
890  even() const -> const multivector_t
891  { // even part of x, sum of the pure(count) with even count
892  auto result = multivector_t();
893  for (auto& this_term : *this)
894  if ((this_term.first.count() % 2) == 0)
895  result.insert(this_term);
896  return result;
897  }
898 
900  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
901  auto
903  odd() const -> const multivector_t
904  { // even part of x, sum of the pure(count) with even count
905  auto result = multivector_t();
906  for (auto& this_term : *this)
907  if ((this_term.first.count() % 2) == 1)
908  result.insert(this_term);
909  return result;
910  }
911 
913  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
914  auto
916  vector_part() const -> const vector_t
917  { return this->vector_part(this->frame(), true); }
918 
920  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
921  auto
923  vector_part(const index_set_t frm, const bool prechecked) const -> const vector_t
924  {
925  if (!prechecked && (this->frame() | frm) != frm)
926  throw error_t("vector_part(frm): value is outside of requested frame");
927  auto result = vector_t();
928  result.reserve(frm.count());
929  const auto frm_end = frm.max()+1;
930  for (auto
931  idx = frm.min();
932  idx != frm_end;
933  ++idx)
934  // Frame may contain indices which do not correspond to a grade 1 term but
935  // frame cannot omit any index corresponding to a grade 1 term
936  if (frm[idx])
937  result.push_back((*this)[index_set_t(idx)]);
938  return result;
939  }
940 
942  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
943  auto
945  involute() const -> const multivector_t
946  {
947  auto result = *this;
948  for (auto& result_term : result)
949  { // for a k-vector u, involute(u) == (-1)^k * u
950  if ((result_term.first.count() % 2) == 1)
951  result_term.second *= Scalar_T(-1);
952  }
953  return result;
954  }
955 
957  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
958  auto
960  reverse() const -> const multivector_t
961  {
962  auto result = *this;
963  for (auto& result_term : result)
964  // For a k-vector u, reverse(u) = { -u, k == 2,3 (mod 4)
965  // { u, k == 0,1 (mod 4)
966  switch (result_term.first.count() % 4)
967  {
968  case 2:
969  case 3:
970  result_term.second *= Scalar_T(-1);
971  break;
972  default:
973  break;
974  }
975  return result;
976  }
977 
979  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
980  auto
982  conj() const -> const multivector_t
983  {
984  auto result = *this;
985  for (auto& result_term : result)
986  // For a k-vector u, conj(u) = { -u, k == 1,2 (mod 4)
987  // { u, k == 0,3 (mod 4)
988  switch (result_term.first.count() % 4)
989  {
990  case 1:
991  case 2:
992  result_term.second *= Scalar_T(-1);
993  break;
994  default:
995  break;
996  }
997  return result;
998  }
999 
1001  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1002  auto
1004  quad() const -> Scalar_T
1005  {
1006  // scalar(conj(x)*x) = 2*quad(even(x)) - quad(x)
1007  // ref: old clical: quadfunction(p:pter):pterm in file compmod.pas
1008  auto result = Scalar_T(0);
1009  for (auto& this_term : *this)
1010  {
1011  const auto sign =
1012  (this_term.first.count_neg() % 2)
1013  ? -Scalar_T(1)
1014  : Scalar_T(1);
1015  result += sign * (this_term.second) * (this_term.second);
1016  }
1017  return result;
1018  }
1019 
1021  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1022  auto
1024  norm() const -> Scalar_T
1025  {
1026  using traits_t = numeric_traits<Scalar_T>;
1027 
1028  auto result = Scalar_T(0);
1029  for (auto& this_term : *this)
1030  {
1031  const auto abs_crd = traits_t::abs(this_term.second);
1032  result += abs_crd * abs_crd;
1033  }
1034  return result;
1035  }
1036 
1038  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1039  auto
1041  max_abs() const -> Scalar_T
1042  {
1043  using traits_t = numeric_traits<Scalar_T>;
1044 
1045  auto result = Scalar_T(0);
1046  for (auto& this_term : *this)
1047  {
1048  const auto abs_crd = traits_t::abs(this_term.second);
1049  if (abs_crd > result)
1050  result = abs_crd;
1051  }
1052  return result;
1053  }
1054 
1056  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1057  auto
1059  random(const index_set_t frm, Scalar_T fill) -> const multivector_t
1060  {
1062  using index_set_t = typename multivector_t::index_set_t;
1063  using term_t = typename multivector_t::term_t;
1064 
1065  using random_generator_t = random_generator<Scalar_T>;
1066  auto& generator = random_generator_t::generator();
1067 
1068  fill =
1069  (fill < Scalar_T(0))
1070  ? Scalar_T(0)
1071  : (fill > Scalar_T(1))
1072  ? Scalar_T(1)
1073  : fill;
1074  const auto algebra_dim = set_value_t(1) << frm.count();
1075  using traits_t = numeric_traits<Scalar_T>;
1076  const auto mean_abs = traits_t::sqrt(Scalar_T(double(algebra_dim)));
1077  auto result = multivector_t();
1078  for (auto
1079  stv = set_value_t(0);
1080  stv != algebra_dim;
1081  ++stv)
1082  if (generator.uniform() < fill)
1083  {
1084  const auto& result_crd = generator.normal() / mean_abs;
1085  result.insert(term_t(index_set_t(stv, frm, true), result_crd));
1086  }
1087  return result;
1088  }
1089 
1091  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1092  inline
1093  void
1095  write(const std::string& msg) const
1096  { std::cout << msg << std::endl << " " << (*this) << std::endl; }
1097 
1099  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1100  inline
1101  void
1103  write(std::ofstream& ofile, const std::string& msg) const
1104  {
1105  if (!ofile)
1106  throw error_t("write(ofile,msg): cannot write to output file");
1107  ofile << msg << std::endl << " " << (*this) << std::endl;
1108  }
1109 
1111  template< typename Map_T,typename Sorted_Map_T >
1113  {
1114  public:
1115  using map_t = Map_T;
1116  using sorted_map_t = Sorted_Map_T;
1117  using sorted_iterator = typename Sorted_Map_T::const_iterator;
1118 
1119  sorted_range (Sorted_Map_T &sorted_val, const Map_T& val)
1120  {
1121  for (auto& val_term : val)
1122  sorted_val.insert(val_term);
1123  sorted_begin = sorted_val.begin();
1124  sorted_end = sorted_val.end();
1125  }
1128  };
1129 
1130  template< typename Sorted_Map_T >
1131  class sorted_range< Sorted_Map_T, Sorted_Map_T >
1132  {
1133  public:
1134  using map_t = Sorted_Map_T;
1135  using sorted_map_t = Sorted_Map_T;
1136  using sorted_iterator = typename Sorted_Map_T::const_iterator;
1137 
1138  sorted_range (Sorted_Map_T &sorted_val, const Sorted_Map_T& val)
1139  : sorted_begin( val.begin() ),
1140  sorted_end( val.end() )
1141  { }
1144  };
1145 
1147  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1148  auto
1149  operator<< (std::ostream& os, const framed_multi<Scalar_T,LO,HI,Tune_P>& val) -> std::ostream&
1150  {
1151  using limits_t = std::numeric_limits<Scalar_T>;
1152  if (val.empty())
1153  os << 0;
1154  else if (val.isnan())
1155  os << limits_t::quiet_NaN();
1156  else if (val.isinf())
1157  {
1158  const Scalar_T& inf = limits_t::infinity();
1159  os << (scalar(val) < 0.0 ? -inf : inf);
1160  }
1161  else
1162  {
1163  using traits_t = numeric_traits<Scalar_T>;
1164  using multivector_t = framed_multi<Scalar_T,LO,HI,Tune_P>;
1165  Scalar_T truncation;
1166  switch (os.flags() & std::ios::floatfield)
1167  {
1168  case std::ios_base::scientific:
1169  truncation = Scalar_T(1) / traits_t::pow(Scalar_T(10), int(os.precision()) + 1);
1170  break;
1171  case std::ios_base::fixed:
1172  truncation = Scalar_T(1) / (traits_t::pow(Scalar_T(10), int(os.precision())) * val.max_abs());
1173  break;
1174  case std::ios_base::fixed | std::ios_base::scientific:
1175  truncation = multivector_t::default_truncation;
1176  break;
1177  default:
1178  truncation = Scalar_T(1) / traits_t::pow(Scalar_T(10), int(os.precision()));
1179  break;
1180  }
1181  auto truncated_val = val.truncated(truncation);
1182  if (truncated_val.empty())
1183  os << 0;
1184  else
1185  {
1186  using map_t = typename multivector_t::map_t;
1187  using sorted_map_t = typename multivector_t::sorted_map_t;
1188  using sorted_iterator = typename sorted_map_t::const_iterator;
1189  auto sorted_val = sorted_map_t();
1190  const auto sorted_val_range = sorted_range< map_t, sorted_map_t >(sorted_val, truncated_val);
1191  auto sorted_it = sorted_val_range.sorted_begin;
1192  os << *sorted_it;
1193  for (++sorted_it;
1194  sorted_it != sorted_val_range.sorted_end;
1195  ++sorted_it)
1196  {
1197  const Scalar_T& scr = sorted_it->second;
1198  if (scr >= 0.0)
1199  os << '+';
1200  os << *sorted_it;
1201  }
1202  }
1203  }
1204  return os;
1205  }
1206 
1208  template< typename Scalar_T, const index_t LO, const index_t HI >
1209  auto
1210  operator<< (std::ostream& os, const std::pair< const index_set<LO,HI>, Scalar_T >& term) -> std::ostream&
1211  {
1212  const auto second_as_double = numeric_traits<Scalar_T>::to_double(term.second);
1213  const auto use_double =
1214  (os.precision() <= std::numeric_limits<double>::digits10) ||
1215  (term.second == Scalar_T(second_as_double));
1216  if (term.first.count() == 0)
1217  if (use_double)
1218  os << second_as_double;
1219  else
1220  os << term.second;
1221  else if (term.second == Scalar_T(-1))
1222  {
1223  os << '-';
1224  os << term.first;
1225  }
1226  else if (term.second != Scalar_T(1))
1227  {
1228  if (use_double)
1229  {
1230  auto tol = std::pow(10.0,-os.precision());
1231  if ( std::fabs(second_as_double + 1.0) < tol )
1232  os << '-';
1233  else if ( std::fabs(second_as_double - 1.0) >= tol )
1234  os << second_as_double;
1235  }
1236  else
1237  os << term.second;
1238  os << term.first;
1239  }
1240  else
1241  os << term.first;
1242  return os;
1243  }
1244 
1246  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1247  auto
1248  operator>> (std::istream& s, framed_multi<Scalar_T,LO,HI,Tune_P> & val) -> std::istream&
1249  { // Input looks like 1.0-2.0{1,2}+3.2{3,4}.
1250  using multivector_t = framed_multi<Scalar_T,LO,HI,Tune_P>;
1251  // Parsing variables.
1252  auto local_val = multivector_t();
1253  auto c = 0;
1254  // Parsing control variables.
1255  auto negative = false;
1256  auto expect_term = true;
1257  // The multivector may begin with '+' or '-'. Check for this.
1258  c = s.peek();
1259  if (s.good() && (c == int('+') || c == int('-')))
1260  { // A '-' here negates the following term.
1261  negative = (c == int('-'));
1262  // Consume the '+' or '-'.
1263  s.get();
1264  }
1265  while (s.good())
1266  { // Parse a term.
1267  // A term consists of an optional scalar, followed by an optional index set.
1268  // At least one of the two must be present.
1269  // Default coordinate is Scalar_T(1).
1270  auto coordinate = Scalar_T(1);
1271  // Default index set is empty.
1272  auto ist = index_set<LO,HI>();
1273  // First, check for an opening brace.
1274  c = s.peek();
1275  if (s.good())
1276  { // If the character is not an opening brace,
1277  // a coordinate value is expected here.
1278  if (c != int('{'))
1279  { // Try to read a coordinate value.
1280  double coordinate_as_double;
1281  s >> coordinate_as_double;
1282  // Reading the coordinate may have resulted in an end of file condition.
1283  // This is not a failure.
1284  if (s)
1285  coordinate = Scalar_T(coordinate_as_double);
1286  }
1287  }
1288  else
1289  { // End of file here ends parsing while a term may still be expected.
1290  break;
1291  }
1292  // Coordinate is now Scalar_T(1) or a Scalar_T value.
1293  // Parse an optional index set.
1294  if (s.good())
1295  {
1296  c = s.peek();
1297  if (s.good() && c == int('{'))
1298  { // Try to read index set.
1299  s >> ist;
1300  }
1301  }
1302  // Reading the term may have resulted in an end of file condition.
1303  // This is not a failure.
1304  if (s)
1305  {
1306  // Immediately after parsing a term, another term is not expected.
1307  expect_term = false;
1308  if (coordinate != Scalar_T(0))
1309  {
1310  // Add the term to the local multivector.
1311  coordinate =
1312  negative
1313  ? -coordinate
1314  : coordinate;
1315  using term_t = typename multivector_t::term_t;
1316  local_val += term_t(ist, coordinate);
1317  }
1318  }
1319  // Check if anything follows the current term.
1320  if (s.good())
1321  {
1322  c = s.peek();
1323  if (s.good())
1324  { // Only '+' and '-' are valid here.
1325  if (c == int('+') || c == int('-'))
1326  { // A '-' here negates the following term.
1327  negative = (c == int('-'));
1328  // Consume the '+' or '-'.
1329  s.get();
1330  // Immediately after '+' or '-',
1331  // expect another term.
1332  expect_term = true;
1333  }
1334  else
1335  { // Any other character here is a not failure,
1336  // but still ends the parsing of the multivector.
1337  break;
1338  }
1339  }
1340  }
1341  }
1342  // If a term is still expected, this is a failure.
1343  if (expect_term)
1344  s.clear(std::istream::failbit);
1345  // End of file is not a failure.
1346  if (s)
1347  { // The multivector has been successfully parsed.
1348  val = local_val;
1349  }
1350  return s;
1351  }
1352 
1354  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1355  auto
1357  nbr_terms () const -> unsigned long
1358  { return this->size(); }
1359 
1361  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1362  inline
1363  auto
1365  operator+= (const term_t& term) -> multivector_t&
1366  { // Do not insert terms with 0 coordinate
1367  if (term.second != Scalar_T(0))
1368  {
1369  const auto& this_it = this->find(term.first);
1370  if (this_it == this->end())
1371  this->insert(term);
1372  else if (this_it->second + term.second == Scalar_T(0))
1373  // Erase term if resulting coordinate is 0
1374  this->erase(this_it);
1375  else
1376  this_it->second += term.second;
1377  }
1378  return *this;
1379  }
1380 
1382  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1383  auto
1385  isinf() const -> bool
1386  {
1387  using traits_t = numeric_traits<Scalar_T>;
1388 
1389  if (std::numeric_limits<Scalar_T>::has_infinity)
1390  for (auto& this_term : *this)
1391  if (traits_t::isInf(this_term.second))
1392  return true;
1393  return false;
1394  }
1395 
1397  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1398  auto
1400  isnan() const -> bool
1401  {
1402  using traits_t = numeric_traits<Scalar_T>;
1403 
1404  if (std::numeric_limits<Scalar_T>::has_quiet_NaN)
1405  for (auto& this_term : *this)
1406  if (traits_t::isNaN(this_term.second))
1407  return true;
1408  return false;
1409  }
1410 
1412  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1413  auto
1415  truncated(const Scalar_T& limit) const -> const multivector_t
1416  {
1417  using traits_t = numeric_traits<Scalar_T>;
1418 
1419  if (this->isnan() || this->isinf())
1420  return *this;
1421  const auto truncation = traits_t::abs(limit);
1422  const auto top = max_abs();
1423  auto result = multivector_t();
1424  if (top != Scalar_T(0))
1425  for (auto& this_term : *this)
1426  if (traits_t::abs(this_term.second) > top * truncation)
1427  result.insert(this_term);
1428  return result;
1429  }
1430 
1432  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1433  auto
1435  fold(const index_set_t frm) const -> multivector_t
1436  {
1437  if (frm.is_contiguous())
1438  return *this;
1439  else
1440  {
1441  auto result = multivector_t();
1442  for (auto& this_term : *this)
1443  result.insert(term_t(this_term.first.fold(frm), this_term.second));
1444  return result;
1445  }
1446  }
1447 
1449  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1450  auto
1452  unfold(const index_set_t frm) const -> multivector_t
1453  {
1454  if (frm.is_contiguous())
1455  return *this;
1456  else
1457  {
1458  auto result = multivector_t();
1459  for (auto& this_term : *this)
1460  result.insert(term_t(this_term.first.unfold(frm), this_term.second));
1461  return result;
1462  }
1463  }
1464 
1466  // Reference: [L] 16.4 Periodicity of 8, p216
1467  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1468  auto
1471  {
1472  // We add 4 to q by subtracting 4 from p
1473  if (q+4 > -LO)
1474  throw error_t("centre_pm4_qp4(p,q): LO is too high to represent this value");
1475  if (this->frame().max() > p-4)
1476  {
1477  using index_pair_t = typename index_set_t::index_pair_t;
1478  const auto pm3210 = index_set_t(index_pair_t(p-3,p), true);
1479  const auto qm4321 = index_set_t(index_pair_t(-q-4,-q-1), true);
1480  const auto& tqm4321 = term_t(qm4321, Scalar_T(1));
1481  auto result = multivector_t();
1482  for (auto& this_term : *this)
1483  {
1484  const auto ist = this_term.first;
1485  if (ist.max() > p-4)
1486  {
1487  auto var_term = var_term_t();
1488  for (auto
1489  n = index_t(0);
1490  n != index_t(4);
1491  ++n)
1492  if (ist[n+p-3])
1493  var_term *= term_t(index_set_t(n-q-4), Scalar_T(1)) * tqm4321;
1494  // Mask out {p-3}..{p}
1495  result.insert(term_t(ist & ~pm3210, this_term.second) *
1496  term_t(var_term.first, var_term.second));
1497  }
1498  else
1499  result.insert(this_term);
1500  }
1501  *this = result;
1502  }
1503  p -=4; q += 4;
1504  return *this;
1505  }
1506 
1508  // Reference: [L] 16.4 Periodicity of 8, p216
1509  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1510  auto
1513  {
1514  // We add 4 to p by subtracting 4 from q
1515  if (p+4 > HI)
1516  throw error_t("centre_pp4_qm4(p,q): HI is too low to represent this value");
1517  if (this->frame().min() < -q+4)
1518  {
1519  using index_pair_t = typename index_set_t::index_pair_t;
1520  const auto qp0123 = index_set_t(index_pair_t(-q,-q+3), true);
1521  const auto pp1234 = index_set_t(index_pair_t(p+1,p+4), true);
1522  const auto& tpp1234 = term_t(pp1234, Scalar_T(1));
1523  auto result = multivector_t();
1524  for (auto& this_term : *this)
1525  {
1526  index_set_t ist = this_term.first;
1527  if (ist.min() < -q+4)
1528  {
1529  auto var_term = var_term_t();
1530  for (auto
1531  n = index_t(0);
1532  n != index_t(4);
1533  ++n)
1534  if (ist[n-q])
1535  var_term *= term_t(index_set_t(n+p+1), Scalar_T(1)) * tpp1234;
1536  // Mask out {-q}..{-q+3}
1537  result.insert(term_t(var_term.first, var_term.second) *
1538  term_t(ist & ~qp0123, this_term.second));
1539  }
1540  else
1541  result.insert(this_term);
1542  }
1543  *this = result;
1544  }
1545  p +=4; q -= 4;
1546  return *this;
1547  }
1548 
1550  // Reference: [P] Proposition 15.20, p 131
1551  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1552  auto
1555  {
1556  if (q+1 > HI)
1557  throw error_t("centre_qp1_pm1(p,q): HI is too low to represent this value");
1558  if (p-1 > -LO)
1559  throw error_t("centre_qp1_pm1(p,q): LO is too high to represent this value");
1560  const auto qp1 = index_set_t(q+1);
1561  const auto& tqp1 = term_t(qp1, Scalar_T(1));
1562  auto result = multivector_t();
1563  for (auto& this_term : *this)
1564  {
1565  const auto ist = this_term.first;
1566  auto var_term = var_term_t(index_set_t(), this_term.second);
1567  for (auto
1568  n = -q;
1569  n != p;
1570  ++n)
1571  if (n != 0 && ist[n])
1572  var_term *= term_t(index_set_t(-n) | qp1, Scalar_T(1));
1573  if (p != 0 && ist[p])
1574  var_term *= tqp1;
1575  result.insert(term_t(var_term.first, var_term.second));
1576  }
1577  index_t orig_p = p;
1578  p = q+1;
1579  q = orig_p-1;
1580  return *this = result;
1581  }
1582 
1584  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1585  auto
1587  divide(const index_set_t ist) const -> const framed_pair_t
1588  {
1589  auto quo = multivector_t();
1590  auto rem = multivector_t();
1591  for (auto& this_term : *this)
1592  if ((this_term.first | ist) == this_term.first)
1593  quo.insert(term_t(this_term.first ^ ist, this_term.second));
1594  else
1595  rem.insert(this_term);
1596  return framed_pair_t(quo, rem);
1597  }
1598 
1600  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1601  auto
1603  fast(const index_t level, const bool odd) const -> const matrix_t
1604  {
1605  // Assume val is already folded and centred
1606  if (this->empty())
1607  {
1608  using matrix_index_t = typename matrix_multi_t::matrix_index_t;
1609  const auto dim = matrix_index_t(1) << level;
1610  auto result =matrix_t(dim, dim);
1611  result.clear();
1612  return result;
1613  }
1614  if (level == 0)
1615  return matrix::unit<matrix_t>(1) * this->scalar();
1616 
1617  using basis_matrix_t = typename matrix_multi_t::basis_matrix_t;
1618  using basis_scalar_t = typename basis_matrix_t::value_type;
1619 
1620  const auto& I = matrix::unit<basis_matrix_t>(2);
1621  auto J = basis_matrix_t(2,2,2);
1622  J.clear();
1623  J(0,1) = basis_scalar_t(-1);
1624  J(1,0) = basis_scalar_t( 1);
1625  auto K = J;
1626  K(0,1) = basis_scalar_t( 1);
1627  auto JK = I;
1628  JK(0,0) = basis_scalar_t(-1);
1629 
1630  const auto ist_mn = index_set_t(-level);
1631  const auto ist_pn = index_set_t(level);
1632  if (level == 1)
1633  {
1634  if (odd)
1635  return matrix_t(J) * (*this)[ist_mn] + matrix_t(K) * (*this)[ist_pn];
1636  else
1637  return matrix_t(I) * this->scalar() + matrix_t(JK) * (*this)[ist_mn ^ ist_pn];
1638  }
1639  else
1640  {
1641  const auto& pair_mn = this->divide(ist_mn);
1642  const auto& quo_mn = pair_mn.first;
1643  const auto& rem_mn = pair_mn.second;
1644  const auto& pair_quo_mnpn = quo_mn.divide(ist_pn);
1645  const auto& val_mnpn = pair_quo_mnpn.first;
1646  const auto& val_mn = pair_quo_mnpn.second;
1647  const auto& pair_rem_mnpn = rem_mn.divide(ist_pn);
1648  const auto& val_pn = pair_rem_mnpn.first;
1649  const auto& val_1 = pair_rem_mnpn.second;
1650  using matrix::kron;
1651  if (odd)
1652  return - kron(JK, val_1.fast (level-1, 1))
1653  + kron(I, val_mnpn.fast(level-1, 1))
1654  + kron(J, val_mn.fast (level-1, 0))
1655  + kron(K, val_pn.fast (level-1, 0));
1656  else
1657  return kron(I, val_1.fast (level-1, 0))
1658  + kron(JK, val_mnpn.fast(level-1, 0))
1659  + kron(K, val_mn.fast (level-1, 1))
1660  - kron(J, val_pn.fast (level-1, 1));
1661  }
1662  }
1663 
1665  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1666  template< typename Other_Scalar_T >
1667  auto
1670  {
1671  // Fold val
1672  auto val = this->fold(frm);
1673  auto p = frm.count_pos();
1674  auto q = frm.count_neg();
1675  const auto bott_offset = gen::offset_to_super[pos_mod(p - q, 8)];
1676  p += std::max(bott_offset,index_t(0));
1677  q -= std::min(bott_offset,index_t(0));
1678  if (p > HI)
1679  throw error_t("fast_matrix_multi(frm): HI is too low to represent this value");
1680  if (q > -LO)
1681  throw error_t("fast_matrix_multi(frm): LO is too high to represent this value");
1682  // Centre val
1683  while (p - q > 4)
1684  val.centre_pm4_qp4(p, q);
1685  while (p - q < -3)
1686  val.centre_pp4_qm4(p, q);
1687  if (p - q > 1)
1688  val.centre_qp1_pm1(p, q);
1689  const index_t level = (p + q)/2;
1690 
1691  // Do the fast transform
1692  const auto& ev_val = val.even();
1693  const auto& od_val = val.odd();
1694  return matrix_multi<Other_Scalar_T,LO,HI,Tune_P>(ev_val.fast(level, 0) + od_val.fast(level, 1), frm);
1695  }
1696 
1697  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1698  inline
1699  auto
1702  { return *this; }
1703 
1705  template< typename Scalar_T, const index_t LO, const index_t HI >
1706  inline
1707  static
1708  auto
1709  crd_of_mult(const std::pair<const index_set<LO,HI>, Scalar_T>& lhs,
1710  const std::pair<const index_set<LO,HI>, Scalar_T>& rhs) -> Scalar_T
1711  { return lhs.first.sign_of_mult(rhs.first) * lhs.second * rhs.second; }
1712 
1714  template< typename Scalar_T, const index_t LO, const index_t HI >
1715  inline
1716  auto
1717  operator* (const std::pair<const index_set<LO,HI>, Scalar_T>& lhs,
1718  const std::pair<const index_set<LO,HI>, Scalar_T>& rhs) -> const std::pair<const index_set<LO,HI>, Scalar_T>
1719  {
1720  using term_t = std::pair<const index_set<LO,HI>, Scalar_T>;
1721  return term_t(lhs.first ^ rhs.first, crd_of_mult(lhs, rhs));
1722  }
1723 
1725  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1726  auto
1728  {
1729  using traits_t = numeric_traits<Scalar_T>;
1730  if (val.isnan())
1731  return traits_t::NaN();
1732 
1733  check_complex(val, i, prechecked);
1734 
1735  const auto realval = val.scalar();
1736  if (val == realval)
1737  {
1738  if (realval < Scalar_T(0))
1739  return i * traits_t::sqrt(-realval);
1740  else
1741  return traits_t::sqrt(realval);
1742  }
1743  using matrix_multi_t = typename framed_multi<Scalar_T,LO,HI,Tune_P>::matrix_multi_t;
1744  return sqrt(matrix_multi_t(val), matrix_multi_t(i), prechecked);
1745  }
1746 
1748  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1749  auto
1751  {
1752  using traits_t = numeric_traits<Scalar_T>;
1753  if (val.isnan())
1754  return traits_t::NaN();
1755 
1756  const auto s = scalar(val);
1757  if (val == s)
1758  return traits_t::exp(s);
1759 
1760  const double size = val.size();
1761  const auto frm_count = val.frame().count();
1762  const auto algebra_dim = set_value_t(1) << frm_count;
1763 
1764  if( (size * size <= double(algebra_dim)) || (frm_count < Tune_P::mult_matrix_threshold))
1765  {
1766  switch (Tune_P::function_precision)
1767  {
1768  case precision_demoted:
1769  {
1770  using demoted_scalar_t = typename traits_t::demoted::type;
1771  using demoted_multivector_t = framed_multi<demoted_scalar_t,LO,HI,Tune_P>;
1772 
1773  const auto& demoted_val = demoted_multivector_t(val);
1774  return clifford_exp(demoted_val);
1775  }
1776  break;
1777  case precision_promoted:
1778  {
1779  using promoted_scalar_t = typename traits_t::promoted::type;
1780  using promoted_multivector_t = framed_multi<promoted_scalar_t,LO,HI,Tune_P>;
1781 
1782  const auto& promoted_val = promoted_multivector_t(val);
1783  return clifford_exp(promoted_val);
1784  }
1785  break;
1786  default:
1787  return clifford_exp(val);
1788  }
1789  }
1790  else
1791  {
1792  using matrix_multi_t = matrix_multi<Scalar_T,LO,HI,Tune_P>;
1793  return exp(matrix_multi_t(val));
1794  }
1795  }
1796 
1798  template< typename Scalar_T, const index_t LO, const index_t HI, typename Tune_P >
1799  auto
1801  {
1802  using traits_t = numeric_traits<Scalar_T>;
1803  if (val == Scalar_T(0) || val.isnan())
1804  return traits_t::NaN();
1805 
1806  check_complex(val, i, prechecked);
1807 
1808  const auto realval = val.scalar();
1809  if (val == realval)
1810  {
1811  if (realval < Scalar_T(0))
1812  return i * traits_t::pi() + traits_t::log(-realval);
1813  else
1814  return traits_t::log(realval);
1815  }
1816  using matrix_multi_t = typename framed_multi<Scalar_T,LO,HI,Tune_P>::matrix_multi_t;
1817  return log(matrix_multi_t(val), matrix_multi_t(i), prechecked);
1818  }
1819 }
1820 #endif // _GLUCAT_FRAMED_MULTI_IMP_H
virtual auto operator*=(const Scalar_T &scr) -> multivector_t &=0
Product of multivector and scalar.
std::pair< const multivector_t, const multivector_t > framed_pair_t
Definition: framed_multi.h:164
auto operator%(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const RHS< Scalar_T, LO, HI, Tune_P > &rhs) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Left contraction.
virtual auto quad() const -> Scalar_T=0
Scalar_T quadratic form == (rev(x)*x)(0)
auto operator^(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const RHS< Scalar_T, LO, HI, Tune_P > &rhs) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Outer product.
auto unfold(const index_set_t frm) const -> multivector_t
Subalgebra isomorphism: unfold each term within the given frame.
Random number generator with single instance per Scalar_T.
Definition: random.h:42
auto exp(const framed_multi< Scalar_T, LO, HI, Tune_P > &val) -> const framed_multi< Scalar_T, LO, HI, Tune_P >
Exponential of multivector.
virtual auto outer_pow(int m) const -> const multivector_t=0
Outer product power.
A framed_multi<Scalar_T,LO,HI,Tune_P> is a framed approximation to a multivector. ...
Definition: framed_multi.h:56
virtual auto operator%=(const multivector_t &rhs) -> multivector_t &=0
Contraction.
#define _GLUCAT_HASH_N(x)
auto odd(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Odd part.
float pi
Definition: PyClical.pyx:1907
virtual auto norm() const -> Scalar_T=0
Scalar_T norm == sum of norm of coordinates.
auto fold(const index_set_t frm) const -> multivector_t
Subalgebra isomorphism: fold each term within the given frame.
auto centre_pm4_qp4(index_t &p, index_t &q) -> multivector_t &
Subalgebra isomorphism: R_{p,q} to R_{p-4,q+4}.
auto scalar(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> Scalar_T
Scalar part.
virtual auto operator[](const index_set_t ist) const -> Scalar_T=0
Subscripting: map from index set to scalar coordinate.
sorted_iterator sorted_end
Extra traits which extend numeric limits.
Definition: scalar.h:47
virtual auto inv() const -> const multivector_t=0
Geometric multiplicative inverse.
index_set< LO, HI > index_set_t
Definition: framed_multi.h:135
auto star(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const RHS< Scalar_T, LO, HI, Tune_P > &rhs) -> Scalar_T
Hestenes scalar product.
virtual auto conj() const -> const multivector_t=0
Conjugation, reverse o involute == involute o reverse.
virtual auto even() const -> const multivector_t=0
Even part of multivector, sum of even grade terms.
virtual auto operator-=(const multivector_t &rhs) -> multivector_t &=0
Geometric difference.
auto pos_mod(LHS_T lhs, RHS_T rhs) -> LHS_T
Modulo function which works reliably for lhs < 0.
Definition: global.h:117
typename matrix_multi_t::matrix_t matrix_t
Definition: framed_multi.h:148
auto clifford_exp(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Exponential of multivector.
auto log(const Multivector< Scalar_T, LO, HI, Tune_P > &val, const Multivector< Scalar_T, LO, HI, Tune_P > &i, const bool prechecked=false) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Natural logarithm of multivector with specified complexifier.
auto isnan(const Matrix_T &m) -> bool
Not a Number.
Definition: matrix_imp.h:292
auto fast_framed_multi() const -> const framed_multi_t
Use inverse generalized FFT to construct a framed_multi_t.
A matrix_multi<Scalar_T,LO,HI,Tune_P> is a matrix approximation to a multivector. ...
Definition: framed_multi.h:59
virtual void write(const std::string &msg="") const=0
Write formatted multivector to output.
virtual auto operator^=(const multivector_t &rhs) -> multivector_t &=0
Outer product.
unsigned long set_value_t
Size of set_value_t should be enough to contain index_set<LO,HI>
Definition: global.h:79
auto operator>>(std::istream &s, framed_multi< Scalar_T, LO, HI, Tune_P > &val) -> std::istream &
Read multivector from input.
Sorted range for use with output.
class var_term var_term_t
Definition: framed_multi.h:147
auto basis_element(const index_set< LO, HI > &ist) const -> const basis_matrix_t
Create a basis element matrix within the current frame.
virtual auto reverse() const -> const multivector_t=0
Reversion, eg. {1}*{2} -> {2}*{1}.
auto kron(const LHS_T &lhs, const RHS_T &rhs) -> const RHS_T
Kronecker tensor product of matrices - as per Matlab kron.
Definition: matrix_imp.h:83
virtual auto operator()(index_t grade) const -> const multivector_t=0
Pure grade-vector part.
virtual auto frame() const -> const index_set_t=0
Subalgebra generated by all generators of terms of given multivector.
static auto 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) -> Scalar_T
Coordinate of product of terms.
virtual auto vector_part() const -> const vector_t=0
Vector part of multivector, as a vector_t with respect to frame()
auto abs(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> Scalar_T
Absolute value == sqrt(norm)
std::vector< Scalar_T > vector_t
Definition: framed_multi.h:137
std::unordered_map< index_set_t, Scalar_T, index_set_hash< LO, HI > > map_t
Definition: framed_multi.h:150
Abstract exception class.
Definition: errors.h:41
auto operator &(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const RHS< Scalar_T, LO, HI, Tune_P > &rhs) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Inner product.
typename Sorted_Map_T::const_iterator sorted_iterator
_GLUCAT_CLIFFORD_ALGEBRA_OPERATIONS auto nbr_terms() const -> unsigned long
Number of terms.
virtual auto truncated(const Scalar_T &limit=default_truncation) const -> const multivector_t=0
Remove all terms with relative size smaller than limit.
auto centre_qp1_pm1(index_t &p, index_t &q) -> multivector_t &
Subalgebra isomorphism: R_{p,q} to R_{q+1,p-1}.
std::pair< index_t, index_t > index_pair_t
Definition: index_set.h:85
static auto to_double(const Scalar_T &val) -> double
Cast to double.
Definition: scalar.h:133
auto operator*(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const Scalar_T &scr) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Product of multivector and scalar.
virtual auto odd() const -> const multivector_t=0
Odd part of multivector, sum of odd grade terms.
auto centre_pp4_qm4(index_t &p, index_t &q) -> multivector_t &
Subalgebra isomorphism: R_{p,q} to R_{p+4,q-4}.
Specific exception class.
Definition: errors.h:56
typename matrix_t::size_type matrix_index_t
Definition: matrix_multi.h:159
auto divide(const index_set_t ist) const -> const framed_pair_t
Divide multivector into part divisible by index_set and remainder.
sorted_iterator sorted_begin
framed_multi()
Default constructor.
auto fast_matrix_multi(const index_set_t frm) const -> const matrix_multi< Other_Scalar_T, LO, HI, Tune_P >
Use generalized FFT to construct a matrix_multi_t.
int index_t
Size of index_t should be enough to represent LO, HI.
Definition: global.h:77
static const std::array< index_t, 8 > offset_to_super
Offsets between the current signature and that of the real superalgebra.
Definition: generation.h:86
auto max() const -> index_t
Maximum member.
sorted_range(Sorted_Map_T &sorted_val, const Map_T &val)
virtual auto grade() const -> index_t=0
Maximum of the grades of each term.
auto sqrt(const Multivector< Scalar_T, LO, HI, Tune_P > &val, const Multivector< Scalar_T, LO, HI, Tune_P > &i, const bool prechecked=false) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Square root of multivector with specified complexifier.
auto nnz(const Matrix_T &m) -> typename Matrix_T::size_type
Number of non-zeros.
Definition: matrix_imp.h:258
Index set class based on std::bitset<> in Gnu standard C++ library.
Definition: index_set.h:45
virtual auto isinf() const -> bool=0
Check if a multivector contains any infinite values.
auto operator/(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const Scalar_T &scr) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Quotient of multivector and scalar.
virtual auto operator==(const multivector_t &val) const -> bool=0
Test for equality of multivectors.
auto vector_part(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const std::vector< Scalar_T >
Vector part of multivector, as a vector_t with respect to frame()
virtual auto operator/=(const Scalar_T &scr) -> multivector_t &=0
Quotient of multivector and scalar.
auto pow(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, int rhs) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Integer power of multivector.
static auto classname() -> const std::string
Class name used in messages.
ublas::compressed_matrix< int, orientation_t > basis_matrix_t
Definition: matrix_multi.h:157
#define _GLUCAT_HASH_SIZE_T(x)
sorted_range(Sorted_Map_T &sorted_val, const Sorted_Map_T &val)
auto fast(const index_t level, const bool odd) const -> const matrix_t
Generalized FFT from multivector_t to matrix_t.
virtual auto max_abs() const -> Scalar_T=0
Maximum of absolute values of components of multivector: multivector infinity norm.
static void check_complex(const Multivector< Scalar_T, LO, HI, Tune_P > &val, const Multivector< Scalar_T, LO, HI, Tune_P > &i, const bool prechecked=false)
Check that i is a valid complexifier for val.
typename Sorted_Map_T::const_iterator sorted_iterator
auto max_abs(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> Scalar_T
Maximum of absolute values of components of multivector: multivector infinity norm.
matrix_t m_matrix
Matrix value representing the multivector within the folded frame.
Definition: matrix_multi.h:280
virtual auto operator|=(const multivector_t &rhs) -> multivector_t &=0
Transformation via twisted adjoint action.
def e(obj)
Definition: PyClical.pyx:1936
virtual auto involute() const -> const multivector_t=0
Main involution, each {i} is replaced by -{i} in each term, eg. {1} -> -{1}.
static auto random(const index_set_t frm, Scalar_T fill=Scalar_T(1)) -> const multivector_t
Random multivector within a frame.
virtual auto isnan() const -> bool=0
Check if a multivector contains any IEEE NaN values.
error< multivector_t > error_t
Definition: framed_multi.h:138
auto min() const -> index_t
Minimum member.
auto isinf(const Matrix_T &m) -> bool
Infinite.
Definition: matrix_imp.h:275
framed_multi multivector_t
Definition: framed_multi.h:131
std::pair< const index_set_t, Scalar_T > term_t
Definition: framed_multi.h:136
virtual auto operator&=(const multivector_t &rhs) -> multivector_t &=0
Inner product.
auto count() const -> index_t
Cardinality: Number of indices included in set.
auto operator|(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const RHS< Scalar_T, LO, HI, Tune_P > &rhs) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Transformation via twisted adjoint action.
auto operator+=(const term_t &term) -> multivector_t &
Add a term, if non-zero.