diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 8abc3de77c916165572b74003474fb2032a516c6..0ff1732dc580acfe22d0831308851efc7005fce7 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -19,14 +19,16 @@ build_docker:
 test_gcc:
   stage: testing
   script:
-    - python3 setup.py install --user -f
-    - pytest-3 -q python/test
+    - python3 setup.py sdist
+    - 'tar xvzf dist/ducc0*.tar.gz'
+    - 'cd ducc0-* && python3 setup.py install --user -f && pytest-3 -q python/test'
 
-#test_clang:
-#  stage: testing
-#  script:
-#    - CC="clang -fsized-deallocation" python3 setup.py install --user -f
-#    - pytest-3 -q python/test
+test_clang:
+  stage: testing
+  script:
+    - python3 setup.py sdist
+    - 'tar xvzf dist/ducc0*.tar.gz'
+    - 'cd ducc0-* && CC="clang -fsized-deallocation" python3 setup.py install --user -f && pytest-3 -q python/test'
 
 release:
   stage: build_tarballs
diff --git a/ChangeLog b/ChangeLog
index fb7dc3e5d84ad8fb30a266aa09413f3ddb0824b7..57386037e1effcd8ad9c072934b04c05a2252b0c 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,12 @@
+0.3.0:
+- general:
+  - The package should now be installable from PyPI via pip even on MacOS.
+    However, MacOS >= 10.14 is required.
+
+- wgridder:
+  - very substantial performance and scaling improvements
+
+
 0.2.0:
 
 - wgridder:
