glucat  0.12.0
index_set_imp.h
Go to the documentation of this file.
1 #ifndef _GLUCAT_INDEX_SET_IMP_H
2 #define _GLUCAT_INDEX_SET_IMP_H
3 /***************************************************************************
4  GluCat : Generic library of universal Clifford algebra templates
5  index_set_imp.h : Implement a class for a set of non-zero integer indices
6  -------------------
7  begin : Sun 2001-12-09
8  copyright : (C) 2001-2016 by Paul C. Leopardi
9  ***************************************************************************
10 
11  This library is free software: you can redistribute it and/or modify
12  it under the terms of the GNU Lesser General Public License as published
13  by the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  This library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19  GNU Lesser General Public License for more details.
20 
21  You should have received a copy of the GNU Lesser General Public License
22  along with this library. If not, see <http://www.gnu.org/licenses/>.
23 
24  ***************************************************************************
25  This library is based on a prototype written by Arvind Raja and was
26  licensed under the LGPL with permission of the author. See Arvind Raja,
27  "Object-oriented implementations of Clifford algebras in C++: a prototype",
28  in Ablamowicz, Lounesto and Parra (eds.)
29  "Clifford algebras with numeric and symbolic computations", Birkhauser, 1996.
30  ***************************************************************************
31  See also Arvind Raja's original header comments in glucat.h
32  ***************************************************************************/
33 
34 #include "glucat/index_set.h"
35 
36 #include <string>
37 #include <sstream>
38 
39 namespace glucat
40 {
41  // References for algorithms:
42  // [JA]: Joerg Arndt, "Algorithms for programmers", http://www.jjj.de/fxt/fxtbook.pdf
43  // Chapter 1, Bit wizardry, http://www.jjj.de/bitwizardry/bitwizardrypage.html
44  // [L]: Pertti Lounesto, "Clifford algebras and spinors", Cambridge UP, 1997.
45 
46  template<const index_t LO, const index_t HI>
47  inline
48  auto
50  classname() -> const std::string
51  { return "index_set"; }
52 
54  template<const index_t LO, const index_t HI>
56  index_set(const index_t idx)
57  { this->set(idx); }
58 
60  template<const index_t LO, const index_t HI>
62  index_set(const bitset_t bst):
63  bitset_t(bst)
64  { }
65 
67  template<const index_t LO, const index_t HI>
69  index_set(const set_value_t folded_val, const index_set_t frm, const bool prechecked)
70  {
71  if (!prechecked && folded_val >= (set_value_t(1) << frm.count()))
72  throw error_t("index_set(val,frm): cannot create: value gives an index set outside of frame");
73  const index_set_t folded_frame = frm.fold();
74  const index_t min_index = folded_frame.min();
75  const index_t skip = min_index > 0 ? 1 : 0;
76  const index_set_t folded_set = index_set_t(bitset_t(folded_val) << (min_index - skip - LO));
77  *this = folded_set.unfold(frm);
78  }
79 
81  template<const index_t LO, const index_t HI>
83  index_set(const index_pair_t& range, const bool prechecked)
84  {
85  if (!prechecked && (range.first < LO || range.second > HI))
86  throw error_t("index_set(range): cannot create: range is too large");
87  const index_t begin_bit = (range.first < 0)
88  ? range.first-LO
89  : range.first-LO-1;
90  const index_t end_bit = (range.second < 0)
91  ? range.second-LO+1
92  : range.second-LO;
93  unsigned long mask = ( (end_bit == _GLUCAT_BITS_PER_ULONG)
94  ? -1UL
95  : (1UL << end_bit)-1UL)
96  & ~((1UL << begin_bit)-1UL);
97  *this = bitset_t(mask);
98  }
99 
101  template<const index_t LO, const index_t HI>
103  index_set(const std::string& str)
104  {
105  std::istringstream ss(str);
106  ss >> *this;
107  if (!ss)
108  throw error_t("index_set_t(str): could not parse string");
109  // Peek to see if the end of the string has been reached.
110  ss.peek();
111  if (!ss.eof())
112  throw error_t("index_set_t(str): could not parse entire string");
113  }
114 
116  template<const index_t LO, const index_t HI>
117  inline
118  auto
120  operator== (const index_set_t rhs) const -> bool
121  {
122  const auto* pthis = static_cast<const bitset_t*>(this);
123  return *pthis == static_cast<bitset_t>(rhs);
124  }
125 
127  template<const index_t LO, const index_t HI>
128  inline
129  auto
131  operator!= (const index_set_t rhs) const -> bool
132  {
133  const auto* pthis = static_cast<const bitset_t*>(this);
134  return *pthis != static_cast<bitset_t>(rhs);
135  }
136 
138  template<const index_t LO, const index_t HI>
139  inline
140  auto
143  { return bitset_t::operator~(); }
144 
146  template<const index_t LO, const index_t HI>
147  inline
148  auto
151  {
152  bitset_t* pthis = this;
153  *pthis ^= static_cast<bitset_t>(rhs);
154  return *this;
155  }
156 
158  template<const index_t LO, const index_t HI>
159  inline
160  auto
162  const index_set<LO,HI>& rhs) -> const
164  {
165  using index_set_t = index_set<LO, HI>;
166  using bitset_t = typename index_set_t::bitset_t;
167  return static_cast<bitset_t>(lhs) ^ static_cast<bitset_t>(rhs);
168  }
169 
171  template<const index_t LO, const index_t HI>
172  inline
173  auto
175  operator&= (const index_set_t rhs) -> index_set_t&
176  {
177  bitset_t* pthis = this;
178  *pthis &= static_cast<bitset_t>(rhs);
179  return *this;
180  }
181 
183  template<const index_t LO, const index_t HI>
184  inline
185  auto
187  const index_set<LO,HI>& rhs) -> const
189  {
190  using index_set_t = index_set<LO, HI>;
191  using bitset_t = typename index_set_t::bitset_t;
192  return static_cast<bitset_t>(lhs) & static_cast<bitset_t>(rhs);
193  }
194 
196  template<const index_t LO, const index_t HI>
197  inline
198  auto
201  {
202  bitset_t* pthis = this;
203  *pthis |= static_cast<bitset_t>(rhs);
204  return *this;
205  }
206 
208  template<const index_t LO, const index_t HI>
209  inline
210  auto
212  const index_set<LO,HI>& rhs) -> const
214  {
215  using index_set_t = index_set<LO, HI>;
216  using bitset_t = typename index_set_t::bitset_t;
217  return static_cast<bitset_t>(lhs) | static_cast<bitset_t>(rhs);
218  }
219 
221  template<const index_t LO, const index_t HI>
222  inline
223  auto
226  { return reference(*this, idx); }
227 
229  template<const index_t LO, const index_t HI>
230  inline
231  auto
233  operator[] (const index_t idx) const -> bool
234  { return this->test(idx); }
235 
237  template<const index_t LO, const index_t HI>
238  inline
239  auto
241  test(const index_t idx) const -> bool
242  {
243  // Reference: [JA], 1.2.1
244  return (idx < 0)
245  ? bool(bitset_t::to_ulong() & (1UL << (idx - LO)))
246  : (idx > 0)
247  ? bool(bitset_t::to_ulong() & (1UL << (idx - LO - 1)))
248  : false;
249  }
250 
252  template<const index_t LO, const index_t HI>
253  inline
254  auto
257  {
258  bitset_t::set();
259  return *this;
260  }
261 
263  template<const index_t LO, const index_t HI>
264  inline
265  auto
268  {
269  if (idx > 0)
270  bitset_t::set(idx-LO-1);
271  else if (idx < 0)
272  bitset_t::set(idx-LO);
273  return *this;
274  }
275 
277  template<const index_t LO, const index_t HI>
278  inline
279  auto
281  set(const index_t idx, const int val) -> index_set_t&
282  {
283  if (idx > 0)
284  bitset_t::set(idx-LO-1, val);
285  else if (idx < 0)
286  bitset_t::set(idx-LO, val);
287  return *this;
288  }
289 
291  template<const index_t LO, const index_t HI>
292  inline
293  auto
296  {
297  bitset_t::reset();
298  return *this;
299  }
300 
302  template<const index_t LO, const index_t HI>
303  inline
304  auto
306  reset(const index_t idx) -> index_set_t&
307  {
308  if (idx > 0)
309  bitset_t::reset(idx-LO-1);
310  else if (idx < 0)
311  bitset_t::reset(idx-LO);
312  return *this;
313  }
314 
316  template<const index_t LO, const index_t HI>
317  inline
318  auto
321  {
322  bitset_t::flip();
323  return *this;
324  }
325 
327  template<const index_t LO, const index_t HI>
328  inline
329  auto
331  flip(const index_t idx) -> index_set_t&
332  {
333  if (idx > 0)
334  bitset_t::flip(idx-LO-1);
335  else if (idx < 0)
336  bitset_t::flip(idx-LO);
337  return *this;
338  }
339 
341  template<const index_t LO, const index_t HI>
342  inline
343  auto
345  count() const -> index_t
346  {
347  unsigned long val = bitset_t::to_ulong();
348  // Reference: [JA], 1.3
349  if (val == 0)
350  return 0;
351  else
352  {
353  index_t result = 1;
354  while (val &= val-1)
355  ++result;
356  return result;
357  }
358  }
359 
361  template<const index_t LO, const index_t HI>
362  inline
363  auto
365  count_neg() const -> index_t
366  {
367  static const index_set_t lo_mask = bitset_t((1UL << -LO) - 1UL);
368  const index_set_t neg_part = *this & lo_mask;
369  return neg_part.count();
370  }
371 
373  template<const index_t LO, const index_t HI>
374  inline
375  auto
377  count_pos() const -> index_t
378  {
379  const auto* pthis = static_cast<const bitset_t*>(this);
380  const index_set_t pos_part = *pthis >> -LO;
381  return pos_part.count();
382  }
383 
384 #if (_GLUCAT_BITS_PER_ULONG == 64)
385  template<const index_t LO, const index_t HI>
387  inline
388  auto
390  min() const -> index_t
391  {
392  // Reference: [JA], 1.3
393  unsigned long val = bitset_t::to_ulong();
394  if (val == 0)
395  return 0;
396  else
397  {
398  val -= val & (val-1); // isolate lowest bit
399 
400  index_t idx = 0;
401  const index_t nbits = HI - LO;
402 
403  if (nbits > 8)
404  {
405  if (val & 0xffffffff00000000ul)
406  idx += 32;
407  if (val & 0xffff0000ffff0000ul)
408  idx += 16;
409  if (val & 0xff00ff00ff00ff00ul)
410  idx += 8;
411  }
412  if (val & 0xf0f0f0f0f0f0f0f0ul)
413  idx += 4;
414  if (val & 0xccccccccccccccccul)
415  idx += 2;
416  if (val & 0xaaaaaaaaaaaaaaaaul)
417  idx += 1;
418 
419  return idx + ((idx < -LO) ? LO : LO+1);
420  }
421  }
422 #elif (_GLUCAT_BITS_PER_ULONG == 32)
423  template<const index_t LO, const index_t HI>
425  inline
426  index_t
428  min() const
429  {
430  // Reference: [JA], 1.3
431  unsigned long val = bitset_t::to_ulong();
432  if (val == 0)
433  return 0;
434  else
435  {
436  val -= val & (val-1); // isolate lowest bit
437 
438  index_t idx = 0;
439  const index_t nbits = HI - LO;
440  if (nbits > 8)
441  {
442  if (val & 0xffff0000ul)
443  idx += 16;
444  if (val & 0xff00ff00ul)
445  idx += 8;
446  }
447  if (val & 0xf0f0f0f0ul)
448  idx += 4;
449  if (val & 0xccccccccul)
450  idx += 2;
451  if (val & 0xaaaaaaaaul)
452  idx += 1;
453 
454  return idx + ((idx < -LO) ? LO : LO+1);
455  }
456  }
457 #else
458  template<const index_t LO, const index_t HI>
460  auto
462  min() const -> index_t
463  {
464  for (auto
465  idx = LO;
466  idx != 0;
467  ++idx)
468  if (this->test(idx))
469  return idx;
470  for (auto
471  idx = index_t(1);
472  idx <= HI;
473  ++idx)
474  if (this->test(idx))
475  return idx;
476  return 0;
477  }
478 #endif
479 
480 #if (_GLUCAT_BITS_PER_ULONG == 64)
481  template<const index_t LO, const index_t HI>
483  inline
484  auto
486  max() const -> index_t
487  {
488  // Reference: [JA], 1.6
489  auto val = bitset_t::to_ulong();
490  if (val == 0)
491  return 0;
492  else
493  {
494  auto idx = index_t(0);
495  const auto nbits = HI - LO;
496  if (nbits > 8)
497  {
498  if (val & 0xffffffff00000000ul)
499  { val >>= 32; idx += 32; }
500  if (val & 0x00000000ffff0000ul)
501  { val >>= 16; idx += 16; }
502  if (val & 0x000000000000ff00ul)
503  { val >>= 8; idx += 8; }
504  }
505  if (val & 0x00000000000000f0ul)
506  { val >>= 4; idx += 4; }
507  if (val & 0x000000000000000cul)
508  { val >>= 2; idx += 2; }
509  if (val & 0x0000000000000002ul)
510  { idx += 1; }
511  return idx + ((idx < -LO) ? LO : LO+1);
512  }
513  }
514 #elif (_GLUCAT_BITS_PER_ULONG == 32)
515  template<const index_t LO, const index_t HI>
517  inline
518  auto
520  max() const -> index_t
521  {
522  // Reference: [JA], 1.6
523  auto val = bitset_t::to_ulong();
524  if (val == 0)
525  return 0;
526  else
527  {
528  auto idx = index_t(0);
529  const auto nbits = HI - LO;
530  if (nbits > 8)
531  {
532  if (val & 0xffff0000ul)
533  { val >>= 16; idx += 16; }
534  if (val & 0x0000ff00ul)
535  { val >>= 8; idx += 8; }
536  }
537  if (val & 0x000000f0ul)
538  { val >>= 4; idx += 4; }
539  if (val & 0x0000000cul)
540  { val >>= 2; idx += 2; }
541  if (val & 0x00000002ul)
542  { idx += 1; }
543  return idx + ((idx < -LO) ? LO : LO+1);
544  }
545  }
546 #else
547  template<const index_t LO, const index_t HI>
549  auto
551  max() const -> index_t
552  {
553  for (auto
554  idx = HI;
555  idx != 0;
556  --idx)
557  if (this->test(idx))
558  return idx;
559  for (auto
560  idx = index_t(-1);
561  idx >= LO;
562  --idx)
563  if (this->test(idx))
564  return idx;
565  return 0;
566  }
567 #endif
568 
570  // eg. {3,4,5} is less than {3,7,8}
571  template<const index_t LO, const index_t HI>
572  inline
573  auto
574  compare(const index_set<LO,HI>& a, const index_set<LO,HI>& b) -> int
575  {
576  return (a == b)
577  ? 0
578  : a.lex_less_than(b)
579  ? -1
580  : 1;
581  }
582 
584  // eg. {3,4,5} is less than {3,7,8}
585  template<const index_t LO, const index_t HI>
586  inline
587  auto
589  lex_less_than(const index_set_t rhs) const -> bool
590  { return bitset_t::to_ulong() < rhs.bitset_t::to_ulong(); }
591 
593  // Order by count, then order lexicographically within the equivalence class of count.
594  template<const index_t LO, const index_t HI>
595  inline
596  auto
598  operator< (const index_set_t rhs) const -> bool
599  {
600  const auto this_grade = this->count();
601  const auto rhs_grade = rhs.count();
602  return (this_grade < rhs_grade)
603  ? true
604  : (this_grade > rhs_grade)
605  ? false
606  : this->lex_less_than(rhs);
607  }
608 
610  template<const index_t LO, const index_t HI>
611  auto
612  operator<< (std::ostream& os, const index_set<LO,HI>& ist) -> std::ostream&
613  {
614  index_t i;
615  os << '{';
616  for (i = LO;
617  (i <= HI) && !(ist[i]);
618  ++i)
619  { }
620  if (i <= HI)
621  os << i;
622  for (++i;
623  i <= HI;
624  ++i)
625  if (ist[i])
626  os << ',' << i;
627  os << '}';
628  return os;
629  }
630 
632  template<const index_t LO, const index_t HI>
633  auto
634  operator>> (std::istream& s, index_set<LO,HI>& ist) -> std::istream&
635  {
636  // Parsing variables.
637  auto i = index_t(0);
638  using index_set_t = index_set<LO,HI>;
639  auto local_ist = index_set_t();
640  // Parsing control variables.
641  auto parse_index_list = true;
642  auto expect_closing_brace = false;
643  auto expect_index = false;
644  // Parse an optional opening brace.
645  auto c = s.peek();
646  // If there is a failure or end of file, this ends parsing.
647  if (!s.good())
648  parse_index_list = false;
649  else
650  { // Check for an opening brace.
651  expect_closing_brace = (c == int('{'));
652  if (expect_closing_brace)
653  { // Consume the opening brace.
654  s.get();
655  // The next character may be a closing brace,
656  // indicating the empty index set.
657  c = s.peek();
658  if (s.good() && (c == int('}')))
659  { // A closing brace has been parsed and is no longer expected.
660  expect_closing_brace = false;
661  // Consume the closing brace.
662  s.get();
663  // This ends parsing.
664  parse_index_list = false;
665  }
666  }
667  }
668  if (s.good() && parse_index_list)
669  { // Parse an optional index list.
670  // The index list starts with a first index.
671  for (s >> i;
672  !s.fail();
673  s >> i)
674  { // An index has been parsed. Check to see if it is in range.
675  if ((i < LO) || (i > HI))
676  { // An index out of range is a failure.
677  s.clear(std::istream::failbit);
678  break;
679  }
680  // Add the index to the index set local_ist.
681  local_ist.set(i);
682  // Immediately after parsing an index, an index is no longer expected.
683  expect_index = false;
684  // Reading the index may have resulted in an end of file condition.
685  // If so, this ends the index list.
686  if (s.eof())
687  break;
688  // The index list continues with a comma, and
689  // may be ended by a closing brace, if it was begun with an opening brace.
690  // Parse a possible comma or closing brace.
691  c = s.peek();
692  if (!s.good())
693  break;
694  // First, test for a closing brace, if expected.
695  if (expect_closing_brace && (c == int('}')))
696  { // Consume the closing brace.
697  s.get();
698  // Immediately after parsing the closing brace, it is no longer expected.
699  expect_closing_brace = false;
700  // A closing brace ends the index list.
701  break;
702  }
703  // Now test for a comma.
704  if (c == int(','))
705  { // Consume the comma.
706  s.get();
707  // A index is expected after the comma.
708  expect_index = true;
709  }
710  else
711  { // Any other character here is a failure.
712  s.clear(std::istream::failbit);
713  break;
714  }
715  }
716  }
717  // If an index or a closing brace is still expected, this is a failure.
718  if (expect_index || expect_closing_brace)
719  s.clear(std::istream::failbit);
720  // End of file is not a failure.
721  if (s)
722  { // The index set has been successfully parsed.
723  ist = local_ist;
724  }
725  return s;
726  }
727 
729  template<const index_t LO, const index_t HI>
730  inline
731  auto
733  is_contiguous () const -> bool
734  {
735  const auto min_index = this->min();
736  const auto max_index = this->max();
737  return (min_index < 0 && max_index > 0)
738  ? max_index - min_index == this->count()
739  : (min_index == 1 || max_index == -1) &&
740  (max_index - min_index == this->count() - 1);
741  }
742 
744  template<const index_t LO, const index_t HI>
745  inline
746  auto
748  fold() const -> const
749  index_set<LO,HI>
750  { return this->fold(*this, true); }
751 
753  template<const index_t LO, const index_t HI>
754  auto
756  fold(const index_set_t frm, const bool prechecked) const -> const
758  {
759  if (!prechecked && ((*this | frm) != frm))
760  throw error_t("fold(frm): cannot fold from outside of frame");
761  const auto frm_min = frm.min();
762  const auto frm_max = frm.max();
763  auto result = index_set_t();
764  auto fold_idx = index_t(-1);
765  for (auto
766  unfold_idx = fold_idx;
767  unfold_idx >= frm_min;
768  --unfold_idx)
769  if (frm.test(unfold_idx))
770  // result.set(fold_idx--, this->test(unfold_idx));
771  {
772  if (this->test(unfold_idx))
773  result.set(fold_idx);
774  --fold_idx;
775  }
776  fold_idx = index_t(1);
777  for (auto
778  unfold_idx = fold_idx;
779  unfold_idx <= frm_max;
780  ++unfold_idx)
781  if (frm.test(unfold_idx))
782  // result.set(fold_idx++, this->test(unfold_idx));
783  {
784  if (this->test(unfold_idx))
785  result.set(fold_idx);
786  ++fold_idx;
787  }
788  return result;
789  }
790 
792  template<const index_t LO, const index_t HI>
793  auto
795  unfold(const index_set_t frm, const bool prechecked) const -> const index_set_t
796  {
797  const char* msg =
798  "unfold(frm): cannot unfold into a smaller frame";
799  const auto frm_min = frm.min();
800  const auto frm_max = frm.max();
801  auto result = index_set_t();
802  auto fold_idx = index_t(-1);
803  for (auto
804  unfold_idx = fold_idx;
805  unfold_idx >= frm_min;
806  --unfold_idx)
807  if (frm.test(unfold_idx))
808  if (this->test(fold_idx--))
809  result.set(unfold_idx);
810  if (!prechecked && ((fold_idx+1) > this->min()))
811  throw error_t(msg);
812  fold_idx = index_t(1);
813  for (auto
814  unfold_idx = fold_idx;
815  unfold_idx <= frm_max;
816  ++unfold_idx)
817  if (frm.test(unfold_idx))
818  if (this->test(fold_idx++))
819  result.set(unfold_idx);
820  if (!prechecked && ((fold_idx-1) < this->max()))
821  throw error_t(msg);
822  return result;
823  }
824 
826  template<const index_t LO, const index_t HI>
827  inline
828  auto
831  {
832  const auto min_index = frm.fold().min();
833  if (min_index == 0)
834  return 0;
835  else
836  {
837  const auto folded_set = this->fold(frm);
838  const auto skip = min_index > 0 ? index_t(1) : index_t(0);
839  return folded_set.bitset_t::to_ulong() >> (min_index-LO-skip);
840  }
841  }
842 
844  inline
845  static
846  auto inverse_reversed_gray(unsigned long x) -> unsigned long
847  {
848  // Reference: [JA]
849 #if (_GLUCAT_BITS_PER_ULONG >= 64)
850  x ^= x << 32; // for 64-bit words
851 #endif
852  x ^= x << 16; // reversed_gray ** 16
853  x ^= x << 8; // reversed_gray ** 8
854  x ^= x << 4; // reversed_gray ** 4
855  x ^= x << 2; // reversed_gray ** 2
856  x ^= x << 1; // reversed_gray ** 1
857  return x;
858  }
859 
861  inline
862  static
863  auto inverse_gray(unsigned long x) -> unsigned long
864  {
865  // Reference: [JA]
866 #if (_GLUCAT_BITS_PER_ULONG >= 64)
867  x ^= x >> 32; // for 64-bit words
868 #endif
869  x ^= x >> 16; // gray ** 16
870  x ^= x >> 8; // gray ** 8
871  x ^= x >> 4; // gray ** 4
872  x ^= x >> 2; // gray ** 2
873  x ^= x >> 1; // gray ** 1
874  return x;
875  }
876 
878  template<const index_t LO, const index_t HI>
879  auto
881  sign_of_mult(const index_set_t rhs) const -> int
882  {
883  // Implemented using Walsh functions and Gray codes.
884  // Reference: [L] Chapter 21, 21.3
885  // Reference: [JA]
886  const auto uthis = this->bitset_t::to_ulong();
887  const auto urhs = rhs.bitset_t::to_ulong();
888  const auto nbits = HI - LO;
889  auto negative = 0UL;
890  if (nbits > 8)
891  {
892  // Set h to be the inverse reversed Gray code of rhs.
893  // This sets each bit of h to be the cumulative ^ of
894  // the same and lower bits of rhs.
895  const auto h = inverse_reversed_gray(urhs);
896  // Set k to be the inverse Gray code of *this & h.
897  // This sets the low bit of k to be parity(*this & h).
898  const auto k = inverse_gray(uthis & h);
899  // Set q to be the inverse Gray code of the positive part of *this & rhs.
900  const auto q = inverse_gray((uthis & urhs) >> -LO);
901  negative = k ^ q;
902  }
903  else
904  {
905  auto h = 0UL;
906  for (auto
907  j = index_t(0);
908  j < -LO;
909  ++j)
910  {
911  h ^= urhs >> j;
912  negative ^= h & (uthis >> j);
913  }
914  for (auto
915  j = index_t(-LO);
916  j < nbits;
917  ++j)
918  {
919  negative ^= h & (uthis >> j);
920  h ^= urhs >> j;
921  }
922  }
923  return 1 - int((negative & 1) << 1);
924  }
925 
927  template<const index_t LO, const index_t HI>
928  inline
929  auto
931  sign_of_square() const -> int
932  {
933  auto result = 1 - int((this->count_neg() % 2) << 1);
934  switch (this->count() % 4)
935  {
936  case 2:
937  case 3:
938  result *= -1;
939  break;
940  default:
941  break;
942  }
943  return result;
944  }
945 
947  template<const index_t LO, const index_t HI>
948  inline
949  auto
951  hash_fn() const -> size_t
952  {
953  static const auto lo_mask = (1UL << -LO) - 1UL;
954  const auto uthis = bitset_t::to_ulong();
955  const auto neg_part = uthis & lo_mask;
956  const auto pos_part = uthis >> -LO;
957  return size_t(neg_part ^ pos_part);
958  }
959 
961  inline
962  auto
964  { return (j < 0) ? -1 : 1; }
965 
967  template<const index_t LO, const index_t HI>
968  inline
969  auto
971  { return std::min(ist.min(), 0); }
972 
974  template<const index_t LO, const index_t HI>
975  inline
976  auto
978  { return std::max(ist.max(), 0); }
979 
980 // index_set reference
981 
983  template<const index_t LO, const index_t HI>
984  inline
987  m_pst(&ist),
988  m_idx(idx)
989  { }
990 
992  template<const index_t LO, const index_t HI>
993  inline
994  auto
996  operator== (const reference& c_j) const -> bool
997  { return m_pst == c_j.m_pst && m_idx == c_j.m_idx; }
998 
1000  template<const index_t LO, const index_t HI>
1001  inline
1002  auto
1005  {
1006  if ( x )
1007  m_pst->set(m_idx);
1008  else
1009  m_pst->reset(m_idx);
1010  return *this;
1011  }
1012 
1014  template<const index_t LO, const index_t HI>
1015  inline
1016  auto
1019  {
1020  if (&c_j != this && c_j != *this)
1021  {
1022  if ( (c_j.m_pst)[c_j.m_idx] )
1023  m_pst->set(m_idx);
1024  else
1025  m_pst->reset(m_idx);
1026  }
1027  return *this;
1028  }
1029 
1031  template<const index_t LO, const index_t HI>
1032  inline
1033  auto
1035  operator~ () const -> bool
1036  { return !(m_pst->test(m_idx)); }
1037 
1039  template<const index_t LO, const index_t HI>
1040  inline
1042  operator bool () const
1043  { return m_pst->test(m_idx); }
1044 
1046  template<const index_t LO, const index_t HI>
1047  inline
1048  auto
1051  {
1052  m_pst->flip(m_idx);
1053  return *this;
1054  }
1055 }
1056 #endif // _GLUCAT_INDEX_SET_IMP_H
auto operator[](const index_t idx) const -> bool
Subscripting: Test idx for membership: test value of bit idx.
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 hash_fn() const -> size_t
Hash function.
auto fold() const -> const index_set_t
Fold this index set within itself as a frame.
auto test(const index_t idx) const -> bool
Test idx for membership: test value of bit idx.
auto operator~() const -> index_set_t
Set complement: not.
auto max_pos(const index_set< LO, HI > &ist) -> index_t
Maximum positive index, or 0 if none.
auto operator~() const -> bool
Flips a bit.
auto flip() -> index_set_t &
Set complement, except 0: flip all bits, except 0.
auto min_neg(const index_set< LO, HI > &ist) -> index_t
Minimum negative index, or 0 if none.
auto unfold(const index_set_t frm, const bool prechecked=false) const -> const index_set_t
Unfold this index set within the given frame.
auto count_pos() const -> index_t
Number of positive indices included in set.
auto operator=(const bool x) -> reference &
for b[i] = x;
auto operator^=(const index_set_t rhs) -> index_set_t &
Symmetric set difference: exclusive or.
auto set() -> index_set_t &
Include all indices except 0: set all bits except 0.
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.
auto sign_of_mult(const index_set_t ist) const -> int
Sign of geometric product of two Clifford basis elements.
std::bitset< HI - LO > bitset_t
Definition: index_set.h:81
auto operator|=(const index_set_t rhs) -> index_set_t &
Set union: or.
auto compare(const index_set< LO, HI > &a, const index_set< LO, HI > &b) -> int
"lexicographic compare" eg. {3,4,5} is less than {3,7,8}
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.
std::pair< index_t, index_t > index_pair_t
Definition: index_set.h:85
auto value_of_fold(const index_set_t frm) const -> set_value_t
The set value of the fold of this index set within the given frame.
Specific exception class.
Definition: errors.h:56
auto count_neg() const -> index_t
Number of negative indices included in set.
auto lex_less_than(const index_set_t rhs) const -> bool
Lexicographic ordering of two sets: *this < rhs.
auto operator==(const reference &c_j) const -> bool
for b[i] == c[j];
int index_t
Size of index_t should be enough to represent LO, HI.
Definition: global.h:77
auto max() const -> index_t
Maximum member.
auto operator &=(const index_set_t rhs) -> index_set_t &
Set intersection: and.
Index set class based on std::bitset<> in Gnu standard C++ library.
Definition: index_set.h:45
Index set member reference.
Definition: index_set.h:177
static auto inverse_reversed_gray(unsigned long x) -> unsigned long
Inverse reversed Gray code.
auto operator<(const index_set_t rhs) const -> bool
Less than operator used for comparisons, map, etc.
auto reset() -> index_set_t &
Make set empty: Set all bits to 0.
auto sign_of_square(index_t j) -> int
Square of generator {j}.
auto operator==(const index_set_t rhs) const -> bool
Equality.
reference()=delete
Default constructor is deleted.
auto sign_of_square() const -> int
Sign of geometric square of a Clifford basis element.
static auto inverse_gray(unsigned long x) -> unsigned long
Inverse Gray code.
index_set()=default
Default constructor creates an empty set.
auto operator!=(const index_set_t rhs) const -> bool
Inequality.
auto min() const -> index_t
Minimum member.
static auto classname() -> const std::string
Definition: index_set_imp.h:50
auto flip() -> reference &
for b[i].flip();
auto is_contiguous() const -> bool
Determine if the index set is contiguous, ie. has no gaps.
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.