diff --git a/.gitignore b/.gitignore
index 3f6ea78e..d30b5fff 100644
--- a/.gitignore
+++ b/.gitignore
@@ -112,3 +112,8 @@ m4/libtool.m4
 buck-out
 BUCK
 make-bin-BUCK.sh
+
+# generated sources #
+#####################
+lib/qcd/spin/gamma-gen/*.h
+lib/qcd/spin/gamma-gen/*.cc
\ No newline at end of file
diff --git a/benchmarks/Benchmark_dwf.cc b/benchmarks/Benchmark_dwf.cc
index 10e4521b..753cd60f 100644
--- a/benchmarks/Benchmark_dwf.cc
+++ b/benchmarks/Benchmark_dwf.cc
@@ -37,11 +37,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 typedef WilsonFermion5D<DomainWallVec5dImplR> WilsonFermion5DR;
@@ -321,7 +321,7 @@ int main (int argc, char ** argv)
     ref = zero;
     for(int mu=0;mu<Nd;mu++){
 
-      //    ref =  src - Gamma(Gamma::GammaX)* src ; // 1+gamma_x
+      //    ref =  src - Gamma(Gamma::Algebra::GammaX)* src ; // 1+gamma_x
       tmp = U[mu]*Cshift(src,mu+1,1);
       for(int i=0;i<ref._odata.size();i++){
   ref._odata[i]+= tmp._odata[i] + Gamma(Gmu[mu])*tmp._odata[i]; ;
diff --git a/benchmarks/Benchmark_dwf_sweep.cc b/benchmarks/Benchmark_dwf_sweep.cc
index 21db823d..37e47062 100644
--- a/benchmarks/Benchmark_dwf_sweep.cc
+++ b/benchmarks/Benchmark_dwf_sweep.cc
@@ -37,11 +37,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 void benchDw(std::vector<int> & L, int Ls, int threads, int report =0 );
diff --git a/benchmarks/Benchmark_wilson.cc b/benchmarks/Benchmark_wilson.cc
index 4930713c..11f1353f 100644
--- a/benchmarks/Benchmark_wilson.cc
+++ b/benchmarks/Benchmark_wilson.cc
@@ -37,11 +37,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 bool overlapComms = false;
@@ -106,7 +106,7 @@ int main (int argc, char ** argv)
   { // Naive wilson implementation
     ref = zero;
     for(int mu=0;mu<Nd;mu++){
-      //    ref =  src + Gamma(Gamma::GammaX)* src ; // 1-gamma_x
+      //    ref =  src + Gamma(Gamma::Algebra::GammaX)* src ; // 1-gamma_x
       tmp = U[mu]*Cshift(src,mu,1);
       for(int i=0;i<ref._odata.size();i++){
 	ref._odata[i]+= tmp._odata[i] - Gamma(Gmu[mu])*tmp._odata[i]; ;
@@ -159,7 +159,7 @@ int main (int argc, char ** argv)
     ref = zero;
     for(int mu=0;mu<Nd;mu++){
 
-      //    ref =  src - Gamma(Gamma::GammaX)* src ; // 1+gamma_x
+      //    ref =  src - Gamma(Gamma::Algebra::GammaX)* src ; // 1+gamma_x
       tmp = U[mu]*Cshift(src,mu,1);
       for(int i=0;i<ref._odata.size();i++){
 	ref._odata[i]+= tmp._odata[i] + Gamma(Gmu[mu])*tmp._odata[i]; ;
diff --git a/benchmarks/Benchmark_wilson_sweep.cc b/benchmarks/Benchmark_wilson_sweep.cc
index 96e5b5e4..aae34f6c 100644
--- a/benchmarks/Benchmark_wilson_sweep.cc
+++ b/benchmarks/Benchmark_wilson_sweep.cc
@@ -30,11 +30,11 @@ struct scal {
   d internal;
 };
 
-Gamma::GammaMatrix Gmu [] = {
-  Gamma::GammaX,
-  Gamma::GammaY,
-  Gamma::GammaZ,
-  Gamma::GammaT
+Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
 };
 
 bool overlapComms = false;
diff --git a/extras/Hadrons/Modules/MContraction/Meson.hpp b/extras/Hadrons/Modules/MContraction/Meson.hpp
index 3a8c2533..270875a6 100644
--- a/extras/Hadrons/Modules/MContraction/Meson.hpp
+++ b/extras/Hadrons/Modules/MContraction/Meson.hpp
@@ -37,17 +37,17 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Grid/Hadrons/ModuleFactory.hpp>
 
 namespace Grid {
-    // Overload >> to extract gamma pair from "[g1 g2]" string.
+    // Overload >> to extract gamma pair from "<g1 g2>" string.
     template <typename T1, typename T2>
     inline std::istringstream &operator>>(std::istringstream &sstr,
                                           std::pair<T1, T2> &buf)
     {
-        T1 buf1;
-        T2 buf2;
+        unsigned int buf1;
+        unsigned int buf2;
         char c;
         sstr >> c >> buf1 >> buf2 >> c;
         sstr.peek();
-        buf = std::make_pair(buf1, buf2);
+        buf = std::make_pair((T1)buf1, (T2)buf2);
         return sstr;
     }
 }
@@ -76,7 +76,7 @@ BEGIN_HADRONS_NAMESPACE
  ******************************************************************************/
 BEGIN_MODULE_NAMESPACE(MContraction)
 
-typedef std::pair<unsigned int, unsigned int> GammaPair;
+typedef std::pair<Gamma::Algebra, Gamma::Algebra> GammaPair;
 
 class MesonPar: Serializable
 {
@@ -99,8 +99,8 @@ public:
     {
     public:
         GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
-                                        unsigned int, gamma_snk,
-                                        unsigned int, gamma_src,
+                                        Gamma::Algebra, gamma_snk,
+                                        Gamma::Algebra, gamma_src,
                                         std::vector<Complex>, corr);
     };
 public:
@@ -151,13 +151,14 @@ void TMeson<FImpl1, FImpl2>::parseGammaString(std::vector<GammaPair> &gammaList)
     if (par().gammas.compare("all") == 0)
     {
         // Do all contractions.
-        unsigned int n_gam = Ns*Ns;
+        unsigned int n_gam = Ns * Ns;
         gammaList.resize(n_gam*n_gam);
-        for (unsigned int i = 0; i < n_gam; ++i)
+        for (unsigned int i = 1; i < Gamma::nGamma; i += 2)
         {
-            for (unsigned int j = 0; j < n_gam; ++j)
+            for (unsigned int j = 1; j < Gamma::nGamma; j += 2)
             {
-                gammaList.push_back(std::make_pair(i, j));
+                gammaList.push_back(std::make_pair((Gamma::Algebra)i, 
+                                                   (Gamma::Algebra)j));
             }
         }
     }
@@ -181,7 +182,7 @@ void TMeson<FImpl1, FImpl2>::execute(void)
     PropagatorField1       &q1 = *env().template getObject<PropagatorField1>(par().q1);
     PropagatorField2       &q2 = *env().template getObject<PropagatorField2>(par().q2);
     LatticeComplex         c(env().getGrid());
-    SpinMatrix             g[Ns*Ns], g5;
+    Gamma                  g5(Gamma::Algebra::Gamma5);
     std::vector<GammaPair> gammaList;
     std::vector<TComplex>  buf;
     std::vector<Result>    result;
@@ -198,17 +199,14 @@ void TMeson<FImpl1, FImpl2>::execute(void)
     }
     ph = exp(-2*M_PI*i*ph);
     
-    g5 = makeGammaProd(Ns*Ns - 1);
-    for (int i = 0; i < Ns*Ns; ++i)
-    {
-        g[i] = makeGammaProd(i);
-    }
     parseGammaString(gammaList);
 
     result.resize(gammaList.size());
     for (unsigned int i = 0; i < result.size(); ++i)
     {
-        c = trace(g[gammaList[i].first]*q1*g[gammaList[i].second]*g5*adj(q2)*g5*ph);
+        Gamma gSnk(gammaList[i].first);
+        Gamma gSrc(gammaList[i].second);
+        c = trace((g5*gSnk)*q1*(gSrc*g5)*adj(q2))*ph;
         sliceSum(c, buf, Tp);
 
         result[i].gamma_snk = gammaList[i].first;
diff --git a/extras/Hadrons/Modules/MSource/SeqGamma.hpp b/extras/Hadrons/Modules/MSource/SeqGamma.hpp
index 814da273..fab3751f 100644
--- a/extras/Hadrons/Modules/MSource/SeqGamma.hpp
+++ b/extras/Hadrons/Modules/MSource/SeqGamma.hpp
@@ -61,11 +61,11 @@ class SeqGammaPar: Serializable
 {
 public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(SeqGammaPar,
-                                    std::string,  q,
-                                    unsigned int, tA,
-                                    unsigned int, tB,
-                                    unsigned int, gamma,
-                                    std::string,  mom);
+                                    std::string,    q,
+                                    unsigned int,   tA,
+                                    unsigned int,   tB,
+                                    Gamma::Algebra, gamma,
+                                    std::string,    mom);
 };
 
 template <typename FImpl>
@@ -141,11 +141,10 @@ void TSeqGamma<FImpl>::execute(void)
     PropagatorField &q   = *env().template getObject<PropagatorField>(par().q);
     Lattice<iScalar<vInteger>> t(env().getGrid());
     LatticeComplex             ph(env().getGrid()), coor(env().getGrid());
-    SpinMatrix                 g;
+    Gamma                      g(par().gamma);
     std::vector<Real>          p;
     Complex                    i(0.0,1.0);
     
-    g  = makeGammaProd(par().gamma);
     p  = strToVec<Real>(par().mom);
     ph = zero;
     for(unsigned int mu = 0; mu < env().getNd(); mu++)
@@ -155,7 +154,7 @@ void TSeqGamma<FImpl>::execute(void)
     }
     ph = exp(2*M_PI*i*ph);
     LatticeCoordinate(t, Tp);
-    src = where((t >= par().tA) and (t <= par().tB), g*ph*q, 0.*q);
+    src = where((t >= par().tA) and (t <= par().tB), ph*(g*q), 0.*q);
 }
 
 END_MODULE_NAMESPACE
diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc
index e2b02290..b447f914 100644
--- a/lib/qcd/action/fermion/WilsonFermion.cc
+++ b/lib/qcd/action/fermion/WilsonFermion.cc
@@ -149,11 +149,11 @@ void WilsonFermion<Impl>::MeooeDag(const FermionField &in, FermionField &out) {
 
     typedef Lattice<iSinglet<vector_type> > LatComplex;
 
-    Gamma::GammaMatrix Gmu [] = {
-      Gamma::GammaX,
-      Gamma::GammaY,
-      Gamma::GammaZ,
-      Gamma::GammaT
+    Gamma::Algebra Gmu [] = {
+      Gamma::Algebra::GammaX,
+      Gamma::Algebra::GammaY,
+      Gamma::Algebra::GammaZ,
+      Gamma::Algebra::GammaT
     };
 
     std::vector<int> latt_size   = _grid->_fdimensions;
diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc
index 252fe368..a82aaa2d 100644
--- a/lib/qcd/action/fermion/WilsonFermion5D.cc
+++ b/lib/qcd/action/fermion/WilsonFermion5D.cc
@@ -503,11 +503,11 @@ void WilsonFermion5D<Impl>::MomentumSpacePropagatorHt(FermionField &out,const Fe
   typedef iSinglet<ScalComplex> Tcomplex;
   typedef Lattice<iSinglet<vector_type> > LatComplex;
   
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
   std::vector<int> latt_size   = _grid->_fdimensions;
@@ -574,11 +574,11 @@ void WilsonFermion5D<Impl>::MomentumSpacePropagatorHt(FermionField &out,const Fe
 template<class Impl>
 void WilsonFermion5D<Impl>::MomentumSpacePropagatorHw(FermionField &out,const FermionField &in,RealD mass) 
 {
-    Gamma::GammaMatrix Gmu [] = {
-      Gamma::GammaX,
-      Gamma::GammaY,
-      Gamma::GammaZ,
-      Gamma::GammaT
+    Gamma::Algebra Gmu [] = {
+      Gamma::Algebra::GammaX,
+      Gamma::Algebra::GammaY,
+      Gamma::Algebra::GammaZ,
+      Gamma::Algebra::GammaT
     };
 
     GridBase *_grid = _FourDimGrid;
diff --git a/lib/qcd/action/fermion/g5HermitianLinop.h b/lib/qcd/action/fermion/g5HermitianLinop.h
index af23c36f..cca7a113 100644
--- a/lib/qcd/action/fermion/g5HermitianLinop.h
+++ b/lib/qcd/action/fermion/g5HermitianLinop.h
@@ -80,7 +80,7 @@ class Gamma5HermitianLinearOperator : public LinearOperatorBase<Field> {
   Matrix &_Mat;
   Gamma g5;
 public:
-    Gamma5HermitianLinearOperator(Matrix &Mat): _Mat(Mat), g5(Gamma::Gamma5) {};
+    Gamma5HermitianLinearOperator(Matrix &Mat): _Mat(Mat), g5(Gamma::Algebra::Gamma5) {};
   void Op     (const Field &in, Field &out){
     HermOp(in,out);
   }
diff --git a/lib/qcd/spin/Dirac.cc b/lib/qcd/spin/Dirac.cc
index 3e0c4b9b..a0e1086c 100644
--- a/lib/qcd/spin/Dirac.cc
+++ b/lib/qcd/spin/Dirac.cc
@@ -1,95 +1,72 @@
-    /*************************************************************************************
+/*************************************************************************************
 
-    Grid physics library, www.github.com/paboyle/Grid 
+Grid physics library, www.github.com/paboyle/Grid 
 
-    Source file: ./lib/qcd/spin/Dirac.cc
+Source file: lib/qcd/spin/Dirac.cc
 
-    Copyright (C) 2015
+Copyright (C) 2015
+Copyright (C) 2016
 
+Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: paboyle <paboyle@ph.ed.ac.uk>
 
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
 
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
 
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
-    See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-    /*  END LEGAL */
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
 #include <Grid/Grid.h>
 
 namespace Grid {
+namespace QCD {
 
-  namespace QCD {
+#include "GammaMulTable.h"
 
-    Gamma::GammaMatrix  Gamma::GammaMatrices [] = {
-      Gamma::Identity,
-      Gamma::GammaX,
-      Gamma::GammaY,
-      Gamma::GammaZ,
-      Gamma::GammaT,
-      Gamma::Gamma5,
-      Gamma::MinusIdentity,
-      Gamma::MinusGammaX,
-      Gamma::MinusGammaY,
-      Gamma::MinusGammaZ,
-      Gamma::MinusGammaT,
-      Gamma::MinusGamma5
-    };
-    const char *Gamma::GammaMatrixNames[] = { 
-      "Identity ",
-      "GammaX   ",
-      "GammaY   ",
-      "GammaZ   ",
-      "GammaT   ",
-      "Gamma5   ",
-      "-Identity",
-      "-GammaX  ",
-      "-GammaY  ",
-      "-GammaZ  ",
-      "-GammaT  ",
-      "-Gamma5  ",
-      "         "
-    };
-    
-    SpinMatrix makeGammaProd(const unsigned int i)
-    {
-      SpinMatrix g;
-      
-      g = 1.;
-      if (i & 0x1)
-      {
-        g = g*Gamma(Gamma::GammaMatrix::GammaX);
-      }
-      if (i & 0x2)
-      {
-        g = g*Gamma(Gamma::GammaMatrix::GammaY);
-      }
-      if (i & 0x4)
-      {
-        g = g*Gamma(Gamma::GammaMatrix::GammaZ);
-      }
-      if (i & 0x8)
-      {
-        g = g*Gamma(Gamma::GammaMatrix::GammaT);
-      }
-      
-      return g;
-    }
+const std::array<const char *, Gamma::nGamma> Gamma::name = {{
+  "-Gamma5      ",
+  "Gamma5       ",
+  "-GammaT      ",
+  "GammaT       ",
+  "-GammaTGamma5",
+  "GammaTGamma5 ",
+  "-GammaX      ",
+  "GammaX       ",
+  "-GammaXGamma5",
+  "GammaXGamma5 ",
+  "-GammaY      ",
+  "GammaY       ",
+  "-GammaYGamma5",
+  "GammaYGamma5 ",
+  "-GammaZ      ",
+  "GammaZ       ",
+  "-GammaZGamma5",
+  "GammaZGamma5 ",
+  "-Identity    ",
+  "Identity     ",
+  "-SigmaXT     ",
+  "SigmaXT      ",
+  "-SigmaXY     ",
+  "SigmaXY      ",
+  "-SigmaXZ     ",
+  "SigmaXZ      ",
+  "-SigmaYT     ",
+  "SigmaYT      ",
+  "-SigmaYZ     ",
+  "SigmaYZ      ",
+  "-SigmaZT     ",
+  "SigmaZT      "}};
 
-    //    void sprojMul( vHalfSpinColourVector &out,vColourMatrix &u, vSpinColourVector &in){
-    //      vHalfSpinColourVector hspin;
-    //      spProjXp(hspin,in);
-    //      mult(&out,&u,&hspin);
-    //    }
-  }
-}
+}}
diff --git a/lib/qcd/spin/Dirac.h b/lib/qcd/spin/Dirac.h
index 00bd519c..6053f935 100644
--- a/lib/qcd/spin/Dirac.h
+++ b/lib/qcd/spin/Dirac.h
@@ -1,628 +1,1452 @@
-    /*************************************************************************************
+/*************************************************************************************
 
-    Grid physics library, www.github.com/paboyle/Grid 
+Grid physics library, www.github.com/paboyle/Grid 
 
-    Source file: ./lib/qcd/spin/Dirac.h
+Source file: lib/qcd/spin/Dirac.h
 
-    Copyright (C) 2015
+Copyright (C) 2015
+Copyright (C) 2016
 
+Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
+Author: paboyle <paboyle@ph.ed.ac.uk>
 
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
 
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
 
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
-    See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-    /*  END LEGAL */
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
 #ifndef GRID_QCD_DIRAC_H
 #define GRID_QCD_DIRAC_H
-namespace Grid{
 
+// code generated by the Mathematica notebook gamma-gen/gamma-gen.nb
+#include <array>
+
+namespace Grid {
 namespace QCD {
 
-
-  class Gamma {
-
+class Gamma {
   public:
+    GRID_SERIALIZABLE_ENUM(Algebra, undef,
+                           MinusGamma5      , 0,
+                           Gamma5           , 1,
+                           MinusGammaT      , 2,
+                           GammaT           , 3,
+                           MinusGammaTGamma5, 4,
+                           GammaTGamma5     , 5,
+                           MinusGammaX      , 6,
+                           GammaX           , 7,
+                           MinusGammaXGamma5, 8,
+                           GammaXGamma5     , 9,
+                           MinusGammaY      , 10,
+                           GammaY           , 11,
+                           MinusGammaYGamma5, 12,
+                           GammaYGamma5     , 13,
+                           MinusGammaZ      , 14,
+                           GammaZ           , 15,
+                           MinusGammaZGamma5, 16,
+                           GammaZGamma5     , 17,
+                           MinusIdentity    , 18,
+                           Identity         , 19,
+                           MinusSigmaXT     , 20,
+                           SigmaXT          , 21,
+                           MinusSigmaXY     , 22,
+                           SigmaXY          , 23,
+                           MinusSigmaXZ     , 24,
+                           SigmaXZ          , 25,
+                           MinusSigmaYT     , 26,
+                           SigmaYT          , 27,
+                           MinusSigmaYZ     , 28,
+                           SigmaYZ          , 29,
+                           MinusSigmaZT     , 30,
+                           SigmaZT          , 31);
+    static constexpr unsigned int nGamma = 32;
+    static const std::array<const char *, nGamma>                name;
+    static const std::array<std::array<Algebra, nGamma>, nGamma> mul;
+    static const std::array<Algebra, nGamma>                     adj;
+    Algebra                                                      g;
+  public:
+    Gamma(Algebra initg): g(initg) {}  
+};
 
-    const int Ns=4;
-    
-    enum GammaMatrix {
-      Identity,                        
-      GammaX, 
-      GammaY, 
-      GammaZ, 
-      GammaT,  
-      Gamma5,
-      MinusIdentity,                        
-      MinusGammaX, 
-      MinusGammaY, 
-      MinusGammaZ, 
-      MinusGammaT,  
-      MinusGamma5
-      //      GammaXGamma5, // Rest are composite (willing to take hit for two calls sequentially)
-      //      GammaYGamma5, // as they are less commonly used.
-      //      GammaZGamma5, 
-      //      GammaTGamma5,  
-      //      SigmaXY,
-      //      SigmaXZ,
-      //      SigmaYZ,
-      //      SigmaXT,
-      //      SigmaYT,
-      //      SigmaZT,
-      //      MinusGammaXGamma5, easiest to form by composition
-      //      MinusGammaYGamma5, as performance is not critical for these
-      //      MinusGammaZGamma5, 
-      //      MinusGammaTGamma5,  
-      //      MinusSigmaXY,
-      //      MinusSigmaXZ,
-      //      MinusSigmaYZ,
-      //      MinusSigmaXT,
-      //      MinusSigmaYT,
-      //      MinusSigmaZT
-    };
-    
-    static GammaMatrix GammaMatrices[];
-    static const char *GammaMatrixNames[];
+inline Gamma operator*(const Gamma &g1, const Gamma &g2)
+{
+  return Gamma(Gamma::mul[g1.g][g2.g]);
+}
 
-    Gamma (GammaMatrix g) { _g=g; }
+template<class vtype>
+inline void multMinusGamma5(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = -rhs(0);
+  ret(1) = -rhs(1);
+  ret(2) = rhs(2);
+  ret(3) = rhs(3);
+};
 
-    GammaMatrix _g;
+template<class vtype>
+inline void lmultMinusGamma5(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = -rhs(0, i);
+    ret(1, i) = -rhs(1, i);
+    ret(2, i) = rhs(2, i);
+    ret(3, i) = rhs(3, i);
+  }
+};
 
-  };
+template<class vtype>
+inline void rmultMinusGamma5(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = -rhs(i, 0);
+    ret(i, 1) = -rhs(i, 1);
+    ret(i, 2) = rhs(i, 2);
+    ret(i, 3) = rhs(i, 3);
+  }
+};
+
+template<class vtype>
+inline void multGamma5(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = rhs(0);
+  ret(1) = rhs(1);
+  ret(2) = -rhs(2);
+  ret(3) = -rhs(3);
+};
+
+template<class vtype>
+inline void lmultGamma5(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = rhs(0, i);
+    ret(1, i) = rhs(1, i);
+    ret(2, i) = -rhs(2, i);
+    ret(3, i) = -rhs(3, i);
+  }
+};
+
+template<class vtype>
+inline void rmultGamma5(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = rhs(i, 0);
+    ret(i, 1) = rhs(i, 1);
+    ret(i, 2) = -rhs(i, 2);
+    ret(i, 3) = -rhs(i, 3);
+  }
+};
+
+template<class vtype>
+inline void multMinusGammaT(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = -rhs(2);
+  ret(1) = -rhs(3);
+  ret(2) = -rhs(0);
+  ret(3) = -rhs(1);
+};
+
+template<class vtype>
+inline void lmultMinusGammaT(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = -rhs(2, i);
+    ret(1, i) = -rhs(3, i);
+    ret(2, i) = -rhs(0, i);
+    ret(3, i) = -rhs(1, i);
+  }
+};
+
+template<class vtype>
+inline void rmultMinusGammaT(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = -rhs(i, 2);
+    ret(i, 1) = -rhs(i, 3);
+    ret(i, 2) = -rhs(i, 0);
+    ret(i, 3) = -rhs(i, 1);
+  }
+};
+
+template<class vtype>
+inline void multGammaT(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = rhs(2);
+  ret(1) = rhs(3);
+  ret(2) = rhs(0);
+  ret(3) = rhs(1);
+};
+
+template<class vtype>
+inline void lmultGammaT(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = rhs(2, i);
+    ret(1, i) = rhs(3, i);
+    ret(2, i) = rhs(0, i);
+    ret(3, i) = rhs(1, i);
+  }
+};
+
+template<class vtype>
+inline void rmultGammaT(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = rhs(i, 2);
+    ret(i, 1) = rhs(i, 3);
+    ret(i, 2) = rhs(i, 0);
+    ret(i, 3) = rhs(i, 1);
+  }
+};
+
+template<class vtype>
+inline void multMinusGammaTGamma5(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = rhs(2);
+  ret(1) = rhs(3);
+  ret(2) = -rhs(0);
+  ret(3) = -rhs(1);
+};
+
+template<class vtype>
+inline void lmultMinusGammaTGamma5(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = rhs(2, i);
+    ret(1, i) = rhs(3, i);
+    ret(2, i) = -rhs(0, i);
+    ret(3, i) = -rhs(1, i);
+  }
+};
+
+template<class vtype>
+inline void rmultMinusGammaTGamma5(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = -rhs(i, 2);
+    ret(i, 1) = -rhs(i, 3);
+    ret(i, 2) = rhs(i, 0);
+    ret(i, 3) = rhs(i, 1);
+  }
+};
+
+template<class vtype>
+inline void multGammaTGamma5(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = -rhs(2);
+  ret(1) = -rhs(3);
+  ret(2) = rhs(0);
+  ret(3) = rhs(1);
+};
+
+template<class vtype>
+inline void lmultGammaTGamma5(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = -rhs(2, i);
+    ret(1, i) = -rhs(3, i);
+    ret(2, i) = rhs(0, i);
+    ret(3, i) = rhs(1, i);
+  }
+};
+
+template<class vtype>
+inline void rmultGammaTGamma5(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = rhs(i, 2);
+    ret(i, 1) = rhs(i, 3);
+    ret(i, 2) = -rhs(i, 0);
+    ret(i, 3) = -rhs(i, 1);
+  }
+};
+
+template<class vtype>
+inline void multMinusGammaX(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = timesMinusI(rhs(3));
+  ret(1) = timesMinusI(rhs(2));
+  ret(2) = timesI(rhs(1));
+  ret(3) = timesI(rhs(0));
+};
+
+template<class vtype>
+inline void lmultMinusGammaX(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = timesMinusI(rhs(3, i));
+    ret(1, i) = timesMinusI(rhs(2, i));
+    ret(2, i) = timesI(rhs(1, i));
+    ret(3, i) = timesI(rhs(0, i));
+  }
+};
+
+template<class vtype>
+inline void rmultMinusGammaX(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = timesI(rhs(i, 3));
+    ret(i, 1) = timesI(rhs(i, 2));
+    ret(i, 2) = timesMinusI(rhs(i, 1));
+    ret(i, 3) = timesMinusI(rhs(i, 0));
+  }
+};
+
+template<class vtype>
+inline void multGammaX(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = timesI(rhs(3));
+  ret(1) = timesI(rhs(2));
+  ret(2) = timesMinusI(rhs(1));
+  ret(3) = timesMinusI(rhs(0));
+};
+
+template<class vtype>
+inline void lmultGammaX(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = timesI(rhs(3, i));
+    ret(1, i) = timesI(rhs(2, i));
+    ret(2, i) = timesMinusI(rhs(1, i));
+    ret(3, i) = timesMinusI(rhs(0, i));
+  }
+};
+
+template<class vtype>
+inline void rmultGammaX(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = timesMinusI(rhs(i, 3));
+    ret(i, 1) = timesMinusI(rhs(i, 2));
+    ret(i, 2) = timesI(rhs(i, 1));
+    ret(i, 3) = timesI(rhs(i, 0));
+  }
+};
+
+template<class vtype>
+inline void multMinusGammaXGamma5(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = timesI(rhs(3));
+  ret(1) = timesI(rhs(2));
+  ret(2) = timesI(rhs(1));
+  ret(3) = timesI(rhs(0));
+};
+
+template<class vtype>
+inline void lmultMinusGammaXGamma5(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = timesI(rhs(3, i));
+    ret(1, i) = timesI(rhs(2, i));
+    ret(2, i) = timesI(rhs(1, i));
+    ret(3, i) = timesI(rhs(0, i));
+  }
+};
+
+template<class vtype>
+inline void rmultMinusGammaXGamma5(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = timesI(rhs(i, 3));
+    ret(i, 1) = timesI(rhs(i, 2));
+    ret(i, 2) = timesI(rhs(i, 1));
+    ret(i, 3) = timesI(rhs(i, 0));
+  }
+};
+
+template<class vtype>
+inline void multGammaXGamma5(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = timesMinusI(rhs(3));
+  ret(1) = timesMinusI(rhs(2));
+  ret(2) = timesMinusI(rhs(1));
+  ret(3) = timesMinusI(rhs(0));
+};
+
+template<class vtype>
+inline void lmultGammaXGamma5(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = timesMinusI(rhs(3, i));
+    ret(1, i) = timesMinusI(rhs(2, i));
+    ret(2, i) = timesMinusI(rhs(1, i));
+    ret(3, i) = timesMinusI(rhs(0, i));
+  }
+};
+
+template<class vtype>
+inline void rmultGammaXGamma5(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = timesMinusI(rhs(i, 3));
+    ret(i, 1) = timesMinusI(rhs(i, 2));
+    ret(i, 2) = timesMinusI(rhs(i, 1));
+    ret(i, 3) = timesMinusI(rhs(i, 0));
+  }
+};
+
+template<class vtype>
+inline void multMinusGammaY(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = rhs(3);
+  ret(1) = -rhs(2);
+  ret(2) = -rhs(1);
+  ret(3) = rhs(0);
+};
+
+template<class vtype>
+inline void lmultMinusGammaY(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = rhs(3, i);
+    ret(1, i) = -rhs(2, i);
+    ret(2, i) = -rhs(1, i);
+    ret(3, i) = rhs(0, i);
+  }
+};
+
+template<class vtype>
+inline void rmultMinusGammaY(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = rhs(i, 3);
+    ret(i, 1) = -rhs(i, 2);
+    ret(i, 2) = -rhs(i, 1);
+    ret(i, 3) = rhs(i, 0);
+  }
+};
+
+template<class vtype>
+inline void multGammaY(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = -rhs(3);
+  ret(1) = rhs(2);
+  ret(2) = rhs(1);
+  ret(3) = -rhs(0);
+};
+
+template<class vtype>
+inline void lmultGammaY(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = -rhs(3, i);
+    ret(1, i) = rhs(2, i);
+    ret(2, i) = rhs(1, i);
+    ret(3, i) = -rhs(0, i);
+  }
+};
+
+template<class vtype>
+inline void rmultGammaY(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = -rhs(i, 3);
+    ret(i, 1) = rhs(i, 2);
+    ret(i, 2) = rhs(i, 1);
+    ret(i, 3) = -rhs(i, 0);
+  }
+};
+
+template<class vtype>
+inline void multMinusGammaYGamma5(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = -rhs(3);
+  ret(1) = rhs(2);
+  ret(2) = -rhs(1);
+  ret(3) = rhs(0);
+};
+
+template<class vtype>
+inline void lmultMinusGammaYGamma5(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = -rhs(3, i);
+    ret(1, i) = rhs(2, i);
+    ret(2, i) = -rhs(1, i);
+    ret(3, i) = rhs(0, i);
+  }
+};
+
+template<class vtype>
+inline void rmultMinusGammaYGamma5(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = rhs(i, 3);
+    ret(i, 1) = -rhs(i, 2);
+    ret(i, 2) = rhs(i, 1);
+    ret(i, 3) = -rhs(i, 0);
+  }
+};
+
+template<class vtype>
+inline void multGammaYGamma5(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = rhs(3);
+  ret(1) = -rhs(2);
+  ret(2) = rhs(1);
+  ret(3) = -rhs(0);
+};
+
+template<class vtype>
+inline void lmultGammaYGamma5(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = rhs(3, i);
+    ret(1, i) = -rhs(2, i);
+    ret(2, i) = rhs(1, i);
+    ret(3, i) = -rhs(0, i);
+  }
+};
+
+template<class vtype>
+inline void rmultGammaYGamma5(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = -rhs(i, 3);
+    ret(i, 1) = rhs(i, 2);
+    ret(i, 2) = -rhs(i, 1);
+    ret(i, 3) = rhs(i, 0);
+  }
+};
+
+template<class vtype>
+inline void multMinusGammaZ(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = timesMinusI(rhs(2));
+  ret(1) = timesI(rhs(3));
+  ret(2) = timesI(rhs(0));
+  ret(3) = timesMinusI(rhs(1));
+};
+
+template<class vtype>
+inline void lmultMinusGammaZ(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = timesMinusI(rhs(2, i));
+    ret(1, i) = timesI(rhs(3, i));
+    ret(2, i) = timesI(rhs(0, i));
+    ret(3, i) = timesMinusI(rhs(1, i));
+  }
+};
+
+template<class vtype>
+inline void rmultMinusGammaZ(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = timesI(rhs(i, 2));
+    ret(i, 1) = timesMinusI(rhs(i, 3));
+    ret(i, 2) = timesMinusI(rhs(i, 0));
+    ret(i, 3) = timesI(rhs(i, 1));
+  }
+};
+
+template<class vtype>
+inline void multGammaZ(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = timesI(rhs(2));
+  ret(1) = timesMinusI(rhs(3));
+  ret(2) = timesMinusI(rhs(0));
+  ret(3) = timesI(rhs(1));
+};
+
+template<class vtype>
+inline void lmultGammaZ(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = timesI(rhs(2, i));
+    ret(1, i) = timesMinusI(rhs(3, i));
+    ret(2, i) = timesMinusI(rhs(0, i));
+    ret(3, i) = timesI(rhs(1, i));
+  }
+};
+
+template<class vtype>
+inline void rmultGammaZ(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = timesMinusI(rhs(i, 2));
+    ret(i, 1) = timesI(rhs(i, 3));
+    ret(i, 2) = timesI(rhs(i, 0));
+    ret(i, 3) = timesMinusI(rhs(i, 1));
+  }
+};
+
+template<class vtype>
+inline void multMinusGammaZGamma5(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = timesI(rhs(2));
+  ret(1) = timesMinusI(rhs(3));
+  ret(2) = timesI(rhs(0));
+  ret(3) = timesMinusI(rhs(1));
+};
+
+template<class vtype>
+inline void lmultMinusGammaZGamma5(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = timesI(rhs(2, i));
+    ret(1, i) = timesMinusI(rhs(3, i));
+    ret(2, i) = timesI(rhs(0, i));
+    ret(3, i) = timesMinusI(rhs(1, i));
+  }
+};
+
+template<class vtype>
+inline void rmultMinusGammaZGamma5(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = timesI(rhs(i, 2));
+    ret(i, 1) = timesMinusI(rhs(i, 3));
+    ret(i, 2) = timesI(rhs(i, 0));
+    ret(i, 3) = timesMinusI(rhs(i, 1));
+  }
+};
+
+template<class vtype>
+inline void multGammaZGamma5(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = timesMinusI(rhs(2));
+  ret(1) = timesI(rhs(3));
+  ret(2) = timesMinusI(rhs(0));
+  ret(3) = timesI(rhs(1));
+};
+
+template<class vtype>
+inline void lmultGammaZGamma5(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = timesMinusI(rhs(2, i));
+    ret(1, i) = timesI(rhs(3, i));
+    ret(2, i) = timesMinusI(rhs(0, i));
+    ret(3, i) = timesI(rhs(1, i));
+  }
+};
+
+template<class vtype>
+inline void rmultGammaZGamma5(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = timesMinusI(rhs(i, 2));
+    ret(i, 1) = timesI(rhs(i, 3));
+    ret(i, 2) = timesMinusI(rhs(i, 0));
+    ret(i, 3) = timesI(rhs(i, 1));
+  }
+};
+
+template<class vtype>
+inline void multMinusIdentity(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = -rhs(0);
+  ret(1) = -rhs(1);
+  ret(2) = -rhs(2);
+  ret(3) = -rhs(3);
+};
+
+template<class vtype>
+inline void lmultMinusIdentity(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = -rhs(0, i);
+    ret(1, i) = -rhs(1, i);
+    ret(2, i) = -rhs(2, i);
+    ret(3, i) = -rhs(3, i);
+  }
+};
+
+template<class vtype>
+inline void rmultMinusIdentity(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = -rhs(i, 0);
+    ret(i, 1) = -rhs(i, 1);
+    ret(i, 2) = -rhs(i, 2);
+    ret(i, 3) = -rhs(i, 3);
+  }
+};
+
+template<class vtype>
+inline void multIdentity(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = rhs(0);
+  ret(1) = rhs(1);
+  ret(2) = rhs(2);
+  ret(3) = rhs(3);
+};
+
+template<class vtype>
+inline void lmultIdentity(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = rhs(0, i);
+    ret(1, i) = rhs(1, i);
+    ret(2, i) = rhs(2, i);
+    ret(3, i) = rhs(3, i);
+  }
+};
+
+template<class vtype>
+inline void rmultIdentity(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = rhs(i, 0);
+    ret(i, 1) = rhs(i, 1);
+    ret(i, 2) = rhs(i, 2);
+    ret(i, 3) = rhs(i, 3);
+  }
+};
+
+template<class vtype>
+inline void multMinusSigmaXT(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = timesMinusI(rhs(1));
+  ret(1) = timesMinusI(rhs(0));
+  ret(2) = timesI(rhs(3));
+  ret(3) = timesI(rhs(2));
+};
+
+template<class vtype>
+inline void lmultMinusSigmaXT(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = timesMinusI(rhs(1, i));
+    ret(1, i) = timesMinusI(rhs(0, i));
+    ret(2, i) = timesI(rhs(3, i));
+    ret(3, i) = timesI(rhs(2, i));
+  }
+};
+
+template<class vtype>
+inline void rmultMinusSigmaXT(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = timesMinusI(rhs(i, 1));
+    ret(i, 1) = timesMinusI(rhs(i, 0));
+    ret(i, 2) = timesI(rhs(i, 3));
+    ret(i, 3) = timesI(rhs(i, 2));
+  }
+};
+
+template<class vtype>
+inline void multSigmaXT(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = timesI(rhs(1));
+  ret(1) = timesI(rhs(0));
+  ret(2) = timesMinusI(rhs(3));
+  ret(3) = timesMinusI(rhs(2));
+};
+
+template<class vtype>
+inline void lmultSigmaXT(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = timesI(rhs(1, i));
+    ret(1, i) = timesI(rhs(0, i));
+    ret(2, i) = timesMinusI(rhs(3, i));
+    ret(3, i) = timesMinusI(rhs(2, i));
+  }
+};
+
+template<class vtype>
+inline void rmultSigmaXT(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = timesI(rhs(i, 1));
+    ret(i, 1) = timesI(rhs(i, 0));
+    ret(i, 2) = timesMinusI(rhs(i, 3));
+    ret(i, 3) = timesMinusI(rhs(i, 2));
+  }
+};
+
+template<class vtype>
+inline void multMinusSigmaXY(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = timesI(rhs(0));
+  ret(1) = timesMinusI(rhs(1));
+  ret(2) = timesI(rhs(2));
+  ret(3) = timesMinusI(rhs(3));
+};
+
+template<class vtype>
+inline void lmultMinusSigmaXY(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = timesI(rhs(0, i));
+    ret(1, i) = timesMinusI(rhs(1, i));
+    ret(2, i) = timesI(rhs(2, i));
+    ret(3, i) = timesMinusI(rhs(3, i));
+  }
+};
+
+template<class vtype>
+inline void rmultMinusSigmaXY(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = timesI(rhs(i, 0));
+    ret(i, 1) = timesMinusI(rhs(i, 1));
+    ret(i, 2) = timesI(rhs(i, 2));
+    ret(i, 3) = timesMinusI(rhs(i, 3));
+  }
+};
+
+template<class vtype>
+inline void multSigmaXY(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = timesMinusI(rhs(0));
+  ret(1) = timesI(rhs(1));
+  ret(2) = timesMinusI(rhs(2));
+  ret(3) = timesI(rhs(3));
+};
+
+template<class vtype>
+inline void lmultSigmaXY(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = timesMinusI(rhs(0, i));
+    ret(1, i) = timesI(rhs(1, i));
+    ret(2, i) = timesMinusI(rhs(2, i));
+    ret(3, i) = timesI(rhs(3, i));
+  }
+};
+
+template<class vtype>
+inline void rmultSigmaXY(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = timesMinusI(rhs(i, 0));
+    ret(i, 1) = timesI(rhs(i, 1));
+    ret(i, 2) = timesMinusI(rhs(i, 2));
+    ret(i, 3) = timesI(rhs(i, 3));
+  }
+};
+
+template<class vtype>
+inline void multMinusSigmaXZ(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = rhs(1);
+  ret(1) = -rhs(0);
+  ret(2) = rhs(3);
+  ret(3) = -rhs(2);
+};
+
+template<class vtype>
+inline void lmultMinusSigmaXZ(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = rhs(1, i);
+    ret(1, i) = -rhs(0, i);
+    ret(2, i) = rhs(3, i);
+    ret(3, i) = -rhs(2, i);
+  }
+};
+
+template<class vtype>
+inline void rmultMinusSigmaXZ(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = -rhs(i, 1);
+    ret(i, 1) = rhs(i, 0);
+    ret(i, 2) = -rhs(i, 3);
+    ret(i, 3) = rhs(i, 2);
+  }
+};
+
+template<class vtype>
+inline void multSigmaXZ(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = -rhs(1);
+  ret(1) = rhs(0);
+  ret(2) = -rhs(3);
+  ret(3) = rhs(2);
+};
+
+template<class vtype>
+inline void lmultSigmaXZ(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = -rhs(1, i);
+    ret(1, i) = rhs(0, i);
+    ret(2, i) = -rhs(3, i);
+    ret(3, i) = rhs(2, i);
+  }
+};
+
+template<class vtype>
+inline void rmultSigmaXZ(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = rhs(i, 1);
+    ret(i, 1) = -rhs(i, 0);
+    ret(i, 2) = rhs(i, 3);
+    ret(i, 3) = -rhs(i, 2);
+  }
+};
+
+template<class vtype>
+inline void multMinusSigmaYT(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = rhs(1);
+  ret(1) = -rhs(0);
+  ret(2) = -rhs(3);
+  ret(3) = rhs(2);
+};
+
+template<class vtype>
+inline void lmultMinusSigmaYT(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = rhs(1, i);
+    ret(1, i) = -rhs(0, i);
+    ret(2, i) = -rhs(3, i);
+    ret(3, i) = rhs(2, i);
+  }
+};
+
+template<class vtype>
+inline void rmultMinusSigmaYT(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = -rhs(i, 1);
+    ret(i, 1) = rhs(i, 0);
+    ret(i, 2) = rhs(i, 3);
+    ret(i, 3) = -rhs(i, 2);
+  }
+};
+
+template<class vtype>
+inline void multSigmaYT(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = -rhs(1);
+  ret(1) = rhs(0);
+  ret(2) = rhs(3);
+  ret(3) = -rhs(2);
+};
+
+template<class vtype>
+inline void lmultSigmaYT(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = -rhs(1, i);
+    ret(1, i) = rhs(0, i);
+    ret(2, i) = rhs(3, i);
+    ret(3, i) = -rhs(2, i);
+  }
+};
+
+template<class vtype>
+inline void rmultSigmaYT(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = rhs(i, 1);
+    ret(i, 1) = -rhs(i, 0);
+    ret(i, 2) = -rhs(i, 3);
+    ret(i, 3) = rhs(i, 2);
+  }
+};
+
+template<class vtype>
+inline void multMinusSigmaYZ(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = timesI(rhs(1));
+  ret(1) = timesI(rhs(0));
+  ret(2) = timesI(rhs(3));
+  ret(3) = timesI(rhs(2));
+};
+
+template<class vtype>
+inline void lmultMinusSigmaYZ(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = timesI(rhs(1, i));
+    ret(1, i) = timesI(rhs(0, i));
+    ret(2, i) = timesI(rhs(3, i));
+    ret(3, i) = timesI(rhs(2, i));
+  }
+};
+
+template<class vtype>
+inline void rmultMinusSigmaYZ(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = timesI(rhs(i, 1));
+    ret(i, 1) = timesI(rhs(i, 0));
+    ret(i, 2) = timesI(rhs(i, 3));
+    ret(i, 3) = timesI(rhs(i, 2));
+  }
+};
+
+template<class vtype>
+inline void multSigmaYZ(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = timesMinusI(rhs(1));
+  ret(1) = timesMinusI(rhs(0));
+  ret(2) = timesMinusI(rhs(3));
+  ret(3) = timesMinusI(rhs(2));
+};
+
+template<class vtype>
+inline void lmultSigmaYZ(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = timesMinusI(rhs(1, i));
+    ret(1, i) = timesMinusI(rhs(0, i));
+    ret(2, i) = timesMinusI(rhs(3, i));
+    ret(3, i) = timesMinusI(rhs(2, i));
+  }
+};
+
+template<class vtype>
+inline void rmultSigmaYZ(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = timesMinusI(rhs(i, 1));
+    ret(i, 1) = timesMinusI(rhs(i, 0));
+    ret(i, 2) = timesMinusI(rhs(i, 3));
+    ret(i, 3) = timesMinusI(rhs(i, 2));
+  }
+};
+
+template<class vtype>
+inline void multMinusSigmaZT(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = timesMinusI(rhs(0));
+  ret(1) = timesI(rhs(1));
+  ret(2) = timesI(rhs(2));
+  ret(3) = timesMinusI(rhs(3));
+};
+
+template<class vtype>
+inline void lmultMinusSigmaZT(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = timesMinusI(rhs(0, i));
+    ret(1, i) = timesI(rhs(1, i));
+    ret(2, i) = timesI(rhs(2, i));
+    ret(3, i) = timesMinusI(rhs(3, i));
+  }
+};
+
+template<class vtype>
+inline void rmultMinusSigmaZT(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = timesMinusI(rhs(i, 0));
+    ret(i, 1) = timesI(rhs(i, 1));
+    ret(i, 2) = timesI(rhs(i, 2));
+    ret(i, 3) = timesMinusI(rhs(i, 3));
+  }
+};
+
+template<class vtype>
+inline void multSigmaZT(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> &rhs)
+{
+  ret(0) = timesI(rhs(0));
+  ret(1) = timesMinusI(rhs(1));
+  ret(2) = timesMinusI(rhs(2));
+  ret(3) = timesI(rhs(3));
+};
+
+template<class vtype>
+inline void lmultSigmaZT(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(0, i) = timesI(rhs(0, i));
+    ret(1, i) = timesMinusI(rhs(1, i));
+    ret(2, i) = timesMinusI(rhs(2, i));
+    ret(3, i) = timesI(rhs(3, i));
+  }
+};
+
+template<class vtype>
+inline void rmultSigmaZT(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)
+{
+  for(unsigned int i = 0; i < Ns; ++i)
+  {
+    ret(i, 0) = timesI(rhs(i, 0));
+    ret(i, 1) = timesMinusI(rhs(i, 1));
+    ret(i, 2) = timesMinusI(rhs(i, 2));
+    ret(i, 3) = timesI(rhs(i, 3));
+  }
+};
+
+
+// FIXME
+//
+// Optimisation; switch over to a "multGammaX(ret._internal,arg._internal)" style early and
+// note that doing so from the lattice operator will avoid copy back and case switch overhead, as
+// was done for the tensor math operator to remove operator * notation early
+//
+
+//left multiply
+template<class vtype> 
+inline auto operator*(const Gamma &G, const iScalar<vtype> &arg)
+->typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype>>::type 
+{
+  iScalar<vtype> ret;
+  ret._internal=G*arg._internal;
+  return ret;
+}
+
+template<class vtype,int N>
+inline auto operator*(const Gamma &G, const iVector<vtype, N> &arg)
+->typename std::enable_if<matchGridTensorIndex<iVector<vtype,N>,SpinorIndex>::notvalue,iVector<vtype,N>>::type 
+{
+  iVector<vtype,N> ret;
+  for(int i=0;i<N;i++){
+    ret._internal[i]=G*arg._internal[i];
+  }
+  return ret;
+}
+
+template<class vtype, int N>
+inline auto operator*(const Gamma &G, const iMatrix<vtype, N> &arg) 
+->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,iMatrix<vtype,N>>::type 
+{
+  iMatrix<vtype,N> ret;
+  for(int i=0;i<N;i++){
+  for(int j=0;j<N;j++){
+    ret._internal[i][j]=G*arg._internal[i][j];
+  }}
+  return ret;
+}
+
+//right multiply
+template<class vtype>
+inline auto operator*(const iScalar<vtype> &arg, const Gamma &G)
+->typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype>>::type 
+{
+  iScalar<vtype> ret;
+  ret._internal=arg._internal*G;
+  return ret;
+}
+
+template<class vtype, int N>
+inline auto operator * (const iMatrix<vtype, N> &arg, const Gamma &G) 
+->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,iMatrix<vtype,N>>::type 
+{
+  iMatrix<vtype,N> ret;
+  for(int i=0;i<N;i++){
+  for(int j=0;j<N;j++){
+    ret._internal[i][j]=arg._internal[i][j]*G;
+  }}
+  return ret;
+}
+
+template<class vtype> 
+inline auto operator*(const Gamma &G, const iVector<vtype, Ns> &arg)
+->typename std::enable_if<matchGridTensorIndex<iVector<vtype, Ns>, SpinorIndex>::value, iVector<vtype, Ns>>::type
+{
+  iVector<vtype, Ns> ret;
+
+  switch (G.g) 
+  {
+  case Gamma::Algebra::MinusGamma5:
+    multMinusGamma5(ret, arg); break;
+  case Gamma::Algebra::Gamma5:
+    multGamma5(ret, arg); break;
+  case Gamma::Algebra::MinusGammaT:
+    multMinusGammaT(ret, arg); break;
+  case Gamma::Algebra::GammaT:
+    multGammaT(ret, arg); break;
+  case Gamma::Algebra::MinusGammaTGamma5:
+    multMinusGammaTGamma5(ret, arg); break;
+  case Gamma::Algebra::GammaTGamma5:
+    multGammaTGamma5(ret, arg); break;
+  case Gamma::Algebra::MinusGammaX:
+    multMinusGammaX(ret, arg); break;
+  case Gamma::Algebra::GammaX:
+    multGammaX(ret, arg); break;
+  case Gamma::Algebra::MinusGammaXGamma5:
+    multMinusGammaXGamma5(ret, arg); break;
+  case Gamma::Algebra::GammaXGamma5:
+    multGammaXGamma5(ret, arg); break;
+  case Gamma::Algebra::MinusGammaY:
+    multMinusGammaY(ret, arg); break;
+  case Gamma::Algebra::GammaY:
+    multGammaY(ret, arg); break;
+  case Gamma::Algebra::MinusGammaYGamma5:
+    multMinusGammaYGamma5(ret, arg); break;
+  case Gamma::Algebra::GammaYGamma5:
+    multGammaYGamma5(ret, arg); break;
+  case Gamma::Algebra::MinusGammaZ:
+    multMinusGammaZ(ret, arg); break;
+  case Gamma::Algebra::GammaZ:
+    multGammaZ(ret, arg); break;
+  case Gamma::Algebra::MinusGammaZGamma5:
+    multMinusGammaZGamma5(ret, arg); break;
+  case Gamma::Algebra::GammaZGamma5:
+    multGammaZGamma5(ret, arg); break;
+  case Gamma::Algebra::MinusIdentity:
+    multMinusIdentity(ret, arg); break;
+  case Gamma::Algebra::Identity:
+    multIdentity(ret, arg); break;
+  case Gamma::Algebra::MinusSigmaXT:
+    multMinusSigmaXT(ret, arg); break;
+  case Gamma::Algebra::SigmaXT:
+    multSigmaXT(ret, arg); break;
+  case Gamma::Algebra::MinusSigmaXY:
+    multMinusSigmaXY(ret, arg); break;
+  case Gamma::Algebra::SigmaXY:
+    multSigmaXY(ret, arg); break;
+  case Gamma::Algebra::MinusSigmaXZ:
+    multMinusSigmaXZ(ret, arg); break;
+  case Gamma::Algebra::SigmaXZ:
+    multSigmaXZ(ret, arg); break;
+  case Gamma::Algebra::MinusSigmaYT:
+    multMinusSigmaYT(ret, arg); break;
+  case Gamma::Algebra::SigmaYT:
+    multSigmaYT(ret, arg); break;
+  case Gamma::Algebra::MinusSigmaYZ:
+    multMinusSigmaYZ(ret, arg); break;
+  case Gamma::Algebra::SigmaYZ:
+    multSigmaYZ(ret, arg); break;
+  case Gamma::Algebra::MinusSigmaZT:
+    multMinusSigmaZT(ret, arg); break;
+  case Gamma::Algebra::SigmaZT:
+    multSigmaZT(ret, arg); break;
+  default: assert(0);
+  }
   
-    // Make gamma products (Chroma convention)
-    SpinMatrix makeGammaProd(const unsigned int i);
-    
-    /* Gx
-     *  0 0  0  i    
-     *  0 0  i  0    
-     *  0 -i 0  0
-     * -i 0  0  0
-     */
-  // right multiplication makes sense for matrix args, not for vector since there is 
-  // no concept of row versus columnar indices
-    template<class vtype> inline void rmultMinusGammaX(iMatrix<vtype,Ns> &ret,const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(i,0) = timesI(rhs(i,3));
-	ret(i,1) = timesI(rhs(i,2));
-	ret(i,2) = timesMinusI(rhs(i,1));
-	ret(i,3) = timesMinusI(rhs(i,0));
-      }
-    };
-    template<class vtype> inline void rmultGammaX(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(i,0) = timesMinusI(rhs(i,3));
-	ret(i,1) = timesMinusI(rhs(i,2));
-	ret(i,2) = timesI(rhs(i,1));
-	ret(i,3) = timesI(rhs(i,0));
-      }
-    };
-    template<class vtype> inline void multGammaX(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(0,i) = timesI(rhs(3,i));
-	ret(1,i) = timesI(rhs(2,i));
-	ret(2,i) = timesMinusI(rhs(1,i));
-	ret(3,i) = timesMinusI(rhs(0,i));
-      }
-    };
-    template<class vtype> inline void multMinusGammaX(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(0,i) = timesMinusI(rhs(3,i));
-	ret(1,i) = timesMinusI(rhs(2,i));
-	ret(2,i) = timesI(rhs(1,i));
-	ret(3,i) = timesI(rhs(0,i));
-      }
-    };
+  return ret;
+}
 