diff --git a/Dockerfile b/Dockerfile
index 00751db6b7026c4feac1d074a5bb0b391405293b..d849a4e7897f8b64c7d23af677cb4b4ddaca4bd7 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -1,2 +1,2 @@
-FROM debian:testing-slim
+FROM debian:buster-slim
 RUN apt-get update && apt-get install -y git python3-pip python3-scipy python3-pytest python3-pybind11 pybind11-dev clang && rm -rf /var/lib/apt/lists/*
diff --git a/MANIFEST.in b/MANIFEST.in
index 8197390d15a28b7b6f88d2c037370b3803bef24e..cac1ab8ab188086c9a2d159f37b37e01b5e90b47 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -5,6 +5,7 @@ include LICENSE
 include src/ducc0/infra/aligned_array.h
 include src/ducc0/infra/error_handling.h
 include src/ducc0/infra/mav.h
+include src/ducc0/infra/misc_utils.h
 include src/ducc0/infra/simd.h
 include src/ducc0/infra/string_utils.cc
 include src/ducc0/infra/string_utils.h
diff --git a/python/gridder_cxx.h b/python/gridder_cxx.h
index 9c48bc31ac25a0ea7b005dbaba65af59dc90a004..b15e30a6bf44937a57dab32d79f9c42dec2f0595 100644
--- a/python/gridder_cxx.h
+++ b/python/gridder_cxx.h
@@ -33,6 +33,7 @@
 #include "ducc0/infra/error_handling.h"
 #include "ducc0/math/fft.h"
 #include "ducc0/infra/threading.h"
+#include "ducc0/infra/misc_utils.h"
 #include "ducc0/infra/useful_macros.h"
 #include "ducc0/infra/mav.h"
 #include "ducc0/infra/simd.h"
@@ -237,15 +238,15 @@ template<typename T> class GridderConfig
     size_t vlim;
     bool uv_side_fast;
 
-    complex<T> wscreen(T x, T y, T w, bool adjoint) const
+    T phase (T x, T y, T w, bool adjoint) const
       {
       constexpr T pi = T(3.141592653589793238462643383279502884197);
       T tmp = 1-x-y;
       if (tmp<=0) return 1; // no phase factor beyond the horizon
       T nm1 = (-x-y)/(sqrt(tmp)+1); // more accurate form of sqrt(1-x-y)-1
-      T phase = 2*pi*w*nm1;
-      if (adjoint) phase *= -1;
-      return complex<T>(cos(phase), sin(phase));
+      T phs = 2*pi*w*nm1;
+      if (adjoint) phs *= -1;
+      return phs;
       }
 
   public:
@@ -307,8 +308,6 @@ template<typename T> class GridderConfig
             if (i2>=nu) i2-=nu;
             size_t j2 = nv-ny_dirty/2+j;
             if (j2>=nv) j2-=nv;
-            // FIXME: for some reason g++ warns about double-to-float conversion
-            // here, even though there is an explicit cast...
             dirty.v(i,j) = tmav(i2,j2)*T(cfu[icfu]*cfv[icfv]);
             }
           }
@@ -322,33 +321,45 @@ template<typename T> class GridderConfig
              y0 = -0.5*ny_dirty*psy;
       execStatic(nx_dirty/2+1, nthreads, 0, [&](Scheduler &sched)
         {
+        using vtype = native_simd<T>;
+        constexpr size_t vlen=vtype::size();
+        size_t nvec = (ny_dirty/2+1+(vlen-1))/vlen;
+        vector<vtype> ph(nvec), sp(nvec), cp(nvec);
         while (auto rng=sched.getNext()) for(auto i=rng.lo; i<rng.hi; ++i)
           {
           T fx = T(x0+i*psx);
           fx *= fx;
+          size_t ix = nu-nx_dirty/2+i;
+          if (ix>=nu) ix-=nu;
+          size_t i2 = nx_dirty-i;
+          size_t ix2 = nu-nx_dirty/2+i2;
+          if (ix2>=nu) ix2-=nu;
           for (size_t j=0; j<=ny_dirty/2; ++j)
             {
             T fy = T(y0+j*psy);
-            auto ws = wscreen(fx, fy*fy, w, true);
-            size_t ix = nu-nx_dirty/2+i;
-            if (ix>=nu) ix-=nu;
-            size_t jx = nv-ny_dirty/2+j;
-            if (jx>=nv) jx-=nv;
-            dirty.v(i,j) += (tmav(ix,jx)*ws).real(); // lower left
-            size_t i2 = nx_dirty-i, j2 = ny_dirty-j;
-            size_t ix2 = nu-nx_dirty/2+i2;
-            if (ix2>=nu) ix2-=nu;
-            size_t jx2 = nv-ny_dirty/2+j2;
-            if (jx2>=nv) jx2-=nv;
-            if ((i>0)&&(i<i2))
+            ph[j/vlen][j%vlen] = phase(fx, fy*fy, w, true);
+            }
+          for (size_t j=0; j<nvec; ++j)
+            for (size_t k=0; k<vlen; ++k)
+               sp[j][k]=sin(ph[j][k]);
+          for (size_t j=0; j<nvec; ++j)
+            for (size_t k=0; k<vlen; ++k)
+              cp[j][k]=cos(ph[j][k]);
+          if ((i>0)&&(i<i2))
+            for (size_t j=0, jx=nv-ny_dirty/2; j<ny_dirty; ++j, jx=(jx+1>=nv)? jx+1-nv : jx+1)
               {
-              dirty.v(i2,j) += (tmav(ix2,jx)*ws).real(); // lower right
-              if ((j>0)&&(j<j2))
-                dirty.v(i2,j2) += (tmav(ix2,jx2)*ws).real(); // upper right
+              size_t j2 = min(j, ny_dirty-j);
+              T re = cp[j2/vlen][j2%vlen], im = sp[j2/vlen][j2%vlen];
+              dirty.v(i,j) += tmav(ix,jx).real()*re - tmav(ix,jx).imag()*im;
+              dirty.v(i2,j) += tmav(ix2,jx).real()*re - tmav(ix2,jx).imag()*im;
+              }
+          else
+            for (size_t j=0, jx=nv-ny_dirty/2; j<ny_dirty; ++j, jx=(jx+1>=nv)? jx+1-nv : jx+1)
+              {
+              size_t j2 = min(j, ny_dirty-j);
+              T re = cp[j2/vlen][j2%vlen], im = sp[j2/vlen][j2%vlen];
+              dirty.v(i,j) += tmav(ix,jx).real()*re - tmav(ix,jx).imag()*im; // lower left
               }
-            if ((j>0)&&(j<j2))
-              dirty.v(i,j2) += (tmav(ix,jx2)*ws).real(); // upper left
-            }
           }
         });
       }
@@ -420,33 +431,45 @@ template<typename T> class GridderConfig
              y0 = -0.5*ny_dirty*psy;
       execStatic(nx_dirty/2+1, nthreads, 0, [&](Scheduler &sched)
         {
+        using vtype = native_simd<T>;
+        constexpr size_t vlen=vtype::size();
+        size_t nvec = (ny_dirty/2+1+(vlen-1))/vlen;
+        vector<vtype> ph(nvec), sp(nvec), cp(nvec);
         while (auto rng=sched.getNext()) for(auto i=rng.lo; i<rng.hi; ++i)
           {
           T fx = T(x0+i*psx);
           fx *= fx;
+          size_t ix = nu-nx_dirty/2+i;
+          if (ix>=nu) ix-=nu;
+          size_t i2 = nx_dirty-i;
+          size_t ix2 = nu-nx_dirty/2+i2;
+          if (ix2>=nu) ix2-=nu;
           for (size_t j=0; j<=ny_dirty/2; ++j)
             {
             T fy = T(y0+j*psy);
-            auto ws = wscreen(fx, fy*fy, w, false);
-            size_t ix = nu-nx_dirty/2+i;
-            if (ix>=nu) ix-=nu;
-            size_t jx = nv-ny_dirty/2+j;
-            if (jx>=nv) jx-=nv;
-            grid.v(ix,jx) = dirty(i,j)*ws; // lower left
-            size_t i2 = nx_dirty-i, j2 = ny_dirty-j;
-            size_t ix2 = nu-nx_dirty/2+i2;
-            if (ix2>=nu) ix2-=nu;
-            size_t jx2 = nv-ny_dirty/2+j2;
-            if (jx2>=nv) jx2-=nv;
-            if ((i>0)&&(i<i2))
+            ph[j/vlen][j%vlen] = phase(fx, fy*fy, w, false);
+            }
+          for (size_t j=0; j<nvec; ++j)
+            for (size_t k=0; k<vlen; ++k)
+               sp[j][k]=sin(ph[j][k]);
+          for (size_t j=0; j<nvec; ++j)
+            for (size_t k=0; k<vlen; ++k)
+              cp[j][k]=cos(ph[j][k]);
+          if ((i>0)&&(i<i2))
+            for (size_t j=0, jx=nv-ny_dirty/2; j<ny_dirty; ++j, jx=(jx+1>=nv)? jx+1-nv : jx+1)
               {
+              size_t j2 = min(j, ny_dirty-j);
+              auto ws = complex<T>(cp[j2/vlen][j2%vlen],sp[j2/vlen][j2%vlen]);
+              grid.v(ix,jx) = dirty(i,j)*ws; // lower left
               grid.v(ix2,jx) = dirty(i2,j)*ws; // lower right
-              if ((j>0)&&(j<j2))
-                grid.v(ix2,jx2) = dirty(i2,j2)*ws; // upper right
               }
-            if ((j>0)&&(j<j2))
-              grid.v(ix,jx2) = dirty(i,j2)*ws; // upper left
-            }
+          else
+            for (size_t j=0, jx=nv-ny_dirty/2; j<ny_dirty; ++j, jx=(jx+1>=nv)? jx+1-nv : jx+1)
+              {
+              size_t j2 = min(j, ny_dirty-j);
+              auto ws = complex<T>(cp[j2/vlen][j2%vlen],sp[j2/vlen][j2%vlen]);
+              grid.v(ix,jx) = dirty(i,j)*ws; // lower left
+              }
           }
         });
       }
@@ -485,56 +508,24 @@ template<typename T> class GridderConfig
       v=fmod1(v_in*psy)*nv;
       iv0 = min(int(v+vshift)-int(nv), maxiv0);
       }
-
-    void apply_wscreen(const mav<complex<T>,2> &dirty,
-      mav<complex<T>,2> &dirty2, double w, bool adjoint) const
-      {
-      checkShape(dirty.shape(), {nx_dirty, ny_dirty});
-      checkShape(dirty2.shape(), {nx_dirty, ny_dirty});
-
-      double x0 = -0.5*nx_dirty*psx,
-             y0 = -0.5*ny_dirty*psy;
-      execStatic(nx_dirty/2+1, nthreads, 0, [&](Scheduler &sched)
-        {
-        while (auto rng=sched.getNext()) for(auto i=rng.lo; i<rng.hi; ++i)
-          {
-          T fx = T(x0+i*psx);
-          fx *= fx;
-          for (size_t j=0; j<=ny_dirty/2; ++j)
-            {
-            T fy = T(y0+j*psy);
-            auto ws = wscreen(fx, fy*fy, T(w), adjoint);
-            dirty2(i,j) = dirty(i,j)*ws; // lower left
-            size_t i2 = nx_dirty-i, j2 = ny_dirty-j;
-            if ((i>0)&&(i<i2))
-              {
-              dirty2(i2,j) = dirty(i2,j)*ws; // lower right
-              if ((j>0)&&(j<j2))
-                dirty2(i2,j2) = dirty(i2,j2)*ws; // upper right
-              }
-            if ((j>0)&&(j<j2))
-              dirty2(i,j2) = dirty(i,j2)*ws; // upper left
-            }
-          }
-        });
-      }
   };
 
 constexpr int logsquare=4;
 
-template<typename T, typename T2=complex<T>> class Helper
+template<typename T> class Helper
   {
   private:
+    using T2=complex<T>;
     const GridderConfig<T> &gconf;
     int nu, nv, nsafe, supp;
     const T2 *grid_r;
     T2 *grid_w;
-    int su, sv;
+    int su, sv, svvec;
     int iu0, iv0; // start index of the current visibility
     int bu0, bv0; // start index of the current buffer
     T wfac;
 
-    vector<T2> rbuf, wbuf;
+    mav<T,2> rbufr, rbufi, wbufr, wbufi;
     bool do_w_gridding;
     double w0, xdw;
     vector<std::mutex> &locks;
@@ -552,7 +543,7 @@ template<typename T, typename T2=complex<T>> class Helper
         std::lock_guard<std::mutex> lock(locks[idxu]);
         for (int iv=0; iv<sv; ++iv)
           {
-          grid_w[idxu*nv + idxv] += wbuf[iu*sv + iv];
+          grid_w[idxu*nv + idxv] += complex<T>(wbufr(iu,iv), wbufi(iu,iv));
           if (++idxv>=nv) idxv=0;
           }
         }
@@ -569,7 +560,8 @@ template<typename T, typename T2=complex<T>> class Helper
         int idxv = idxv0;
         for (int iv=0; iv<sv; ++iv)
           {
-          rbuf[iu*sv + iv] = grid_r[idxu*nv + idxv];
+          rbufr.v(iu,iv) = grid_r[idxu*nv + idxv].real();
+          rbufi.v(iu,iv) = grid_r[idxu*nv + idxv].imag();
           if (++idxv>=nv) idxv=0;
           }
         if (++idxu>=nu) idxu=0;
@@ -578,8 +570,8 @@ template<typename T, typename T2=complex<T>> class Helper
 
   public:
     size_t nvec;
-    const T2 *p0r;
-    T2 *p0w;
+    const T *p0rr, *p0ri;
+    T *p0wr, *p0wi;
     static constexpr size_t vlen=native_simd<T>::size();
     union kbuf {
       T scalar[64];
@@ -591,10 +583,15 @@ template<typename T, typename T2=complex<T>> class Helper
       vector<std::mutex> &locks_, double w0_=-1, double dw_=-1)
       : gconf(gconf_), nu(gconf.Nu()), nv(gconf.Nv()), nsafe(gconf.Nsafe()),
         supp(gconf.Supp()), grid_r(grid_r_),
-        grid_w(grid_w_), su(2*nsafe+(1<<logsquare)), sv(2*nsafe+(1<<logsquare)),
+        grid_w(grid_w_),
+        su(2*nsafe+(1<<logsquare)),
+        sv(2*nsafe+(1<<logsquare)),
+        svvec(((sv+vlen-1)/vlen)*vlen),
         bu0(-1000000), bv0(-1000000),
-        rbuf(su*sv*(grid_r!=nullptr),T(0)),
-        wbuf(su*sv*(grid_w!=nullptr),T(0)),
+        rbufr({size_t(su),size_t(svvec)*(grid_r!=nullptr)}),
+        rbufi({size_t(su),size_t(svvec)*(grid_r!=nullptr)}),
+        wbufr({size_t(su),size_t(svvec)*(grid_w!=nullptr)}),
+        wbufi({size_t(su),size_t(svvec)*(grid_w!=nullptr)}),
         do_w_gridding(dw_>0),
         w0(w0_),
         xdw(T(1)/dw_),
@@ -605,7 +602,7 @@ template<typename T, typename T2=complex<T>> class Helper
       }
     ~Helper() { if (grid_w) dump(); }
 
-    int lineJump() const { return sv; }
+    int lineJump() const { return svvec; }
     T Wfac() const { return wfac; }
     void prep(const UVW &in)
       {
@@ -621,13 +618,20 @@ template<typename T, typename T2=complex<T>> class Helper
         wfac = krn.eval_single(T(xdw*xsupp*abs(w0-in.w)));
       if ((iu0<bu0) || (iv0<bv0) || (iu0+supp>bu0+su) || (iv0+supp>bv0+sv))
         {
-        if (grid_w) { dump(); fill(wbuf.begin(), wbuf.end(), T(0)); }
+        if (grid_w)
+          {
+          dump();
+          wbufr.apply([](T &v){v=0;});
+          wbufi.apply([](T &v){v=0;});
+          }
         bu0=((((iu0+nsafe)>>logsquare)<<logsquare))-nsafe;
         bv0=((((iv0+nsafe)>>logsquare)<<logsquare))-nsafe;
         if (grid_r) load();
         }
-      p0r = grid_r ? rbuf.data() + sv*(iu0-bu0) + iv0-bv0 : nullptr;
-      p0w = grid_w ? wbuf.data() + sv*(iu0-bu0) + iv0-bv0 : nullptr;
+      p0rr = grid_r ? &rbufr(iu0-bu0, iv0-bv0) : nullptr;
+      p0ri = grid_r ? &rbufi(iu0-bu0, iv0-bv0) : nullptr;
+      p0wr = grid_w ? &wbufr.v(iu0-bu0, iv0-bv0) : nullptr;
+      p0wi = grid_w ? &wbufi.v(iu0-bu0, iv0-bv0) : nullptr;
       }
   };
 
@@ -703,60 +707,88 @@ template<class T, class T2> MsServ<T, T2> makeMsServ
    const mav<idx_t,1> &idx, T2 &ms, const mav<T,2> &wgt)
   { return MsServ<T, T2>(baselines, idx, ms, wgt); }
 
-template<typename T, typename Serv> void x2grid_c
+
+template<size_t SUPP, typename T, typename Serv> [[gnu::hot]] void x2grid_c_helper
   (const GridderConfig<T> &gconf, Serv &srv, mav<complex<T>,2> &grid,
   double w0=-1, double dw=-1)
   {
-  checkShape(grid.shape(), {gconf.Nu(), gconf.Nv()});
-  MR_assert(grid.contiguous(), "grid is not contiguous");
-  size_t supp = gconf.Supp();
+  constexpr size_t vlen=native_simd<T>::size();
+  constexpr size_t NVEC((SUPP+vlen-1)/vlen);
   size_t nthreads = gconf.Nthreads();
   bool do_w_gridding = dw>0;
   vector<std::mutex> locks(gconf.Nu());
-
   size_t np = srv.Nvis();
   execGuided(np, nthreads, 100, 0.2, [&](Scheduler &sched)
     {
     Helper<T> hlp(gconf, nullptr, grid.vdata(), locks, w0, dw);
     int jump = hlp.lineJump();
     const T * DUCC0_RESTRICT ku = hlp.buf.scalar;
-    const T * DUCC0_RESTRICT kv = hlp.buf.scalar+hlp.vlen*hlp.nvec;
+    const auto * DUCC0_RESTRICT kv = hlp.buf.simd+NVEC;
 
     while (auto rng=sched.getNext()) for(auto ipart=rng.lo; ipart<rng.hi; ++ipart)
       {
       UVW coord = srv.getCoord(ipart);
       auto flip = coord.FixW();
       hlp.prep(coord);
-      auto * DUCC0_RESTRICT ptr = hlp.p0w;
+      auto * DUCC0_RESTRICT ptrr = hlp.p0wr;
+      auto * DUCC0_RESTRICT ptri = hlp.p0wi;
       auto v(srv.getVis(ipart));
       if (do_w_gridding) v*=hlp.Wfac();
       if (flip) v=conj(v);
-      for (size_t cu=0; cu<supp; ++cu)
+      for (size_t cu=0; cu<SUPP; ++cu)
         {
         complex<T> tmp(v*ku[cu]);
-        size_t cv=0;
-        for (; cv+3<supp; cv+=4)
+        for (size_t cv=0; cv<NVEC; ++cv)
           {
-          ptr[cv  ] += tmp*kv[cv  ];
-          ptr[cv+1] += tmp*kv[cv+1];
-          ptr[cv+2] += tmp*kv[cv+2];
-          ptr[cv+3] += tmp*kv[cv+3];
+          auto tr = native_simd<T>::loadu(ptrr+cv*hlp.vlen);
+          tr += tmp.real()*kv[cv];
+          tr.storeu(ptrr+cv*hlp.vlen);
+          auto ti = native_simd<T>::loadu(ptri+cv*hlp.vlen);
+          ti += tmp.imag()*kv[cv];
+          ti.storeu(ptri+cv*hlp.vlen);
           }
-        for (; cv<supp; ++cv)
-          ptr[cv] += tmp*kv[cv];
-        ptr+=jump;
+        ptrr+=jump;
+        ptri+=jump;
         }
       }
     });
   }
 
-template<typename T, typename Serv> void grid2x_c
-  (const GridderConfig<T> &gconf, const mav<complex<T>,2> &grid,
-  Serv &srv, double w0=-1, double dw=-1)
+template<typename T, typename Serv> void x2grid_c
+  (const GridderConfig<T> &gconf, Serv &srv, mav<complex<T>,2> &grid,
+  double w0=-1, double dw=-1)
   {
   checkShape(grid.shape(), {gconf.Nu(), gconf.Nv()});
   MR_assert(grid.contiguous(), "grid is not contiguous");
-  size_t supp = gconf.Supp();
+
+  switch(gconf.Supp())
+    {
+    case  1: x2grid_c_helper< 1>(gconf, srv, grid, w0, dw); break;
+    case  2: x2grid_c_helper< 2>(gconf, srv, grid, w0, dw); break;
+    case  3: x2grid_c_helper< 3>(gconf, srv, grid, w0, dw); break;
+    case  4: x2grid_c_helper< 4>(gconf, srv, grid, w0, dw); break;
+    case  5: x2grid_c_helper< 5>(gconf, srv, grid, w0, dw); break;
+    case  6: x2grid_c_helper< 6>(gconf, srv, grid, w0, dw); break;
+    case  7: x2grid_c_helper< 7>(gconf, srv, grid, w0, dw); break;
+    case  8: x2grid_c_helper< 8>(gconf, srv, grid, w0, dw); break;
+    case  9: x2grid_c_helper< 9>(gconf, srv, grid, w0, dw); break;
+    case 10: x2grid_c_helper<10>(gconf, srv, grid, w0, dw); break;
+    case 11: x2grid_c_helper<11>(gconf, srv, grid, w0, dw); break;
+    case 12: x2grid_c_helper<12>(gconf, srv, grid, w0, dw); break;
+    case 13: x2grid_c_helper<13>(gconf, srv, grid, w0, dw); break;
+    case 14: x2grid_c_helper<14>(gconf, srv, grid, w0, dw); break;
+    case 15: x2grid_c_helper<15>(gconf, srv, grid, w0, dw); break;
+    case 16: x2grid_c_helper<16>(gconf, srv, grid, w0, dw); break;
+    default: MR_fail("must not happen");
+    }
+  }
+
+template<size_t SUPP, typename T, typename Serv> [[gnu::hot]] void grid2x_c_helper
+  (const GridderConfig<T> &gconf, const mav<complex<T>,2> &grid,
+  Serv &srv, double w0=-1, double dw=-1)
+  {
+  constexpr size_t vlen=native_simd<T>::size();
+  constexpr size_t NVEC((SUPP+vlen-1)/vlen);
   size_t nthreads = gconf.Nthreads();
   bool do_w_gridding = dw>0;
   vector<std::mutex> locks(gconf.Nu());
@@ -768,29 +800,30 @@ template<typename T, typename Serv> void grid2x_c
     Helper<T> hlp(gconf, grid.data(), nullptr, locks, w0, dw);
     int jump = hlp.lineJump();
     const T * DUCC0_RESTRICT ku = hlp.buf.scalar;
-    const T * DUCC0_RESTRICT kv = hlp.buf.scalar+hlp.vlen*hlp.nvec;
+    const auto * DUCC0_RESTRICT kv = hlp.buf.simd+NVEC;
 
     while (auto rng=sched.getNext()) for(auto ipart=rng.lo; ipart<rng.hi; ++ipart)
       {
       UVW coord = srv.getCoord(ipart);
       auto flip = coord.FixW();
       hlp.prep(coord);
-      complex<T> r = 0;
-      const auto * DUCC0_RESTRICT ptr = hlp.p0r;
-      for (size_t cu=0; cu<supp; ++cu)
+      native_simd<T> rr=0, ri=0;
+      const auto * DUCC0_RESTRICT ptrr = hlp.p0rr;
+      const auto * DUCC0_RESTRICT ptri = hlp.p0ri;
+      for (size_t cu=0; cu<SUPP; ++cu)
         {
-        complex<T> tmp(0);
-        size_t cv=0;
-        for (; cv+3<supp; cv+=4)
-          tmp += ptr[cv  ]*kv[cv  ]
-               + ptr[cv+1]*kv[cv+1]
-               + ptr[cv+2]*kv[cv+2]
-               + ptr[cv+3]*kv[cv+3];
-        for (; cv<supp; ++cv)
-          tmp += ptr[cv] * kv[cv];
-        r += tmp*ku[cu];
-        ptr += jump;
+        native_simd<T> tmpr(0), tmpi(0);
+        for (size_t cv=0; cv<NVEC; ++cv)
+          {
+          tmpr += kv[cv]*native_simd<T>::loadu(ptrr+hlp.vlen*cv);
+          tmpi += kv[cv]*native_simd<T>::loadu(ptri+hlp.vlen*cv);
+          }
+        rr += ku[cu]*tmpr;
+        ri += ku[cu]*tmpi;
+        ptrr += jump;
+        ptri += jump;
         }
+      auto r = complex<T>(reduce(rr, std::plus<>()), reduce(ri, std::plus<>()));
       if (flip) r=conj(r);
       if (do_w_gridding) r*=hlp.Wfac();
       srv.addVis(ipart, r);
@@ -798,6 +831,35 @@ template<typename T, typename Serv> void grid2x_c
     });
   }
 
+template<typename T, typename Serv> void grid2x_c
+  (const GridderConfig<T> &gconf, const mav<complex<T>,2> &grid,
+  Serv &srv, double w0=-1, double dw=-1)
+  {
+  checkShape(grid.shape(), {gconf.Nu(), gconf.Nv()});
+  MR_assert(grid.contiguous(), "grid is not contiguous");
+
+  switch(gconf.Supp())
+    {
+    case  1: grid2x_c_helper< 1>(gconf, grid, srv, w0, dw); break;
+    case  2: grid2x_c_helper< 2>(gconf, grid, srv, w0, dw); break;
+    case  3: grid2x_c_helper< 3>(gconf, grid, srv, w0, dw); break;
+    case  4: grid2x_c_helper< 4>(gconf, grid, srv, w0, dw); break;
+    case  5: grid2x_c_helper< 5>(gconf, grid, srv, w0, dw); break;
+    case  6: grid2x_c_helper< 6>(gconf, grid, srv, w0, dw); break;
+    case  7: grid2x_c_helper< 7>(gconf, grid, srv, w0, dw); break;
+    case  8: grid2x_c_helper< 8>(gconf, grid, srv, w0, dw); break;
+    case  9: grid2x_c_helper< 9>(gconf, grid, srv, w0, dw); break;
+    case 10: grid2x_c_helper<10>(gconf, grid, srv, w0, dw); break;
+    case 11: grid2x_c_helper<11>(gconf, grid, srv, w0, dw); break;
+    case 12: grid2x_c_helper<12>(gconf, grid, srv, w0, dw); break;
+    case 13: grid2x_c_helper<13>(gconf, grid, srv, w0, dw); break;
+    case 14: grid2x_c_helper<14>(gconf, grid, srv, w0, dw); break;
+    case 15: grid2x_c_helper<15>(gconf, grid, srv, w0, dw); break;
+    case 16: grid2x_c_helper<16>(gconf, grid, srv, w0, dw); break;
+    default: MR_fail("must not happen");
+    }
+  }
+
 template<typename T> void apply_global_corrections(const GridderConfig<T> &gconf,
   mav<T,2> &dirty, double dw, bool divide_by_n)
   {
@@ -860,7 +922,7 @@ template<typename T, typename Serv> class WgridHelper
   private:
     Serv &srv;
     double wmin, dw;
-    size_t nplanes, supp;
+    size_t nplanes, supp, nthreads;
     vector<vector<idx_t>> minplane;
     size_t verbosity;
 
@@ -883,11 +945,12 @@ template<typename T, typename Serv> class WgridHelper
       }
 
     template<typename T2> static void update_idx(vector<T2> &v, const vector<T2> &add,
-      const vector<T2> &del)
+      const vector<T2> &del, size_t nthreads)
       {
       MR_assert(v.size()>=del.size(), "must not happen");
       vector<T2> res;
       res.reserve((v.size()+add.size())-del.size());
+#if 0
       auto iin=v.begin(), ein=v.end();
       auto iadd=add.begin(), eadd=add.end();
       auto irem=del.begin(), erem=del.end();
@@ -904,16 +967,51 @@ template<typename T, typename Serv> class WgridHelper
       MR_assert(irem==erem, "must not happen");
       while(iadd!=eadd)
         res.push_back(*(iadd++));
+#else
+      if (v.empty()) //special case
+        {
+        MR_assert(del.empty(), "must not happen");
+        for (auto x: add)
+          res.push_back(x);
+        }
+      else
+        {
+        res.resize((v.size()+add.size())-del.size());
+        execParallel(nthreads, [&](Scheduler &sched) {
+          auto tid = sched.thread_num();
+          auto [lo, hi] = calcShare(nthreads, tid, v.size());
+          if (lo==hi) return; // if interval is empty, do nothing
+          auto iin=v.begin()+lo, ein=v.begin()+hi;
+          auto iadd = (iin==v.begin()) ? add.begin() : lower_bound(add.begin(), add.end(), *iin);
+          auto eadd = (ein==v.end()) ? add.end() : lower_bound(add.begin(), add.end(), *ein);
+          auto irem = (iin==v.begin()) ? del.begin() : lower_bound(del.begin(), del.end(), *iin);
+          auto erem = (ein==v.end()) ? del.end() : lower_bound(del.begin(), del.end(), *ein);
+          auto iout = res.begin()+lo-(irem-del.begin())+(iadd-add.begin());
+          while(iin!=ein)
+            {
+            if ((irem!=erem) && (*iin==*irem))
+              { ++irem; ++iin; } // skip removed entry
+            else if ((iadd!=eadd) && (*iadd<*iin))
+              *(iout++) = *(iadd++); // add new entry
+            else
+              *(iout++) = *(iin++);
+            }
+          MR_assert(irem==erem, "must not happen");
+          while(iadd!=eadd)
+            *(iout++) = *(iadd++);
+          });
+        }
+#endif
       MR_assert(res.size()==(v.size()+add.size())-del.size(), "must not happen");
       v.swap(res);
       }
 
   public:
     WgridHelper(const GridderConfig<T> &gconf, Serv &srv_, size_t verbosity_)
-      : srv(srv_), supp(gconf.Supp()), verbosity(verbosity_), curplane(-1)
+      : srv(srv_), supp(gconf.Supp()), nthreads(gconf.Nthreads()),
+        verbosity(verbosity_), curplane(-1)
       {
       size_t nvis = srv.Nvis();
-      size_t nthreads = gconf.Nthreads();
       double wmax;
 
       wminmax(srv, wmin, wmax);
@@ -942,22 +1040,37 @@ template<typename T, typename Serv> class WgridHelper
         }
 #else
       // more efficient: precalculate final vector sizes and avoid reallocations
-      vector<size_t> cnt(nplanes,0);
-      for(size_t ipart=0; ipart<nvis; ++ipart)
+      vector<int> p0(nvis);
+      mav<size_t,2> cnt({nthreads, nplanes+16}); // safety distance against false sharing
+      execStatic(nvis, nthreads, 0, [&](Scheduler &sched)
         {
-        int plane0 = max(0,int(1+(abs(srv.getCoord(ipart).w)-(0.5*supp*dw)-wmin)/dw));
-        ++cnt[plane0];
+        auto tid=sched.thread_num();
+        while (auto rng=sched.getNext()) for(auto i=rng.lo; i<rng.hi; ++i)
+          {
+          p0[i] = max(0,int(1+(abs(srv.getCoord(i).w)-(0.5*supp*dw)-wmin)/dw));
+          ++cnt.v(tid, p0[i]);
+          }
+        });
+
+      for (size_t p=0; p<nplanes; ++p)
+        {
+        size_t offset=0;
+        for (idx_t tid=0; tid<nthreads; ++tid)
+          {
+          auto tmp = cnt(tid, p);
+          cnt.v(tid, p) = offset;
+          offset += tmp;
+          }
+        minplane[p].resize(offset);
         }
 
       // fill minplane
-      for (size_t j=0; j<nplanes; ++j)
-        minplane[j].resize(cnt[j]);
-      vector<size_t> ofs(nplanes, 0);
-      for (size_t ipart=0; ipart<nvis; ++ipart)
+      execStatic(nvis, nthreads, 0, [&](Scheduler &sched)
         {
-        int plane0 = max(0,int(1+(abs(srv.getCoord(ipart).w)-(0.5*supp*dw)-wmin)/dw));
-        minplane[plane0][ofs[plane0]++]=idx_t(ipart);
-        }
+        auto tid=sched.thread_num();
+        while (auto rng=sched.getNext()) for(auto i=rng.lo; i<rng.hi; ++i)
+          minplane[p0[i]][cnt.v(tid,p0[i])++]=idx_t(i);
+        });
 #endif
       }
 
@@ -972,7 +1085,7 @@ template<typename T, typename Serv> class WgridHelper
     bool advance()
       {
       if (++curplane>=int(nplanes)) return false;
-      update_idx(subidx, minplane[curplane], curplane>=int(supp) ? minplane[curplane-supp] : vector<idx_t>());
+      update_idx(subidx, minplane[curplane], curplane>=int(supp) ? minplane[curplane-supp] : vector<idx_t>(), nthreads);
       if (verbosity>1)
         cout << "Working on plane " << curplane << " containing " << subidx.size()
              << " visibilities" << endl;
@@ -1058,23 +1171,14 @@ template<typename T, typename Serv> void dirty2x(
     }
   }
 
-void calc_share(size_t nshares, size_t myshare, size_t nwork, size_t &lo,
-  size_t &hi)
-  {
-  size_t nbase = nwork/nshares;
-  size_t additional = nwork%nshares;
-  lo = myshare*nbase + ((myshare<additional) ? myshare : additional);
-  hi = lo+nbase+(myshare<additional);
-  }
-
-
 template<typename T> vector<idx_t> getWgtIndices(const Baselines &baselines,
   const GridderConfig<T> &gconf, const mav<T,2> &wgt,
   const mav<complex<T>,2> &ms)
   {
   size_t nrow=baselines.Nrows(),
          nchan=baselines.Nchannels(),
-         nsafe=gconf.Nsafe();
+         nsafe=gconf.Nsafe(),
+         nthreads=gconf.Nthreads();
   bool have_wgt=wgt.size()!=0;
   if (have_wgt) checkShape(wgt.shape(),{nrow,nchan});
   bool have_ms=ms.size()!=0;
@@ -1082,13 +1186,18 @@ template<typename T> vector<idx_t> getWgtIndices(const Baselines &baselines,
   constexpr int side=1<<logsquare;
   size_t nbu = (gconf.Nu()+1+side-1) >> logsquare,
          nbv = (gconf.Nv()+1+side-1) >> logsquare;
-  vector<idx_t> acc(nbu*nbv+1,0);
+  mav<idx_t,2> acc({nthreads, (nbu*nbv+16)}); // the 16 is safety distance to avoid false sharing
   vector<idx_t> tmp(nrow*nchan);
 
-  for (idx_t irow=0, idx=0; irow<nrow; ++irow)
-    for (idx_t ichan=0; ichan<nchan; ++ichan, ++idx)
-      if (((!have_ms ) || (norm(ms(irow,ichan))!=0)) &&
-          ((!have_wgt) || (wgt(irow,ichan)!=0)))
+  execParallel(nthreads, [&](Scheduler &sched)
+    {
+    idx_t tid = sched.thread_num();
+    auto [lo, hi] = calcShare(nthreads, tid, nrow);
+    for(auto irow=idx_t(lo); irow<idx_t(hi); ++irow)
+      {
+      for (idx_t ichan=0, idx=irow*nchan; ichan<nchan; ++ichan, ++idx)
+        if (((!have_ms ) || (norm(ms(irow,ichan))!=0)) &&
+            ((!have_wgt) || (wgt(irow,ichan)!=0)))
         {
         auto uvw = baselines.effectiveCoord(RowChan{irow,idx_t(ichan)});
         if (uvw.w<0) uvw.Flip();
@@ -1097,20 +1206,32 @@ template<typename T> vector<idx_t> getWgtIndices(const Baselines &baselines,
         gconf.getpix(uvw.u, uvw.v, u, v, iu0, iv0);
         iu0 = (iu0+nsafe)>>logsquare;
         iv0 = (iv0+nsafe)>>logsquare;
-        ++acc[nbv*iu0 + iv0 + 1];
         tmp[idx] = nbv*iu0 + iv0;
+        ++acc.v(tid, tmp[idx]);
         }
       else
         tmp[idx] = ~idx_t(0);
+      }
+    });
 
-  for (size_t i=1; i<acc.size(); ++i)
-    acc[i] += acc[i-1];
+  idx_t offset=0;
+  for (idx_t idx=0; idx<nbu*nbv; ++idx)
+    for (idx_t tid=0; tid<nthreads; ++tid)
+      {
+      auto tmp = acc(tid, idx);
+      acc.v(tid, idx) = offset;
+      offset += tmp;
+      }
 
-  vector<idx_t> res(acc.back());
-  for (size_t irow=0, idx=0; irow<nrow; ++irow)
-    for (size_t ichan=0; ichan<nchan; ++ichan, ++idx)
-      if (tmp[idx]!=(~idx_t(0)))
-        res[acc[tmp[idx]]++] = baselines.getIdx(irow, ichan);
+  vector<idx_t> res(offset);
+  execStatic(nrow, nthreads, 0, [&](Scheduler &sched)
+    {
+    idx_t tid = sched.thread_num();
+    while (auto rng=sched.getNext()) for(auto irow=idx_t(rng.lo); irow<idx_t(rng.hi); ++irow)
+      for (size_t ichan=0, idx=irow*nchan; ichan<nchan; ++ichan, ++idx)
+        if (tmp[idx]!=(~idx_t(0)))
+          res[acc.v(tid, tmp[idx])++] = baselines.getIdx(irow, ichan);
+    });
   return res;
   }
 
diff --git a/python/misc.cc b/python/misc.cc
index 7ef2c975d582a525368e3ce0bd39db6817b157ca..74791c0e525b508d31b69a1077143d146cae89cc 100644
--- a/python/misc.cc
+++ b/python/misc.cc
@@ -203,6 +203,8 @@ py::array py_transpose(const py::array &in, py::array &out)
   MR_fail("unsupported datatype");
   }
 
+void py_set_thread_pool_size(size_t new_pool_size)
+  { set_pool_size(new_pool_size); }
 
 const char *misc_DS = R"""(
 Various unsorted utilities
@@ -225,6 +227,8 @@ void add_misc(py::module &msup)
 
   m.def("ascontiguousarray",&py_ascontiguousarray, "in"_a);
   m.def("transpose",&py_transpose, "in"_a, "out"_a);
+
+  m.def("set_thread_pool_size",&py_set_thread_pool_size, "new_pool_size"_a);
   }
 
 }
diff --git a/python/test/test_pointing.py b/python/test/test_pointing.py
index 55286917bcfaaec1325db0982b7938ae08ce4cf5..ce5218aeb7e61ae792090a76eb1774525146b728 100644
--- a/python/test/test_pointing.py
+++ b/python/test/test_pointing.py
@@ -64,8 +64,11 @@ def testp2():
     assert_((quat2==quat3).all(), "problem")
     quat2 *= np.sign(quat2[:,0]).reshape((-1,1))
 
-    from scipy.spatial.transform import Rotation as R
-    from scipy.spatial.transform import Slerp
+    try:
+        from scipy.spatial.transform import Rotation as R
+        from scipy.spatial.transform import Slerp
+    except:
+        pytest.skip()
     times1 = t01 + 1./f1*np.arange(size1)
     r1 = R.from_quat(quat1)
     rrquat = R.from_quat(rquat)
diff --git a/python/test/test_wgridder.py b/python/test/test_wgridder.py
index 6e05d53d3f2d90a0a06b3631721be1c1bf545999..1078e475a871dcc0ba3b2e192cbab8d51e8582f4 100644
--- a/python/test/test_wgridder.py
+++ b/python/test/test_wgridder.py
@@ -108,7 +108,7 @@ def test_adjointness_ms2dirty(nxdirty, nydirty, ofactor, nrow, nchan, epsilon,
 @pmp("fov", (1., 20.))
 def test_ms2dirty_against_wdft2(nxdirty, nydirty, ofactor, nrow, nchan, epsilon, singleprec, wstacking, use_wgt, fov, nthreads):
     if singleprec and epsilon < 5e-5:
-        return
+        pytest.skip()
     rng = np.random.default_rng(42)
     pixsizex = fov*np.pi/180/nxdirty
     pixsizey = fov*np.pi/180/nydirty*1.1
diff --git a/src/ducc0/infra/misc_utils.h b/src/ducc0/infra/misc_utils.h
new file mode 100644
index 0000000000000000000000000000000000000000..9c862f669b93732e3f8fc0e086f45bee0651390b
--- /dev/null
+++ b/src/ducc0/infra/misc_utils.h
@@ -0,0 +1,53 @@
+/*
+ *  This file is part of the MR utility library.
+ *
+ *  This code 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 code 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 code; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+/* Copyright (C) 2019-2020 Max-Planck-Society
+   Author: Martin Reinecke */
+
+#ifndef DUCC0_MISC_UTILS_H
+#define DUCC0_MISC_UTILS_H
+
+#include <tuple>
+
+namespace ducc0 {
+
+namespace detail_misc_utils {
+
+using namespace std;
+
+template<typename T> auto calcShare(size_t nshares, size_t myshare,
+  const T &begin, const T &end)
+  {
+  auto nwork = end-begin;
+  auto nbase = nwork/nshares;
+  auto additional = nwork%nshares;
+  auto lo = begin + (myshare*nbase + ((myshare<additional) ? myshare : additional));
+  auto hi = lo+nbase+(myshare<additional);
+  return make_tuple(lo, hi);
+  }
+
+template<typename T> auto calcShare(size_t nshares, size_t myshare, const T &end)
+  { return calcShare(nshares, myshare, T(0), end); }
+
+}
+
+using detail_misc_utils::calcShare;
+
+}
+
+#endif
diff --git a/src/ducc0/infra/threading.cc b/src/ducc0/infra/threading.cc
index 8948fb3ebafaa3dc68275a1eafd821e41eabbd62..39d4b20265d380b0dd716f3bd365f5985b651320 100644
--- a/src/ducc0/infra/threading.cc
+++ b/src/ducc0/infra/threading.cc
@@ -29,6 +29,7 @@
 #include <queue>
 #include <atomic>
 #include <vector>
