diff --git a/ql/math/integrals/discreteintegrals.cpp b/ql/math/integrals/discreteintegrals.cpp
index d3381c83e0b..0b7c583725e 100644
--- a/ql/math/integrals/discreteintegrals.cpp
+++ b/ql/math/integrals/discreteintegrals.cpp
@@ -52,13 +52,13 @@ namespace QuantLib {
             const Real dxj   = x[j+1]-x[j];
             const Real dxjp1 = x[j+2]-x[j+1];
 
-            const Real alpha = -dxjp1*(2*x[j]-3*x[j+1]+x[j+2]);
-            const Real dd = x[j+2]-x[j];
+            const Real alpha = dxjp1*(2*dxj-dxjp1);
+            const Real dd = dxj+dxjp1;
             const Real k = dd/(6*dxjp1*dxj);
             const Real beta = dd*dd;
-            const Real gamma = dxj*(x[j]-3*x[j+1]+2*x[j+2]);
+            const Real gamma = dxj*(2*dxjp1-dxj);
 
-            acc(k*alpha*f[j]+k*beta*f[j+1]+k*gamma*f[j+2]);
+            acc(k*(alpha*f[j]+beta*f[j+1]+gamma*f[j+2]));
         }
         if ((n & 1) == 0U) {
             acc(0.5*(x[n-1]-x[n-2])*(f[n-1]+f[n-2]));
@@ -67,24 +67,52 @@ namespace QuantLib {
         return sum(acc);
     }
 
-
     Real DiscreteTrapezoidIntegrator::integrate(
         const ext::function<Real (Real)>& f, Real a, Real b) const {
-            const Array x(maxEvaluations(), a, (b-a)/(maxEvaluations()-1));
-            Array fv(x.size());
-            std::transform(x.begin(), x.end(), fv.begin(), f);
-
+            const Size n=maxEvaluations()-1;
+            const Real d=(b-a)/n;
+            
+            accumulator_set<Real, features<tag::sum> > acc;
+            
+            acc(f(a)/2);
+            for(Size i=0;i<n-1;++i) {
+                a+=d;
+                acc(f(a));
+            }
+            acc(f(b)/2);
+            
             increaseNumberOfEvaluations(maxEvaluations());
-            return DiscreteTrapezoidIntegral()(x, fv);
+            
+            return d*sum(acc);
     }
 
     Real DiscreteSimpsonIntegrator::integrate(
         const ext::function<Real (Real)>& f, Real a, Real b) const {
-            const Array x(maxEvaluations(), a, (b-a)/(maxEvaluations()-1));
-            Array fv(x.size());
-            std::transform(x.begin(), x.end(), fv.begin(), f);
-
+            const Size n=maxEvaluations()-1;
+            const Real d=(b-a)/n, d2=d*2;
+            
+            Real sum=0.0, x=a+d;
+            for(Size i=1;i<n;i+=2) {//to time 4
+                sum+=f(x);
+                x+=d2;
+            }
+            sum*=2;
+
+            x=a+d2;
+            for(Size i=2;i<n-1;i+=2) {//to time 2
+                sum+=f(x);
+                x+=d2;
+            }
+            sum*=2;
+
+            sum+=f(a);
+            if(n&1)
+                sum+=1.5*f(b)+2.5*f(b-d);
+            else
+                sum+=f(b);
+            
             increaseNumberOfEvaluations(maxEvaluations());
-            return DiscreteSimpsonIntegral()(x, fv);
+            
+            return d/3*sum;
     }
 }