-    template<class vtype> inline void multGammaX(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
-      ret._internal[0] = timesI(rhs._internal[3]);
-      ret._internal[1] = timesI(rhs._internal[2]);
-      ret._internal[2] = timesMinusI(rhs._internal[1]);
-      ret._internal[3] = timesMinusI(rhs._internal[0]);
-    };
-    template<class vtype> inline void multMinusGammaX(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
-	ret(0) = timesMinusI(rhs(3));
-	ret(1) = timesMinusI(rhs(2));
-	ret(2) = timesI(rhs(1));
-	ret(3) = timesI(rhs(0));
-    };
+template<class vtype> 
+inline auto operator*(const Gamma &G, const iMatrix<vtype, Ns> &arg)
+->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype, Ns>, SpinorIndex>::value, iMatrix<vtype, Ns>>::type
+{
+  iMatrix<vtype, Ns> ret;
 
+  switch (G.g) 
+  {
+  case Gamma::Algebra::MinusGamma5:
+    lmultMinusGamma5(ret, arg); break;
+  case Gamma::Algebra::Gamma5:
+    lmultGamma5(ret, arg); break;
+  case Gamma::Algebra::MinusGammaT:
+    lmultMinusGammaT(ret, arg); break;
+  case Gamma::Algebra::GammaT:
+    lmultGammaT(ret, arg); break;
+  case Gamma::Algebra::MinusGammaTGamma5:
+    lmultMinusGammaTGamma5(ret, arg); break;
+  case Gamma::Algebra::GammaTGamma5:
+    lmultGammaTGamma5(ret, arg); break;
+  case Gamma::Algebra::MinusGammaX:
+    lmultMinusGammaX(ret, arg); break;
+  case Gamma::Algebra::GammaX:
+    lmultGammaX(ret, arg); break;
+  case Gamma::Algebra::MinusGammaXGamma5:
+    lmultMinusGammaXGamma5(ret, arg); break;
+  case Gamma::Algebra::GammaXGamma5:
+    lmultGammaXGamma5(ret, arg); break;
+  case Gamma::Algebra::MinusGammaY:
+    lmultMinusGammaY(ret, arg); break;
+  case Gamma::Algebra::GammaY:
+    lmultGammaY(ret, arg); break;
+  case Gamma::Algebra::MinusGammaYGamma5:
+    lmultMinusGammaYGamma5(ret, arg); break;
+  case Gamma::Algebra::GammaYGamma5:
+    lmultGammaYGamma5(ret, arg); break;
+  case Gamma::Algebra::MinusGammaZ:
+    lmultMinusGammaZ(ret, arg); break;
+  case Gamma::Algebra::GammaZ:
+    lmultGammaZ(ret, arg); break;
+  case Gamma::Algebra::MinusGammaZGamma5:
+    lmultMinusGammaZGamma5(ret, arg); break;
+  case Gamma::Algebra::GammaZGamma5:
+    lmultGammaZGamma5(ret, arg); break;
+  case Gamma::Algebra::MinusIdentity:
+    lmultMinusIdentity(ret, arg); break;
+  case Gamma::Algebra::Identity:
+    lmultIdentity(ret, arg); break;
+  case Gamma::Algebra::MinusSigmaXT:
+    lmultMinusSigmaXT(ret, arg); break;
+  case Gamma::Algebra::SigmaXT:
+    lmultSigmaXT(ret, arg); break;
+  case Gamma::Algebra::MinusSigmaXY:
+    lmultMinusSigmaXY(ret, arg); break;
+  case Gamma::Algebra::SigmaXY:
+    lmultSigmaXY(ret, arg); break;
+  case Gamma::Algebra::MinusSigmaXZ:
+    lmultMinusSigmaXZ(ret, arg); break;
+  case Gamma::Algebra::SigmaXZ:
+    lmultSigmaXZ(ret, arg); break;
+  case Gamma::Algebra::MinusSigmaYT:
+    lmultMinusSigmaYT(ret, arg); break;
+  case Gamma::Algebra::SigmaYT:
+    lmultSigmaYT(ret, arg); break;
+  case Gamma::Algebra::MinusSigmaYZ:
+    lmultMinusSigmaYZ(ret, arg); break;
+  case Gamma::Algebra::SigmaYZ:
+    lmultSigmaYZ(ret, arg); break;
+  case Gamma::Algebra::MinusSigmaZT:
+    lmultMinusSigmaZT(ret, arg); break;
+  case Gamma::Algebra::SigmaZT:
+    lmultSigmaZT(ret, arg); break;
+  default: assert(0);
+  }
+  
+  return ret;
+}
 
-    /*Gy
-     *  0 0  0  -1  [0] -+ [3]
-     *  0 0  1  0   [1] +- [2]
-     *  0 1  0  0
-     * -1 0  0  0
-     */
-    template<class vtype> inline void rmultGammaY(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(i,0) = -rhs(i,3);
-	ret(i,1) =  rhs(i,2);
-	ret(i,2) =  rhs(i,1);
-	ret(i,3) = -rhs(i,0);
-      }
-    };
-    template<class vtype> inline void rmultMinusGammaY(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(i,0) =  rhs(i,3);
-	ret(i,1) = -rhs(i,2);
-	ret(i,2) = -rhs(i,1);
-	ret(i,3) =  rhs(i,0);
-      }
-    };
-    template<class vtype> inline void multGammaY(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(0,i) = -rhs(3,i);
-	ret(1,i) =  rhs(2,i);
-	ret(2,i) =  rhs(1,i);
-	ret(3,i) = -rhs(0,i);
-      }
-    };
-    template<class vtype> inline void multMinusGammaY(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(0,i) =  rhs(3,i);
-	ret(1,i) = -rhs(2,i);
-	ret(2,i) = -rhs(1,i);
-	ret(3,i) =  rhs(0,i);
-      }
-    };
-    template<class vtype> inline void multGammaY(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
-      ret(0) = -rhs(3);
-      ret(1) =  rhs(2);
-      ret(2) =  rhs(1);
-      ret(3) = -rhs(0);
-    };
-    template<class vtype> inline void multMinusGammaY(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
-      ret(0) =  rhs(3);
-      ret(1) = -rhs(2);
-      ret(2) = -rhs(1);
-      ret(3) =  rhs(0);
-    };
-    /*Gz
-     *  0 0  i  0   [0]+-i[2]
-     *  0 0  0 -i   [1]-+i[3]
-     * -i 0  0  0
-     *  0 i  0  0
-     */
-    template<class vtype> inline void rmultGammaZ(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(i,0) = timesMinusI(rhs(i,2));
-	ret(i,1) =      timesI(rhs(i,3));
-	ret(i,2) =      timesI(rhs(i,0));
-	ret(i,3) = timesMinusI(rhs(i,1));
-      }
-    };
-    template<class vtype> inline void rmultMinusGammaZ(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(i,0) =      timesI(rhs(i,2));
-	ret(i,1) = timesMinusI(rhs(i,3));
-	ret(i,2) = timesMinusI(rhs(i,0));
-	ret(i,3) =      timesI(rhs(i,1));
-      }
-    };
-    template<class vtype> inline void multGammaZ(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(0,i) = timesI(rhs(2,i));
-	ret(1,i) =timesMinusI(rhs(3,i));
-	ret(2,i) =timesMinusI(rhs(0,i));
-	ret(3,i) = timesI(rhs(1,i));
-      }
-    };
-    template<class vtype> inline void multMinusGammaZ(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(0,i) = timesMinusI(rhs(2,i));
-	ret(1,i) = timesI(rhs(3,i));
-	ret(2,i) = timesI(rhs(0,i));
-	ret(3,i) = timesMinusI(rhs(1,i));
-      }
-    };
-    template<class vtype> inline void multGammaZ(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
-      ret(0) = timesI(rhs(2));
-      ret(1) =timesMinusI(rhs(3));
-      ret(2) =timesMinusI(rhs(0));
-      ret(3) = timesI(rhs(1));
-    };
-    template<class vtype> inline void multMinusGammaZ(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
-      ret(0) = timesMinusI(rhs(2));
-      ret(1) = timesI(rhs(3));
-      ret(2) = timesI(rhs(0));
-      ret(3) = timesMinusI(rhs(1));
-    };
-    /*Gt
-     *  0 0  1  0 [0]+-[2]
-     *  0 0  0  1 [1]+-[3]
-     *  1 0  0  0
-     *  0 1  0  0
-     */
-    template<class vtype> inline void rmultGammaT(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(i,0) = rhs(i,2);
-	ret(i,1) = rhs(i,3);
-	ret(i,2) = rhs(i,0);
-	ret(i,3) = rhs(i,1);
-      }
-    };
-    template<class vtype> inline void rmultMinusGammaT(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(i,0) =- rhs(i,2);
-	ret(i,1) =- rhs(i,3);
-	ret(i,2) =- rhs(i,0);
-	ret(i,3) =- rhs(i,1);
-      }
-    };
-    template<class vtype> inline void multGammaT(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(0,i) = rhs(2,i);
-	ret(1,i) = rhs(3,i);
-	ret(2,i) = rhs(0,i);
-	ret(3,i) = rhs(1,i);
-      }
-    };
-    template<class vtype> inline void multMinusGammaT(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(0,i) =-rhs(2,i);
-	ret(1,i) =-rhs(3,i);
-	ret(2,i) =-rhs(0,i);
-	ret(3,i) =-rhs(1,i);
-      }
-    };
-    template<class vtype> inline void multGammaT(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
-      ret(0) = rhs(2);
-      ret(1) = rhs(3);
-      ret(2) = rhs(0);
-      ret(3) = rhs(1);
-    };
-    template<class vtype> inline void multMinusGammaT(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
-      ret(0) =-rhs(2);
-      ret(1) =-rhs(3);
-      ret(2) =-rhs(0);
-      ret(3) =-rhs(1);
-    };
-    /*G5
-     *  1 0  0  0 [0]+-[2]
-     *  0 1  0  0 [1]+-[3]
-     *  0 0 -1  0
-     *  0 0  0 -1
-     */
-    template<class vtype> inline void rmultGamma5(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(i,0) = rhs(i,0);
-	ret(i,1) = rhs(i,1);
-	ret(i,2) =-rhs(i,2);
-	ret(i,3) =-rhs(i,3);
-      }
-    };
-    template<class vtype> inline void rmultMinusGamma5(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(i,0) =-rhs(i,0);
-	ret(i,1) =-rhs(i,1);
-	ret(i,2) = rhs(i,2);
-	ret(i,3) = rhs(i,3);
-      }
-    };
+template<class vtype> 
+inline auto operator*(const iMatrix<vtype, Ns> &arg, const Gamma &G)
+->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype, Ns>, SpinorIndex>::value, iMatrix<vtype, Ns>>::type
+{
+  iMatrix<vtype, Ns> ret;
 
-    template<class vtype> inline void multGamma5(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(0,i) = rhs(0,i);
-	ret(1,i) = rhs(1,i);
-	ret(2,i) =-rhs(2,i);
-	ret(3,i) =-rhs(3,i);
-      }
-    };
-    template<class vtype> inline void multMinusGamma5(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
-      for(int i=0;i<Ns;i++){
-	ret(0,i) =-rhs(0,i);
-	ret(1,i) =-rhs(1,i);
-	ret(2,i) = rhs(2,i);
-	ret(3,i) = rhs(3,i);
-      }
-    };
+  switch (G.g) 
+  {
+  case Gamma::Algebra::MinusGamma5:
+    rmultMinusGamma5(ret, arg); break;
+  case Gamma::Algebra::Gamma5:
+    rmultGamma5(ret, arg); break;
+  case Gamma::Algebra::MinusGammaT:
+    rmultMinusGammaT(ret, arg); break;
+  case Gamma::Algebra::GammaT:
+    rmultGammaT(ret, arg); break;
+  case Gamma::Algebra::MinusGammaTGamma5:
+    rmultMinusGammaTGamma5(ret, arg); break;
+  case Gamma::Algebra::GammaTGamma5:
+    rmultGammaTGamma5(ret, arg); break;
+  case Gamma::Algebra::MinusGammaX:
+    rmultMinusGammaX(ret, arg); break;
+  case Gamma::Algebra::GammaX:
+    rmultGammaX(ret, arg); break;
+  case Gamma::Algebra::MinusGammaXGamma5:
+    rmultMinusGammaXGamma5(ret, arg); break;
+  case Gamma::Algebra::GammaXGamma5:
+    rmultGammaXGamma5(ret, arg); break;
+  case Gamma::Algebra::MinusGammaY:
+    rmultMinusGammaY(ret, arg); break;
+  case Gamma::Algebra::GammaY:
+    rmultGammaY(ret, arg); break;
+  case Gamma::Algebra::MinusGammaYGamma5:
+    rmultMinusGammaYGamma5(ret, arg); break;
+  case Gamma::Algebra::GammaYGamma5:
+    rmultGammaYGamma5(ret, arg); break;
+  case Gamma::Algebra::MinusGammaZ:
+    rmultMinusGammaZ(ret, arg); break;
+  case Gamma::Algebra::GammaZ:
+    rmultGammaZ(ret, arg); break;
+  case Gamma::Algebra::MinusGammaZGamma5:
+    rmultMinusGammaZGamma5(ret, arg); break;
+  case Gamma::Algebra::GammaZGamma5:
+    rmultGammaZGamma5(ret, arg); break;
+  case Gamma::Algebra::MinusIdentity:
+    rmultMinusIdentity(ret, arg); break;
+  case Gamma::Algebra::Identity:
+    rmultIdentity(ret, arg); break;
+  case Gamma::Algebra::MinusSigmaXT:
+    rmultMinusSigmaXT(ret, arg); break;
+  case Gamma::Algebra::SigmaXT:
+    rmultSigmaXT(ret, arg); break;
+  case Gamma::Algebra::MinusSigmaXY:
+    rmultMinusSigmaXY(ret, arg); break;
+  case Gamma::Algebra::SigmaXY:
+    rmultSigmaXY(ret, arg); break;
+  case Gamma::Algebra::MinusSigmaXZ:
+    rmultMinusSigmaXZ(ret, arg); break;
+  case Gamma::Algebra::SigmaXZ:
+    rmultSigmaXZ(ret, arg); break;
+  case Gamma::Algebra::MinusSigmaYT:
+    rmultMinusSigmaYT(ret, arg); break;
+  case Gamma::Algebra::SigmaYT:
+    rmultSigmaYT(ret, arg); break;
+  case Gamma::Algebra::MinusSigmaYZ:
+    rmultMinusSigmaYZ(ret, arg); break;
+  case Gamma::Algebra::SigmaYZ:
+    rmultSigmaYZ(ret, arg); break;
+  case Gamma::Algebra::MinusSigmaZT:
+    rmultMinusSigmaZT(ret, arg); break;
+  case Gamma::Algebra::SigmaZT:
+    rmultSigmaZT(ret, arg); break;
+  default: assert(0);
+  }
+  
+  return ret;
+}
 
-    template<class vtype> inline void multGamma5(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
-      ret(0) = rhs(0);
-      ret(1) = rhs(1);
-      ret(2) =-rhs(2);
-      ret(3) =-rhs(3);
-    };
-    template<class vtype> inline void multMinusGamma5(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
-      ret(0) =-rhs(0);
-      ret(1) =-rhs(1);
-      ret(2) = rhs(2);
-      ret(3) = rhs(3);
-    };
+}
 
+inline QCD::Gamma adj(const QCD::Gamma &g)
+{
+  return QCD::Gamma(QCD::Gamma::adj[g.g]);
+}
 
+}
 
-#ifdef GRID_WARN_SUBOPTIMAL
-#warning "Optimisation alert switch over to multGammaX early "
-#endif     
-
-    ///////////////////////////////////////////////////////////////////////////////////////////////////
-    // Operator * : first case this is not a spin index, so recurse
-    ///////////////////////////////////////////////////////////////////////////////////////////////////
-    
-    // FIXME
-    //
-    // Optimisation; switch over to a "multGammaX(ret._internal,arg._internal)" style early and
-    // note that doing so from the lattice operator will avoid copy back and case switch overhead, as
-    // was done for the tensor math operator to remove operator * notation early
-    //
-
-    //left multiply
-    template<class vtype> inline auto operator * ( const Gamma &G,const iScalar<vtype> &arg) ->
-      typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype> >::type 
-
-	{
-	  iScalar<vtype> ret;
-	  ret._internal=G*arg._internal;
-	  return ret;
-	}
-    template<class vtype,int N> inline auto operator * ( const Gamma &G,const iVector<vtype,N> &arg) ->
-      typename std::enable_if<matchGridTensorIndex<iVector<vtype,N>,SpinorIndex>::notvalue,iVector<vtype,N> >::type 
-	{
-	  iVector<vtype,N> ret;
-	  for(int i=0;i<N;i++){
-	    ret._internal[i]=G*arg._internal[i];
-	  }
-	  return ret;
-	}
-    template<class vtype,int N> inline auto operator * ( const Gamma &G,const iMatrix<vtype,N> &arg) ->
-      typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,iMatrix<vtype,N> >::type 
-	{
-	  iMatrix<vtype,N> ret;
-	  for(int i=0;i<N;i++){
-	  for(int j=0;j<N;j++){
-	    ret._internal[i][j]=G*arg._internal[i][j];
-	  }}
-	  return ret;
-	}
-
-
-    //right multiply
-    template<class vtype> inline auto operator * (const iScalar<vtype> &arg, const Gamma &G) ->
-      typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype> >::type 
-
-	{
-	  iScalar<vtype> ret;
-	  ret._internal=arg._internal*G;
-	  return ret;
-	}
-    template<class vtype,int N> inline auto operator * (const iVector<vtype,N> &arg, const Gamma &G) ->
-      typename std::enable_if<matchGridTensorIndex<iVector<vtype,N>,SpinorIndex>::notvalue,iVector<vtype,N> >::type 
-	{
-	  iVector<vtype,N> ret;
-	  for(int i=0;i<N;i++){
-	    ret._internal=arg._internal[i]*G;
-	  }
-	  return ret;
-	}
-    template<class vtype,int N> inline auto operator * (const iMatrix<vtype,N> &arg, const Gamma &G) ->
-      typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,iMatrix<vtype,N> >::type 
-	{
-	  iMatrix<vtype,N> ret;
-	  for(int i=0;i<N;i++){
-	  for(int j=0;j<N;j++){
-	    ret._internal[i][j]=arg._internal[i][j]*G;
-	  }}
-	  return ret;
-	}
-
-    ////////////////////////////////////////////////////////
-    // When we hit the spin index this matches and we stop
-    ////////////////////////////////////////////////////////
-    template<class vtype> inline auto operator * ( const Gamma &G,const iMatrix<vtype,Ns> &arg) ->
-      typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,Ns>,SpinorIndex>::value,iMatrix<vtype,Ns> >::type 
-      {
-	iMatrix<vtype,Ns> ret;
-	switch (G._g) {
-	case Gamma::Identity:
-	  ret = arg;
-	  break;
-	case Gamma::MinusIdentity:
-	  ret = -arg;
-	  break;
-	case Gamma::GammaX:
-	  multGammaX(ret,arg);
-	  break;
-	case Gamma::MinusGammaX:
-	  multMinusGammaX(ret,arg);
-	  break;
-	case Gamma::GammaY:
-	  multGammaY(ret,arg);
-	  break;
-	case Gamma::MinusGammaY:
-	  multMinusGammaY(ret,arg);
-	  break;
-	case Gamma::GammaZ:
-	  multGammaZ(ret,arg);
-	  break;
-	case Gamma::MinusGammaZ:
-	  multMinusGammaZ(ret,arg);
-	  break;
-	case Gamma::GammaT:
-	  multGammaT(ret,arg);
-	  break;
-	case Gamma::MinusGammaT:
-	  multMinusGammaT(ret,arg);
-	  break;
-	case Gamma::Gamma5:
-	  multGamma5(ret,arg);
-	  break;
-	case Gamma::MinusGamma5:
-	  multMinusGamma5(ret,arg);
-	  break;
-	default:
-	  assert(0);
-	  break;
-	}
-	return ret;
-      }
-    // Could have used type trait for Matrix/vector and then an enable if to share code
-    template<class vtype> inline auto operator * ( const Gamma &G,const iVector<vtype,Ns> &arg) ->
-      typename std::enable_if<matchGridTensorIndex<iVector<vtype,Ns>,SpinorIndex>::value,iVector<vtype,Ns> >::type 
-      {
-	iVector<vtype,Ns> ret;
-	switch (G._g) {
-	case Gamma::Identity:
-	  ret = arg;
-	  break;
-	case Gamma::MinusIdentity:
-	  ret = -arg;
-	  break;
-	case Gamma::GammaX:
-	  multGammaX(ret,arg);
-	  break;
-	case Gamma::MinusGammaX:
-	  multMinusGammaX(ret,arg);
-	  break;
-	case Gamma::GammaY:
-	  multGammaY(ret,arg);
-	  break;
-	case Gamma::MinusGammaY:
-	  multMinusGammaY(ret,arg);
-	  break;
-	case Gamma::GammaZ:
-	  multGammaZ(ret,arg);
-	  break;
-	case Gamma::MinusGammaZ:
-	  multMinusGammaZ(ret,arg);
-	  break;
-	case Gamma::GammaT:
-	  multGammaT(ret,arg);
-	  break;
-	case Gamma::MinusGammaT:
-	  multMinusGammaT(ret,arg);
-	  break;
-	case Gamma::Gamma5:
-	  multGamma5(ret,arg);
-	  break;
-	case Gamma::MinusGamma5:
-	  multMinusGamma5(ret,arg);
-	  break;
-	default:
-	  assert(0);
-	  break;
-	}
-	return ret;
-      }
-
-    template<class vtype> inline auto operator * (const iMatrix<vtype,Ns> &arg, const Gamma &G) ->
-      typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,Ns>,SpinorIndex>::value,iMatrix<vtype,Ns> >::type 
-      {
-	iMatrix<vtype,Ns> ret;
-	switch (G._g) {
-	case Gamma::Identity:
-	  ret = arg;
-	  break;
-	case Gamma::MinusIdentity:
-	  ret = -arg;
-	  break;
-	case Gamma::GammaX:
-	  rmultGammaX(ret,arg);
-	  break;
-	case Gamma::MinusGammaX:
-	  rmultMinusGammaX(ret,arg);
-	  break;
-	case Gamma::GammaY:
-	  rmultGammaY(ret,arg);
-	  break;
-	case Gamma::MinusGammaY:
-	  rmultMinusGammaY(ret,arg);
-	  break;
-	case Gamma::GammaZ:
-	  rmultGammaZ(ret,arg);
-	  break;
-	case Gamma::MinusGammaZ:
-	  rmultMinusGammaZ(ret,arg);
-	  break;
-	case Gamma::GammaT:
-	  rmultGammaT(ret,arg);
-	  break;
-	case Gamma::MinusGammaT:
-	  rmultMinusGammaT(ret,arg);
-	  break;
-	case Gamma::Gamma5:
-	  rmultGamma5(ret,arg);
-	  break;
-	case Gamma::MinusGamma5:
-	  rmultMinusGamma5(ret,arg);
-	  break;
-	default:
-	  assert(0);
-	  break;
-	}
-	return ret;
-      }
-
-
-    /* Output from test
-./Grid_gamma 
-Identity((1,0),(0,0),(0,0),(0,0))
-        ((0,0),(1,0),(0,0),(0,0))
-        ((0,0),(0,0),(1,0),(0,0))
-        ((0,0),(0,0),(0,0),(1,0))   OK
-
-GammaX  ((0,0),(0,0),(0,0),(0,1))
-        ((0,0),(0,0),(0,1),(0,0))
-        ((0,0),(0,-1),(0,0),(0,0))
-        ((0,-1),(0,0),(0,0),(0,0))  OK
-    * Gx
-    *  0 0  0  i    
-    *  0 0  i  0    
-    *  0 -i 0  0
-    * -i 0  0  0
-
-GammaY  ((-0,-0),(-0,-0),(-0,-0),(-1,-0))
-        ((0,0),(0,0),(1,0),(0,0))
-        ((0,0),(1,0),(0,0),(0,0))          OK
-        ((-1,-0),(-0,-0),(-0,-0),(-0,-0))
-     *Gy
-     *  0 0  0  -1  [0] -+ [3]
-     *  0 0  1  0   [1] +- [2]
-     *  0 1  0  0
-     * -1 0  0  0
-
-GammaZ  ((0,0),(0,0),(0,1),(0,0))
-        ((0,0),(0,0),(0,0),(0,-1))
-        ((0,-1),(0,0),(0,0),(0,0))
-        ((0,0),(0,1),(0,0),(0,0))   OK
-     *  0 0  i  0   [0]+-i[2]
-     *  0 0  0 -i   [1]-+i[3]
-     * -i 0  0  0
-     *  0 i  0  0
-
-GammaT  ((0,0),(0,0),(1,0),(0,0))
-        ((0,0),(0,0),(0,0),(1,0))  OK
-        ((1,0),(0,0),(0,0),(0,0))
-        ((0,0),(1,0),(0,0),(0,0))
-     *  0 0  1  0 [0]+-[2]
-     *  0 0  0  1 [1]+-[3]
-     *  1 0  0  0
-     *  0 1  0  0
-
-Gamma5  ((1,0),(0,0),(0,0),(0,0))
-        ((0,0),(1,0),(0,0),(0,0))
-        ((-0,-0),(-0,-0),(-1,-0),(-0,-0))
-        ((-0,-0),(-0,-0),(-0,-0),(-1,-0))
-     *  1 0  0  0 [0]+-[2]
-     *  0 1  0  0 [1]+-[3]  OK
-     *  0 0 -1  0
-     *  0 0  0 -1
-     */
-
-}   //namespace QCD
-} // Grid
 #endif