+#include <memory>
 #if __has_include(<pthread.h>)
 #include <pthread.h>
 #endif
@@ -40,9 +41,7 @@ namespace detail_threading {
 
 #ifndef DUCC0_NO_THREADING
 
-static const size_t max_threads_ = std::max<size_t>(1, std::thread::hardware_concurrency());
-
-std::atomic<size_t> default_nthreads_(max_threads_);
+std::atomic<size_t> default_nthreads_(std::max<size_t>(1, std::thread::hardware_concurrency()));
 
 size_t get_default_nthreads()
   { return default_nthreads_; }
@@ -50,8 +49,6 @@ size_t get_default_nthreads()
 void set_default_nthreads(size_t new_default_nthreads)
   { default_nthreads_ = std::max<size_t>(1, new_default_nthreads); }
 
-size_t max_threads() { return max_threads_; }
-
 class latch
   {
     std::atomic<size_t> num_left_;
@@ -155,10 +152,10 @@ class thread_pool
       threads_(nthreads)
       { create_threads(); }
 
-    thread_pool(): thread_pool(max_threads_) {}
-
     ~thread_pool() { shutdown(); }
 
+    size_t size() const { return threads_.size(); }
+
     void submit(std::function<void()> work)
       {
       work_queue_.push(move(work));
@@ -179,9 +176,12 @@ class thread_pool
       }
   };
 
-inline thread_pool &get_pool()
+thread_pool &get_pool(size_t nthreads=0)
   {
-  static thread_pool pool;
+  static std::unique_ptr<thread_pool> pool;
+  if ((!pool) && (nthreads==0)) nthreads=default_nthreads_;
+  if ((!pool) || ((nthreads!=0) && (nthreads!=pool->size()))) // resize
+    pool = std::make_unique<thread_pool>(nthreads);
 #if __has_include(<pthread.h>)
   static std::once_flag f;
   call_once(f,
@@ -194,9 +194,12 @@ inline thread_pool &get_pool()
     });
 #endif
 
-  return pool;
+  return *pool;
   }
 
