SymmetryGroup.h

#ifndef UTIL_SYMMETRY_GROUP_H
#define UTIL_SYMMETRY_GROUP_H

#include <iostream>
#include <vector>

namespace Util {


   template <class Symmetry>    
   class SymmetryGroup 
   {

   public:

      SymmetryGroup();

      SymmetryGroup(const SymmetryGroup& other);

      ~SymmetryGroup();

      bool add(Symmetry& symmetry);

      void makeCompleteGroup();

      const Symmetry* find(const Symmetry& symmetry) const;

      const Symmetry& identity() const;

      int size() const;

      SymmetryGroup& operator = (const SymmetryGroup& other);
 
      Symmetry& operator [] (int i);
 
      const Symmetry& operator [] (int i) const;

      bool isValid() const;
 
   private:

      std::vector<Symmetry> elements_;

      Symmetry identity_;

   };


   // Method definitions

   /*
   * Default constructor
   */
   template <class Symmetry>
   SymmetryGroup<Symmetry>::SymmetryGroup()
   {
      identity_ = Symmetry::identity();
      elements_.push_back(identity_);
   }

   /*
   * Copy constructor
   */
   template <class Symmetry>
   SymmetryGroup<Symmetry>::SymmetryGroup(const SymmetryGroup& other)
   {
      identity_ = other.identity();
      for (int i = 0; i < other.size(); ++i) {
         elements_.push_back(other.elements_[i]);
      }
   }

   /*
   * Destructor
   */
   template <class Symmetry>
   SymmetryGroup<Symmetry>::~SymmetryGroup()
   {}

   /*
   * Assignment operator.
   */
   template <class Symmetry>
   SymmetryGroup<Symmetry>& 
   SymmetryGroup<Symmetry>::operator = (const SymmetryGroup& other)
   {
      if (this != &other) {
         identity_ = other.identity();
         elements_.clear();
         for (int i = 0; i < other.size(); ++i) {
            elements_.push_back(other.elements_[i]);
         }
      }
      return *this;
   }

   /*
   * Find an element in the group, return const pointer to its address.
   *
   * Return a null pointer if symmetry is not a member of the group.
   */
   template <class Symmetry>
   const Symmetry* SymmetryGroup<Symmetry>::find(const Symmetry& symmetry) const
   {
      for (int i=0; i < size(); ++i) {
         if (symmetry == elements_[i]) {
            return &(elements_[i]);
         }
      }
      // Return null pointer if not found
      return 0;
   }
 
   /*
   * Add a new element to the group.
   */
   template <class Symmetry>
   bool SymmetryGroup<Symmetry>::add(Symmetry& symmetry)
   {
      const Symmetry* ptr = find(symmetry);
      bool added; 

      if (ptr == 0) {
         elements_.push_back(symmetry);
         added = true;
      } else {
         added = false;
      }
      return added;
   }

   /*
   * Create a complete group.
   */
   template <class Symmetry>
   void SymmetryGroup<Symmetry>::makeCompleteGroup()
   {
      Symmetry a, b, c;
      int      i, j, n;
      bool added, complete;

      // Add all inverses
      n = size();
      for (i = 0; i < n; ++i) {
         a = elements_[i].inverse();
         add(a);
      }

      // Add all products of existing elements, and their inverses
      complete = false;
      while (!complete) {
         complete = true;
         n = size();
         for (i = 0; i < n; ++i) {
            a = elements_[i];
            for (j = 0; j < n; ++j) {
               b = elements_[j];
               c = a*b;
               added = add(c);
               if (added) {
                  complete = false;
               }
               b = c.inverse(); 
               added = add(b);
               if (added) { 
                  complete = false;
               }
            } 
         } 
      } 

   }


   /*
   * Check validity of this group.
   */
   template <class Symmetry>
   bool SymmetryGroup<Symmetry>::isValid() const
   {
      Symmetry a, b, c;
      int      i, j, n;
 
      // Check for inclusion of identity element
      c = Symmetry::identity();
      if (find(c) == 0) {
         UTIL_THROW("Identity element is not in group");
      }

      // Check inverses, uniqueness, and completeness
      n = size();
      for (i = 0; i < n; ++i) {
         a = elements_[i];
         c = a.inverse();
         if (find(c) == 0) {
            UTIL_THROW("Inverse of element not in group");
         }
         for (j = 0; j < n; ++j) {
            b = elements_[j];
            if (i != j) {
               if (a == b) {
                  UTIL_THROW("An element of the group is not unique");
               }
            }
            c = a*b;
            if (find(c) == 0) {
               UTIL_THROW("Product of two element is not in group");
            }
         }
      }

      // If no Exceptions have been thrown, return true
      return true; 

   }

   /*
   * Return the current size of the group.
   */
   template <class Symmetry>
   inline int SymmetryGroup<Symmetry>::size() const
   {  return elements_.size(); }
 
   /*
   * Return identity element.
   */
   template <class Symmetry>
   inline const Symmetry& SymmetryGroup<Symmetry>::identity() const
   {  return identity_; }
 
   /*
   * Element access operator (by reference).
   */
   template <class Symmetry> inline 
   Symmetry& SymmetryGroup<Symmetry>::operator [] (int i)
   { return elements_[i]; }
 
   /*
   * Element access operator (by reference).
   */
   template <class Symmetry> inline 
   const Symmetry& SymmetryGroup<Symmetry>::operator [] (int i) const
   { return elements_[i]; }
 
}
#endif