diff --git a/lib/qcd/spin/GammaMulTable.h b/lib/qcd/spin/GammaMulTable.h
new file mode 100644
index 00000000..a71cab5b
--- /dev/null
+++ b/lib/qcd/spin/GammaMulTable.h
@@ -0,0 +1,1121 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: lib/qcd/spin/GammaMulTable.h
+
+Copyright (C) 2015
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+const std::array<std::array<Gamma::Algebra, Gamma::nGamma>, Gamma::nGamma> 
+Gamma::mul = {{
+  {Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::SigmaXY}
+  ,
+  {Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::MinusSigmaXY}
+  ,
+  {Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaZ}
+  ,
+  {Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaZ}
+  ,
+  {Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaZGamma5}
+  ,
+  {Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaZGamma5}
+  ,
+  {Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaYGamma5}
+  ,
+  {Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaYGamma5}
+  ,
+  {Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaY}
+  ,
+  {Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaY}
+  ,
+  {Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaXGamma5}
+  ,
+  {Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaXGamma5}
+  ,
+  {Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaX}
+  ,
+  {Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaX}
+  ,
+  {Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaT}
+  ,
+  {Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaT}
+  ,
+  {Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaTGamma5}
+  ,
+  {Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaTGamma5}
+  ,
+  {Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::MinusSigmaZT}
+  ,
+  {Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::SigmaZT}
+  ,
+  {Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::SigmaXZ}
+  ,
+  {Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusSigmaXZ}
+  ,
+  {Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusGamma5}
+  ,
+  {Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::Gamma5}
+  ,
+  {Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusSigmaXT}
+  ,
+  {Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::SigmaXT}
+  ,
+  {Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::SigmaYZ}
+  ,
+  {Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::MinusSigmaYZ}
+  ,
+  {Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaYT}
+  ,
+  {Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaYT}
+  ,
+  {Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusIdentity,
+   Gamma::Algebra::Identity}
+  ,
+  {Gamma::Algebra::SigmaXY,
+   Gamma::Algebra::MinusSigmaXY,
+   Gamma::Algebra::MinusGammaZ,
+   Gamma::Algebra::GammaZ,
+   Gamma::Algebra::MinusGammaZGamma5,
+   Gamma::Algebra::GammaZGamma5,
+   Gamma::Algebra::GammaYGamma5,
+   Gamma::Algebra::MinusGammaYGamma5,
+   Gamma::Algebra::GammaY,
+   Gamma::Algebra::MinusGammaY,
+   Gamma::Algebra::MinusGammaXGamma5,
+   Gamma::Algebra::GammaXGamma5,
+   Gamma::Algebra::MinusGammaX,
+   Gamma::Algebra::GammaX,
+   Gamma::Algebra::GammaT,
+   Gamma::Algebra::MinusGammaT,
+   Gamma::Algebra::GammaTGamma5,
+   Gamma::Algebra::MinusGammaTGamma5,
+   Gamma::Algebra::MinusSigmaZT,
+   Gamma::Algebra::SigmaZT,
+   Gamma::Algebra::MinusSigmaXZ,
+   Gamma::Algebra::SigmaXZ,
+   Gamma::Algebra::MinusGamma5,
+   Gamma::Algebra::Gamma5,
+   Gamma::Algebra::SigmaXT,
+   Gamma::Algebra::MinusSigmaXT,
+   Gamma::Algebra::MinusSigmaYZ,
+   Gamma::Algebra::SigmaYZ,
+   Gamma::Algebra::SigmaYT,
+   Gamma::Algebra::MinusSigmaYT,
+   Gamma::Algebra::Identity,
+   Gamma::Algebra::MinusIdentity}
+}};
+
+const std::array<Gamma::Algebra, Gamma::nGamma> Gamma::adj = {
+  Gamma::Algebra::MinusGamma5,
+  Gamma::Algebra::Gamma5,
+  Gamma::Algebra::MinusGammaT,
+  Gamma::Algebra::GammaT,
+  Gamma::Algebra::GammaTGamma5,
+  Gamma::Algebra::MinusGammaTGamma5,
+  Gamma::Algebra::MinusGammaX,
+  Gamma::Algebra::GammaX,
+  Gamma::Algebra::GammaXGamma5,
+  Gamma::Algebra::MinusGammaXGamma5,
+  Gamma::Algebra::MinusGammaY,
+  Gamma::Algebra::GammaY,
+  Gamma::Algebra::GammaYGamma5,
+  Gamma::Algebra::MinusGammaYGamma5,
+  Gamma::Algebra::MinusGammaZ,
+  Gamma::Algebra::GammaZ,
+  Gamma::Algebra::GammaZGamma5,
+  Gamma::Algebra::MinusGammaZGamma5,
+  Gamma::Algebra::MinusIdentity,
+  Gamma::Algebra::Identity,
+  Gamma::Algebra::SigmaXT,
+  Gamma::Algebra::MinusSigmaXT,
+  Gamma::Algebra::SigmaXY,
+  Gamma::Algebra::MinusSigmaXY,
+  Gamma::Algebra::SigmaXZ,
+  Gamma::Algebra::MinusSigmaXZ,
+  Gamma::Algebra::SigmaYT,
+  Gamma::Algebra::MinusSigmaYT,
+  Gamma::Algebra::SigmaYZ,
+  Gamma::Algebra::MinusSigmaYZ,
+  Gamma::Algebra::SigmaZT,
+  Gamma::Algebra::MinusSigmaZT
+};
diff --git a/lib/qcd/spin/gamma-gen/gamma-gen.nb b/lib/qcd/spin/gamma-gen/gamma-gen.nb
new file mode 100644
index 00000000..d41b3e2e
--- /dev/null
+++ b/lib/qcd/spin/gamma-gen/gamma-gen.nb
@@ -0,0 +1,1520 @@
+(* Content-type: application/vnd.wolfram.mathematica *)
+
+(*** Wolfram Notebook File ***)
+(* http://www.wolfram.com/nb *)
+
+(* CreatedBy='Mathematica 11.0' *)
+
+(*CacheID: 234*)
+(* Internal cache information:
+NotebookFileLineBreakTest
+NotebookFileLineBreakTest
+NotebookDataPosition[       158,          7]
+NotebookDataLength[     58006,       1512]
+NotebookOptionsPosition[     56426,       1458]
+NotebookOutlinePosition[     56786,       1474]
+CellTagsIndexPosition[     56743,       1471]
+WindowFrame->Normal*)
+
+(* Beginning of Notebook Content *)
+Notebook[{
+Cell[TextData[StyleBox["Grid physics library, www.github.com/paboyle/Grid \n\n\
+Source file: lib/qcd/spin/gamma-gen/gamma-gen.nb\n\nCopyright (C) 2015\n\
+Copyright (C) 2016\nCopyright (C) 2017\n\nAuthor: Antonin Portelli \
+<antonin.portelli@me.com>\n\nThis program is free software; you can \
+redistribute it and/or modify\nit under the terms of the GNU General Public \
+License as published by\nthe Free Software Foundation; either version 2 of \
+the License, or\n(at your option) any later version.\n\nThis program is \
+distributed in the hope that it will be useful,\nbut WITHOUT ANY WARRANTY; \
+without even the implied warranty of\nMERCHANTABILITY or FITNESS FOR A \
+PARTICULAR PURPOSE.  See the\nGNU General Public License for more details.\n\n\
+You should have received a copy of the GNU General Public License along\nwith \
+this program; if not, write to the Free Software Foundation, Inc.,\n51 \
+Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.\n\nSee the full \
+license in the file \[OpenCurlyDoubleQuote]LICENSE\[CloseCurlyDoubleQuote] in \
+the top level distribution directory",
+ Background->RGBColor[1, 1, 0.85]]], "Text",
+ CellChangeTimes->{{3.694192702097444*^9, 3.694192776106186*^9}, {
+  3.6942089129622507`*^9, 3.694208916624515*^9}},
+ Background->RGBColor[1, 1, 0.85]],
+
+Cell[CellGroupData[{
+
+Cell[BoxData[{
+ RowBox[{"ClearAll", "[", "\"\<Global`*\>\"", "]"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"Unprotect", "[", "\"\<Dot\>\"", "]"}], ";", 
+  RowBox[{"Unprotect", "[", "\"\<ConjugateTranspose\>\"", "]"}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{"SetDirectory", "[", 
+  RowBox[{"NotebookDirectory", "[", "]"}], "]"}]}], "Input",
+ CellChangeTimes->{{3.69418610909842*^9, 3.694186122331771*^9}, {
+  3.694189694542165*^9, 3.6941897146300087`*^9}, {3.694297413625847*^9, 
+  3.694297419003489*^9}}],
+
+Cell[BoxData["\<\"/Users/antonin/Development/Grid/lib/qcd/spin/gamma-gen\"\>"]\
+, "Output",
+ CellChangeTimes->{{3.694189715067813*^9, 3.694189718130962*^9}, {
+   3.694207343654113*^9, 3.694207368715312*^9}, 3.694208353902673*^9, 
+   3.694208422714191*^9, 3.694208720092059*^9, 3.694208887255764*^9, 
+   3.6942115023319197`*^9, 3.694211674672401*^9, 3.694211769773765*^9, 
+   3.694212832251072*^9, 3.694212924703574*^9, 3.694213477890645*^9, 
+   3.694214215876555*^9, 3.694214744582322*^9, 3.694214868607993*^9, {
+   3.694214930774279*^9, 3.694214956189962*^9}, 3.69421520097348*^9, 
+   3.6942153028002653`*^9, 3.694215388363636*^9, 3.694215531793765*^9, 
+   3.694293467607703*^9, 3.694294096138891*^9, 3.694297419616928*^9, 
+   3.694297655191285*^9, 3.69429782738804*^9, 3.694297901657515*^9, 
+   3.6942984168973837`*^9}]
+}, Open  ]],
+
+Cell[CellGroupData[{
+
+Cell["Clifford algebra generation", "Section",
+ CellChangeTimes->{{3.6942089434583883`*^9, 3.694208978559093*^9}}],
+
+Cell[BoxData[{
+ RowBox[{
+  RowBox[{"com", "[", 
+   RowBox[{"a_", ",", "b_"}], "]"}], ":=", 
+  RowBox[{
+   RowBox[{"a", ".", "b"}], "-", 
+   RowBox[{"b", ".", "a"}]}]}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"acom", "[", 
+   RowBox[{"a_", ",", "b_"}], "]"}], ":=", 
+  RowBox[{
+   RowBox[{"a", ".", "b"}], "+", 
+   RowBox[{"b", ".", "a"}]}]}]}], "Input",
+ CellChangeTimes->{{3.694184330267939*^9, 3.694184337479828*^9}, {
+  3.694184821238667*^9, 3.6941848260602217`*^9}}],
+
+Cell["Definition of the matrix representation of the algebra:", "Text",
+ CellChangeTimes->{{3.6942090405172586`*^9, 3.694209073962101*^9}, {
+  3.6942974330697393`*^9, 3.694297433821431*^9}}],
+
+Cell[BoxData[{
+ RowBox[{
+  RowBox[{
+   RowBox[{"mat", "[", "id", "]"}], "=", 
+   RowBox[{"IdentityMatrix", "[", "4", "]"}]}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"mat", "[", "gx", "]"}], "=", 
+   RowBox[{"(", GridBox[{
+      {"0", "0", "0", "\[ImaginaryI]"},
+      {"0", "0", "\[ImaginaryI]", "0"},
+      {"0", 
+       RowBox[{"-", "\[ImaginaryI]"}], "0", "0"},
+      {
+       RowBox[{"-", "\[ImaginaryI]"}], "0", "0", "0"}
+     }], ")"}]}], ";", 
+  RowBox[{
+   RowBox[{"mat", "[", "gy", "]"}], "=", 
+   RowBox[{"(", GridBox[{
+      {"0", "0", "0", 
+       RowBox[{"-", "1"}]},
+      {"0", "0", "1", "0"},
+      {"0", "1", "0", "0"},
+      {
+       RowBox[{"-", "1"}], "0", "0", "0"}
+     }], ")"}]}], ";", 
+  RowBox[{
+   RowBox[{"mat", "[", "gz", "]"}], "=", 
+   RowBox[{"(", GridBox[{
+      {"0", "0", "\[ImaginaryI]", "0"},
+      {"0", "0", "0", 
+       RowBox[{"-", "\[ImaginaryI]"}]},
+      {
+       RowBox[{"-", "\[ImaginaryI]"}], "0", "0", "0"},
+      {"0", "\[ImaginaryI]", "0", "0"}
+     }], ")"}]}], ";", 
+  RowBox[{
+   RowBox[{"mat", "[", "gt", "]"}], "=", 
+   RowBox[{"(", GridBox[{
+      {"0", "0", "1", "0"},
+      {"0", "0", "0", "1"},
+      {"1", "0", "0", "0"},
+      {"0", "1", "0", "0"}
+     }], ")"}]}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"mat", "[", "g5", "]"}], "=", 
+   RowBox[{
+    RowBox[{"mat", "[", "gx", "]"}], ".", 
+    RowBox[{"mat", "[", "gy", "]"}], ".", 
+    RowBox[{"mat", "[", "gz", "]"}], ".", 
+    RowBox[{"mat", "[", "gt", "]"}]}]}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"mat", "[", "gx5", "]"}], "=", 
+   RowBox[{
+    RowBox[{"mat", "[", "gx", "]"}], ".", 
+    RowBox[{"mat", "[", "g5", "]"}]}]}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"mat", "[", "gy5", "]"}], "=", 
+   RowBox[{
+    RowBox[{"mat", "[", "gy", "]"}], ".", 
+    RowBox[{"mat", "[", "g5", "]"}]}]}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"mat", "[", "gz5", "]"}], "=", 
+   RowBox[{
+    RowBox[{"mat", "[", "gz", "]"}], ".", 
+    RowBox[{"mat", "[", "g5", "]"}]}]}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"mat", "[", "gt5", "]"}], "=", 
+   RowBox[{
+    RowBox[{"mat", "[", "gt", "]"}], ".", 
+    RowBox[{"mat", "[", "g5", "]"}]}]}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"mat", "[", "\[Sigma]xy", "]"}], "=", 
+   RowBox[{
+    RowBox[{"com", "[", 
+     RowBox[{
+      RowBox[{"mat", "[", "gx", "]"}], ",", 
+      RowBox[{"mat", "[", "gy", "]"}]}], "]"}], "/", "2"}]}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"mat", "[", "\[Sigma]xz", "]"}], "=", 
+   RowBox[{
+    RowBox[{"com", "[", 
+     RowBox[{
+      RowBox[{"mat", "[", "gx", "]"}], ",", 
+      RowBox[{"mat", "[", "gz", "]"}]}], "]"}], "/", "2"}]}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"mat", "[", "\[Sigma]xt", "]"}], "=", 
+   RowBox[{
+    RowBox[{"com", "[", 
+     RowBox[{
+      RowBox[{"mat", "[", "gx", "]"}], ",", 
+      RowBox[{"mat", "[", "gt", "]"}]}], "]"}], "/", "2"}]}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"mat", "[", "\[Sigma]yz", "]"}], "=", 
+   RowBox[{
+    RowBox[{"com", "[", 
+     RowBox[{
+      RowBox[{"mat", "[", "gy", "]"}], ",", 
+      RowBox[{"mat", "[", "gz", "]"}]}], "]"}], "/", "2"}]}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"mat", "[", "\[Sigma]yt", "]"}], "=", 
+   RowBox[{
+    RowBox[{"com", "[", 
+     RowBox[{
+      RowBox[{"mat", "[", "gy", "]"}], ",", 
+      RowBox[{"mat", "[", "gt", "]"}]}], "]"}], "/", "2"}]}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"mat", "[", "\[Sigma]zt", "]"}], "=", 
+   RowBox[{
+    RowBox[{"com", "[", 
+     RowBox[{
+      RowBox[{"mat", "[", "gz", "]"}], ",", 
+      RowBox[{"mat", "[", "gt", "]"}]}], "]"}], "/", "2"}]}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"palg", "=", 
+   RowBox[{"{", 
+    RowBox[{
+    "id", ",", "gx", ",", "gy", ",", "gz", ",", "gt", ",", "g5", ",", "gx5", 
+     ",", "gy5", ",", "gz5", ",", "gt5", ",", "\[Sigma]xy", ",", "\[Sigma]xz",
+      ",", "\[Sigma]xt", ",", "\[Sigma]yz", ",", "\[Sigma]yt", ",", 
+     "\[Sigma]zt"}], "}"}]}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"malg", "=", 
+   RowBox[{"Table", "[", 
+    RowBox[{
+     RowBox[{"-", "i"}], ",", 
+     RowBox[{"{", 
+      RowBox[{"i", ",", "palg"}], "}"}]}], "]"}]}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"alg", "=", 
+   RowBox[{"Union", "[", 
+    RowBox[{"palg", ",", "malg"}], "]"}]}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{"Do", "[", 
+  RowBox[{
+   RowBox[{
+    RowBox[{"mat", "[", 
+     RowBox[{"-", "i"}], "]"}], "=", 
+    RowBox[{"-", 
+     RowBox[{"mat", "[", "i", "]"}]}]}], ",", 
+   RowBox[{"{", 
+    RowBox[{"i", ",", "palg"}], "}"}]}], "]"}]}], "Input",
+ CellChangeTimes->{{3.694182703426999*^9, 3.694182706401226*^9}, {
+   3.6941838790187283`*^9, 3.694184129323387*^9}, {3.694184186220387*^9, 
+   3.6941842073992023`*^9}, {3.694184241024517*^9, 3.6941843050620823`*^9}, {
+   3.694184341303878*^9, 3.694184431096826*^9}, {3.694184469241919*^9, 
+   3.694184532366125*^9}, {3.694184830238468*^9, 3.694184950533866*^9}, {
+   3.694185044355978*^9, 3.694185099415689*^9}, {3.694185502749824*^9, 
+   3.694185675128971*^9}, {3.694185728773429*^9, 3.694185729056695*^9}, 
+   3.694185780274218*^9, 3.6941858224264593`*^9, {3.694185941282981*^9, 
+   3.694185950262871*^9}}],
+
+Cell["Generation of the abstract algebra:", "Text",
+ CellChangeTimes->{{3.6942090658330803`*^9, 3.694209076132119*^9}}],
+
+Cell[BoxData[{
+ RowBox[{"Do", "[", 
+  RowBox[{
+   RowBox[{
+    RowBox[{"Dot", "[", 
+     RowBox[{"i1", ",", "i2"}], "]"}], "=", 
+    RowBox[{
+     RowBox[{"Select", "[", 
+      RowBox[{"alg", ",", 
+       RowBox[{
+        RowBox[{
+         RowBox[{
+          RowBox[{"mat", "[", "i1", "]"}], ".", 
+          RowBox[{"mat", "[", "i2", "]"}]}], "\[Equal]", 
+         RowBox[{"mat", "[", "#", "]"}]}], "&"}]}], "]"}], "[", 
+     RowBox[{"[", "1", "]"}], "]"}]}], ",", 
+   RowBox[{"{", 
+    RowBox[{"i1", ",", "palg"}], "}"}], ",", 
+   RowBox[{"{", 
+    RowBox[{"i2", ",", "palg"}], "}"}]}], "]"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"Dot", "[", 
+    RowBox[{
+     RowBox[{"-", "i1_"}], ",", "i2_"}], "]"}], ":=", 
+   RowBox[{"-", 
+    RowBox[{"i1", ".", "i2"}]}]}], ";", 
+  RowBox[{
+   RowBox[{"Dot", "[", 
+    RowBox[{"i1_", ",", 
+     RowBox[{"-", "i2_"}]}], "]"}], ":=", 
+   RowBox[{"-", 
+    RowBox[{"i1", ".", "i2"}]}]}], ";", 
+  RowBox[{
+   RowBox[{"Dot", "[", 
+    RowBox[{
+     RowBox[{"-", "i1_"}], ",", 
+     RowBox[{"-", "i2_"}]}], "]"}], ":=", 
+   RowBox[{"i1", ".", "i2"}]}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"Do", "[", 
+   RowBox[{
+    RowBox[{
+     RowBox[{"ConjugateTranspose", "[", "i", "]"}], "=", 
+     RowBox[{
+      RowBox[{"Select", "[", 
+       RowBox[{"alg", ",", 
+        RowBox[{
+         RowBox[{
+          RowBox[{"ConjugateTranspose", "[", 
+           RowBox[{"mat", "[", "i", "]"}], "]"}], "\[Equal]", 
+          RowBox[{"mat", "[", "#", "]"}]}], "&"}]}], "]"}], "[", 
+      RowBox[{"[", "1", "]"}], "]"}]}], ",", 
+    RowBox[{"{", 
+     RowBox[{"i", ",", "palg"}], "}"}]}], "]"}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"ConjugateTranspose", "[", 
+   RowBox[{"-", "i_"}], "]"}], ":=", 
+  RowBox[{"-", 
+   RowBox[{"ConjugateTranspose", "[", "i", "]"}]}]}]}], "Input",
+ CellChangeTimes->{{3.6941860329437103`*^9, 3.6941860343133917`*^9}, {
+  3.694186163571176*^9, 3.6941862016761427`*^9}, {3.69418700219066*^9, 
+  3.6941870425469627`*^9}, {3.694297326197534*^9, 3.6942974062629423`*^9}, {
+  3.694297634175386*^9, 3.6942976496897383`*^9}}],
+
+Cell["Check that we can reconstruct the Euclidean metric:", "Text",
+ CellChangeTimes->{{3.694209080190936*^9, 3.694209096585559*^9}}],
+
+Cell[CellGroupData[{
+
+Cell[BoxData[
+ RowBox[{
+  RowBox[{"Table", "[", 
+   RowBox[{
+    RowBox[{
+     RowBox[{"acom", "[", 
+      RowBox[{"i1", ",", "i2"}], "]"}], "/", "2"}], ",", 
+    RowBox[{"{", 
+     RowBox[{"i1", ",", 
+      RowBox[{"{", 
+       RowBox[{"gx", ",", "gy", ",", "gz", ",", "gt"}], "}"}]}], "}"}], ",", 
+    RowBox[{"{", 
+     RowBox[{"i2", ",", 
+      RowBox[{"{", 
+       RowBox[{"gx", ",", "gy", ",", "gz", ",", "gt"}], "}"}]}], "}"}]}], 
+   "]"}], "//", "MatrixForm"}]], "Input"],
+
+Cell[BoxData[
+ TagBox[
+  RowBox[{"(", "\[NoBreak]", GridBox[{
+     {"id", "0", "0", "0"},
+     {"0", "id", "0", "0"},
+     {"0", "0", "id", "0"},
+     {"0", "0", "0", "id"}
+    },
+    GridBoxAlignment->{
+     "Columns" -> {{Center}}, "ColumnsIndexed" -> {}, "Rows" -> {{Baseline}}, 
+      "RowsIndexed" -> {}},
+    GridBoxSpacings->{"Columns" -> {
+        Offset[0.27999999999999997`], {
+         Offset[0.7]}, 
+        Offset[0.27999999999999997`]}, "ColumnsIndexed" -> {}, "Rows" -> {
+        Offset[0.2], {
+         Offset[0.4]}, 
+        Offset[0.2]}, "RowsIndexed" -> {}}], "\[NoBreak]", ")"}],
+  Function[BoxForm`e$, 
+   MatrixForm[BoxForm`e$]]]], "Output",
+ CellChangeTimes->{
+  3.694186559657666*^9, 3.6941870440388947`*^9, 3.694189718434745*^9, {
+   3.694207344033985*^9, 3.694207369106125*^9}, 3.694208354345052*^9, 
+   3.694208423089445*^9, 3.6942087204812307`*^9, 3.6942088876269073`*^9, 
+   3.694211502659554*^9, 3.6942116750282307`*^9, 3.694211770126875*^9, 
+   3.6942128326290417`*^9, 3.6942129251958523`*^9, 3.694213478493619*^9, 
+   3.694214216283064*^9, 3.694214744924881*^9, 3.694214869005068*^9, {
+   3.694214931160738*^9, 3.6942149565175247`*^9}, 3.69421520140513*^9, 
+   3.694215303184371*^9, 3.694215388782826*^9, 3.694215532274271*^9, 
+   3.694293468120837*^9, 3.69429409699356*^9, 3.694297655587861*^9, 
+   3.694297827778974*^9, 3.694297902050129*^9, 3.694298417447421*^9}]
+}, Open  ]],
+
+Cell["Full multiplication table:", "Text",
+ CellChangeTimes->{{3.694209113187169*^9, 3.6942091210767593`*^9}}],
+
+Cell[CellGroupData[{
+
+Cell[BoxData[
+ RowBox[{
+  RowBox[{"Table", "[", 
+   RowBox[{
+    RowBox[{"i1", ".", "i2"}], ",", 
+    RowBox[{"{", 
+     RowBox[{"i1", ",", "palg"}], "}"}], ",", 
+    RowBox[{"{", 
+     RowBox[{"i2", ",", "palg"}], "}"}]}], "]"}], "//", 
+  "MatrixForm"}]], "Input",
+ CellChangeTimes->{{3.6941862426584797`*^9, 3.694186256858178*^9}, {
+  3.694186605271886*^9, 3.694186617894228*^9}, {3.694186972131384*^9, 
+  3.69418697419895*^9}, {3.694192885918524*^9, 3.694192888888296*^9}}],
+
+Cell[BoxData[
+ TagBox[
+  RowBox[{"(", "\[NoBreak]", GridBox[{
+     {"id", "gx", "gy", "gz", "gt", "g5", "gx5", "gy5", "gz5", "gt5", 
+      "\[Sigma]xy", "\[Sigma]xz", "\[Sigma]xt", "\[Sigma]yz", "\[Sigma]yt", 
+      "\[Sigma]zt"},
+     {"gx", "id", "\[Sigma]xy", "\[Sigma]xz", "\[Sigma]xt", "gx5", "g5", 
+      RowBox[{"-", "\[Sigma]zt"}], "\[Sigma]yt", 
+      RowBox[{"-", "\[Sigma]yz"}], "gy", "gz", "gt", 
+      RowBox[{"-", "gt5"}], "gz5", 
+      RowBox[{"-", "gy5"}]},
+     {"gy", 
+      RowBox[{"-", "\[Sigma]xy"}], "id", "\[Sigma]yz", "\[Sigma]yt", "gy5", 
+      "\[Sigma]zt", "g5", 
+      RowBox[{"-", "\[Sigma]xt"}], "\[Sigma]xz", 
+      RowBox[{"-", "gx"}], "gt5", 
+      RowBox[{"-", "gz5"}], "gz", "gt", "gx5"},
+     {"gz", 
+      RowBox[{"-", "\[Sigma]xz"}], 
+      RowBox[{"-", "\[Sigma]yz"}], "id", "\[Sigma]zt", "gz5", 
+      RowBox[{"-", "\[Sigma]yt"}], "\[Sigma]xt", "g5", 
+      RowBox[{"-", "\[Sigma]xy"}], 
+      RowBox[{"-", "gt5"}], 
+      RowBox[{"-", "gx"}], "gy5", 
+      RowBox[{"-", "gy"}], 
+      RowBox[{"-", "gx5"}], "gt"},
+     {"gt", 
+      RowBox[{"-", "\[Sigma]xt"}], 
+      RowBox[{"-", "\[Sigma]yt"}], 
+      RowBox[{"-", "\[Sigma]zt"}], "id", "gt5", "\[Sigma]yz", 
+      RowBox[{"-", "\[Sigma]xz"}], "\[Sigma]xy", "g5", "gz5", 
+      RowBox[{"-", "gy5"}], 
+      RowBox[{"-", "gx"}], "gx5", 
+      RowBox[{"-", "gy"}], 
+      RowBox[{"-", "gz"}]},
+     {"g5", 
+      RowBox[{"-", "gx5"}], 
+      RowBox[{"-", "gy5"}], 
+      RowBox[{"-", "gz5"}], 
+      RowBox[{"-", "gt5"}], "id", 
+      RowBox[{"-", "gx"}], 
+      RowBox[{"-", "gy"}], 
+      RowBox[{"-", "gz"}], 
+      RowBox[{"-", "gt"}], 
+      RowBox[{"-", "\[Sigma]zt"}], "\[Sigma]yt", 
+      RowBox[{"-", "\[Sigma]yz"}], 
+      RowBox[{"-", "\[Sigma]xt"}], "\[Sigma]xz", 
+      RowBox[{"-", "\[Sigma]xy"}]},
+     {"gx5", 
+      RowBox[{"-", "g5"}], "\[Sigma]zt", 
+      RowBox[{"-", "\[Sigma]yt"}], "\[Sigma]yz", "gx", 
+      RowBox[{"-", "id"}], 
+      RowBox[{"-", "\[Sigma]xy"}], 
+      RowBox[{"-", "\[Sigma]xz"}], 
+      RowBox[{"-", "\[Sigma]xt"}], "gy5", "gz5", "gt5", 
+      RowBox[{"-", "gt"}], "gz", 
+      RowBox[{"-", "gy"}]},
+     {"gy5", 
+      RowBox[{"-", "\[Sigma]zt"}], 
+      RowBox[{"-", "g5"}], "\[Sigma]xt", 
+      RowBox[{"-", "\[Sigma]xz"}], "gy", "\[Sigma]xy", 
+      RowBox[{"-", "id"}], 
+      RowBox[{"-", "\[Sigma]yz"}], 
+      RowBox[{"-", "\[Sigma]yt"}], 
+      RowBox[{"-", "gx5"}], "gt", 
+      RowBox[{"-", "gz"}], "gz5", "gt5", "gx"},
+     {"gz5", "\[Sigma]yt", 
+      RowBox[{"-", "\[Sigma]xt"}], 
+      RowBox[{"-", "g5"}], "\[Sigma]xy", "gz", "\[Sigma]xz", "\[Sigma]yz", 
+      RowBox[{"-", "id"}], 
+      RowBox[{"-", "\[Sigma]zt"}], 
+      RowBox[{"-", "gt"}], 
+      RowBox[{"-", "gx5"}], "gy", 
+      RowBox[{"-", "gy5"}], 
+      RowBox[{"-", "gx"}], "gt5"},
+     {"gt5", 
+      RowBox[{"-", "\[Sigma]yz"}], "\[Sigma]xz", 
+      RowBox[{"-", "\[Sigma]xy"}], 
+      RowBox[{"-", "g5"}], "gt", "\[Sigma]xt", "\[Sigma]yt", "\[Sigma]zt", 
+      RowBox[{"-", "id"}], "gz", 
+      RowBox[{"-", "gy"}], 
+      RowBox[{"-", "gx5"}], "gx", 
+      RowBox[{"-", "gy5"}], 
+      RowBox[{"-", "gz5"}]},
+     {"\[Sigma]xy", 
+      RowBox[{"-", "gy"}], "gx", 
+      RowBox[{"-", "gt5"}], "gz5", 
+      RowBox[{"-", "\[Sigma]zt"}], 
+      RowBox[{"-", "gy5"}], "gx5", 
+      RowBox[{"-", "gt"}], "gz", 
+      RowBox[{"-", "id"}], 
+      RowBox[{"-", "\[Sigma]yz"}], 
+      RowBox[{"-", "\[Sigma]yt"}], "\[Sigma]xz", "\[Sigma]xt", "g5"},
+     {"\[Sigma]xz", 
+      RowBox[{"-", "gz"}], "gt5", "gx", 
+      RowBox[{"-", "gy5"}], "\[Sigma]yt", 
+      RowBox[{"-", "gz5"}], "gt", "gx5", 
+      RowBox[{"-", "gy"}], "\[Sigma]yz", 
+      RowBox[{"-", "id"}], 
+      RowBox[{"-", "\[Sigma]zt"}], 
+      RowBox[{"-", "\[Sigma]xy"}], 
+      RowBox[{"-", "g5"}], "\[Sigma]xt"},
+     {"\[Sigma]xt", 
+      RowBox[{"-", "gt"}], 
+      RowBox[{"-", "gz5"}], "gy5", "gx", 
+      RowBox[{"-", "\[Sigma]yz"}], 
+      RowBox[{"-", "gt5"}], 
+      RowBox[{"-", "gz"}], "gy", "gx5", "\[Sigma]yt", "\[Sigma]zt", 
+      RowBox[{"-", "id"}], "g5", 
+      RowBox[{"-", "\[Sigma]xy"}], 
+      RowBox[{"-", "\[Sigma]xz"}]},
+     {"\[Sigma]yz", 
+      RowBox[{"-", "gt5"}], 
+      RowBox[{"-", "gz"}], "gy", "gx5", 
+      RowBox[{"-", "\[Sigma]xt"}], 
+      RowBox[{"-", "gt"}], 
+      RowBox[{"-", "gz5"}], "gy5", "gx", 
+      RowBox[{"-", "\[Sigma]xz"}], "\[Sigma]xy", "g5", 
+      RowBox[{"-", "id"}], 
+      RowBox[{"-", "\[Sigma]zt"}], "\[Sigma]yt"},
+     {"\[Sigma]yt", "gz5", 
+      RowBox[{"-", "gt"}], 
+      RowBox[{"-", "gx5"}], "gy", "\[Sigma]xz", "gz", 
+      RowBox[{"-", "gt5"}], 
+      RowBox[{"-", "gx"}], "gy5", 
+      RowBox[{"-", "\[Sigma]xt"}], 
+      RowBox[{"-", "g5"}], "\[Sigma]xy", "\[Sigma]zt", 
+      RowBox[{"-", "id"}], 
+      RowBox[{"-", "\[Sigma]yz"}]},
+     {"\[Sigma]zt", 
+      RowBox[{"-", "gy5"}], "gx5", 
+      RowBox[{"-", "gt"}], "gz", 
+      RowBox[{"-", "\[Sigma]xy"}], 
+      RowBox[{"-", "gy"}], "gx", 
+      RowBox[{"-", "gt5"}], "gz5", "g5", 
+      RowBox[{"-", "\[Sigma]xt"}], "\[Sigma]xz", 
+      RowBox[{"-", "\[Sigma]yt"}], "\[Sigma]yz", 
+      RowBox[{"-", "id"}]}
+    },
+    GridBoxAlignment->{
+     "Columns" -> {{Center}}, "ColumnsIndexed" -> {}, "Rows" -> {{Baseline}}, 
+      "RowsIndexed" -> {}},
+    GridBoxSpacings->{"Columns" -> {
+        Offset[0.27999999999999997`], {
+         Offset[0.7]}, 
+        Offset[0.27999999999999997`]}, "ColumnsIndexed" -> {}, "Rows" -> {
+        Offset[0.2], {
+         Offset[0.4]}, 
+        Offset[0.2]}, "RowsIndexed" -> {}}], "\[NoBreak]", ")"}],
+  Function[BoxForm`e$, 
+   MatrixForm[BoxForm`e$]]]], "Output",
+ CellChangeTimes->{{3.694186251703266*^9, 3.6941862573246403`*^9}, {
+   3.69418660949746*^9, 3.694186618329227*^9}, 3.6941869746472473`*^9, 
+   3.694187045160287*^9, 3.69418971854458*^9, 3.694192889450261*^9, {
+   3.694207344114737*^9, 3.694207369223531*^9}, 3.694208354449048*^9, 
+   3.694208423183125*^9, 3.6942087205954943`*^9, 3.694208887707542*^9, 
+   3.694211502726407*^9, 3.694211675104418*^9, 3.694211770202587*^9, 
+   3.694212832700527*^9, 3.69421292529856*^9, 3.694213478604137*^9, 
+   3.694214216368701*^9, 3.694214744992138*^9, 3.694214869090289*^9, {
+   3.694214931237383*^9, 3.6942149565897913`*^9}, 3.6942152015279818`*^9, 
+   3.694215303248364*^9, 3.694215388866436*^9, 3.694215532359406*^9, 
+   3.694293468220346*^9, 3.694294097121203*^9, 3.694297655690278*^9, 
+   3.694297827858612*^9, 3.694297902125073*^9, 3.694298417530912*^9}]
+}, Open  ]]
+}, Open  ]],
+
+Cell[CellGroupData[{
+
+Cell["Header file Dirac.h generation", "Section",
+ CellChangeTimes->{{3.694208986784461*^9, 3.6942090005062523`*^9}}],
+
+Cell["File skeleton:", "Text",
+ CellFrame->{{0, 0}, {0, 0.5}},
+ CellChangeTimes->{{3.694209131604498*^9, 3.694209133792495*^9}}],
+
+Cell[BoxData[
+ RowBox[{
+  RowBox[{"makeHeaderCode", "[", 
+   RowBox[{"enumCode_", ",", "funcCode_"}], "]"}], ":=", 
+  RowBox[{"Module", "[", 
+   RowBox[{
+    RowBox[{"{", 
+     RowBox[{"out", "=", "\"\<\>\""}], "}"}], ",", "\[IndentingNewLine]", 
+    RowBox[{
+     RowBox[{"out", " ", "=", " ", 
+      RowBox[{
+      "out", "<>", "\[IndentingNewLine]", 
+       "\"\<#ifndef GRID_QCD_DIRAC_H\n#define GRID_QCD_DIRAC_H\n\n// code \
+generated by the Mathematica notebook gamma-gen/gamma-gen.nb\n\n#include \
+<array>\n\nnamespace Grid {\nnamespace QCD {\n\nclass Gamma {\n  public:\n\
+\>\""}]}], ";", "\[IndentingNewLine]", 
+     RowBox[{"out", " ", "=", " ", 
+      RowBox[{"out", "<>", "enumCode"}]}], ";", "\[IndentingNewLine]", 
+     RowBox[{"out", "=", 
+      RowBox[{
+      "out", "<>", 
+       "\"\<    static const std::array<const char *, nGamma>                \
+name;\n    static const std::array<std::array<Algebra, nGamma>, nGamma> mul;\n\
+    static const std::array<Algebra, nGamma>                     adj;\n    \
+Algebra                                                      g;\n  public:\n  \
+  Gamma(Algebra initg): g(initg) {}  \n};\n\ninline Gamma operator*(const \
+Gamma &g1, const Gamma &g2)\n{\n  return Gamma(Gamma::mul[g1.g][g2.g]);\n}\n\n\
+inline Gamma adj(const Gamma &g)\n{\n  return Gamma(Gamma::adj[g.g]);\n}\n\n\
+\>\""}]}], ";", "\[IndentingNewLine]", 
+     RowBox[{"out", " ", "=", 
+      RowBox[{"out", "<>", "funcCode"}]}], ";", "\[IndentingNewLine]", 
+     RowBox[{"out", "=", 
+      RowBox[{"out", "<>", "\"\<\n}}\n\n#endif\n\>\""}]}], ";", 
+     "\[IndentingNewLine]", 
+     RowBox[{"Return", "[", "out", "]"}]}]}], "]"}]}]], "Input",
+ CellChangeTimes->{{3.69419645529875*^9, 3.694196727399803*^9}, {
+   3.694196786125751*^9, 3.694196882528638*^9}, {3.694197240621915*^9, 
+   3.694197350313327*^9}, {3.694206801774742*^9, 3.6942068152492313`*^9}, {
+   3.694206953106125*^9, 3.694206955712624*^9}, 3.694207329784376*^9, {
+   3.694207361301278*^9, 3.694207363264554*^9}, {3.694207717624036*^9, 
+   3.694207756242198*^9}, {3.694208382782*^9, 3.694208417198256*^9}, {
+   3.694208503080744*^9, 3.6942086978552227`*^9}, {3.6942113930386953`*^9, 
+   3.694211409741232*^9}, {3.6942115367346992`*^9, 3.694211695577471*^9}, {
+   3.694211732616888*^9, 3.69421173594759*^9}, {3.694212824802087*^9, 
+   3.69421282832565*^9}, {3.69421325181595*^9, 3.6942133005791903`*^9}, {
+   3.694214209845763*^9, 3.6942142108217773`*^9}, {3.694214621189205*^9, 
+   3.694214647829077*^9}, {3.694214730808668*^9, 3.6942148593149223`*^9}, {
+   3.694214900404429*^9, 3.694214905595294*^9}, {3.6942155281289682`*^9, 
+   3.6942155289376802`*^9}, {3.694297793776082*^9, 3.694297819122456*^9}, {
+   3.694297854809973*^9, 3.694297889420939*^9}}],
+
+Cell["Gamma enum generation:", "Text",
+ CellFrame->{{0, 0}, {0, 0.5}},
+ CellChangeTimes->{{3.694209168488991*^9, 3.6942091715073423`*^9}, 
+   3.694209215969149*^9}],
+
+Cell[BoxData[{
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", "id", "]"}], "=", "\"\<Algebra::Identity\>\""}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", "gx", "]"}], "=", "\"\<Algebra::GammaX\>\""}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", "gy", "]"}], "=", "\"\<Algebra::GammaY\>\""}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", "gz", "]"}], "=", "\"\<Algebra::GammaZ\>\""}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", "gt", "]"}], "=", "\"\<Algebra::GammaT\>\""}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", "g5", "]"}], "=", "\"\<Algebra::Gamma5\>\""}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", "gx5", "]"}], "=", "\"\<Algebra::GammaXGamma5\>\""}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", "gy5", "]"}], "=", "\"\<Algebra::GammaYGamma5\>\""}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", "gz5", "]"}], "=", "\"\<Algebra::GammaZGamma5\>\""}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", "gt5", "]"}], "=", "\"\<Algebra::GammaTGamma5\>\""}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", "\[Sigma]xy", "]"}], "=", 
+   "\"\<Algebra::SigmaXY\>\""}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", "\[Sigma]xz", "]"}], "=", 
+   "\"\<Algebra::SigmaXZ\>\""}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", "\[Sigma]xt", "]"}], "=", 
+   "\"\<Algebra::SigmaXT\>\""}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", "\[Sigma]yz", "]"}], "=", 
+   "\"\<Algebra::SigmaYZ\>\""}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", "\[Sigma]yt", "]"}], "=", 
+   "\"\<Algebra::SigmaYT\>\""}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", "\[Sigma]zt", "]"}], "=", 
+   "\"\<Algebra::SigmaZT\>\""}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", 
+    RowBox[{"-", "id"}], "]"}], "=", "\"\<Algebra::MinusIdentity\>\""}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", 
+    RowBox[{"-", "gx"}], "]"}], "=", "\"\<Algebra::MinusGammaX\>\""}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", 
+    RowBox[{"-", "gy"}], "]"}], "=", "\"\<Algebra::MinusGammaY\>\""}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", 
+    RowBox[{"-", "gz"}], "]"}], "=", "\"\<Algebra::MinusGammaZ\>\""}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", 
+    RowBox[{"-", "gt"}], "]"}], "=", "\"\<Algebra::MinusGammaT\>\""}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", 
+    RowBox[{"-", "g5"}], "]"}], "=", "\"\<Algebra::MinusGamma5\>\""}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", 
+    RowBox[{"-", "gx5"}], "]"}], "=", "\"\<Algebra::MinusGammaXGamma5\>\""}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", 
+    RowBox[{"-", "gy5"}], "]"}], "=", "\"\<Algebra::MinusGammaYGamma5\>\""}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", 
+    RowBox[{"-", "gz5"}], "]"}], "=", "\"\<Algebra::MinusGammaZGamma5\>\""}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", 
+    RowBox[{"-", "gt5"}], "]"}], "=", "\"\<Algebra::MinusGammaTGamma5\>\""}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", 
+    RowBox[{"-", "\[Sigma]xy"}], "]"}], "=", 
+   "\"\<Algebra::MinusSigmaXY\>\""}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", 
+    RowBox[{"-", "\[Sigma]xz"}], "]"}], "=", 
+   "\"\<Algebra::MinusSigmaXZ\>\""}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", 
+    RowBox[{"-", "\[Sigma]xt"}], "]"}], "=", 
+   "\"\<Algebra::MinusSigmaXT\>\""}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", 
+    RowBox[{"-", "\[Sigma]yz"}], "]"}], "=", 
+   "\"\<Algebra::MinusSigmaYZ\>\""}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", 
+    RowBox[{"-", "\[Sigma]yt"}], "]"}], "=", 
+   "\"\<Algebra::MinusSigmaYT\>\""}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"enum", "[", 
+    RowBox[{"-", "\[Sigma]zt"}], "]"}], "=", 
+   "\"\<Algebra::MinusSigmaZT\>\""}], ";"}]}], "Input",
+ CellChangeTimes->{{3.69418665896658*^9, 3.6941867305497723`*^9}, {
+  3.694186782865391*^9, 3.694186840513199*^9}, {3.694186889568404*^9, 
+  3.694186968177154*^9}, {3.6941870767730503`*^9, 3.69418716300373*^9}, {
+  3.694213209628356*^9, 3.6942132459364033`*^9}}],
+
+Cell[BoxData[
+ RowBox[{
+  RowBox[{"makeEnumCode", "[", "indent_", "]"}], ":=", 
+  RowBox[{"Module", "[", 
+   RowBox[{
+    RowBox[{"{", "out", "}"}], ",", "\[IndentingNewLine]", 
+    RowBox[{
+     RowBox[{"out", "=", 
+      RowBox[{
+      "indent", "<>", "\"\<GRID_SERIALIZABLE_ENUM(Algebra, undef,\>\""}]}], 
+     ";", "\[IndentingNewLine]", 
+     RowBox[{"ind", "=", "0"}], ";", "\[IndentingNewLine]", 
+     RowBox[{"Do", "[", "\[IndentingNewLine]", 
+      RowBox[{
+       RowBox[{
+        RowBox[{"out", "=", 
+         RowBox[{
+         "out", "<>", "\"\<\\n\>\"", "<>", "indent", "<>", " ", 
+          "\"\<                       \>\"", "<>", 
+          RowBox[{"StringPadRight", "[", 
+           RowBox[{
+            RowBox[{"StringTrim", "[", 
+             RowBox[{
+              RowBox[{"enum", "[", "i", "]"}], ",", "\"\<Algebra::\>\""}], 
+             "]"}], ",", "17"}], "]"}], "<>", "\"\<, \>\"", "<>", 
+          RowBox[{"ToString", "[", "ind", "]"}]}]}], ";", 
+        "\[IndentingNewLine]", 
+        RowBox[{"out", " ", "=", " ", 
+         RowBox[{"out", "<>", 
+          RowBox[{"If", "[", 
+           RowBox[{
+            RowBox[{"i", "\[Equal]", 
+             RowBox[{"Last", "[", "alg", "]"}]}], ",", 
+            RowBox[{
+            "\"\<);\\n\>\"", "<>", "indent", "<>", 
+             "\"\<static constexpr unsigned int nGamma = \>\"", "<>", 
+             RowBox[{"ToString", "[", 
+              RowBox[{"ind", "+", "1"}], "]"}], "<>", "\"\<;\\n\>\""}], ",", 
+            "\"\<,\>\"", ",", "\"\<,\>\""}], "]"}]}]}], ";", 
+        "\[IndentingNewLine]", 
+        RowBox[{"ind", "=", 
+         RowBox[{"ind", "+", "1"}]}]}], ",", "\[IndentingNewLine]", 
+       RowBox[{"{", 
+        RowBox[{"i", ",", "alg"}], "}"}]}], "]"}], ";", "\[IndentingNewLine]", 
+     RowBox[{"Return", "[", "out", "]"}]}]}], "]"}]}]], "Input",
+ CellChangeTimes->{{3.694187393082726*^9, 3.694187405147518*^9}, {
+  3.694187542837305*^9, 3.694187612925947*^9}, {3.694187643751446*^9, 
+  3.694187659177355*^9}, {3.69418773039787*^9, 3.694187848150207*^9}, {
+  3.6941879877586403`*^9, 3.6941880928686857`*^9}, {3.694188156471807*^9, 
+  3.694188204997113*^9}, {3.694188290260517*^9, 3.694188291135356*^9}, {
+  3.6941933492591*^9, 3.6941934762397633`*^9}, {3.694206335439877*^9, 
+  3.694206371730096*^9}, {3.694206667692852*^9, 3.6942066758217916`*^9}, {
+  3.694206752732321*^9, 3.694206753090602*^9}, {3.6942071072527027`*^9, 
+  3.694207214318696*^9}, {3.694211442308366*^9, 3.694211490100521*^9}, {
+  3.6942115668101377`*^9, 3.694211571070611*^9}, {3.6942133172135267`*^9, 
+  3.694213325783718*^9}}],
+
+Cell["Multiplication functions generation:", "Text",
+ CellFrame->{{0, 0}, {0, 0.5}},
+ CellChangeTimes->{{3.69420919761381*^9, 3.694209206431526*^9}}],
+
+Cell[BoxData[{
+ RowBox[{
+  RowBox[{"rhsv", "=", 
+   RowBox[{"{", 
+    RowBox[{
+     RowBox[{"rhs", "[", "0", "]"}], ",", 
+     RowBox[{"rhs", "[", "1", "]"}], ",", 
+     RowBox[{"rhs", "[", "2", "]"}], ",", 
+     RowBox[{"rhs", "[", "3", "]"}]}], "}"}]}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"vecCoef", "[", 
+   RowBox[{"v_", "[", "i_", "]"}], "]"}], ":=", 
+  RowBox[{
+   RowBox[{"ToString", "[", "v", "]"}], "<>", "\"\<(\>\"", "<>", 
+   RowBox[{"ToString", "[", "i", "]"}], "<>", 
+   "\"\<)\>\""}]}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"lmatCoef", "[", 
+   RowBox[{"v_", "[", "i_", "]"}], "]"}], ":=", 
+  RowBox[{
+   RowBox[{"ToString", "[", "v", "]"}], "<>", "\"\<(\>\"", "<>", 
+   RowBox[{"ToString", "[", "i", "]"}], "<>", 
+   "\"\<, i)\>\""}]}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"rmatCoef", "[", 
+   RowBox[{"v_", "[", "i_", "]"}], "]"}], ":=", 
+  RowBox[{
+   RowBox[{"ToString", "[", "v", "]"}], "<>", "\"\<(i, \>\"", "<>", 
+   RowBox[{"ToString", "[", "i", "]"}], "<>", 
+   "\"\<)\>\""}]}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"expr", "[", 
+   RowBox[{"ir_", ",", 
+    RowBox[{"rhs", "[", "i_", "]"}], ",", "strFn_"}], "]"}], ":=", 
+  RowBox[{
+   RowBox[{"strFn", "[", 
+    RowBox[{"ret", "[", "ir", "]"}], "]"}], "<>", "\"\< = \>\"", "<>", " ", 
+   RowBox[{"strFn", "[", 
+    RowBox[{"rhs", "[", "i", "]"}], "]"}], "<>", 
+   "\"\<;\>\""}]}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"expr", "[", 
+   RowBox[{"ir_", ",", 
+    RowBox[{"Times", "[", 
+     RowBox[{
+      RowBox[{"-", "1"}], ",", 
+      RowBox[{"rhs", "[", "i_", "]"}]}], "]"}], ",", "strFn_"}], "]"}], ":=", 
+  
+  RowBox[{
+   RowBox[{"strFn", "[", 
+    RowBox[{"ret", "[", "ir", "]"}], "]"}], "<>", "\"\< = -\>\"", "<>", " ", 
+   RowBox[{"strFn", "[", 
+    RowBox[{"rhs", "[", "i", "]"}], "]"}], "<>", 
+   "\"\<;\>\""}]}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"expr", "[", 
+   RowBox[{"ir_", ",", 
+    RowBox[{"Times", "[", 
+     RowBox[{
+      RowBox[{"Complex", " ", "[", 
+       RowBox[{"0", ",", "1"}], "]"}], ",", 
+      RowBox[{"rhs", "[", "i_", "]"}]}], "]"}], ",", "strFn_"}], "]"}], ":=", 
+  
+  RowBox[{
+   RowBox[{"strFn", "[", 
+    RowBox[{"ret", "[", "ir", "]"}], "]"}], "<>", "\"\< = timesI(\>\"", "<>", 
+   " ", 
+   RowBox[{"strFn", "[", 
+    RowBox[{"rhs", "[", "i", "]"}], "]"}], "<>", 
+   "\"\<);\>\""}]}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"expr", "[", 
+   RowBox[{"ir_", ",", 
+    RowBox[{"Times", "[", 
+     RowBox[{
+      RowBox[{"Complex", " ", "[", 
+       RowBox[{"0", ",", 
+        RowBox[{"-", "1"}]}], "]"}], ",", 
+      RowBox[{"rhs", "[", "i_", "]"}]}], "]"}], ",", "strFn_"}], "]"}], ":=", 
+  
+  RowBox[{
+   RowBox[{"strFn", "[", 
+    RowBox[{"ret", "[", "ir", "]"}], "]"}], "<>", "\"\< = timesMinusI(\>\"", "<>",
+    " ", 
+   RowBox[{"strFn", "[", 
+    RowBox[{"rhs", "[", "i", "]"}], "]"}], "<>", 
+   "\"\<);\>\""}]}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{
+   RowBox[{"vecFunc", "[", "g_", "]"}], ":=", 
+   RowBox[{"Module", "[", 
+    RowBox[{
+     RowBox[{"{", 
+      RowBox[{
+       RowBox[{"out", "=", "\"\<\>\""}], ",", 
+       RowBox[{"ir", "=", "0"}]}], "}"}], ",", "\[IndentingNewLine]", 
+     RowBox[{
+      RowBox[{"out", "=", 
+       RowBox[{
+       "\"\<template<class vtype>\ninline void mult\>\"", "<>", 
+        "\[IndentingNewLine]", 
+        RowBox[{"StringTrim", "[", 
+         RowBox[{
+          RowBox[{"enum", "[", "g", "]"}], ",", "\"\<Algebra::\>\""}], "]"}], 
+        "<>", "\[IndentingNewLine]", 
+        "\"\<(iVector<vtype, Ns> &ret, const iVector<vtype, Ns> \
+&rhs)\\n{\\n\>\""}]}], ";", "\[IndentingNewLine]", 
+      RowBox[{"Do", "[", 
+       RowBox[{
+        RowBox[{
+         RowBox[{"out", "=", 
+          RowBox[{"out", "<>", "\"\<  \>\"", "<>", 
+           RowBox[{"expr", "[", 
+            RowBox[{"ir", ",", "c", ",", "vecCoef"}], "]"}], "<>", 
+           "\"\<\\n\>\""}]}], ";", 
+         RowBox[{"ir", "=", 
+          RowBox[{"ir", "+", "1"}]}]}], ",", 
+        RowBox[{"{", 
+         RowBox[{"c", ",", 
+          RowBox[{
+           RowBox[{"mat", "[", "g", "]"}], ".", "rhsv"}]}], "}"}]}], "]"}], 
+      ";", "\[IndentingNewLine]", 
+      RowBox[{"out", " ", "=", 
+       RowBox[{"out", "<>", "\"\<};\\n\>\""}]}], ";", "\[IndentingNewLine]", 
+      RowBox[{"Return", "[", "out", "]"}]}]}], "]"}]}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"lmatFunc", "[", "g_", "]"}], ":=", 
+  RowBox[{"Module", "[", 
+   RowBox[{
+    RowBox[{"{", 
+     RowBox[{
+      RowBox[{"out", "=", "\"\<\>\""}], ",", 
+      RowBox[{"ir", "=", "0"}]}], "}"}], ",", "\[IndentingNewLine]", 
+    RowBox[{
+     RowBox[{"out", "=", 
+      RowBox[{
+      "\"\<template<class vtype>\ninline void lmult\>\"", "<>", 
+       "\[IndentingNewLine]", 
+       RowBox[{"StringTrim", "[", 
+        RowBox[{
+         RowBox[{"enum", "[", "g", "]"}], ",", "\"\<Algebra::\>\""}], "]"}], "<>",
+        "\[IndentingNewLine]", 
+       "\"\<(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)\\n{\\n\>\
+\"", "<>", "\[IndentingNewLine]", 
+       "\"\<  for(unsigned int i = 0; i < Ns; ++i)\\n  {\\n\>\""}]}], ";", 
+     "\[IndentingNewLine]", 
+     RowBox[{"Do", "[", 
+      RowBox[{
+       RowBox[{
+        RowBox[{"out", "=", 
+         RowBox[{"out", "<>", "\"\<    \>\"", "<>", 
+          RowBox[{"expr", "[", 
+           RowBox[{"ir", ",", "c", ",", "lmatCoef"}], "]"}], "<>", 
+          "\"\<\\n\>\""}]}], ";", 
+        RowBox[{"ir", "=", 
+         RowBox[{"ir", "+", "1"}]}]}], ",", 
+       RowBox[{"{", 
+        RowBox[{"c", ",", 
+         RowBox[{
+          RowBox[{"mat", "[", "g", "]"}], ".", "rhsv"}]}], "}"}]}], "]"}], 
+     ";", "\[IndentingNewLine]", 
+     RowBox[{"out", " ", "=", 
+      RowBox[{"out", "<>", "\"\<  }\\n};\\n\>\""}]}], ";", 
+     "\[IndentingNewLine]", 
+     RowBox[{"Return", "[", "out", "]"}]}]}], 
+   "]"}]}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"rmatFunc", "[", "g_", "]"}], ":=", 
+  RowBox[{"Module", "[", 
+   RowBox[{
+    RowBox[{"{", 
+     RowBox[{
+      RowBox[{"out", "=", "\"\<\>\""}], ",", 
+      RowBox[{"ir", "=", "0"}]}], "}"}], ",", "\[IndentingNewLine]", 
+    RowBox[{
+     RowBox[{"out", "=", 
+      RowBox[{
+      "\"\<template<class vtype>\ninline void rmult\>\"", "<>", 
+       "\[IndentingNewLine]", 
+       RowBox[{"StringTrim", "[", 
+        RowBox[{
+         RowBox[{"enum", "[", "g", "]"}], ",", "\"\<Algebra::\>\""}], "]"}], "<>",
+        "\[IndentingNewLine]", 
+       "\"\<(iMatrix<vtype, Ns> &ret, const iMatrix<vtype, Ns> &rhs)\\n{\\n\>\
+\"", "<>", "\[IndentingNewLine]", 
+       "\"\<  for(unsigned int i = 0; i < Ns; ++i)\\n  {\\n\>\""}]}], ";", 
+     "\[IndentingNewLine]", 
+     RowBox[{"Do", "[", 
+      RowBox[{
+       RowBox[{
+        RowBox[{"out", "=", 
+         RowBox[{"out", "<>", "\"\<    \>\"", "<>", 
+          RowBox[{"expr", "[", 
+           RowBox[{"ir", ",", "c", ",", "rmatCoef"}], "]"}], "<>", 
+          "\"\<\\n\>\""}]}], ";", 
+        RowBox[{"ir", "=", 
+         RowBox[{"ir", "+", "1"}]}]}], ",", 
+       RowBox[{"{", 
+        RowBox[{"c", ",", 
+         RowBox[{"rhsv", ".", 
+          RowBox[{"mat", "[", "g", "]"}]}]}], "}"}]}], "]"}], ";", 
+     "\[IndentingNewLine]", 
+     RowBox[{"out", " ", "=", 
+      RowBox[{"out", "<>", "\"\<  }\\n};\\n\>\""}]}], ";", 
+     "\[IndentingNewLine]", 
+     RowBox[{"Return", "[", "out", "]"}]}]}], 
+   "]"}]}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"allFunc", "[", "g_", "]"}], ":=", 
+  RowBox[{"Module", "[", 
+   RowBox[{
+    RowBox[{"{", 
+     RowBox[{"out", "=", "\"\<\>\""}], "}"}], ",", "\[IndentingNewLine]", 
+    RowBox[{
+     RowBox[{"out", "=", 
+      RowBox[{"out", "<>", 
+       RowBox[{"vecFunc", "[", "g", "]"}], "<>", "\"\<\\n\>\""}]}], ";", 
+     "\[IndentingNewLine]", 
+     RowBox[{"out", "=", 
+      RowBox[{"out", "<>", 
+       RowBox[{"lmatFunc", "[", "g", "]"}], "<>", "\"\<\\n\>\""}]}], ";", 
+     "\[IndentingNewLine]", 
+     RowBox[{"out", "=", 
+      RowBox[{"out", "<>", 
+       RowBox[{"rmatFunc", "[", "g", "]"}], "<>", "\"\<\\n\>\""}]}], ";", 
+     "\[IndentingNewLine]", 
+     RowBox[{"Return", "[", "out", "]"}]}]}], 
+   "]"}]}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"vecFinal", "[", "]"}], ":=", 
+  RowBox[{"Module", "[", 
+   RowBox[{
+    RowBox[{"{", 
+     RowBox[{"out", "=", "\"\<\>\""}], "}"}], ",", "\[IndentingNewLine]", 
+    RowBox[{
+     RowBox[{
+     "out", "=", 
+      "\"\<template<class vtype> \ninline auto operator*(const Gamma &G, \
+const iVector<vtype, Ns> &arg)\n->typename \
+std::enable_if<matchGridTensorIndex<iVector<vtype, Ns>, SpinorIndex>::value, \
+iVector<vtype, Ns>>::type\n{\n  iVector<vtype, Ns> ret;\n\n  switch (G.g) \n  \
+{\n\>\""}], ";", "\[IndentingNewLine]", 
+     RowBox[{"Do", "[", 
+      RowBox[{
+       RowBox[{
+        RowBox[{"out", " ", "=", 
+         RowBox[{"out", "<>", "\"\<  case Gamma::\>\"", "<>", 
+          RowBox[{"enum", "[", "g", "]"}], "<>", "\"\<:\\n\>\""}]}], ";", 
+        "\[IndentingNewLine]", 
+        RowBox[{"out", "=", 
+         RowBox[{"out", "<>", "\"\<    mult\>\"", "<>", 
+          RowBox[{"StringTrim", "[", 
+           RowBox[{
+            RowBox[{"enum", "[", "g", "]"}], ",", "\"\<Algebra::\>\""}], 
+           "]"}], "<>", "\"\<(ret, arg); break;\\n\>\""}]}]}], ",", 
+       "\[IndentingNewLine]", 
+       RowBox[{"{", 
+        RowBox[{"g", ",", "alg"}], "}"}]}], "]"}], ";", "\[IndentingNewLine]", 
+     RowBox[{"out", "=", 
+      RowBox[{
+      "out", "<>", 
+       "\"\<  default: assert(0);\n  }\n  \n  return ret;\n}\\n\>\""}]}], ";",
+      "\[IndentingNewLine]", 
+     RowBox[{"Return", "[", "out", "]"}]}]}], 
+   "]"}]}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"lmatFinal", "[", "]"}], ":=", 
+  RowBox[{"Module", "[", 
+   RowBox[{
+    RowBox[{"{", 
+     RowBox[{"out", "=", "\"\<\>\""}], "}"}], ",", "\[IndentingNewLine]", 
+    RowBox[{
+     RowBox[{
+     "out", "=", 
+      "\"\<template<class vtype> \ninline auto operator*(const Gamma &G, \
+const iMatrix<vtype, Ns> &arg)\n->typename \
+std::enable_if<matchGridTensorIndex<iMatrix<vtype, Ns>, SpinorIndex>::value, \
+iMatrix<vtype, Ns>>::type\n{\n  iMatrix<vtype, Ns> ret;\n\n  switch (G.g) \n  \
+{\n\>\""}], ";", "\[IndentingNewLine]", 
+     RowBox[{"Do", "[", 
+      RowBox[{
+       RowBox[{
+        RowBox[{"out", " ", "=", 
+         RowBox[{"out", "<>", "\"\<  case Gamma::\>\"", "<>", 
+          RowBox[{"enum", "[", "g", "]"}], "<>", "\"\<:\\n\>\""}]}], ";", 
+        "\[IndentingNewLine]", 
+        RowBox[{"out", "=", 
+         RowBox[{"out", "<>", "\"\<    lmult\>\"", "<>", 
+          RowBox[{"StringTrim", "[", 
+           RowBox[{
+            RowBox[{"enum", "[", "g", "]"}], ",", "\"\<Algebra::\>\""}], 
+           "]"}], "<>", "\"\<(ret, arg); break;\\n\>\""}]}]}], ",", 
+       "\[IndentingNewLine]", 
+       RowBox[{"{", 
+        RowBox[{"g", ",", "alg"}], "}"}]}], "]"}], ";", "\[IndentingNewLine]", 
+     RowBox[{"out", "=", 
+      RowBox[{
+      "out", "<>", 
+       "\"\<  default: assert(0);\n  }\n  \n  return ret;\n}\\n\>\""}]}], ";",
+      "\[IndentingNewLine]", 
+     RowBox[{"Return", "[", "out", "]"}]}]}], 
+   "]"}]}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"makeFuncCode", "[", "]"}], ":=", 
+  RowBox[{"Module", "[", 
+   RowBox[{
+    RowBox[{"{", 
+     RowBox[{"out", "=", "\"\<\>\""}], "}"}], ",", 
+    RowBox[{
+     RowBox[{"Do", "[", 
+      RowBox[{
+       RowBox[{"out", "=", 
+        RowBox[{"out", "<>", 
+         RowBox[{"allFunc", "[", "g", "]"}]}]}], ",", 
+       RowBox[{"{", 
+        RowBox[{"g", ",", "alg"}], "}"}]}], "]"}], ";", 
+     RowBox[{"Return", "[", "out", "]"}]}]}], 
+   "]"}]}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"rmatFinal", "[", "]"}], ":=", 
+  RowBox[{"Module", "[", 
+   RowBox[{
+    RowBox[{"{", 
+     RowBox[{"out", "=", "\"\<\>\""}], "}"}], ",", "\[IndentingNewLine]", 
+    RowBox[{
+     RowBox[{
+     "out", "=", 
+      "\"\<template<class vtype> \ninline auto operator*(const iMatrix<vtype, \
+Ns> &arg, const Gamma &G)\n->typename \
+std::enable_if<matchGridTensorIndex<iMatrix<vtype, Ns>, SpinorIndex>::value, \
+iMatrix<vtype, Ns>>::type\n{\n  iMatrix<vtype, Ns> ret;\n\n  switch (G.g) \n  \
+{\n\>\""}], ";", "\[IndentingNewLine]", 
+     RowBox[{"Do", "[", 
+      RowBox[{
+       RowBox[{
+        RowBox[{"out", " ", "=", 
+         RowBox[{"out", "<>", "\"\<  case Gamma::\>\"", "<>", 
+          RowBox[{"enum", "[", "g", "]"}], "<>", "\"\<:\\n\>\""}]}], ";", 
+        "\[IndentingNewLine]", 
+        RowBox[{"out", "=", 
+         RowBox[{"out", "<>", "\"\<    rmult\>\"", "<>", 
+          RowBox[{"StringTrim", "[", 
+           RowBox[{
+            RowBox[{"enum", "[", "g", "]"}], ",", "\"\<Algebra::\>\""}], 
+           "]"}], "<>", "\"\<(ret, arg); break;\\n\>\""}]}]}], ",", 
+       "\[IndentingNewLine]", 
+       RowBox[{"{", 
+        RowBox[{"g", ",", "alg"}], "}"}]}], "]"}], ";", "\[IndentingNewLine]", 
+     RowBox[{"out", "=", 
+      RowBox[{
+      "out", "<>", 
+       "\"\<  default: assert(0);\n  }\n  \n  return ret;\n}\\n\>\""}]}], ";",
+      "\[IndentingNewLine]", 
+     RowBox[{"Return", "[", "out", "]"}]}]}], 
+   "]"}]}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"makeFuncCode", "[", "]"}], ":=", 
+  RowBox[{"Module", "[", 
+   RowBox[{
+    RowBox[{"{", 
+     RowBox[{"out", "=", "\"\<\>\""}], "}"}], ",", "\[IndentingNewLine]", 
+    RowBox[{
+     RowBox[{"Do", "[", 
+      RowBox[{
+       RowBox[{"out", "=", 
+        RowBox[{"out", "<>", 
+         RowBox[{"allFunc", "[", "g", "]"}]}]}], ",", 
+       RowBox[{"{", 
+        RowBox[{"g", ",", "alg"}], "}"}]}], "]"}], ";", "\[IndentingNewLine]", 
+     RowBox[{"out", "=", 
+      RowBox[{
+      "out", "<>", 
+       "\"\<\n// FIXME\n//\n// Optimisation; switch over to a \
+\\\"multGammaX(ret._internal,arg._internal)\\\" style early and\n// note that \
+doing so from the lattice operator will avoid copy back and case switch \
+overhead, as\n// was done for the tensor math operator to remove operator * \
+notation early\n//\n\n//left multiply\ntemplate<class vtype> \ninline auto \
+operator*(const Gamma &G, const iScalar<vtype> &arg)\n->typename \
+std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,\
+iScalar<vtype>>::type \n{\n  iScalar<vtype> ret;\n  \
+ret._internal=G*arg._internal;\n  return ret;\n}\n\ntemplate<class vtype,int \
+N>\ninline auto operator*(const Gamma &G, const iVector<vtype, N> &arg)\n\
+->typename \
+std::enable_if<matchGridTensorIndex<iVector<vtype,N>,SpinorIndex>::notvalue,\
+iVector<vtype,N>>::type \n{\n  iVector<vtype,N> ret;\n  for(int i=0;i<N;i++){\
+\n    ret._internal[i]=G*arg._internal[i];\n  }\n  return ret;\n}\n\n\
+template<class vtype, int N>\ninline auto operator*(const Gamma &G, const \
+iMatrix<vtype, N> &arg) \n->typename \
+std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,\
+iMatrix<vtype,N>>::type \n{\n  iMatrix<vtype,N> ret;\n  for(int i=0;i<N;i++){\
+\n  for(int j=0;j<N;j++){\n    ret._internal[i][j]=G*arg._internal[i][j];\n  \
+}}\n  return ret;\n}\n\n//right multiply\ntemplate<class vtype>\ninline auto \
+operator*(const iScalar<vtype> &arg, const Gamma &G)\n->typename \
+std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,\
+iScalar<vtype>>::type \n{\n  iScalar<vtype> ret;\n  \
+ret._internal=arg._internal*G;\n  return ret;\n}\n\ntemplate<class vtype, int \
+N>\ninline auto operator * (const iMatrix<vtype, N> &arg, const Gamma &G) \n\
+->typename \
+std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,\
+iMatrix<vtype,N>>::type \n{\n  iMatrix<vtype,N> ret;\n  for(int i=0;i<N;i++){\
+\n  for(int j=0;j<N;j++){\n    ret._internal[i][j]=arg._internal[i][j]*G;\n  \
+}}\n  return ret;\n}\n\n\>\"", "<>", 
+       RowBox[{"vecFinal", "[", "]"}], "<>", "\"\<\\n\>\"", "<>", 
+       RowBox[{"lmatFinal", "[", "]"}], "<>", "\"\<\\n\>\"", "<>", 
+       RowBox[{"rmatFinal", "[", "]"}]}]}], ";", "\[IndentingNewLine]", 
+     RowBox[{"Return", "[", "out", "]"}]}]}], "]"}]}]}], "Input",
+ CellChangeTimes->{{3.6941903338959*^9, 3.694190334958905*^9}, {
+   3.694190447518065*^9, 3.694190478462919*^9}, {3.69419063973578*^9, 
+   3.694190653669788*^9}, {3.694190716258309*^9, 3.694190763163622*^9}, {
+   3.6941907963765583`*^9, 3.694190805046734*^9}, {3.69419084029659*^9, 
+   3.694190865238916*^9}, {3.694190901904236*^9, 3.694191193971095*^9}, {
+   3.694191319457264*^9, 3.694191355674004*^9}, {3.694191450019418*^9, 
+   3.694191492801243*^9}, {3.694191537908826*^9, 3.6941918755887957`*^9}, {
+   3.694191978131227*^9, 3.694191986666445*^9}, {3.694192020425993*^9, 
+   3.694192167484473*^9}, {3.6941922686760483`*^9, 3.69419237821404*^9}, {
+   3.6941924316822777`*^9, 3.694192445274864*^9}, 3.694193521866951*^9, {
+   3.694206465375606*^9, 3.6942065075728807`*^9}, {3.694206554520149*^9, 
+   3.694206587955097*^9}, {3.694206881090027*^9, 3.6942069067224293`*^9}, {
+   3.694207497315682*^9, 3.6942076418716516`*^9}, {3.694207790223007*^9, 
+   3.694207989637043*^9}, {3.694208028110401*^9, 3.694208057057311*^9}, {
+   3.694208089863509*^9, 3.6942082975753603`*^9}, {3.694212881031795*^9, 
+   3.694212885985661*^9}, {3.694213330227892*^9, 3.694213417940996*^9}, {
+   3.694214921431739*^9, 3.694214951876449*^9}, {3.6942151046483088`*^9, 
+   3.694215118335286*^9}, {3.6942151634191313`*^9, 3.694215188429871*^9}, {
+   3.6942940839999113`*^9, 3.694294090686364*^9}}],
+
+Cell["Header file generation:", "Text",
+ CellFrame->{{0, 0}, {0, 0.5}},
+ CellChangeTimes->{{3.694209229931735*^9, 3.694209240365596*^9}}],
+
+Cell[BoxData[
+ RowBox[{
+  RowBox[{"file", "=", 
+   RowBox[{"OpenWrite", "[", "\"\<Dirac.h\>\"", "]"}]}], ";", 
+  RowBox[{"WriteString", "[", 
+   RowBox[{"file", ",", 
+    RowBox[{"makeHeaderCode", "[", 
+     RowBox[{
+      RowBox[{"makeEnumCode", "[", "\"\<    \>\"", "]"}], ",", 
+      RowBox[{"makeFuncCode", "[", "]"}]}], "]"}]}], "]"}], ";", 
+  RowBox[{"Close", "[", "file", "]"}], ";"}]], "Input",
+ CellChangeTimes->{{3.694207263508134*^9, 3.694207312537463*^9}}]
+}, Open  ]],
+
+Cell[CellGroupData[{
+
+Cell["Dirac.cc source file generation", "Section",
+ CellChangeTimes->{{3.694209250542284*^9, 3.694209266269907*^9}}],
+
+Cell[BoxData[{
+ RowBox[{
+  RowBox[{
+  "out", " ", "=", " ", 
+   "\"\<const std::array<std::array<Gamma::Algebra, Gamma::nGamma>, \
+Gamma::nGamma> \nGamma::mul = {{\\n\>\""}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"Do", "[", "\[IndentingNewLine]", 
+   RowBox[{
+    RowBox[{"Do", "[", "\[IndentingNewLine]", 
+     RowBox[{
+      RowBox[{"out", "=", 
+       RowBox[{"out", "<>", 
+        RowBox[{"If", "[", 
+         RowBox[{
+          RowBox[{"i2", "\[Equal]", 
+           RowBox[{"First", "[", "alg", "]"}]}], ",", "\"\<  {\>\"", ",", 
+          "\"\<   \>\"", ",", "\"\<   \>\""}], "]"}], "\[IndentingNewLine]", 
+        "<>", "\"\<Gamma::\>\"", "<>", 
+        RowBox[{"enum", "[", 
+         RowBox[{"i1", ".", "i2"}], "]"}], "<>", 
+        RowBox[{"If", "[", 
+         RowBox[{
+          RowBox[{"i2", "\[Equal]", 
+           RowBox[{"Last", "[", "alg", "]"}]}], ",", "\"\<}\>\"", ",", 
+          "\"\<\>\"", ",", "\"\<\>\""}], "]"}], "\[IndentingNewLine]", "<>", 
+        RowBox[{"If", "[", 
+         RowBox[{
+          RowBox[{"i2", "\[Equal]", 
+           RowBox[{"Last", "[", "alg", "]"}]}], ",", 
+          RowBox[{"If", "[", 
+           RowBox[{
+            RowBox[{"i1", "\[Equal]", 
+             RowBox[{"Last", "[", "alg", "]"}]}], ",", "\"\<\>\"", ",", 
+            "\"\<\\n  ,\\n\>\"", ",", "\"\<\\n  ,\\n\>\""}], "]"}], ",", 
+          "\"\<,\\n\>\"", ",", "\"\<,\\n\>\""}], "]"}]}]}], ",", 
+      "\[IndentingNewLine]", 
+      RowBox[{"{", 
+       RowBox[{"i2", ",", "alg"}], "}"}]}], "]"}], ",", "\[IndentingNewLine]", 
+    RowBox[{"{", 
+     RowBox[{"i1", ",", "alg"}], "}"}]}], "]"}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"out", "=", 
+   RowBox[{"out", "<>", "\"\<\\n}};\\n\\n\>\""}]}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"out", " ", "=", " ", 
+   RowBox[{
+   "out", "<>", 
+    "\"\<const std::array<Gamma::Algebra, Gamma::nGamma> Gamma::adj = \
+{\\n\>\""}]}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"Do", "[", "\[IndentingNewLine]", 
+   RowBox[{
+    RowBox[{"out", "=", 
+     RowBox[{"out", "<>", "\"\<  Gamma::\>\"", "<>", 
+      RowBox[{"enum", "[", 
+       RowBox[{"i", "\[ConjugateTranspose]"}], "]"}], "<>", 
+      RowBox[{"If", "[", 
+       RowBox[{
+        RowBox[{"i", "\[Equal]", 
+         RowBox[{"Last", "[", "alg", "]"}]}], ",", "\"\<\>\"", ",", 
+        "\"\<,\\n\>\"", ",", "\"\<,\\n\>\""}], "]"}]}]}], ",", 
+    "\[IndentingNewLine]", 
+    RowBox[{"{", 
+     RowBox[{"i", ",", "alg"}], "}"}]}], "]"}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"out", "=", 
+   RowBox[{"out", "<>", "\"\<\\n};\\n\>\""}]}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"file", "=", 
+   RowBox[{"OpenWrite", "[", "\"\<GammaMulTable.h\>\"", "]"}]}], 
+  ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"WriteString", "[", 
+   RowBox[{"file", ",", "out"}], "]"}], ";"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"Close", "[", "file", "]"}], ";"}]}], "Input",
+ CellChangeTimes->{{3.694188457712205*^9, 3.694188459081934*^9}, {
+   3.694188583199066*^9, 3.694188690803248*^9}, {3.6941887378127832`*^9, 
+   3.6941889049946423`*^9}, {3.694188939559833*^9, 3.6941889419840307`*^9}, {
+   3.694189004894498*^9, 3.694189029210237*^9}, 3.694189123631489*^9, {
+   3.694189184854504*^9, 3.694189212219471*^9}, {3.694189278607349*^9, 
+   3.694189315273435*^9}, {3.694189367213715*^9, 3.694189508102169*^9}, {
+   3.694189592594305*^9, 3.69418967222689*^9}, {3.69418976988986*^9, 
+   3.694189793142088*^9}, {3.6941898310221663`*^9, 3.694189833288706*^9}, {
+   3.694189898235135*^9, 3.694189929506019*^9}, 3.6941900145565357`*^9, {
+   3.694190119264579*^9, 3.694190161822473*^9}, {3.694209485052578*^9, 
+   3.6942094865108747`*^9}, {3.694211744604581*^9, 3.6942117644193983`*^9}, {
+   3.694211794987585*^9, 3.694211795663188*^9}, 3.694213436760853*^9, {
+   3.6942134707645893`*^9, 3.69421347411766*^9}, {3.694215381435295*^9, 
+   3.694215385261436*^9}, {3.694297237092742*^9, 3.6942972934131393`*^9}, {
+   3.6942974934555683`*^9, 3.694297566898333*^9}, {3.694297598113819*^9, 
+   3.694297609176114*^9}, {3.694297683998802*^9, 3.694297711493322*^9}}],
+
+Cell[BoxData[
+ RowBox[{
+  RowBox[{
+   RowBox[{"makeSourceCode", "[", "]"}], ":=", 
+   RowBox[{"Module", "[", 
+    RowBox[{
+     RowBox[{"{", 
+      RowBox[{"out", "=", "\"\<\>\""}], "}"}], ",", "\[IndentingNewLine]", 
+     RowBox[{
+      RowBox[{"out", " ", "=", " ", 
+       RowBox[{
+       "out", "<>", "\[IndentingNewLine]", 
+        "\"\<#include <Grid.h>\n\nnamespace Grid {\nnamespace QCD {\n\n\
+#include \\\"GammaMulTable.h\\\"\n\nconst std::array<const char *, \
+Gamma::nGamma> Gamma::name = {{\n\>\""}]}], ";", "\[IndentingNewLine]", 
+      RowBox[{"Do", "[", 
+       RowBox[{
+        RowBox[{"out", " ", "=", " ", 
+         RowBox[{"out", "<>", "\"\<  \\\"\>\"", "<>", "\[IndentingNewLine]", 
+          RowBox[{"StringPadRight", "[", 
+           RowBox[{
+            RowBox[{"StringTrim", "[", 
+             RowBox[{
+              RowBox[{"StringReplace", "[", 
+               RowBox[{
+                RowBox[{"enum", "[", "g", "]"}], ",", 
+                RowBox[{"\"\<Minus\>\"", "\[Rule]", "\"\<-\>\""}]}], "]"}], 
+              ",", "\"\<Algebra::\>\""}], "]"}], ",", "13"}], "]"}], "<>", 
+          "\[IndentingNewLine]", "\"\<\\\"\>\"", "<>", "\[IndentingNewLine]", 
+          
+          RowBox[{"If", "[", 
+           RowBox[{
+            RowBox[{"g", "\[Equal]", 
+             RowBox[{"Last", "[", "alg", "]"}]}], ",", "\"\<}};\\n\>\"", ",", 
+            "\"\<,\\n\>\"", ",", "\"\<,\\n\>\""}], "]"}]}]}], ",", 
+        RowBox[{"{", 
+         RowBox[{"g", ",", "alg"}], "}"}]}], "]"}], ";", 
+      "\[IndentingNewLine]", 
+      RowBox[{"out", " ", "=", 
+       RowBox[{"out", "<>", "\"\<\n}}\n\>\""}]}], ";", "\[IndentingNewLine]", 
+      
+      RowBox[{"Return", "[", "out", "]"}]}]}], "]"}]}], 
+  "\[IndentingNewLine]"}]], "Input",
+ CellChangeTimes->{{3.69421121373662*^9, 3.694211248473432*^9}, {
+   3.6942118524105463`*^9, 3.694211892664316*^9}, {3.69421196260581*^9, 
+   3.69421202215168*^9}, {3.69421216516626*^9, 3.694212368624585*^9}, {
+   3.694212488334324*^9, 3.694212500692437*^9}, 3.6942134310746183`*^9, {
+   3.694215284859365*^9, 3.694215285515815*^9}, {3.694215368213806*^9, 
+   3.69421537607312*^9}, {3.694215508990799*^9, 3.694215509814066*^9}}],
+
+Cell[BoxData[
+ RowBox[{
+  RowBox[{"file", "=", 
+   RowBox[{"OpenWrite", "[", "\"\<Dirac.cc\>\"", "]"}]}], ";", 
+  RowBox[{"WriteString", "[", 
+   RowBox[{"file", ",", 
+    RowBox[{"makeSourceCode", "[", "]"}]}], "]"}], ";", 
+  RowBox[{"Close", "[", "file", "]"}], ";"}]], "Input",
+ CellChangeTimes->{{3.694212294084733*^9, 3.694212304244363*^9}, {
+  3.6942124200138683`*^9, 3.694212434840803*^9}}]
+}, Open  ]]
+},
+WindowSize->{1246, 1005},
+WindowMargins->{{64, Automatic}, {Automatic, 0}},
+FrontEndVersion->"11.0 for Mac OS X x86 (32-bit, 64-bit Kernel) (September \
+21, 2016)",
+StyleDefinitions->"Default.nb"
+]
+(* End of Notebook Content *)
+
+(* Internal cache information *)
+(*CellTagsOutline
+CellTagsIndex->{}
+*)
+(*CellTagsIndex
+CellTagsIndex->{}
+*)
+(*NotebookFileOutline
+Notebook[{
+Cell[558, 20, 1295, 18, 502, "Text"],
+Cell[CellGroupData[{
+Cell[1878, 42, 513, 10, 75, "Input"],
+Cell[2394, 54, 821, 12, 32, "Output"]
+}, Open  ]],
+Cell[CellGroupData[{
+Cell[3252, 71, 114, 1, 64, "Section"],
+Cell[3369, 74, 475, 14, 54, "Input"],
+Cell[3847, 90, 190, 2, 30, "Text"],
+Cell[4040, 94, 5454, 168, 427, "Input"],
+Cell[9497, 264, 119, 1, 30, "Text"],
+Cell[9619, 267, 2115, 63, 96, "Input"],
+Cell[11737, 332, 133, 1, 30, "Text"],
+Cell[CellGroupData[{
+Cell[11895, 337, 479, 15, 32, "Input"],
+Cell[12377, 354, 1398, 30, 96, "Output"]
+}, Open  ]],
+Cell[13790, 387, 110, 1, 30, "Text"],
+Cell[CellGroupData[{
+Cell[13925, 392, 476, 12, 32, "Input"],
+Cell[14404, 406, 6453, 162, 312, "Output"]
+}, Open  ]]
+}, Open  ]],
+Cell[CellGroupData[{
+Cell[20906, 574, 117, 1, 64, "Section"],
+Cell[21026, 577, 128, 2, 38, "Text"],
+Cell[21157, 581, 2742, 49, 894, "Input"],
+Cell[23902, 632, 164, 3, 38, "Text"],
+Cell[24069, 637, 4876, 148, 684, "Input"],
+Cell[28948, 787, 2588, 55, 201, "Input"],
+Cell[31539, 844, 149, 2, 38, "Text"],
+Cell[31691, 848, 17242, 426, 3456, "Input"],
+Cell[48936, 1276, 137, 2, 38, "Text"],
+Cell[49076, 1280, 468, 11, 32, "Input"]
+}, Open  ]],
+Cell[CellGroupData[{
+Cell[49581, 1296, 116, 1, 64, "Section"],
+Cell[49700, 1299, 4128, 95, 390, "Input"],
+Cell[53831, 1396, 2179, 48, 432, "Input"],
+Cell[56013, 1446, 397, 9, 32, "Input"]
+}, Open  ]]
+}
+]
+*)
+
diff --git a/lib/qcd/utils/LinalgUtils.h b/lib/qcd/utils/LinalgUtils.h
index e7e6a794..754cdd2c 100644
--- a/lib/qcd/utils/LinalgUtils.h
+++ b/lib/qcd/utils/LinalgUtils.h
@@ -47,7 +47,7 @@ void axpibg5x(Lattice<vobj> &z,const Lattice<vobj> &x,Coeff a,Coeff b)
 
   GridBase *grid=x._grid;
 