+void set_pool_size(size_t new_pool_size)
+  { get_pool(new_pool_size); }
+
 class Distribution
   {
   private:
@@ -317,7 +320,7 @@ void Distribution::thread_map(std::function<void(Scheduler &)> f)
     return;
     }
 
-  auto & pool = get_pool();
+  auto &pool = get_pool();
   latch counter(nthreads_);
   std::exception_ptr ex;
   std::mutex ex_mut;
diff --git a/src/ducc0/infra/threading.h b/src/ducc0/infra/threading.h
index 043de513a8056a626796da981e1acb376ecdcbda..a4b454577e3d1b9c697be7bfdcffbb1199d53a1b 100644
--- a/src/ducc0/infra/threading.h
+++ b/src/ducc0/infra/threading.h
@@ -47,10 +47,11 @@ class Scheduler
     virtual Range getNext() = 0;
   };
 
-size_t max_threads();
 void set_default_nthreads(size_t new_default_nthreads);
 size_t get_default_nthreads();
 
+void set_pool_size(size_t new_pool_size);
+
 void execSingle(size_t nwork,
   std::function<void(Scheduler &)> func);
 void execStatic(size_t nwork, size_t nthreads, size_t chunksize,
@@ -63,7 +64,7 @@ void execParallel(size_t nthreads, std::function<void(Scheduler &)> func);
 
 } // end of namespace detail_threading
 
