diff --git a/dune/elasticity/assemblers/localhyperdualstiffness.hh b/dune/elasticity/assemblers/localhyperdualstiffness.hh
index a91b68ef0d31a6c0bd7da8507a20c301d9d528e0..e3918385e15e2cd8d53cb8df61944011a2e472a4 100644
--- a/dune/elasticity/assemblers/localhyperdualstiffness.hh
+++ b/dune/elasticity/assemblers/localhyperdualstiffness.hh
@@ -56,7 +56,7 @@ energy(const LocalView& localView,
         
     energy = localEnergy_->energy(localView,localHyperDualConfiguration);
 
-    pureEnergy = energy.real();
+    pureEnergy = energy.value();
 
     return pureEnergy;
 }
@@ -78,30 +78,21 @@ assembleGradientAndHessian(const LocalView& localView,
     // Create hyperdual vector localHyperDualConfiguration = double vector localConfiguration
     std::vector<hyperdual> localHyperDualConfiguration(nDoubles);
     for (size_t i=0; i<nDoubles; i++)
-       localHyperDualConfiguration[i] = hyperdual(localConfiguration[i]);
-    
-    // Compute gradient and hessian (symmetry is used) using hyperdual function evaluation
-    // localHyperDualConfiguration puts Ones in the eps1 and eps2 component to evaluate the different partial derivatives
-    for (size_t i=0; i<nDoubles; i++) 
+       localHyperDualConfiguration[i] = hyperdual(localConfiguration[i],i);
+
+    hyperdual temp = localEnergy_->energy(localView, localHyperDualConfiguration);
+
+
+    for (size_t i=0; i<nDoubles; i++)
     {
-      localHyperDualConfiguration[i].seteps1(1.0);
-      localHyperDualConfiguration[i].seteps2(1.0);
-    
-      hyperdual temp = localEnergy_->energy(localView, localHyperDualConfiguration);
-      localGradient[i] = temp.eps1();                   
-      localHessian[i][i] = temp.eps1eps2();
-      
-      localHyperDualConfiguration[i].seteps2(0.0);
-
-      for (size_t j=i+1; j<nDoubles; j++) 
-      { 
-        localHyperDualConfiguration[j].seteps2(1.0);
-        localHessian[i][j] = localEnergy_->energy(localView, localHyperDualConfiguration).eps1eps2();
-        localHessian[j][i] = localHessian[i][j];        // Use symmetry of hessian 
-        localHyperDualConfiguration[j].seteps2(0.0);
+      localGradient[i] = temp.gradient()[i];
+      localHessian[i][i] = temp.hessian()[i][i];
+
+      for (size_t j=i+1; j<nDoubles; j++)
+      {
+        localHessian[i][j] = temp.hessian()[i][j];
+        localHessian[j][i] = localHessian[i][j];        // Use symmetry of hessian
       }
-        
-      localHyperDualConfiguration[i].seteps1(0.0);
     }
 
 }
diff --git a/dune/elasticity/common/hyperdual.hh b/dune/elasticity/common/hyperdual.hh
index 95ab3de7bdbf37e1f49f88f4435ff20565b0041a..6c6df1140860f6a29acb002095ef1145ffc24c83 100644
--- a/dune/elasticity/common/hyperdual.hh
+++ b/dune/elasticity/common/hyperdual.hh
@@ -35,25 +35,34 @@
 #include <iostream>
 #include <math.h>
 
+#include<dune/common/fmatrix.hh>
+#include <boost/hana/size.hpp>
+
 class hyperdual{
-	double f0,f1,f2,f12;
+
 public:
+
+  const static int size = 29;
+
+  using Matrix = Dune::FieldMatrix<double,size,size>;
+  using Vector = Dune::FieldVector<double,size>;
+
 	//creation operators and function to manually set values
 	hyperdual();
-	hyperdual(double x1,double x2,double x3,double x4);
 	hyperdual(double x1);
-	void setvalues(double x1,double x2,double x3,double x4);
-	void setreal(double x1);
-	void seteps1(double x2);
-	void seteps2(double x3);
-	void seteps12(double x4);
-
-	//examine values
-	void view(void);
-	double real(void) const;
-	double eps1(void) const;
-	double eps2(void) const;
-	double eps1eps2(void) const;
+	hyperdual(double x1, int index);
+
+  double real() const;
+
+  double& value();
+  double value() const;
+
+  Vector& gradient();
+  Vector gradient() const;
+
+  Matrix& hessian();
+  Matrix hessian() const;
+
 	friend std::ostream& operator<<(std::ostream& output, const hyperdual& rhs);
 
 	//basic manipulation
@@ -115,76 +124,78 @@ public:
 	friend bool operator!= (hyperdual lhs, hyperdual rhs);
 	friend bool operator!= (double lhs, hyperdual rhs);
 	friend bool operator!= (hyperdual lhs, double rhs);
+
+private:
+
+	double val;
+  Vector grad;
+
+  // TODO we need only one half of the hessian
+  Matrix hess;
+
 };
 
 
 hyperdual::hyperdual()
-{
-	f0 = 0.0;
-	f1 = 0.0;
-	f2 = 0.0;
-	f12 = 0.0;
-}
-hyperdual::hyperdual(double x1,double x2,double x3,double x4)
-{
-	f0 = x1;
-	f1 = x2;
-	f2 = x3;
-	f12 = x4;
-}
+: val(0.)
+, grad(0.)
+, hess(0.)
+{}
+
+
 hyperdual::hyperdual(double x1)
+: val(x1)
+, grad(0.)
+, hess(0.)
+{}
+
+// active the i-th entry
+hyperdual::hyperdual(double x1, int index)
+: val(x1)
+, grad(0.)
+, hess(0.)
 {
-	f0 = x1;
-	f1 = 0.0;
-	f2 = 0.0;
-	f12 = 0.0;
-}
-void hyperdual::setvalues(double x1,double x2,double x3,double x4)
-{
-	f0 = x1;
-	f1 = x2;
-	f2 = x3;
-	f12 = x4;
-}
-void hyperdual::setreal(double x1)
-{
-    f0 = x1;
-}
-void hyperdual::seteps1(double x2)
-{
-    f1 = x2;
+  grad[index] = 1.;
 }
-void hyperdual::seteps2(double x3)
+
+double hyperdual::real() const
 {
-    f2 = x3;
+  return value();
 }
-void hyperdual::seteps12(double x4)
+
+double& hyperdual::value()
 {
-    f12 = x4;
+  return val;
 }
-void hyperdual::view(void)
+
+double hyperdual::value() const
 {
-	printf("%g  +  %g epsilon1  +  %g epsilon2  +  %g epsilon1 epsilon2\n",f0,f1,f2,f12);
+  return val;
 }
-double hyperdual::real(void) const
+
+hyperdual::Vector& hyperdual::gradient()
 {
-	return f0;
+  return grad;
 }
-double hyperdual::eps1(void) const
+
+hyperdual::Vector hyperdual::gradient() const
 {
-	return f1;
+  return grad;
 }
-double hyperdual::eps2(void) const
+
+hyperdual::Matrix& hyperdual::hessian()
 {
-	return f2;
+  return hess;
 }
-double hyperdual::eps1eps2(void) const
+
+hyperdual::Matrix hyperdual::hessian() const
 {
-	return f12;
+  return hess;
 }
+
 std::ostream& operator<<(std::ostream& output, const hyperdual& rhs)
 {
-        output << "(" << rhs.f0 << ","<< rhs.f1 << ","<< rhs.f2 << ","<< rhs.f12 << ")";
+        output << "value: " << rhs.value() << "\ngradient: " << rhs.gradient() << "\nhessian: " << rhs.hessian();
         return output;
 }
 
@@ -192,421 +203,488 @@ hyperdual hyperdual::operator+ () const
 {
 	return *this;
 }
+
 hyperdual hyperdual::operator+ (const hyperdual rhs) const
 {
-	hyperdual temp;
-	temp.f0 = f0 + rhs.f0;
-	temp.f1 = f1 + rhs.f1;
-	temp.f2 = f2 + rhs.f2;
-	temp.f12 = f12 + rhs.f12;
+	hyperdual temp = *this;
+  temp += rhs;
 	return temp;
 }
+
 hyperdual operator+ (const double lhs, const hyperdual rhs)
 {
-	hyperdual temp;
-	temp.f0 = lhs + rhs.f0;
-	temp.f1 = rhs.f1;
-	temp.f2 = rhs.f2;
-	temp.f12 = rhs.f12;
+	hyperdual temp = rhs;
+  temp.value() += lhs;
 	return temp;
 }
+
 hyperdual hyperdual::operator- () const
 {
-	hyperdual temp;
-	temp.f0 = -f0;
-	temp.f1 = -f1;
-	temp.f2 = -f2;
-	temp.f12 = -f12;
-	return temp;
+	hyperdual temp = *this;
+  temp *= -1.;
+  return temp;
 }
+
 hyperdual hyperdual::operator- (const hyperdual rhs) const
 {
-	hyperdual temp;
-	temp.f0 = f0 - rhs.f0;
-	temp.f1 = f1 - rhs.f1;
-	temp.f2 = f2 - rhs.f2;
-	temp.f12 = f12 - rhs.f12;
+	hyperdual temp = *this;
+	temp -= rhs;
 	return temp;
 }
+
 hyperdual operator- (const double lhs, const hyperdual rhs)
 {
-	hyperdual temp;
-	temp.f0 = lhs - rhs.f0;
-	temp.f1 = -rhs.f1;
-	temp.f2 = -rhs.f2;
-	temp.f12 = -rhs.f12;
+	hyperdual temp = -rhs;
+	temp.val += lhs;
 	return temp;
 }
+
 hyperdual hyperdual::operator* (const hyperdual rhs) const
 {
-	hyperdual temp;
-	temp.f0 = f0*rhs.f0;
-	temp.f1 = f0*rhs.f1 + f1*rhs.f0;
-	temp.f2 = f0*rhs.f2 + f2*rhs.f0;
-	temp.f12 = f0*rhs.f12 + f1*rhs.f2 + f2*rhs.f1 + f12*rhs.f0;
+	hyperdual temp = *this;
+  temp *= rhs;
 	return temp;
 }
+
 hyperdual operator* (const double lhs, const hyperdual rhs)
 {
-	hyperdual temp;
-	temp.f0 = lhs*rhs.f0;
-	temp.f1 = lhs*rhs.f1;
-	temp.f2 = lhs*rhs.f2;
-	temp.f12 = lhs*rhs.f12;
+	hyperdual temp = rhs;
+  temp *= lhs;
 	return temp;
 }
+
 hyperdual operator/ (const hyperdual lhs, const hyperdual rhs)
 {
-	hyperdual inv;
-	double prod;
+	hyperdual temp = lhs;
+  temp /= rhs;
+  return temp;
 
-	prod = rhs.f0*rhs.f0;
-	inv.f0 = 1.0/rhs.f0;
-	inv.f1 = -rhs.f1/prod;
-	inv.f2 = -rhs.f2/prod;
-	inv.f12 = 2.0*rhs.f1*rhs.f2/(prod*rhs.f0) - rhs.f12/prod;
-	return lhs*inv;
 }
+
 hyperdual operator/ (const double lhs, const hyperdual rhs)
 {
-	hyperdual temp,inv;
-	temp = hyperdual(lhs)/rhs;
-	return temp;
+	hyperdual temp(lhs);
+  temp /= rhs;
+  return temp;
 }
+
 hyperdual operator/ (const hyperdual lhs, const double rhs)
 {
-	hyperdual temp;
-	double inv;
-	inv = 1.0/rhs;
-	temp.f0 = inv*lhs.f0;
-	temp.f1 = inv*lhs.f1;
-	temp.f2 = inv*lhs.f2;
-	temp.f12 = inv*lhs.f12;
-	return temp;
+	hyperdual temp = lhs;
+  temp /= rhs;
+  return temp;
 }
+
 hyperdual& hyperdual::operator+= (hyperdual rhs)
 {
-	f0 += rhs.f0;
-	f1 += rhs.f1;
-	f2 += rhs.f2;
-	f12 += rhs.f12;
-	return *this;
+	val += rhs.val;
+  grad += rhs.grad;
+  hess += rhs.hess;
+  return *this;
 }
+
 hyperdual& hyperdual::operator-= (hyperdual rhs)
 {
-	f0 -= rhs.f0;
-	f1 -= rhs.f1;
-	f2 -= rhs.f2;
-	f12 -= rhs.f12;
-	return *this;
+	val -= rhs.val;
+  grad -= rhs.grad;
+  hess -= rhs.hess;
+  return *this;
 }
+
 hyperdual& hyperdual::operator*= (hyperdual rhs)
 {
-	double tf0,tf1,tf2,tf12;
-	tf0=f0;
-	tf1=f1;
-	tf2=f2;
-	tf12=f12;
-	f0 = tf0*rhs.f0;
-	f1 = tf0*rhs.f1 + tf1*rhs.f0;
-	f2 = tf0*rhs.f2 + tf2*rhs.f0;
-	f12 = tf0*rhs.f12 + tf1*rhs.f2 + tf2*rhs.f1 + tf12*rhs.f0;
+  // order is important!
+  for(int i=0; i<size; i++)
+  {
+    for(int j=i; j<size; j++)
+      hess[i][j] = rhs.grad[j]*grad[i] + rhs.val*hess[i][j] + grad[j]*rhs.grad[i] + val*rhs.hess[i][j];
+
+    grad[i] = val*rhs.grad[i] + rhs.val*grad[i];
+  }
+
+  val *= rhs.val;
+
 	return *this;
 }
+
 hyperdual& hyperdual::operator*= (double rhs)
 {
-	f0 *= rhs;
-	f1 *= rhs;
-	f2 *= rhs;
-	f12 *= rhs;
-	return *this;
+  val *= rhs;
+  grad *= rhs;
+  hess *= rhs;
+
+  return *this;
 }
+
 hyperdual& hyperdual::operator/= (hyperdual rhs) 
 {
-    *this = *this / rhs;
-    return *this;
+  // order is important!
+  for(int i=0; i<size; i++)
+  {
+    for(int j=i; j<size; j++)
+      hess[i][j] = -(rhs.grad[j]*grad[i] + rhs.grad[i]*grad[j] + val*rhs.hess[i][j])/(rhs.val*rhs.val)
+                   + hess[i][j]/rhs.val
+                   +2.*val*rhs.grad[i]*rhs.grad[j]/(rhs.val*rhs.val*rhs.val);
+
+    grad[i] = grad[i]/rhs.val - val*rhs.grad[i]/(rhs.val*rhs.val);
+  }
+
+  val /= rhs.val;
+
+	return *this;
 }
+
 hyperdual& hyperdual::operator/= (double rhs)
 {
-	f0 /= rhs;
-	f1 /= rhs;
-	f2 /= rhs;
-	f12 /= rhs;
-	return *this;
+  val /= rhs;
+  grad /= rhs;
+  hess /= rhs;
+
+  return *this;
 }
+
 hyperdual pow (hyperdual x, double a)
 {
-	hyperdual temp;
-	double deriv,xval,tol;
-	xval = x.f0;
-	tol = 1e-15;
-	if (fabs(xval) < tol)
-	{
-		if (xval >= 0)
-			xval = tol;
-		if (xval < 0)
-			xval = -tol;
-	}
-	deriv = a*std::pow(xval,(a-1));
-	//temp.f0 = pow(xval,a);
-	temp.f0 = std::pow(x.f0,a);  //Use actual x value, only use tol for derivs
-	temp.f1 = x.f1*deriv;
-	temp.f2 = x.f2*deriv;
-	temp.f12 = x.f12*deriv + a*(a-1)*x.f1*x.f2*std::pow(xval,(a-2));
+  // pow is expensive, use it just once
+	double xPowMinusTwo = std::pow(x.val,a-2);
+  double xPowMinusOne = xPowMinusTwo*x.val;
+  double xPow = xPowMinusOne*x.val;
+
+  double aAMinusOne = a*(a-1.);
+  double aXPowMinusOne = a*xPowMinusOne;
+  double aAMinusOneXPowMinusTwo = aAMinusOne*xPowMinusTwo;
+
+  hyperdual temp;
+  temp.val = xPow;
+
+  for(int i=0; i<hyperdual::size; i++)
+  {
+    for(int j=i; j<hyperdual::size; j++)
+      temp.hess[i][j] = aAMinusOneXPowMinusTwo*x.grad[i]*x.grad[j] + aXPowMinusOne*x.hess[i][j];
+
+    temp.grad[i] = aXPowMinusOne*x.grad[i];
+  }
 
 	return temp;
 }
+
 hyperdual pow (hyperdual x, hyperdual a)
 {
 	return exp(a*log(x));
 }
+
 hyperdual exp(hyperdual x)
 {
+  // exp is expensive
+  double eX = std::exp(x.val);
+
 	hyperdual temp;
-	double deriv;
-	deriv = std::exp(x.f0);
-	temp.f0 = deriv;
-	temp.f1 = deriv*x.f1;
-	temp.f2 = deriv*x.f2;
-	temp.f12 = deriv*(x.f12 + x.f1*x.f2);
-	return temp;
+	temp.val = eX;
+  for(int i=0; i<hyperdual::size; i++)
+  {
+    for(int j=i; j<hyperdual::size; j++)
+      temp.hess[i][j] = eX*(x.grad[i]*x.grad[j] + x.hess[i][j]);
+
+    temp.grad[i] = eX*x.grad[i];
+  }
+
+  return temp;
 }
+
 hyperdual log(hyperdual x)
 {
+  // log is expensive
+  double logX = std::log(x.val);
+  double inv = 1./x.val;
+
 	hyperdual temp;
-	double deriv1,deriv2;
-	deriv1 = x.f1/x.f0;
-	deriv2 = x.f2/x.f0;
-	temp.f0 = std::log(x.f0);
-	temp.f1 = deriv1;
-	temp.f2 = deriv2;
-	temp.f12 = x.f12/x.f0 - (deriv1*deriv2);
-	return temp;
+	temp.val = logX;
+  for(int i=0; i<hyperdual::size; i++)
+  {
+    for(int j=i; j<hyperdual::size; j++)
+      temp.hess[i][j] = -inv*inv*(x.grad[i]*x.grad[j]) + inv*x.hess[i][j];
+
+    temp.grad[i] = inv*x.grad[i];
+  }
+
+  return temp;
 }
+
 hyperdual sin(hyperdual x)
 {
+  double sinX = std::sin(x.val);
+  double cosX = std::cos(x.val);
+
 	hyperdual temp;
-	double funval,deriv;
-	funval = std::sin(x.f0);
-	deriv = std::cos(x.f0);
-	temp.f0 = funval;
-	temp.f1 = deriv*x.f1;
-	temp.f2 = deriv*x.f2;
-	temp.f12 = deriv*x.f12 - funval*x.f1*x.f2;
-	return temp;
+	temp.val = sinX;
+  for(int i=0; i<hyperdual::size; i++)
+  {
+    for(int j=i; j<hyperdual::size; j++)
+      temp.hess[i][j] = -sinX*(x.grad[i]*x.grad[j]) + cosX*x.hess[i][j];
+
+    temp.grad[i] = cosX*x.grad[i];
+  }
+
+  return temp;
 }
+
 hyperdual cos(hyperdual x)
 {
+  double sinX = std::sin(x.val);
+  double cosX = std::cos(x.val);
+
 	hyperdual temp;
-	double funval,deriv;
-	funval = std::cos(x.f0);
-	deriv = -std::sin(x.f0);
-	temp.f0 = funval;
-	temp.f1 = deriv*x.f1;
-	temp.f2 = deriv*x.f2;
-	temp.f12 = deriv*x.f12 - funval*x.f1*x.f2;
-	return temp;
+	temp.val = cosX;
+  for(int i=0; i<hyperdual::size; i++)
+  {
+    for(int j=i; j<hyperdual::size; j++)
+      temp.hess[i][j] = -cosX*(x.grad[i]*x.grad[j]) - sinX*x.hess[i][j];
+
+    temp.grad[i] = -sinX*x.grad[i];
+  }
+
+  return temp;
 }
 hyperdual tan(hyperdual x)
 {
+  double tanX = std::tan(x.val);
+  double deriv = 1.+tanX*tanX;
+
 	hyperdual temp;
-	double funval,deriv;
-	funval = std::tan(x.f0);
-	deriv  = funval*funval + 1.0;
-	temp.f0 = funval;
-	temp.f1 = deriv*x.f1;
-	temp.f2 = deriv*x.f2;
-	temp.f12 = deriv*x.f12 + x.f1*x.f2*(2*funval*deriv);
-	return temp;
+	temp.val = tanX;
+  for(int i=0; i<hyperdual::size; i++)
+  {
+    for(int j=i; j<hyperdual::size; j++)
+      temp.hess[i][j] = 2.*tanX*deriv*(x.grad[i]*x.grad[j]) + deriv*x.hess[i][j];
+
+    temp.grad[i] = deriv*x.grad[i];
+  }
+
+  return temp;
 }
+
 hyperdual asin(hyperdual x)
 {
+  double asinX = std::asin(x.val);
+  double deriv = 1./std::sqrt(1.-x.val*x.val);
+  double deriv2 = 2.*x.val*deriv*deriv*deriv;
+
 	hyperdual temp;
-	double funval,deriv1,deriv;
-	funval = std::asin(x.f0);
-	deriv1 = 1.0-x.f0*x.f0;
-	deriv = 1.0/std::sqrt(deriv1);
-	temp.f0 = funval;
-	temp.f1 = deriv*x.f1;
-	temp.f2 = deriv*x.f2;
-	temp.f12 = deriv*x.f12 + x.f1*x.f2*(x.f0*std::pow(deriv1,-1.5));
-	return temp;
+	temp.val = asinX;
+  for(int i=0; i<hyperdual::size; i++)
+  {
+    for(int j=i; j<hyperdual::size; j++)
+      temp.hess[i][j] = deriv2*(x.grad[i]*x.grad[j]) + deriv*x.hess[i][j];
+
+    temp.grad[i] = deriv*x.grad[i];
+  }
+
+  return temp;
 }
+
 hyperdual acos(hyperdual x)
 {
+  double acosX = std::acos(x.val);
+  double deriv = -1./std::sqrt(1.-x.val*x.val);
+  double deriv2 = -2.*x.val*deriv*deriv*deriv;
+
 	hyperdual temp;
-	double funval,deriv1,deriv;
-	funval = std::acos(x.f0);
-	deriv1 = 1.0-x.f0*x.f0;
-	deriv = -1.0/std::sqrt(deriv1);
-	temp.f0 = funval;
-	temp.f1 = deriv*x.f1;
-	temp.f2 = deriv*x.f2;
-	temp.f12 = deriv*x.f12 + x.f1*x.f2*(-x.f0*std::pow(deriv1,-1.5));
-	return temp;
+	temp.val = acosX;
+  for(int i=0; i<hyperdual::size; i++)
+  {
+    for(int j=i; j<hyperdual::size; j++)
+      temp.hess[i][j] = deriv2*(x.grad[i]*x.grad[j]) + deriv*x.hess[i][j];
+
+    temp.grad[i] = deriv*x.grad[i];
+  }
+
+  return temp;
 }
+
 hyperdual atan(hyperdual x)
 {
+  double atanX = std::atan(x.val);
+  double deriv = 1./(1.+x.val*x.val);
+  double deriv2 = -2.*x.val*deriv*deriv;
+
 	hyperdual temp;
-	double funval,deriv1,deriv;
-	funval = std::atan(x.f0);
-	deriv1 = 1.0+x.f0*x.f0;
-	deriv = 1.0/deriv1;
-	temp.f0 = funval;
-	temp.f1 = deriv*x.f1;
-	temp.f2 = deriv*x.f2;
-	temp.f12 = deriv*x.f12 + x.f1*x.f2*(-2*x.f0/(deriv1*deriv1));
-	return temp;
+	temp.val = atanX;
+  for(int i=0; i<hyperdual::size; i++)
+  {
+    for(int j=i; j<hyperdual::size; j++)
+      temp.hess[i][j] = deriv2*(x.grad[i]*x.grad[j]) + deriv*x.hess[i][j];
+
+    temp.grad[i] = deriv*x.grad[i];
+  }
+
+  return temp;
 }
+
 hyperdual sqrt(hyperdual x)
 {
+  double sqrtX = std::sqrt(x.val);
+  double deriv = 0.5/sqrtX;
+  double deriv2 = -2.0*deriv*deriv*deriv;
+
 	hyperdual temp;
-	temp.f0 = std::sqrt(x.f0);
-	temp.f1 = 0.5*x.f1/temp.f0;
-	temp.f2 = 0.5*x.f2/temp.f0;
-	temp.f12 = 0.5*(x.f12 - 2.*temp.f1*temp.f2)/temp.f0;
-	return temp;
+	temp.val = sqrtX;
+  for(int i=0; i<hyperdual::size; i++)
+  {
+    for(int j=i; j<hyperdual::size; j++)
+      temp.hess[i][j] = deriv2*(x.grad[i]*x.grad[j]) + deriv*x.hess[i][j];
+
+    temp.grad[i] = deriv*x.grad[i];
+  }
+
+  return temp;
 }
+
 hyperdual fabs(hyperdual x)
 {
-	hyperdual temp;
 	if (x < 0.0)
-		temp = -x;
+		return -x;
 	else
-		temp = x;
-	return temp;
+		return x;
 }
 
 hyperdual abs(hyperdual x)
 {
 	return fabs(x);
 }
+
 hyperdual max(hyperdual x1, hyperdual x2)
 {
-	hyperdual temp;
 	if (x1>x2)
-		temp = x1;
+		return x1;
 	else
-		temp = x2;
-	return temp;
+		return x2;
 }
+
 hyperdual max(hyperdual x1, double x2)
 {
-	hyperdual temp;
 	if (x1>x2)
-		temp = x1;
+    return x1;
 	else
-		temp = x2;
-	return temp;
+		return x2;
 }
+
 hyperdual max(double x1, hyperdual x2)
 {
-	hyperdual temp;
 	if (x1>x2)
-		temp = x1;
+    return x1;
 	else
-		temp = x2;
-	return temp;
+		return x2;
 }
+
 hyperdual min(hyperdual x1, hyperdual x2)
 {
-	hyperdual temp;
 	if (x1<x2)
-		temp = x1;
+    return x1;
 	else
-		temp = x2;
-	return temp;
+		return x2;
 }
+
 hyperdual min(hyperdual x1, double x2)
 {
-	hyperdual temp;
 	if (x1<x2)
-		temp = x1;
+    return x1;
 	else
-		temp = x2;
-	return temp;
+		return x2;
 }
+
 hyperdual min(double x1, hyperdual x2)
 {
-	hyperdual temp;
 	if (x1<x2)
-		temp = x1;
+    return x1;
 	else
-		temp = x2;
-	return temp;
+		return x2;
 }
 
 bool operator> (hyperdual lhs, hyperdual rhs)
 {
-	return (lhs.f0 > rhs.f0);
+	return (lhs.val > rhs.val);
 }
+
 bool operator> (double lhs, hyperdual rhs)
 {
-	return (lhs > rhs.f0);
+	return (lhs > rhs.val);
 }
+
 bool operator> (hyperdual lhs, double rhs)
 {
-	return (lhs.f0 > rhs);
+	return (lhs.val > rhs);
 }
+
 bool operator>= (hyperdual lhs, hyperdual rhs)
 {
-	return (lhs.f0 >= rhs.f0);
+	return (lhs.val >= rhs.val);
 }
 bool operator>= (double lhs, hyperdual rhs)
 {
-	return (lhs >= rhs.f0);
+	return (lhs >= rhs.val);
 }
+
 bool operator>= (hyperdual lhs, double rhs)
 {
-	return (lhs.f0 >= rhs);
+	return (lhs.val >= rhs);
 }
+
 bool operator< (hyperdual lhs, hyperdual rhs)
 {
-	return (lhs.f0 < rhs.f0);
+	return (lhs.val < rhs.val);
 }
+
 bool operator< (double lhs, hyperdual rhs)
 {
-	return (lhs < rhs.f0);
+	return (lhs < rhs.val);
 }
+
 bool operator< (hyperdual lhs, double rhs)
 {
-	return (lhs.f0 < rhs);
+	return (lhs.val < rhs);
 }
+
 bool operator<= (hyperdual lhs, hyperdual rhs)
 {
-	return (lhs.f0 <= rhs.f0);
+	return (lhs.val <= rhs.val);
 }
+
 bool operator<= (double lhs, hyperdual rhs)
 {
-	return (lhs <= rhs.f0);
+	return (lhs <= rhs.val);
 }
 bool operator<= (hyperdual lhs, double rhs)
 {
-	return (lhs.f0 <= rhs);
+	return (lhs.val <= rhs);
 }
+
 bool operator== (hyperdual lhs, hyperdual rhs)
 {
-	return (lhs.f0 == rhs.f0);
+	return (lhs.val == rhs.val);
 }
+
 bool operator== (double lhs, hyperdual rhs)
 {
-	return (lhs == rhs.f0);
+	return (lhs == rhs.val);
 }
 bool operator== (hyperdual lhs, double rhs)
 {
-	return (lhs.f0 == rhs);
+	return (lhs.val == rhs);
 }
+
 bool operator!= (hyperdual lhs, hyperdual rhs)
 {
-	return (lhs.f0 != rhs.f0);
+	return (lhs.val != rhs.val);
 }
+
 bool operator!= (double lhs, hyperdual rhs)
 {
-	return (lhs != rhs.f0);
+	return (lhs != rhs.val);
 }
+
 bool operator!= (hyperdual lhs, double rhs)
 {
-	return (lhs.f0 != rhs);
+	return (lhs.val != rhs);
 }