-  Gamma G5(Gamma::Gamma5);
+  Gamma G5(Gamma::Algebra::Gamma5);
 PARALLEL_FOR_LOOP
   for(int ss=0;ss<grid->oSites();ss++){
     vobj tmp;
@@ -80,7 +80,7 @@ void ag5xpby_ssp(Lattice<vobj> &z,Coeff a,const Lattice<vobj> &x,Coeff b,const L
   conformable(x,z);
   GridBase *grid=x._grid;
   int Ls = grid->_rdimensions[0];
-  Gamma G5(Gamma::Gamma5);
+  Gamma G5(Gamma::Algebra::Gamma5);
 PARALLEL_FOR_LOOP
   for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
     vobj tmp;
@@ -98,7 +98,7 @@ void axpbg5y_ssp(Lattice<vobj> &z,Coeff a,const Lattice<vobj> &x,Coeff b,const L
   conformable(x,z);
   GridBase *grid=x._grid;
   int Ls = grid->_rdimensions[0];
-  Gamma G5(Gamma::Gamma5);
+  Gamma G5(Gamma::Algebra::Gamma5);
 PARALLEL_FOR_LOOP
   for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
     vobj tmp;
@@ -116,7 +116,7 @@ void ag5xpbg5y_ssp(Lattice<vobj> &z,Coeff a,const Lattice<vobj> &x,Coeff b,const
   conformable(x,z);
   GridBase *grid=x._grid;
   int Ls = grid->_rdimensions[0];
-  Gamma G5(Gamma::Gamma5);
+  Gamma G5(Gamma::Algebra::Gamma5);
 PARALLEL_FOR_LOOP
   for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
     vobj tmp1;
@@ -168,7 +168,7 @@ void G5R5(Lattice<vobj> &z,const Lattice<vobj> &x)
   z.checkerboard = x.checkerboard;
   conformable(x,z);
   int Ls = grid->_rdimensions[0];
-  Gamma G5(Gamma::Gamma5);
+  Gamma G5(Gamma::Algebra::Gamma5);
 PARALLEL_FOR_LOOP
   for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
     vobj tmp;
diff --git a/lib/serialisation/MacroMagic.h b/lib/serialisation/MacroMagic.h
index 8b027f30..855014e0 100644
--- a/lib/serialisation/MacroMagic.h
+++ b/lib/serialisation/MacroMagic.h
@@ -150,14 +150,14 @@ friend inline bool operator==(const cname &lhs, const cname &rhs) {\
 class name: public Grid::Serializable\
 {\
 public:\
-  enum EnumType\
+  enum\
   {\
     GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMVAL,__VA_ARGS__))\
     undefname = -1\
   };\
 public:\
   name(void): value_(undefname) {};\
-  name(EnumType value): value_(value) {};\
+  name(int value): value_(value) {};\
   template <typename T>\
   static inline void write(Grid::Writer<T> &WR,const std::string &s, const name &obj)\
   {\
@@ -177,7 +177,7 @@ public:\
     GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMTEST,__VA_ARGS__))\
     else {obj = name::undefname;}\
   }\
-  inline operator EnumType(void) const\
+  inline operator int(void) const\
   {\
     return value_;\
   }\
@@ -190,7 +190,7 @@ public:\
     return os;\
   }\
 private:\
-  EnumType value_;\
+  int value_;\
 };
 
 