-using detail_threading::max_threads;
+using detail_threading::set_pool_size;
 using detail_threading::get_default_nthreads;
 using detail_threading::set_default_nthreads;
 using detail_threading::Scheduler;
diff --git a/src/ducc0/math/fft.h b/src/ducc0/math/fft.h
index fbfb97e017834abd565962516fa91e53c2449510..18b6337de421953538f562daa3083a3d5979483b 100644
--- a/src/ducc0/math/fft.h
+++ b/src/ducc0/math/fft.h
@@ -49,6 +49,7 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #include <complex>
 #include <algorithm>
 #include "ducc0/infra/threading.h"
+#include "ducc0/infra/misc_utils.h"
 #include "ducc0/infra/simd.h"
 #include "ducc0/infra/mav.h"
 #ifndef DUCC0_NO_THREADING
@@ -117,7 +118,7 @@ struct util // hack to avoid duplicate symbols
     size_t parallel = size / (info.shape(axis) * vlen);
     if (info.shape(axis) < 1000)
       parallel /= 4;
-    size_t max_threads = (nthreads==0) ? ducc0::max_threads() : nthreads;
+    size_t max_threads = (nthreads==0) ? ducc0::get_default_nthreads() : nthreads;
     return std::max(size_t(1), std::min(parallel, max_threads));
     }
 #endif
