diff --git a/bootstrap.sh b/bootstrap.sh
index 7075888c..505fb2b8 100755
--- a/bootstrap.sh
+++ b/bootstrap.sh
@@ -1,4 +1,4 @@
-]#!/usr/bin/env bash
+#!/usr/bin/env bash
 
 EIGEN_URL='http://bitbucket.org/eigen/eigen/get/3.3.3.tar.bz2'
 
diff --git a/lib/lattice/Lattice_unary.h b/lib/lattice/Lattice_unary.h
index 0afeaa32..44b7b4f1 100644
--- a/lib/lattice/Lattice_unary.h
+++ b/lib/lattice/Lattice_unary.h
@@ -62,17 +62,12 @@ namespace Grid {
     return ret;
   }
 
-  template<class obj> Lattice<obj> expMat(const Lattice<obj> &rhs, ComplexD alpha, Integer Nexp = DEFAULT_MAT_EXP){
+  template<class obj> Lattice<obj> expMat(const Lattice<obj> &rhs, RealD alpha, Integer Nexp = DEFAULT_MAT_EXP){
     Lattice<obj> ret(rhs._grid);
     ret.checkerboard = rhs.checkerboard;
     conformable(ret,rhs);
-    obj unit(1.0);
     parallel_for(int ss=0;ss<rhs._grid->oSites();ss++){
-      //ret._odata[ss]=Exponentiate(rhs._odata[ss],alpha, Nexp);
-      ret._odata[ss] = unit;
-      for(int i=Nexp; i>=1;--i)
-	ret._odata[ss] = unit + ret._odata[ss]*rhs._odata[ss]*(alpha/RealD(i));
-
+      ret._odata[ss]=Exponentiate(rhs._odata[ss],alpha, Nexp);
     }
 
     return ret;
diff --git a/lib/qcd/action/gauge/GaugeImplTypes.h b/lib/qcd/action/gauge/GaugeImplTypes.h
index 9d36eead..29e79581 100644
--- a/lib/qcd/action/gauge/GaugeImplTypes.h
+++ b/lib/qcd/action/gauge/GaugeImplTypes.h
@@ -59,7 +59,7 @@ public:
   typedef iImplGaugeLink<Simd>  SiteLink;
   typedef iImplGaugeField<Simd> SiteField;
 
-  typedef Lattice<SiteLink>  LinkField; 
+  typedef Lattice<SiteLink>  LinkField;
   typedef Lattice<SiteField> Field;
 
   // Guido: we can probably separate the types from the HMC functions
@@ -80,7 +80,7 @@ public:
 
   ///////////////////////////////////////////////////////////
   // Move these to another class
-  // HMC auxiliary functions 
+  // HMC auxiliary functions
   static inline void generate_momenta(Field &P, GridParallelRNG &pRNG) {
     // specific for SU gauge fields
     LinkField Pmu(P._grid);
@@ -92,14 +92,19 @@ public:
   }
 
   static inline Field projectForce(Field &P) { return Ta(P); }
-  
+
   static inline void update_field(Field& P, Field& U, double ep){
-    for (int mu = 0; mu < Nd; mu++) {
-      auto Umu = PeekIndex<LorentzIndex>(U, mu);
-      auto Pmu = PeekIndex<LorentzIndex>(P, mu);
-      Umu = expMat(Pmu, ep, Nexp) * Umu;
-      PokeIndex<LorentzIndex>(U, ProjectOnGroup(Umu), mu);
+    //static std::chrono::duration<double> diff;
+
+    //auto start = std::chrono::high_resolution_clock::now();
+    parallel_for(int ss=0;ss<P._grid->oSites();ss++){
+      for (int mu = 0; mu < Nd; mu++) 
+        U[ss]._internal[mu] = ProjectOnGroup(Exponentiate(P[ss]._internal[mu], ep, Nexp) * U[ss]._internal[mu]);
     }
+    
+    //auto end = std::chrono::high_resolution_clock::now();
+   // diff += end - start;
+   // std::cout << "Time to exponentiate matrix " << diff.count() << " s\n";
   }
 
   static inline RealD FieldSquareNorm(Field& U){
diff --git a/lib/qcd/smearing/StoutSmearing.h b/lib/qcd/smearing/StoutSmearing.h
index b50ffe21..bfc37d0d 100644
--- a/lib/qcd/smearing/StoutSmearing.h
+++ b/lib/qcd/smearing/StoutSmearing.h
@@ -58,6 +58,8 @@ class Smear_Stout : public Smear<Gimpl> {
     SmearBase->smear(C, U);
   };
 
+
+  // Repetion of code here (use the Tensor_exp.h function)
   void exponentiate_iQ(GaugeLinkField& e_iQ, const GaugeLinkField& iQ) const {
     // Put this outside
     // only valid for SU(3) matrices
diff --git a/lib/tensors/Tensor_exp.h b/lib/tensors/Tensor_exp.h
index ad0d2071..e18fed70 100644
--- a/lib/tensors/Tensor_exp.h
+++ b/lib/tensors/Tensor_exp.h
@@ -37,30 +37,105 @@ namespace Grid {
   /////////////////////////////////////////////// 
 
 
-  template<class vtype> inline iScalar<vtype> Exponentiate(const iScalar<vtype>&r, ComplexD alpha ,  Integer Nexp = DEFAULT_MAT_EXP)
+  template<class vtype> inline iScalar<vtype> Exponentiate(const iScalar<vtype>&r, RealD alpha ,  Integer Nexp = DEFAULT_MAT_EXP)
     {
       iScalar<vtype> ret;
       ret._internal = Exponentiate(r._internal, alpha, Nexp);
       return ret;
     }
 
-
-  template<class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr> 
-    inline iMatrix<vtype,N> Exponentiate(const iMatrix<vtype,N> &arg, ComplexD alpha  , Integer Nexp = DEFAULT_MAT_EXP )
+template<class vtype, int N> inline iVector<vtype, N> Exponentiate(const iVector<vtype,N>&r, RealD alpha ,  Integer Nexp = DEFAULT_MAT_EXP)
     {
-      iMatrix<vtype,N> unit(1.0);
-      iMatrix<vtype,N> temp(unit);
-      
-      for(int i=Nexp; i>=1;--i){
-	temp *= alpha/ComplexD(i);
-	temp = unit + temp*arg;
-      }
-      
-      return temp;
-      
+      iVector<vtype, N> ret;
+      for (int i = 0; i < N; i++)
+        ret._internal[i] = Exponentiate(r._internal[i], alpha, Nexp);
+      return ret;
     }
 
 
 
+    // Specialisation: Cayley-Hamilton exponential for SU(3)
+    template<class vtype, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0>::type * =nullptr> 
+    inline iMatrix<vtype,3> Exponentiate(const iMatrix<vtype,3> &arg, RealD alpha  , Integer Nexp = DEFAULT_MAT_EXP )
+    {
+    // for SU(3) 2x faster than the std implementation using Nexp=12
+    // notice that it actually computes
+    // exp ( input matrix )
+    // the i sign is coming from outside
+    // input matrix is anti-hermitian NOT hermitian
+      typedef iMatrix<vtype,3> mat;
+      typedef iScalar<vtype> scalar;
+      mat unit(1.0);
+      mat temp(unit);
+      const Complex one_over_three = 1.0 / 3.0;
+      const Complex one_over_two = 1.0 / 2.0;
+
+      scalar c0, c1, tmp, c0max, theta, u, w;
+      scalar xi0, u2, w2, cosw;
+      scalar fden, h0, h1, h2;
+      scalar e2iu, emiu, ixi0, qt;
+      scalar f0, f1, f2;
+      scalar unity(1.0);
+      
+      mat iQ2 = arg*arg*alpha*alpha;
+      mat iQ3 = arg*iQ2*alpha;   
+      // sign in c0 from the conventions on the Ta
+      c0 = -imag( trace(iQ3) ) * one_over_three;  
+      c1 = -real( trace(iQ2) ) * one_over_two;
+
+      // Cayley Hamilton checks to machine precision, tested
+      tmp = c1 * one_over_three;
+      c0max = 2.0 * pow(tmp, 1.5);
+
+      theta = acos(c0 / c0max) * one_over_three;
+      u = sqrt(tmp) * cos(theta);
+      w = sqrt(c1) * sin(theta);
+
+      xi0 = sin(w) / w;
+      u2 = u * u;
+      w2 = w * w;
+      cosw = cos(w);
+
+      ixi0 = timesI(xi0);
+      emiu = cos(u) - timesI(sin(u));
+      e2iu = cos(2.0 * u) + timesI(sin(2.0 * u));
+
+      h0 = e2iu * (u2 - w2) +
+           emiu * ((8.0 * u2 * cosw) + (2.0 * u * (3.0 * u2 + w2) * ixi0));
+      h1 = e2iu * (2.0 * u) - emiu * ((2.0 * u * cosw) - (3.0 * u2 - w2) * ixi0);
+      h2 = e2iu - emiu * (cosw + (3.0 * u) * ixi0);
+
+      fden = unity / (9.0 * u2 - w2);  // reals
+      f0 = h0 * fden;
+      f1 = h1 * fden;
+      f2 = h2 * fden;
+
+      return (f0 * unit + timesMinusI(f1) * arg*alpha - f2 * iQ2);
+    }
+
+
+
+// General exponential
+template<class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr> 
+    inline iMatrix<vtype,N> Exponentiate(const iMatrix<vtype,N> &arg, RealD alpha  , Integer Nexp = DEFAULT_MAT_EXP )
+    {
+    // notice that it actually computes
+    // exp ( input matrix )
+    // the i sign is coming from outside
+    // input matrix is anti-hermitian NOT hermitian
+      typedef iMatrix<vtype,N> mat;
+      mat unit(1.0);
+      mat temp(unit);
+      for(int i=Nexp; i>=1;--i){
+	      temp *= alpha/RealD(i);
+	      temp = unit + temp*arg;
+      }
+      return temp;
+
+    }
+
+
+
+
 }
 #endif