diff --git a/scripts/filelist b/scripts/filelist
index bf2fbc41..de8351f8 100755
--- a/scripts/filelist
+++ b/scripts/filelist
@@ -4,8 +4,9 @@ home=`pwd`
 
 # library Make.inc
 cd $home/lib
-HFILES=`find . -type f -name '*.h' -not -name '*Hdf5*' -not -path '*/Old/*' -not -path '*/Eigen/*'`
-CCFILES=`find . -type f -name '*.cc' -not -name '*Communicator*.cc' -not -name '*Hdf5*'`
+HFILES=`find . -type f -name '*.h' -not -name '*Hdf5*' -not -path '*/gamma-gen/*' -not -path '*/Old/*' -not -path '*/Eigen/*'`
+HFILES="$HFILES"
+CCFILES=`find . -type f -name '*.cc' -not -path '*/gamma-gen/*' -not -name '*Communicator*.cc' -not -name '*Hdf5*'`
 echo HFILES=$HFILES > Make.inc
 echo >> Make.inc
 echo CCFILES=$CCFILES >> Make.inc
diff --git a/tests/Test_cayley_even_odd_vec.cc b/tests/Test_cayley_even_odd_vec.cc
index f8600782..f39f4ac1 100644
--- a/tests/Test_cayley_even_odd_vec.cc
+++ b/tests/Test_cayley_even_odd_vec.cc
@@ -36,12 +36,12 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
-  };
+    Gamma::Algebra Gmu [] = {
+        Gamma::Algebra::GammaX,
+        Gamma::Algebra::GammaY,
+        Gamma::Algebra::GammaZ,
+        Gamma::Algebra::GammaT
+    };
 
 typedef DomainWallFermion<DomainWallVec5dImplR>                      DomainWallVecFermionR;
 typedef ZMobiusFermion<ZDomainWallVec5dImplR>                        ZMobiusVecFermionR;