@@ -437,10 +438,7 @@ template<size_t N> class multi_iter
       if (nshares==1) return;
       if (nshares==0) throw std::runtime_error("can't run with zero threads");
       if (myshare>=nshares) throw std::runtime_error("impossible share requested");
-      size_t nbase = rem/nshares;
-      size_t additional = rem%nshares;
-      size_t lo = myshare*nbase + ((myshare<additional) ? myshare : additional);
-      size_t hi = lo+nbase+(myshare<additional);
+      auto [lo, hi] = calcShare(nshares, myshare, rem);
       size_t todo = hi-lo;
 
       size_t chunk = rem;
diff --git a/src/ducc0/math/gridding_kernel.h b/src/ducc0/math/gridding_kernel.h
index 883b45348ee166c7f41f9a67082f47be6711f033..e44a2c4df420bcc0634bafa86c0739a5ff84de19 100644
--- a/src/ducc0/math/gridding_kernel.h
+++ b/src/ducc0/math/gridding_kernel.h
@@ -161,6 +161,8 @@ template<typename T> class HornerKernel: public GriddingKernel<T>
     size_t W, D, nvec;
 
     vector<Tsimd> coeff;
+    const T *scoeff;
+    size_t sstride;
     void (HornerKernel<T>::* evalfunc) (T, native_simd<T> *) const;
     T (HornerKernel<T>::* evalsinglefunc) (T) const;
 
