diff --git a/pocketfft_hdronly.h b/pocketfft_hdronly.h
index 7515f79413ff99516511fd384148fea3c30fb8b9..f507572f753f996ee46db3f672c7c14a5034f91c 100644
--- a/pocketfft_hdronly.h
+++ b/pocketfft_hdronly.h
@@ -196,6 +196,13 @@ template<typename T> struct cmplx {
     { r+=other.r; i+=other.i; return *this; }
   template<typename T2>cmplx &operator*= (T2 other)
     { r*=other; i*=other; return *this; }
+  template<typename T2>cmplx &operator*= (const cmplx<T2> &other)
+    {
+    T tmp = r*other.r - i*other.i;
+    i = r*other.i + i*other.r;
+    r = tmp;
+    return *this;
+    }
   cmplx operator+ (const cmplx &other) const
     { return cmplx(r+other.r, i+other.i); }
   cmplx operator- (const cmplx &other) const
@@ -2221,10 +2228,10 @@ template<typename T0> class pocketfft_c
         packplan=unique_ptr<cfftp<T0>>(new cfftp<T0>(length));
       }
 
-    template<typename T> POCKETFFT_NOINLINE void backward(cmplx<T> c[], T0 fct)
+    template<typename T> POCKETFFT_NOINLINE void backward(cmplx<T> c[], T0 fct) const
       { packplan ? packplan->backward(c,fct) : blueplan->backward(c,fct); }
 
-    template<typename T> POCKETFFT_NOINLINE void forward(cmplx<T> c[], T0 fct)
+    template<typename T> POCKETFFT_NOINLINE void forward(cmplx<T> c[], T0 fct) const
       { packplan ? packplan->forward(c,fct) : blueplan->forward(c,fct); }
 
     size_t length() const { return len; }
@@ -2419,31 +2426,65 @@ template<typename T0> class T_dct3
 template<typename T0> class T_dct4
   {
   private:
+    size_t N;
     T_dct2<T0> dct2;
     arr<T0> C;
+    pocketfft_c<T0> fft;
+    arr<cmplx<T0>> C2;
 
   public:
     POCKETFFT_NOINLINE T_dct4(size_t length)
-      : dct2(length), C(length)
+      : N(length), dct2((N&1)?N:1), C((N&1)?N:0),
+        fft((N&1)?1:N/2), C2((N&1) ? 0:N/2)
       {
       constexpr T0 pi = T0(3.141592653589793238462643383279502884197L);
-      for (size_t i=0; i<length; ++i)
-        C[i] = cos(T0(0.5)*pi*(T0(i)+T0(0.5))/T0(length));
+      if (N&1)
+        {
+        for (size_t i=0; i<N; ++i)
+          C[i] = cos(T0(0.5)*pi*(T0(i)+T0(0.5))/T0(N));
+        }
+      else
+        {
+        for (size_t i=0; i<N/2; ++i)
+          {
+          T0 ang = -pi/T0(N)*(T0(i)+T0(0.125));
+          C2[i].Set(cos(ang), sin(ang));
+          }
+        }
       }
 
     template<typename T> POCKETFFT_NOINLINE void exec(T c[], T0 fct) const
       {
       // FIXME: due to the recurrence below, the algorithm is not very accurate
       // Trying to find a better one...
-      size_t N = length();
-      for (size_t i=0; i<N; ++i)
-        c[i] *= C[i];
-      dct2.exec(c, fct);
-      for (size_t i=1; i<N; ++i)
-        c[i] = 2*c[i] - c[i-1];
+      if (N&1)
+        {
+        for (size_t i=0; i<N; ++i)
+          c[i] *= C[i];
+        dct2.exec(c, fct);
+        for (size_t i=1; i<N; ++i)
+          c[i] = 2*c[i] - c[i-1];
+        }
+      else
+        {
+        arr<cmplx<T>> y(N/2);
+        for(size_t i=0; i<N/2; ++i)
+          {
+          y[i].Set(c[2*i],c[N-1-2*i]);
+          y[i] *= C2[i];
+          }
+        fft.forward(y.data(), fct);
+        for(size_t i=0; i<N/2; ++i)
+          y[i] *= C2[i];
+        for(size_t i=0; i<N/2; ++i)
+          {
+          c[2*i] = 2*y[i].r;
+          c[2*i+1] = -2*y[N/2-1-i].i;
+          }
+        }
       }
 
-    size_t length() const { return dct2.length(); }
+    size_t length() const { return N; }
   };
 
 template<typename T0> class T_dst1