@@ -339,7 +339,7 @@ void  TestMoo(This & Dw, That &sDw)
   LatticeFermion ndiff(ngrid);
   LatticeFermion sdiff(sgrid);
 
-  Gamma g5( Gamma::Gamma5 );
+    Gamma g5( Gamma::Algebra::Gamma5 );
 
   std::vector<int> seeds({1,2,3,4,5,7,8});
   GridParallelRNG    RNG5(ngrid);  
diff --git a/tests/Test_dwf_mixedcg_prec.cc b/tests/Test_dwf_mixedcg_prec.cc
index 56394a5e..2601b76c 100644
--- a/tests/Test_dwf_mixedcg_prec.cc
+++ b/tests/Test_dwf_mixedcg_prec.cc
@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)
diff --git a/tests/core/Test_cf_coarsen_support.cc b/tests/core/Test_cf_coarsen_support.cc
index 325cc3a5..fab8c97c 100644
--- a/tests/core/Test_cf_coarsen_support.cc
+++ b/tests/core/Test_cf_coarsen_support.cc
@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)
diff --git a/tests/core/Test_checker.cc b/tests/core/Test_checker.cc
index a974965b..56b0eb60 100644
--- a/tests/core/Test_checker.cc
+++ b/tests/core/Test_checker.cc
@@ -37,11 +37,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 int toint(const char* str){
diff --git a/tests/core/Test_contfrac_even_odd.cc b/tests/core/Test_contfrac_even_odd.cc
index 28f56797..4698f2b9 100644
--- a/tests/core/Test_contfrac_even_odd.cc
+++ b/tests/core/Test_contfrac_even_odd.cc
@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 
diff --git a/tests/core/Test_dwf_even_odd.cc b/tests/core/Test_dwf_even_odd.cc
index 5c60d84b..8fa3f1d1 100644
--- a/tests/core/Test_dwf_even_odd.cc
+++ b/tests/core/Test_dwf_even_odd.cc
@@ -37,11 +37,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 
diff --git a/tests/core/Test_dwf_rb5d.cc b/tests/core/Test_dwf_rb5d.cc
index ab9e4d6e..a3099c7c 100644
--- a/tests/core/Test_dwf_rb5d.cc
+++ b/tests/core/Test_dwf_rb5d.cc
@@ -37,11 +37,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 typedef WilsonFermion5D<DomainWallVec5dImplR> WilsonFermion5DR;
diff --git a/tests/core/Test_fft.cc b/tests/core/Test_fft.cc
index 03c8bc52..877683f0 100644
--- a/tests/core/Test_fft.cc
+++ b/tests/core/Test_fft.cc
@@ -177,7 +177,7 @@ int main (int argc, char ** argv)
     const int sdir=0;
     RealD mass=0.01;
     RealD M5  =1.0;
-    Gamma G5(Gamma::Gamma5);
+    Gamma G5(Gamma::Algebra::Gamma5);
 
     GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,&GRID);
     GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,&GRID);
@@ -218,12 +218,12 @@ int main (int argc, char ** argv)
     /////////////////////////////////////////////////////////////////
     // work out the predicted from Fourier
     /////////////////////////////////////////////////////////////////