@@ -182,11 +184,10 @@ template<typename T> class HornerKernel: public GriddingKernel<T>
       {
       auto nth = min(W-1, size_t(max(T(0), (x+1)*W*T(0.5))));
       x = (x+1)*W-2*nth-1;
-      auto i = nth/vlen;
-      auto imod = nth%vlen;
-      auto tval = coeff[i][imod];
+      auto ptr = scoeff+nth;
+      auto tval = *ptr;
       for (size_t j=1; j<=DEG; ++j)
-        tval = tval*x + coeff[j*nvec+i][imod];
+        tval = tval*x + ptr[j*sstride];
       return tval;
       }
 
@@ -206,11 +207,10 @@ template<typename T> class HornerKernel: public GriddingKernel<T>
       {
       auto nth = min(W-1, size_t(max(T(0), (x+1)*W*T(0.5))));
       x = (x+1)*W-2*nth-1;
-      auto i = nth/vlen;
-      auto imod = nth%vlen;
-      auto tval = coeff[i][imod];
+      auto ptr = scoeff+nth;
+      auto tval = *ptr;
       for (size_t j=1; j<=D; ++j)
-        tval = tval*x + coeff[j*nvec+i][imod];
+        tval = tval*x + ptr[j*sstride];
       return tval;
       }
 
@@ -274,7 +274,8 @@ template<typename T> class HornerKernel: public GriddingKernel<T>
     HornerKernel(size_t W_, size_t D_, const function<double(double)> &func,
       const KernelCorrection &corr_)
       : W(W_), D(D_), nvec((W+vlen-1)/vlen),
-        coeff(makeCoeff(W_, D_, func)), corr(corr_)
+        coeff(makeCoeff(W_, D_, func)), scoeff(reinterpret_cast<T *>(&coeff[0])),
+        sstride(vlen*nvec), corr(corr_)
       { wire_eval(); }
 
     virtual size_t support() const { return W; }