diff --git a/pub/src/expr.cpp b/pub/src/expr.cpp
index 23718d04cd22f65afb5889fbe0035818e570b8d3..3845a53ca3a40f57f4320619e143df9c8526d126 100644
--- a/pub/src/expr.cpp
+++ b/pub/src/expr.cpp
@@ -491,10 +491,27 @@ void CTree::tree_val( CTOut *ret, const CContext *ctx ) const
         uint kconv = setConvK( ctx->k );
         uint jconv = kconv == -1 ? -1 : setConvJ( ctx->k, ctx->j, kconv );
         double conv_norm; // normalization constant (integral of resolution)
+        double conv_step; // step in x/t, if equidistant
         COld *fv;
         CScan *sv;
         uint nv;
         if( kconv!=-1 ){
+            // Check whether return grid is equidistant
+            if( n>1 ){
+                conv_step = ((*(ctx->fargs))[n-1]-(*(ctx->fargs))[0]) / (n-1);
+                if( conv_step<0 )
+                    conv_step = 0;
+                else {
+                    for( uint i=0; i<n-1; ++i ) {
+                        if( fabs((*(ctx->fargs))[i+1]-(*(ctx->fargs))[i]
+                                 -conv_step)>1e-8*conv_step ){
+                            conv_step = 0; // not equidistant
+                            break;
+                        }
+                    }
+                }
+            } else
+                conv_step = 0;
             // Preset fv,sv,nv.
             fv = NOlm::MOM[kconv].D();
             if ( !fv )
@@ -504,32 +521,42 @@ void CTree::tree_val( CTOut *ret, const CContext *ctx ) const
             if ( nv < 2 )
                 throw "Resolution " + strg(kconv) + ", scan " + strg(jconv) +
                     ", has only " + strg(nv) + " points";
-            // Calculate normalization integral. At this occasion, check order.
+            // Calculate normalization integral.
+            // At this occasion, assert that x is ordered,
             conv_norm = 0;
-            double step;
+            double dx;
             // special treatment for first and last step:
-            step = sv->x[1] - sv->x[0];
-            if ( step <= 0 )
+            dx = sv->x[1] - sv->x[0];
+            if ( dx <= 0 )
                 throw "Resolution " + strg(kconv) + ", scan " + strg(jconv) +
                     ", not ordered in x";
-            conv_norm += sv->y[0] * step;
-            step = sv->x[nv-1] - sv->x[nv-2];
-            if ( step <= 0 )
+            conv_norm += sv->y[0] * dx;
+            dx = sv->x[nv-1] - sv->x[nv-2];
+            if ( dx <= 0 )
                 throw "Resolution " + strg(kconv) + ", scan " + strg(jconv) +
                     ", not ordered in x";
-            conv_norm += sv->y[nv-1] * step;
+            conv_norm += sv->y[nv-1] * dx;
             // the other steps:
             for( uint i=1; i<nv-1; ++i ) {
-                step = ( sv->x[i+1] - sv->x[i-1] ) / 2;
-                if ( step <= 0 )
+                dx = ( sv->x[i+1] - sv->x[i-1] ) / 2;
+                if ( dx <= 0 )
                     throw "Resolution " + strg(kconv) + ", scan " + 
                         strg(jconv) + ", not ordered in x";
-                conv_norm += sv->y[i] * step;
+                conv_norm += sv->y[i] * dx;
             }
             // check result:
             if ( conv_norm <= 0 )
                 throw "Resolution " + strg(kconv) + ", scan " + 
                     strg(jconv) + ", has non-positive norm " + strg(conv_norm);
+            // Check whether x is equidistant
+            if( conv_step ){
+                for( uint i=0; i<nv-1; ++i ) {
+                    if( fabs(sv->x[i+1]-sv->x[i]-conv_step)>1e-8*conv_step ){
+                        conv_step = 0; // not equidistant
+                        break;
+                    }
+                }
+            }
         }
         // The actual computation is rather different for _CONV and _DIRAC:
         if        ( typ == _CONV ) {
@@ -548,27 +575,43 @@ void CTree::tree_val( CTOut *ret, const CContext *ctx ) const
                 for ( uint i=0; i<n; ++i )
                     ret->vd[i] = r.to_d(i) - shift;
             } else if ( NRead::iofmac("\\cv_mode")==-1 ) {
-                // Convolute theory with resolution.
-                // Simplest version: non-equidistant x_out and x_res.
                 CTOut r;
                 CContext cvctx(*ctx);
-		vector<double> tt(n);
-		cvctx.fargs = &tt;
+                vector<double> tt;
+                cvctx.fargs = &tt;
                 for ( uint i=0; i<n; ++i )
                     ret->vd[i] = 0;
-                for ( uint iv=0; iv<nv; ++iv ) {
-                    double step;
-                    if     ( iv==0 )
-                        step = sv->x[1] - sv->x[0];
-                    else if( iv==nv-1 )
-                        step = sv->x[nv-1] - sv->x[nv-2];
-                    else
-                        step = ( sv->x[iv+1] - sv->x[iv-1] ) / 2;
-                    for ( uint i=0; i<n; ++i )
-                        tt[i] = (*(ctx->fargs))[i] - sv->x[iv] - shift;
+                // Convolute theory with resolution.
+                if( conv_step ){
+                    // Accelerated version with equidistant grids 
+                    uint ntt = n + nv - 1;
+                    tt.resize(ntt);
+                    for ( uint i=0; i<ntt; ++i )
+                        tt[i] = (*(ctx->fargs))[0] - sv->x[0] - shift +
+                            (i-(nv-1))*conv_step;
                     arg[0]->tree_val( &r, &cvctx );
-                    for ( uint i=0; i<n; ++i )
-                        ret->vd[i] += r.to_d(i) * step * sv->y[iv] / conv_norm;
+                    for ( uint iv=0; iv<nv; ++iv )
+                        for ( uint i=0; i<n; ++i )
+                            ret->vd[i] += r.to_d(nv-1+i-iv) * sv->y[iv] *
+                                conv_step / conv_norm;
+                } else {
+                    // Simplest version: non-equidistant x_out and x_res.
+                    tt.resize(n);
+                    for ( uint iv=0; iv<nv; ++iv ) {
+                        double dx;
+                        if     ( iv==0 )
+                            dx = sv->x[1] - sv->x[0];
+                        else if( iv==nv-1 )
+                            dx = sv->x[nv-1] - sv->x[nv-2];
+                        else
+                            dx = ( sv->x[iv+1] - sv->x[iv-1] ) / 2;
+                        for ( uint i=0; i<n; ++i )
+                            tt[i] = (*(ctx->fargs))[i] - sv->x[iv] - shift;
+                        arg[0]->tree_val( &r, &cvctx );
+                        for ( uint i=0; i<n; ++i )
+                            ret->vd[i] += r.to_d(i) * dx * sv->y[iv] /
+                                conv_norm;
+                    }
                 }
             } else {
                 // Convolute theory with resolution.