-    Gamma::GammaMatrix Gmu [] = {
-      Gamma::GammaX,
-      Gamma::GammaY,
-      Gamma::GammaZ,
-      Gamma::GammaT,
-      Gamma::Gamma5
+    Gamma::Algebra Gmu [] = {
+      Gamma::Algebra::GammaX,
+      Gamma::Algebra::GammaY,
+      Gamma::Algebra::GammaZ,
+      Gamma::Algebra::GammaT,
+      Gamma::Algebra::Gamma5
     };
     LatticeFermionD    Kinetic(FGrid); Kinetic = zero;
     LatticeComplexD    kmu(FGrid); 
@@ -311,7 +311,7 @@ int main (int argc, char ** argv)
     std::cout << " Solving by FFT and Feynman rules" <<std::endl;
     Ddwf.FreePropagator(src,ref,mass) ;
 
-    Gamma G5(Gamma::Gamma5);
+    Gamma G5(Gamma::Algebra::Gamma5);
 
     LatticeFermionD    src5(FGrid); src5=zero;
     LatticeFermionD    tmp5(FGrid); 
@@ -391,7 +391,7 @@ int main (int argc, char ** argv)
     std::cout << " Solving by FFT and Feynman rules" <<std::endl;
     Dov.FreePropagator(src,ref,mass) ;
 
-    Gamma G5(Gamma::Gamma5);
+    Gamma G5(Gamma::Algebra::Gamma5);
 
     LatticeFermionD    src5(FGrid); src5=zero;
     LatticeFermionD    tmp5(FGrid); 
diff --git a/tests/core/Test_gamma.cc b/tests/core/Test_gamma.cc
index 1d3ab28c..2bd6a9b2 100644
--- a/tests/core/Test_gamma.cc
+++ b/tests/core/Test_gamma.cc
@@ -1,223 +1,256 @@
-    /*************************************************************************************
+/*************************************************************************************
 
-    Grid physics library, www.github.com/paboyle/Grid 
+Grid physics library, www.github.com/paboyle/Grid 
 
-    Source file: ./tests/Test_gamma.cc
+Source file: ./tests/Test_gamma.cc
 
-    Copyright (C) 2015
+Copyright (C) 2015-2017
 
 Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Antonin Portelli <antonin.portelli@ed.ac.uk>
 
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
 
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
 
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
-    See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
     /*  END LEGAL */
 #include <Grid/Grid.h>
 
 using namespace std;
 using namespace Grid;
-using namespace Grid::QCD;
+using namespace QCD;
 
-//template<class vobj> class is_pod< iScalar<vobj> >
-//{
-//
-//};
+static constexpr double                      tolerance = 1.0e-6;
+static std::array<SpinMatrix, Gamma::nGamma> testAlgebra;
 
-int main (int argc, char ** argv)
+void print(const SpinMatrix &g)
+{
+  for(int i = 0; i < Ns; i++)
+  {
+    std::cout << GridLogMessage << "(";
+    for(int j=0;j<Ns;j++){
+      if ( abs(g()(i,j)()) == 0 ) {
+        std::cout<< " 0";
+      } else if ( abs(g()(i,j)() - Complex(0,1)) == 0){
+        std::cout<< " i";
+      } else if ( abs(g()(i,j)() + Complex(0,1)) == 0){
+        std::cout<< "-i";
+      } else if ( abs(g()(i,j)() - Complex(1,0)) == 0){
+        std::cout<< " 1";
+      } else if ( abs(g()(i,j)() + Complex(1,0)) == 0){
+        std::cout<< "-1";
+      }
+      std::cout<<((j == Ns-1) ? ")" : "," );
+    }
+    std::cout << std::endl;
+  }
+  std::cout << GridLogMessage << std::endl;
+}
+
+void createTestAlgebra(void)
+{
+  std::array<SpinMatrix, 4> testg;
+  SpinMatrix                testg5;
+  const Complex             I(0., 1.), mI(0., -1.);
+  
+  testg[0] = zero;
+  testg[0]()(0, 3) = I;
+  testg[0]()(1, 2) = I;
+  testg[0]()(2, 1) = mI;
+  testg[0]()(3, 0) = mI;
+  std::cout << GridLogMessage << "test GammaX= " << std::endl;
+  print(testg[0]);
+  testg[1] = zero;
+  testg[1]()(0, 3) = -1.;
+  testg[1]()(1, 2) = 1.;
+  testg[1]()(2, 1) = 1.;
+  testg[1]()(3, 0) = -1.;
+  std::cout << GridLogMessage << "test GammaY= " << std::endl;
+  print(testg[1]);
+  testg[2] = zero;
+  testg[2]()(0, 2) = I;
+  testg[2]()(1, 3) = mI;
+  testg[2]()(2, 0) = mI;
+  testg[2]()(3, 1) = I;
+  std::cout << GridLogMessage << "test GammaZ= " << std::endl;
+  print(testg[2]);
+  testg[3] = zero;
+  testg[3]()(0, 2) = 1.;
+  testg[3]()(1, 3) = 1.;
+  testg[3]()(2, 0) = 1.;
+  testg[3]()(3, 1) = 1.;
+  std::cout << GridLogMessage << "test GammaT= " << std::endl;
+  print(testg[3]);
+  testg5 = testg[0]*testg[1]*testg[2]*testg[3];
+  
+#define DEFINE_TEST_G(g, exp)\
+testAlgebra[Gamma::Algebra::g]        = exp;\
+testAlgebra[Gamma::Algebra::Minus##g] = -exp;\
+  
+  DEFINE_TEST_G(Identity    , 1.);
+  DEFINE_TEST_G(Gamma5      , testg5);
+  DEFINE_TEST_G(GammaX      , testg[0]);
+  DEFINE_TEST_G(GammaY      , testg[1]);
+  DEFINE_TEST_G(GammaZ      , testg[2]);
+  DEFINE_TEST_G(GammaT      , testg[3]);
+  DEFINE_TEST_G(GammaXGamma5, testg[0]*testg5);
+  DEFINE_TEST_G(GammaYGamma5, testg[1]*testg5);
+  DEFINE_TEST_G(GammaZGamma5, testg[2]*testg5);
+  DEFINE_TEST_G(GammaTGamma5, testg[3]*testg5);
+  DEFINE_TEST_G(SigmaXY     , .5*(testg[0]*testg[1] - testg[1]*testg[0]));
+  DEFINE_TEST_G(SigmaXZ     , .5*(testg[0]*testg[2] - testg[2]*testg[0]));
+  DEFINE_TEST_G(SigmaXT     , .5*(testg[0]*testg[3] - testg[3]*testg[0]));
+  DEFINE_TEST_G(SigmaYZ     , .5*(testg[1]*testg[2] - testg[2]*testg[1]));
+  DEFINE_TEST_G(SigmaYT     , .5*(testg[1]*testg[3] - testg[3]*testg[1]));
+  DEFINE_TEST_G(SigmaZT     , .5*(testg[2]*testg[3] - testg[3]*testg[2]));
+  
+#undef DEFINE_TEST_G
+}
+
+template <typename Expr>
+void test(const Expr &a, const Expr &b)
+{
+  if (norm2(a - b) < tolerance)
+  {
+    std::cout << "[OK] ";
+  }
+  else
+  {
+    std::cout << "[fail]" << std::endl;
+    std::cout << GridLogError << "a= " << a << std::endl;
+    std::cout << GridLogError << "is different (tolerance= " << tolerance << ") from " << std::endl;
+    std::cout << GridLogError << "b= " << b << std::endl;
+    exit(EXIT_FAILURE);
+  }
+}
+
+void checkMat(const Gamma::Algebra a, GridSerialRNG &rng)
+{
+  SpinVector v;
+  SpinMatrix m, &testg = testAlgebra[a];
+  Gamma      g(a);
+  bool       pass = true;
+  
+  random(rng, v);
+  random(rng, m);
+  
+  std::cout << GridLogMessage << "Checking " << Gamma::name[a] << ": ";
+  std::cout << "vecmul ";
+  test(g*v, testg*v);
+  std::cout << "matlmul ";
+  test(g*m, testg*m);
+  std::cout << "matrmul ";
+  test(m*g, m*testg);
+  std::cout << std::endl;
+}
+
+void checkProd(const Gamma::Algebra a, const Gamma::Algebra b)
+{
+  SpinMatrix gm, testg = testAlgebra[a]*testAlgebra[b];
+  Gamma      g = Gamma(a)*Gamma(b);
+  bool       pass = true;
+  
+  std::cout << GridLogMessage << "Checking " << Gamma::name[a] << " * "
+            << Gamma::name[b] << ": ";
+  gm = 1.0;
+  gm = g*gm;
+  test(gm, testg);
+  std::cout << "(= " << Gamma::name[g.g] << ")" << std::endl;
+}
+
+void checkAdj(const Gamma::Algebra a)
+{
+  SpinMatrix gm, testg = adj(testAlgebra[a]);
+  Gamma      g(adj(Gamma(a)));
+  bool       pass = true;
+  
+  std::cout << GridLogMessage << "Checking adj(" << Gamma::name[a] << "): ";
+  gm = 1.0;
+  gm = g*gm;
+  test(gm, testg);
+  std::cout << "(= " << Gamma::name[g.g] << ")" << std::endl;
+}
+
+void checkProject(GridSerialRNG &rng)
+{
+  SpinVector     rv, recon, full;
+  HalfSpinVector hsp, hsm;
+  
+  random(rng, rv);
+  
+#define CHECK_PROJ(dir, gamma)\
+std::cout << GridLogMessage << "Checking " << #dir << " projector: ";\
+spProj##dir(hsm,rv);\
+spRecon##dir(recon,hsm);\
+test(recon, rv + Gamma(Gamma::Algebra::gamma)*rv);\
+std::cout << std::endl;
+  
+  CHECK_PROJ(Xp, GammaX);
+  CHECK_PROJ(Yp, GammaY);
+  CHECK_PROJ(Zp, GammaZ);
+  CHECK_PROJ(Tp, GammaT);
+  CHECK_PROJ(5p, Gamma5);
+  CHECK_PROJ(Xm, MinusGammaX);
+  CHECK_PROJ(Ym, MinusGammaY);
+  CHECK_PROJ(Zm, MinusGammaZ);
+  CHECK_PROJ(Tm, MinusGammaT);
+  CHECK_PROJ(5m, MinusGamma5);
+  
+#undef CHECK_PROJ
+}
+
+int main(int argc, char *argv[])
 {
   Grid_init(&argc,&argv);
-
+  
   std::vector<int> latt_size   = GridDefaultLatt();
   std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
   std::vector<int> mpi_layout  = GridDefaultMpi();
-     
-  GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
-
-  GridParallelRNG          pRNG(&Grid);
-  pRNG.SeedRandomDevice();
-
-  GridSerialRNG            sRNG;
+  
+  GridCartesian Grid(latt_size,simd_layout,mpi_layout);
+  GridSerialRNG sRNG;
+  
   sRNG.SeedRandomDevice();
-
-  SpinMatrix ident; ident=zero;
-  SpinMatrix rnd  ; random(sRNG,rnd);
-
-  SpinMatrix ll; ll=zero;
-  SpinMatrix rr; rr=zero;
-  SpinMatrix result;
-
-  SpinVector lv; random(sRNG,lv);
-  SpinVector rv; random(sRNG,rv);
-
-  //  std::cout<<GridLogMessage << " Is pod " << std::is_pod<SpinVector>::value  << std::endl;
-  //  std::cout<<GridLogMessage << " Is pod double   " << std::is_pod<double>::value  << std::endl;
-  //  std::cout<<GridLogMessage << " Is pod ComplexF " << std::is_pod<ComplexF>::value  << std::endl;
-  //  std::cout<<GridLogMessage << " Is triv double " << std::has_trivial_default_constructor<double>::value  << std::endl;
-  //  std::cout<<GridLogMessage << " Is triv ComplexF " << std::has_trivial_default_constructor<ComplexF>::value  << std::endl;
-  //  std::cout<<GridLogMessage << " Is pod Scalar<double> " << std::is_pod<iScalar<double> >::value  << std::endl;
-  //  std::cout<<GridLogMessage << " Is pod Scalar<ComplexF> " << std::is_pod<iScalar<ComplexF> >::value  << std::endl;
-  //  std::cout<<GridLogMessage << " Is pod Scalar<vComplexF> " << std::is_pod<iScalar<vComplexF> >::value  << std::endl;
-  //  std::cout<<GridLogMessage << " Is pod Scalar<vComplexD> " << std::is_pod<iScalar<vComplexD> >::value  << std::endl;
-  //  std::cout<<GridLogMessage << " Is pod Scalar<vRealF> " << std::is_pod<iScalar<vRealF> >::value  << std::endl;
-  //  std::cout<<GridLogMessage << " Is pod Scalar<vRealD> " << std::is_pod<iScalar<vRealD> >::value  << std::endl;
-  //  std::cout<<GridLogMessage << " Is triv Scalar<double> " <<std::has_trivial_default_constructor<iScalar<double> >::value << std::endl;
-  //  std::cout<<GridLogMessage << " Is triv Scalar<vComplexD> "<<std::has_trivial_default_constructor<iScalar<vComplexD> >::value  << std::endl;
-
-  for(int a=0;a<Ns;a++){
-    ident()(a,a) = ComplexF(1.0);
+  
+  std::cout << GridLogMessage << "======== Test algebra" << std::endl;
+  createTestAlgebra();
+  std::cout << GridLogMessage << "======== Multiplication operators check" << std::endl;
+  for (int i = 0; i < Gamma::nGamma; ++i)
+  {
+    checkMat(i, sRNG);
   }
-
-  const Gamma::GammaMatrix *g = Gamma::GammaMatrices;
-  const char **list           = Gamma::GammaMatrixNames;
-
-  result =ll*Gamma(g[0])*rr;
-  result =ll*Gamma(g[0]);
-  rv = Gamma(g[0])*lv;
-
-  for(int mu=0;mu<12;mu++){
-
-    result = Gamma(g[mu])* ident;
-
-    for(int i=0;i<Ns;i++){
-
-      if(i==0) std::cout<<GridLogMessage << list[mu];
-      else     std::cout<<GridLogMessage << list[12];
-
-      std::cout<<"(";
-      for(int j=0;j<Ns;j++){
-	if ( abs(result()(i,j)())==0 ) {
-	  std::cout<< " 0";
-	} else if ( abs(result()(i,j)() - Complex(0,1))==0){
-	  std::cout<< " i";
-	} else if ( abs(result()(i,j)() + Complex(0,1))==0){
-	  std::cout<< "-i";
-	} else if ( abs(result()(i,j)() - Complex(1,0))==0){
-	  std::cout<< " 1";
-	} else if ( abs(result()(i,j)() + Complex(1,0))==0){
-	  std::cout<< "-1";
-	}
-	std::cout<<((j==Ns-1) ? ")" : "," );
-      }
-      std::cout << std::endl;
-    }
-
-    std::cout << std::endl;
-
+  std::cout << GridLogMessage << std::endl;
+  std::cout << GridLogMessage << "======== Algebra multiplication table check" << std::endl;
+  for (int i = 0; i < Gamma::nGamma; ++i)
+  for (int j = 0; j < Gamma::nGamma; ++j)
+  {
+    checkProd(i, j);
   }
-
-  std::cout << "Testing Gamma^2 - 1 = 0"<<std::endl;
-  for(int mu=0;mu<6;mu++){
-    result =  Gamma(g[mu])* ident * Gamma(g[mu]);
-    result = result - ident;
-    RealD mag = norm2(result);
-    std::cout << list[mu]<<" " << mag<<std::endl;
+  std::cout << GridLogMessage << std::endl;
+  std::cout << GridLogMessage << "======== Adjoints check" << std::endl;
+  for (int i = 0; i < Gamma::nGamma; ++i)
+  {
+    checkAdj(i);
   }
-
-  std::cout << "Testing (MinusGamma + G )M = 0"<<std::endl;
-  for(int mu=0;mu<6;mu++){
-    result =          rnd * Gamma(g[mu]);
-    result = result + rnd * Gamma(g[mu+6]);
-    RealD mag = norm2(result);
-    std::cout << list[mu]<<" " << mag<<std::endl;
-  }
-
-  std::cout << "Testing M(MinusGamma + G )  = 0"<<std::endl;
-  for(int mu=0;mu<6;mu++){
-    result =           Gamma(g[mu])  *rnd;
-    result = result +  Gamma(g[mu+6])*rnd;
-    RealD mag = norm2(result);
-    std::cout << list[mu]<<" " << mag<<std::endl;
-  }
-
-  // Testing spins and reconstructs
-  SpinVector     recon; random(sRNG,rv);
-  SpinVector      full;
-  HalfSpinVector hsp,hsm;
-
-  // Xp
-  double mag;
-  spProjXp(hsm,rv);
-  spReconXp(recon,hsm);
-  full = rv + Gamma(Gamma::GammaX) *rv;
-  mag = TensorRemove(norm2(full-recon));
-  std::cout << "Xp "<< mag<<std::endl;
-
-  // Xm
-  spProjXm(hsm,rv);
-  spReconXm(recon,hsm);
-  full = rv - Gamma(Gamma::GammaX) *rv;
-  mag = TensorRemove(norm2(full-recon));
-  std::cout << "Xm "<< mag<<std::endl;
-
-  // Yp
-  spProjYp(hsm,rv);
-  spReconYp(recon,hsm);
-  full = rv + Gamma(Gamma::GammaY) *rv;
-  mag = TensorRemove(norm2(full-recon));
-  std::cout << "Yp "<< mag<<std::endl;
-
-  // Ym
-  spProjYm(hsm,rv);
-  spReconYm(recon,hsm);
-  full = rv - Gamma(Gamma::GammaY) *rv;
-  mag = TensorRemove(norm2(full-recon));
-  std::cout << "Ym "<< mag<<std::endl;
-
-  // Zp
-  spProjZp(hsm,rv);
-  spReconZp(recon,hsm);
-  full = rv + Gamma(Gamma::GammaZ) *rv;
-  mag = TensorRemove(norm2(full-recon));
-  std::cout << "Zp "<< mag<<std::endl;
-
-  // Zm
-  spProjZm(hsm,rv);
-  spReconZm(recon,hsm);
-  full = rv - Gamma(Gamma::GammaZ) *rv;
-  mag = TensorRemove(norm2(full-recon));
-  std::cout << "Zm "<< mag<<std::endl;
-
-  // Tp
-  spProjTp(hsm,rv);
-  spReconTp(recon,hsm);
-  full = rv + Gamma(Gamma::GammaT) *rv;
-  mag = TensorRemove(norm2(full-recon));
-  std::cout << "Tp "<< mag<<std::endl;
-
-  // Tm
-  spProjTm(hsm,rv);
-  spReconTm(recon,hsm);
-  full = rv - Gamma(Gamma::GammaT) *rv;
-  mag = TensorRemove(norm2(full-recon));
-  std::cout << "Tm "<< mag<<std::endl;
-
-  // 5p
-  spProj5p(hsm,rv);
-  spRecon5p(recon,hsm);
-  full = rv + Gamma(Gamma::Gamma5) *rv;
-  mag = TensorRemove(norm2(full-recon));
-  std::cout << "5p "<< mag<<std::endl;
-
-  // 5m
-  spProj5m(hsm,rv);
-  spRecon5m(recon,hsm);
-  full = rv - Gamma(Gamma::Gamma5) *rv;
-  mag = TensorRemove(norm2(full-recon));
-  std::cout << "5m "<< mag<<std::endl;
-
+  std::cout << GridLogMessage << std::endl;
+  std::cout << GridLogMessage << "======== Spin projectors check" << std::endl;
+  checkProject(sRNG);
+  std::cout << GridLogMessage << std::endl;
+  
   Grid_finalize();
+  
+  return EXIT_SUCCESS;
 }
diff --git a/tests/core/Test_wilson_even_odd.cc b/tests/core/Test_wilson_even_odd.cc
index 8028e6de..b5b25058 100644
--- a/tests/core/Test_wilson_even_odd.cc
+++ b/tests/core/Test_wilson_even_odd.cc
@@ -37,11 +37,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)
diff --git a/tests/core/Test_wilson_tm_even_odd.cc b/tests/core/Test_wilson_tm_even_odd.cc
index 2d28240c..2e66438e 100644
--- a/tests/core/Test_wilson_tm_even_odd.cc
+++ b/tests/core/Test_wilson_tm_even_odd.cc
@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)
diff --git a/tests/debug/Test_cayley_cg.cc b/tests/debug/Test_cayley_cg.cc
index 8856ecdf..e0e89ea1 100644
--- a/tests/debug/Test_cayley_cg.cc
+++ b/tests/debug/Test_cayley_cg.cc
@@ -37,11 +37,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 template<class What> 
diff --git a/tests/debug/Test_cayley_coarsen_support.cc b/tests/debug/Test_cayley_coarsen_support.cc
index 9c740b8b..2d503df8 100644
--- a/tests/debug/Test_cayley_coarsen_support.cc
+++ b/tests/debug/Test_cayley_coarsen_support.cc
@@ -38,11 +38,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)
diff --git a/tests/debug/Test_cayley_even_odd.cc b/tests/debug/Test_cayley_even_odd.cc
index 5c15450c..7d8d2a12 100644
--- a/tests/debug/Test_cayley_even_odd.cc
+++ b/tests/debug/Test_cayley_even_odd.cc
@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 
diff --git a/tests/solver/Test_cf_cr_unprec.cc b/tests/solver/Test_cf_cr_unprec.cc
index 59031ba8..d7d8f51a 100644
--- a/tests/solver/Test_cf_cr_unprec.cc
+++ b/tests/solver/Test_cf_cr_unprec.cc
@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)
diff --git a/tests/solver/Test_contfrac_cg.cc b/tests/solver/Test_contfrac_cg.cc
index 5e495071..13f80f69 100644
--- a/tests/solver/Test_contfrac_cg.cc
+++ b/tests/solver/Test_contfrac_cg.cc
@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 
diff --git a/tests/solver/Test_dwf_cg_prec.cc b/tests/solver/Test_dwf_cg_prec.cc
index d248c560..30436e36 100644
--- a/tests/solver/Test_dwf_cg_prec.cc
+++ b/tests/solver/Test_dwf_cg_prec.cc
@@ -37,8 +37,8 @@ struct scal {
   d internal;
 };
 
-Gamma::GammaMatrix Gmu[] = {Gamma::GammaX, Gamma::GammaY, Gamma::GammaZ,
-                            Gamma::GammaT};
+Gamma::Algebra Gmu[] = {Gamma::Algebra::GammaX, Gamma::Algebra::GammaY, Gamma::Algebra::GammaZ,
+                            Gamma::Algebra::GammaT};
 
 int main(int argc, char** argv) {
   Grid_init(&argc, &argv);
diff --git a/tests/solver/Test_dwf_cg_schur.cc b/tests/solver/Test_dwf_cg_schur.cc
index 75ca3521..62a95235 100644
--- a/tests/solver/Test_dwf_cg_schur.cc
+++ b/tests/solver/Test_dwf_cg_schur.cc
@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)
diff --git a/tests/solver/Test_dwf_cg_unprec.cc b/tests/solver/Test_dwf_cg_unprec.cc
index 41a05042..89e9088e 100644
--- a/tests/solver/Test_dwf_cg_unprec.cc
+++ b/tests/solver/Test_dwf_cg_unprec.cc
@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)
diff --git a/tests/solver/Test_dwf_cr_unprec.cc b/tests/solver/Test_dwf_cr_unprec.cc
index 9ee38d35..7549137a 100644
--- a/tests/solver/Test_dwf_cr_unprec.cc
+++ b/tests/solver/Test_dwf_cr_unprec.cc
@@ -37,11 +37,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)
diff --git a/tests/solver/Test_dwf_fpgcr.cc b/tests/solver/Test_dwf_fpgcr.cc
index 3b0f99d0..aaf622b9 100644
--- a/tests/solver/Test_dwf_fpgcr.cc
+++ b/tests/solver/Test_dwf_fpgcr.cc
@@ -38,11 +38,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)
diff --git a/tests/solver/Test_dwf_hdcr.cc b/tests/solver/Test_dwf_hdcr.cc
index 313b2077..64ca0b33 100644
--- a/tests/solver/Test_dwf_hdcr.cc
+++ b/tests/solver/Test_dwf_hdcr.cc
@@ -504,7 +504,7 @@ int main (int argc, char ** argv)
   GridParallelRNG          RNG4(UGrid);   RNG4.SeedFixedIntegers(seeds4);
   GridParallelRNG          CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
 
-  Gamma g5(Gamma::Gamma5);
+  Gamma g5(Gamma::Algebra::Gamma5);
 
   LatticeFermion    src(FGrid); gaussian(RNG5,src);// src=src+g5*src;
   LatticeFermion result(FGrid); result=zero;
diff --git a/tests/solver/Test_wilson_cg_prec.cc b/tests/solver/Test_wilson_cg_prec.cc
index 7cc9d574..e336161f 100644
--- a/tests/solver/Test_wilson_cg_prec.cc
+++ b/tests/solver/Test_wilson_cg_prec.cc
@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)
diff --git a/tests/solver/Test_wilson_cg_schur.cc b/tests/solver/Test_wilson_cg_schur.cc
index 1ea6a07c..8b62b040 100644
--- a/tests/solver/Test_wilson_cg_schur.cc
+++ b/tests/solver/Test_wilson_cg_schur.cc
@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)
diff --git a/tests/solver/Test_wilson_cg_unprec.cc b/tests/solver/Test_wilson_cg_unprec.cc
index 34b0a687..bb429468 100644
--- a/tests/solver/Test_wilson_cg_unprec.cc
+++ b/tests/solver/Test_wilson_cg_unprec.cc
@@ -37,11 +37,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)
diff --git a/tests/solver/Test_wilson_cr_unprec.cc b/tests/solver/Test_wilson_cr_unprec.cc
index 4a517461..add5c2ba 100644
--- a/tests/solver/Test_wilson_cr_unprec.cc
+++ b/tests/solver/Test_wilson_cr_unprec.cc
@@ -36,11 +36,11 @@ struct scal {
   d internal;
 };
 
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
   };
 
 int main (int argc, char ** argv)