From d556ae211ce1da3ff7e2c460cf67e37490ced917 Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Mon, 22 Jul 2019 09:37:19 +0200
Subject: [PATCH 1/4] cleanup DCT/DST, part 1

---
 pocketfft_hdronly.h | 133 +++++++++++++++++---------------------------
 1 file changed, 50 insertions(+), 83 deletions(-)

diff --git a/pocketfft_hdronly.h b/pocketfft_hdronly.h
index eba3398..b875db8 100644
--- a/pocketfft_hdronly.h
+++ b/pocketfft_hdronly.h
@@ -2342,40 +2342,29 @@ template<typename T0> class T_dct2
       {
       constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
       size_t N=length();
-      if (N==1)
-        c[0]*=2*fct;
-      else if (N==2)
+      size_t NS2 = (N+1)/2;
+      for (size_t i=2; i<N; i+=2)
         {
-        T x1 = 2*fct*(c[0]+c[1]);
-        c[1] = sqrt2*fct*(c[0]-c[1]);
-        c[0] = x1;
+        T xim1 = T0(0.5)*(c[i-1]+c[i]);
+        c[i] = T0(0.5)*(c[i]-c[i-1]);
+        c[i-1] = xim1;
         }
-      else
+      fftplan.backward(c, fct);
+      for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
         {
-        size_t NS2 = (N+1)/2;
-        for (size_t i=2; i<N; i+=2)
-          {
-          T xim1 = T0(0.5)*(c[i-1]+c[i]);
-          c[i] = T0(0.5)*(c[i]-c[i-1]);
-          c[i-1] = xim1;
-          }
-        fftplan.backward(c, fct);
-        for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
-          {
-          T tmp = twiddle[k-1]*c[kc]+twiddle[kc-1]*c[k];
-          c[kc] = twiddle[k-1]*c[k]-twiddle[kc-1]*c[kc];
-          c[k] = tmp;
-          }
-        if ((N&1)==0)
-          c[NS2] = twiddle[NS2-1]*(c[NS2]+c[NS2]);
-        for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
-          {
-          T tmp = c[k]+c[kc];
-          c[kc] = c[k]-c[kc];
-          c[k] = tmp;
-          }
-        c[0] *= 2;
+        T tmp = twiddle[k-1]*c[kc]+twiddle[kc-1]*c[k];
+        c[kc] = twiddle[k-1]*c[k]-twiddle[kc-1]*c[kc];
+        c[k] = tmp;
+        }
+      if ((N&1)==0)
+        c[NS2] = twiddle[NS2-1]*(c[NS2]+c[NS2]);
+      for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
+        {
+        T tmp = c[k]+c[kc];
+        c[kc] = c[k]-c[kc];
+        c[k] = tmp;
         }
+      c[0] *= 2;
       if (ortho) c[0]/=sqrt2;
       }
 
@@ -2402,40 +2391,29 @@ template<typename T0> class T_dct3
       constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
       size_t N=length();
       if (ortho) c[0]*=sqrt2;
-      if (N==1)
-        c[0]*=fct;
-      else if (N==2)
+      size_t NS2 = (N+1)/2;
+      for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
         {
-        T TSQX = sqrt2*c[1];
-        c[1] = fct*(c[0]-TSQX);
-        c[0] = fct*(c[0]+TSQX);
+        T tmp = c[k]-c[kc];
+        c[k] = c[k]+c[kc];
+        c[kc] = tmp;
         }
-      else
+      if ((N&1)==0)
+        c[NS2] = c[NS2]+c[NS2];
+      for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
         {
-        size_t NS2 = (N+1)/2;
-        for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
-          {
-          T tmp = c[k]-c[kc];
-          c[k] = c[k]+c[kc];
-          c[kc] = tmp;
-          }
-        if ((N&1)==0)
-          c[NS2] = c[NS2]+c[NS2];
-        for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
-          {
-          T tmp = twiddle[k-1]*c[k]-twiddle[kc-1]*c[kc];
-          c[k] = twiddle[k-1]*c[kc]+twiddle[kc-1]*c[k];
-          c[kc] = tmp;
-          }
-        if ((N&1)==0)
-          c[NS2] = twiddle[NS2-1]*c[NS2];
-        fftplan.forward(c, fct);
-        for (size_t i=2; i<N; i+=2)
-          {
-          T xim1 = c[i-1]-c[i];
-          c[i] += c[i-1];
-          c[i-1] = xim1;
-          }
+        T tmp = twiddle[k-1]*c[k]-twiddle[kc-1]*c[kc];
+        c[k] = twiddle[k-1]*c[kc]+twiddle[kc-1]*c[k];
+        c[kc] = tmp;
+        }
+      if ((N&1)==0)
+        c[NS2] = twiddle[NS2-1]*c[NS2];
+      fftplan.forward(c, fct);
+      for (size_t i=2; i<N; i+=2)
+        {
+        T xim1 = c[i-1]-c[i];
+        c[i] += c[i-1];
+        c[i-1] = xim1;
         }
       }
 
@@ -2577,16 +2555,11 @@ template<typename T0> class T_dst2
       {
       constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
       size_t N=length();
-      if (N==1)
-        c[0]*=2*fct;
-      else
-        {
-        for (size_t k=1; k<N; k+=2)
-          c[k] = -c[k];
-        dct.exec(c, fct, false);
-        for (size_t k=0, kc=N-1; k<kc; ++k, --kc)
-          swap(c[k], c[kc]);
-        }
+      for (size_t k=1; k<N; k+=2)
+        c[k] = -c[k];
+      dct.exec(c, fct, false);
+      for (size_t k=0, kc=N-1; k<kc; ++k, --kc)
+        swap(c[k], c[kc]);
       if (ortho) c[0]/=sqrt2;
       }
 
@@ -2607,17 +2580,12 @@ template<typename T0> class T_dst3
       constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
       size_t N=length();
       if (ortho) c[0]*=sqrt2;
-      if (N==1)
-        c[0]*=fct;
-      else
-        {
-        size_t NS2 = N/2;
-        for (size_t k=0, kc=N-1; k<NS2; ++k, --kc)
-          swap(c[k], c[kc]);
-        dct.exec(c, fct, false);
-        for (size_t k=1; k<N; k+=2)
-          c[k] = -c[k];
-        }
+      size_t NS2 = N/2;
+      for (size_t k=0, kc=N-1; k<NS2; ++k, --kc)
+        swap(c[k], c[kc]);
+      dct.exec(c, fct, false);
+      for (size_t k=1; k<N; k+=2)
+        c[k] = -c[k];
       }
 
     size_t length() const { return dct.length(); }
@@ -2635,7 +2603,6 @@ template<typename T0> class T_dst4
     template<typename T> POCKETFFT_NOINLINE void exec(T c[], T0 fct, bool /*ortho*/)
       {
       size_t N=length();
-      //if (N==1) { c[0]*=fct; return; }
       size_t NS2 = N/2;
       for (size_t k=0, kc=N-1; k<NS2; ++k, --kc)
         swap(c[k], c[kc]);
-- 
GitLab


From c2da68d45734bd3fc22b99cc4b5f3d1603172836 Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Mon, 22 Jul 2019 10:22:43 +0200
Subject: [PATCH 2/4] consolidate trigonometric transforms to save cache space
 and precomputation time

---
 pocketfft_hdronly.h | 258 +++++++++++++++++++-------------------------
 1 file changed, 112 insertions(+), 146 deletions(-)

diff --git a/pocketfft_hdronly.h b/pocketfft_hdronly.h
index b875db8..ac2ac7f 100644
--- a/pocketfft_hdronly.h
+++ b/pocketfft_hdronly.h
@@ -2293,6 +2293,7 @@ template<typename T0> class pocketfft_r
 //
 // sine/cosine transforms
 //
+
 template<typename T0> class T_dct1
   {
   private:
@@ -2302,7 +2303,8 @@ template<typename T0> class T_dct1
     POCKETFFT_NOINLINE T_dct1(size_t length)
       : fftplan(2*(length-1)) {}
 
-    template<typename T> POCKETFFT_NOINLINE void exec(T c[], T0 fct, bool ortho) const
+    template<typename T> POCKETFFT_NOINLINE void exec(T c[], T0 fct, bool ortho,
+      int /*type*/, bool /*cosine*/) const
       {
       constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
       size_t N=fftplan.length(), n=N/2+1;
@@ -2323,22 +2325,42 @@ template<typename T0> class T_dct1
     size_t length() const { return fftplan.length()/2+1; }
   };
 
-template<typename T0> class T_dct2
+template<typename T0> class T_dst1
   {
   private:
     pocketfft_r<T0> fftplan;
-    vector<T0> twiddle;
 
   public:
-    POCKETFFT_NOINLINE T_dct2(size_t length)
-      : fftplan(length), twiddle(length)
+    POCKETFFT_NOINLINE T_dst1(size_t length)
+      : fftplan(2*(length+1)) {}
+
+    template<typename T> POCKETFFT_NOINLINE void exec(T c[], T0 fct,
+      bool /*ortho*/, int /*type*/, bool /*cosine*/) const
       {
-      constexpr T0 pi = T0(3.141592653589793238462643383279502884197L);
-      for (size_t i=0; i<length; ++i)
-        twiddle[i] = T0(cos(0.5*pi*T0(i+1)/T0(length)));
+      size_t N=fftplan.length(), n=N/2-1;
+      arr<T> tmp(N);
+      tmp[0] = tmp[n+1] = c[0]*0;
+      for (size_t i=0; i<n; ++i)
+        {
+        tmp[i+1] = c[i];
+        tmp[N-1-i] = -c[i];
+        }
+      fftplan.forward(tmp.data(), fct);
+      for (size_t i=0; i<n; ++i)
+        c[i] = -tmp[2*i+2];
       }
 
-    template<typename T> POCKETFFT_NOINLINE void exec(T c[], T0 fct, bool ortho) const
+    size_t length() const { return fftplan.length()/2-1; }
+  };
+
+template<typename T0> class T_dcst23
+  {
+  private:
+    pocketfft_r<T0> fftplan;
+    vector<T0> twiddle;
+
+    template<typename T> POCKETFFT_NOINLINE void exec2c(T c[], T0 fct,
+      bool ortho) const
       {
       constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
       size_t N=length();
@@ -2368,25 +2390,8 @@ template<typename T0> class T_dct2
       if (ortho) c[0]/=sqrt2;
       }
 
-    size_t length() const { return fftplan.length(); }
-  };
-
-template<typename T0> class T_dct3
-  {
-  private:
-    pocketfft_r<T0> fftplan;
-    vector<T0> twiddle;
-
-  public:
-    POCKETFFT_NOINLINE T_dct3(size_t length)
-      : fftplan(length), twiddle(length)
-      {
-      constexpr T0 pi = T0(3.141592653589793238462643383279502884197L);
-      for (size_t i=0; i<length; ++i)
-        twiddle[i] = T0(cos(0.5*pi*T0(i+1)/T0(length)));
-      }
-
-    template<typename T> POCKETFFT_NOINLINE void exec(T c[], T0 fct, bool ortho) const
+    template<typename T> POCKETFFT_NOINLINE void exec3c(T c[], T0 fct,
+      bool ortho) const
       {
       constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
       size_t N=length();
@@ -2417,10 +2422,55 @@ template<typename T0> class T_dct3
         }
       }
 
+    template<typename T> POCKETFFT_NOINLINE void exec2s(T c[], T0 fct,
+      bool ortho) const
+      {
+      constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
+      size_t N=length();
+      for (size_t k=1; k<N; k+=2)
+        c[k] = -c[k];
+      exec2c(c, fct, false);
+      for (size_t k=0, kc=N-1; k<kc; ++k, --kc)
+        swap(c[k], c[kc]);
+      if (ortho) c[0]/=sqrt2;
+      }
+
+    template<typename T> POCKETFFT_NOINLINE void exec3s(T c[], T0 fct,
+      bool ortho) const
+      {
+      constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
+      size_t N=length();
+      if (ortho) c[0]*=sqrt2;
+      size_t NS2 = N/2;
+      for (size_t k=0, kc=N-1; k<NS2; ++k, --kc)
+        swap(c[k], c[kc]);
+      exec3c(c, fct, false);
+      for (size_t k=1; k<N; k+=2)
+        c[k] = -c[k];
+      }
+
+  public:
+    POCKETFFT_NOINLINE T_dcst23(size_t length)
+      : fftplan(length), twiddle(length)
+      {
+      constexpr T0 pi = T0(3.141592653589793238462643383279502884197L);
+      for (size_t i=0; i<length; ++i)
+        twiddle[i] = T0(cos(0.5*pi*T0(i+1)/T0(length)));
+      }
+
+    template<typename T> POCKETFFT_NOINLINE void exec(T c[], T0 fct, bool ortho,
+      int type, bool cosine) const
+      {
+      if (type==2)
+        cosine ? exec2c(c, fct, ortho) : exec2s(c, fct, ortho);
+      else
+        cosine ? exec3c(c, fct, ortho) : exec3s(c, fct, ortho);
+      }
+
     size_t length() const { return fftplan.length(); }
   };
 
-template<typename T0> class T_dct4
+template<typename T0> class T_dcst4
   {
   // even length algorithm from
   // https://www.appletonaudio.com/blog/2013/derivation-of-fast-dct-4-algorithm-based-on-dft/
@@ -2430,23 +2480,7 @@ template<typename T0> class T_dct4
     unique_ptr<pocketfft_r<T0>> rfft;
     arr<cmplx<T0>> C2;
 
-  public:
-    POCKETFFT_NOINLINE T_dct4(size_t length)
-      : N(length),
-        fft((N&1) ? nullptr : new pocketfft_c<T0>(N/2)),
-        rfft((N&1)? new pocketfft_r<T0>(N) : nullptr),
-        C2((N&1) ? 0 : N/2)
-      {
-      constexpr T0 pi = T0(3.141592653589793238462643383279502884197L);
-      if ((N&1)==0)
-        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, bool /*ortho*/) const
+    template<typename T> POCKETFFT_NOINLINE void exec4c(T c[], T0 fct) const
       {
       constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
       if (N&1)
@@ -2512,106 +2546,38 @@ template<typename T0> class T_dct4
         }
       }
 
-    size_t length() const { return N; }
-  };
-
-template<typename T0> class T_dst1
-  {
-  private:
-    pocketfft_r<T0> fftplan;
-
-  public:
-    POCKETFFT_NOINLINE T_dst1(size_t length)
-      : fftplan(2*(length+1)) {}
-
-    template<typename T> POCKETFFT_NOINLINE void exec(T c[], T0 fct, bool /*ortho*/) const
+    template<typename T> POCKETFFT_NOINLINE void exec4s(T c[], T0 fct) const
       {
-      size_t N=fftplan.length(), n=N/2-1;
-      arr<T> tmp(N);
-      tmp[0] = tmp[n+1] = c[0]*0;
-      for (size_t i=0; i<n; ++i)
-        {
-        tmp[i+1] = c[i];
-        tmp[N-1-i] = -c[i];
-        }
-      fftplan.forward(tmp.data(), fct);
-      for (size_t i=0; i<n; ++i)
-        c[i] = -tmp[2*i+2];
-      }
-
-    size_t length() const { return fftplan.length()/2-1; }
-  };
-
-template<typename T0> class T_dst2
-  {
-  private:
-    T_dct2<T0> dct;
-
-  public:
-    POCKETFFT_NOINLINE T_dst2(size_t length)
-      : dct(length) {}
-
-    template<typename T> POCKETFFT_NOINLINE void exec(T c[], T0 fct, bool ortho) const
-      {
-      constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
       size_t N=length();
-      for (size_t k=1; k<N; k+=2)
-        c[k] = -c[k];
-      dct.exec(c, fct, false);
-      for (size_t k=0, kc=N-1; k<kc; ++k, --kc)
-        swap(c[k], c[kc]);
-      if (ortho) c[0]/=sqrt2;
-      }
-
-    size_t length() const { return dct.length(); }
-  };
-
-template<typename T0> class T_dst3
-  {
-  private:
-    T_dct3<T0> dct;
-
-  public:
-    POCKETFFT_NOINLINE T_dst3(size_t length)
-      : dct(length) {}
-
-    template<typename T> POCKETFFT_NOINLINE void exec(T c[], T0 fct, bool ortho)
-      {
-      constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
-      size_t N=length();
-      if (ortho) c[0]*=sqrt2;
       size_t NS2 = N/2;
       for (size_t k=0, kc=N-1; k<NS2; ++k, --kc)
         swap(c[k], c[kc]);
-      dct.exec(c, fct, false);
+      exec4c(c, fct);
       for (size_t k=1; k<N; k+=2)
         c[k] = -c[k];
       }
 
-    size_t length() const { return dct.length(); }
-  };
-
-template<typename T0> class T_dst4
-  {
-  private:
-    T_dct4<T0> dct;
-
   public:
-    POCKETFFT_NOINLINE T_dst4(size_t length)
-      : dct(length) {}
-
-    template<typename T> POCKETFFT_NOINLINE void exec(T c[], T0 fct, bool /*ortho*/)
+    POCKETFFT_NOINLINE T_dcst4(size_t length)
+      : N(length),
+        fft((N&1) ? nullptr : new pocketfft_c<T0>(N/2)),
+        rfft((N&1)? new pocketfft_r<T0>(N) : nullptr),
+        C2((N&1) ? 0 : N/2)
       {
-      size_t N=length();
-      size_t NS2 = N/2;
-      for (size_t k=0, kc=N-1; k<NS2; ++k, --kc)
-        swap(c[k], c[kc]);
-      dct.exec(c, fct, false);
-      for (size_t k=1; k<N; k+=2)
-        c[k] = -c[k];
+      constexpr T0 pi = T0(3.141592653589793238462643383279502884197L);
+      if ((N&1)==0)
+        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));
+          }
       }
 
-    size_t length() const { return dct.length(); }
+    template<typename T> POCKETFFT_NOINLINE void exec(T c[], T0 fct,
+      bool /*ortho*/, int /*type*/, bool cosine) const
+      { cosine ? exec4c(c, fct) : exec4s(c, fct); }
+
+    size_t length() const { return N; }
   };
 
 
@@ -3060,7 +3026,7 @@ template<typename T> POCKETFFT_NOINLINE void general_hartley(
 
 template<typename Trafo, typename T> POCKETFFT_NOINLINE void general_dcst(
   const cndarr<T> &in, ndarr<T> &out, const shape_t &axes,
-  T fct, bool ortho, size_t POCKETFFT_NTHREADS)
+  T fct, bool ortho, int type, bool cosine, size_t POCKETFFT_NTHREADS)
   {
   shared_ptr<Trafo> plan;
 
@@ -3088,7 +3054,7 @@ template<typename Trafo, typename T> POCKETFFT_NOINLINE void general_dcst(
         for (size_t i=0; i<len; ++i)
           for (size_t j=0; j<vlen; ++j)
             tdatav[i][j] = tin[it.iofs(j,i)];
-        plan->exec(tdatav, fct, ortho);
+        plan->exec(tdatav, fct, ortho, type, cosine);
         for (size_t i=0; i<len; ++i)
           for (size_t j=0; j<vlen; ++j)
             out[it.oofs(j,i)] = tdatav[i][j];
@@ -3099,18 +3065,18 @@ template<typename Trafo, typename T> POCKETFFT_NOINLINE void general_dcst(
       it.advance(1);
       auto tdata = reinterpret_cast<T *>(storage.data());
       if ((&tin[0]==&out[0]) && (it.stride_out()==sizeof(T))) // fully in-place
-        plan->exec(&out[it.oofs(0)], fct, ortho);
+        plan->exec(&out[it.oofs(0)], fct, ortho, type, cosine);
       else if (it.stride_out()==sizeof(T)) // compute FFT in output location
         {
         for (size_t i=0; i<len; ++i)
           out[it.oofs(i)] = tin[it.iofs(i)];
-        plan->exec(&out[it.oofs(0)], fct, ortho);
+        plan->exec(&out[it.oofs(0)], fct, ortho, type, cosine);
         }
       else
         {
         for (size_t i=0; i<len; ++i)
           tdata[i] = tin[it.iofs(i)];
-        plan->exec(tdata, fct, ortho);
+        plan->exec(tdata, fct, ortho, type, cosine);
         for (size_t i=0; i<len; ++i)
           out[it.oofs(i)] = tdata[i];
         }
@@ -3374,13 +3340,13 @@ template<typename T> void dct(const shape_t &shape,
   cndarr<T> ain(data_in, shape, stride_in);
   ndarr<T> aout(data_out, shape, stride_out);
   if (type==1)
-    general_dcst<T_dct1<T>>(ain, aout, axes, fct, ortho, nthreads);
+    general_dcst<T_dct1<T>>(ain, aout, axes, fct, ortho, type, true, nthreads);
   else if (type==2)
-    general_dcst<T_dct2<T>>(ain, aout, axes, fct, ortho, nthreads);
+    general_dcst<T_dcst23<T>>(ain, aout, axes, fct, ortho, type, true, nthreads);
   else if (type==3)
-    general_dcst<T_dct3<T>>(ain, aout, axes, fct, ortho, nthreads);
+    general_dcst<T_dcst23<T>>(ain, aout, axes, fct, ortho, type, true, nthreads);
   else if (type==4)
-    general_dcst<T_dct4<T>>(ain, aout, axes, fct, ortho, nthreads);
+    general_dcst<T_dcst4<T>>(ain, aout, axes, fct, ortho, type, true, nthreads);
   else
     throw runtime_error("unsupported DCT type");
   }
@@ -3395,13 +3361,13 @@ template<typename T> void dst(const shape_t &shape,
   cndarr<T> ain(data_in, shape, stride_in);
   ndarr<T> aout(data_out, shape, stride_out);
   if (type==1)
-    general_dcst<T_dst1<T>>(ain, aout, axes, fct, ortho, nthreads);
+    general_dcst<T_dst1<T>>(ain, aout, axes, fct, ortho, type, false, nthreads);
   else if (type==2)
-    general_dcst<T_dst2<T>>(ain, aout, axes, fct, ortho, nthreads);
+    general_dcst<T_dcst23<T>>(ain, aout, axes, fct, ortho, type, false, nthreads);
   else if (type==3)
-    general_dcst<T_dst3<T>>(ain, aout, axes, fct, ortho, nthreads);
+    general_dcst<T_dcst23<T>>(ain, aout, axes, fct, ortho, type, false, nthreads);
   else if (type==4)
-    general_dcst<T_dst4<T>>(ain, aout, axes, fct, ortho, nthreads);
+    general_dcst<T_dcst4<T>>(ain, aout, axes, fct, ortho, type, false, nthreads);
   else
     throw runtime_error("unsupported DST type");
   }
-- 
GitLab


From c9dc81da35a351b76657f89917e9ef02ebfe7608 Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Mon, 22 Jul 2019 13:15:59 +0200
Subject: [PATCH 3/4] more polishing

---
 pocketfft_hdronly.h | 24 ++++++++++--------------
 1 file changed, 10 insertions(+), 14 deletions(-)

diff --git a/pocketfft_hdronly.h b/pocketfft_hdronly.h
index ac2ac7f..96cdf05 100644
--- a/pocketfft_hdronly.h
+++ b/pocketfft_hdronly.h
@@ -2492,10 +2492,9 @@ template<typename T0> class T_dcst4
         auto SGN_SET = [](T x, size_t i) {return (i%2) ? -x : x;};
         arr<T> y(N);
         size_t n2 = N/2;
-        size_t i;
         {
-        size_t m;
-        for (i=0, m=n2; m<N; ++i, m+=4)
+        size_t i=0, m=n2;
+        for (; m<N; ++i, m+=4)
           y[i] = c[m];
         for (; m<2*N; ++i, m+=4)
           y[i] = -c[2*N-m-1];
@@ -2508,7 +2507,9 @@ template<typename T0> class T_dcst4
           y[i] = c[m];
         }
         rfft->forward(y.data(), fct);
-        for (i=0; i+i+1<n2; ++i)
+        {
+        size_t i=0;
+        for (; i+i+1<n2; ++i)
           {
           size_t k = i+i+1;
           T c1=y[2*k-1], s1=y[2*k], c2=y[2*k+1], s2=y[2*k+2];
@@ -2523,12 +2524,15 @@ template<typename T0> class T_dcst4
           c[i] = sqrt2 * (SGN_SET(cx, (i+1)/2) + SGN_SET(sx, i/2));
           c[N-(i+1)] = sqrt2 * (SGN_SET(cx, (i+2)/2) + SGN_SET(sx, (i+1)/2));
           }
+        }
         c[n2] = sqrt2 * SGN_SET(y[0], (n2+1)/2);
 
         // FFTW-derived code ends here
         }
       else
         {
+        // even length algorithm from
+        // https://www.appletonaudio.com/blog/2013/derivation-of-fast-dct-4-algorithm-based-on-dft/
         arr<cmplx<T>> y(N/2);
         for(size_t i=0; i<N/2; ++i)
           {
@@ -3341,14 +3345,10 @@ template<typename T> void dct(const shape_t &shape,
   ndarr<T> aout(data_out, shape, stride_out);
   if (type==1)
     general_dcst<T_dct1<T>>(ain, aout, axes, fct, ortho, type, true, nthreads);
-  else if (type==2)
-    general_dcst<T_dcst23<T>>(ain, aout, axes, fct, ortho, type, true, nthreads);
-  else if (type==3)
-    general_dcst<T_dcst23<T>>(ain, aout, axes, fct, ortho, type, true, nthreads);
   else if (type==4)
     general_dcst<T_dcst4<T>>(ain, aout, axes, fct, ortho, type, true, nthreads);
   else
-    throw runtime_error("unsupported DCT type");
+    general_dcst<T_dcst23<T>>(ain, aout, axes, fct, ortho, type, true, nthreads);
   }
 
 template<typename T> void dst(const shape_t &shape,
@@ -3362,14 +3362,10 @@ template<typename T> void dst(const shape_t &shape,
   ndarr<T> aout(data_out, shape, stride_out);
   if (type==1)
     general_dcst<T_dst1<T>>(ain, aout, axes, fct, ortho, type, false, nthreads);
-  else if (type==2)
-    general_dcst<T_dcst23<T>>(ain, aout, axes, fct, ortho, type, false, nthreads);
-  else if (type==3)
-    general_dcst<T_dcst23<T>>(ain, aout, axes, fct, ortho, type, false, nthreads);
   else if (type==4)
     general_dcst<T_dcst4<T>>(ain, aout, axes, fct, ortho, type, false, nthreads);
   else
-    throw runtime_error("unsupported DST type");
+    general_dcst<T_dcst23<T>>(ain, aout, axes, fct, ortho, type, false, nthreads);
   }
 
 template<typename T> void r2c(const shape_t &shape_in,
-- 
GitLab


From 88f64fcf1d8ef54025a19c619f5ce8b3e8c58a02 Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Mon, 22 Jul 2019 15:04:17 +0200
Subject: [PATCH 4/4] more compactification

---
 pocketfft_hdronly.h | 157 +++++++++++++++++---------------------------
 1 file changed, 62 insertions(+), 95 deletions(-)

diff --git a/pocketfft_hdronly.h b/pocketfft_hdronly.h
index 96cdf05..b74b1c6 100644
--- a/pocketfft_hdronly.h
+++ b/pocketfft_hdronly.h
@@ -209,13 +209,19 @@ template<typename T> struct cmplx {
     r = tmp;
     return *this;
     }
-  cmplx operator+ (const cmplx &other) const
-    { return cmplx(r+other.r, i+other.i); }
-  cmplx operator- (const cmplx &other) const
-    { return cmplx(r-other.r, i-other.i); }
+  template<typename T2>cmplx &operator+= (const cmplx<T2> &other)
+    { r+=other.r; i+=other.i; return *this; }
+  template<typename T2>cmplx &operator-= (const cmplx<T2> &other)
+    { r-=other.r; i-=other.i; return *this; }
   template<typename T2> auto operator* (const T2 &other) const
     -> cmplx<decltype(r*other)>
     { return {r*other, i*other}; }
+  template<typename T2> auto operator+ (const cmplx<T2> &other) const
+    -> cmplx<decltype(r+other.r)>
+    { return {r+other.r, i+other.i}; }
+  template<typename T2> auto operator- (const cmplx<T2> &other) const
+    -> cmplx<decltype(r+other.r)>
+    { return {r-other.r, i-other.i}; }
   template<typename T2> auto operator* (const cmplx<T2> &other) const
     -> cmplx<decltype(r+other.r)>
     { return {r*other.r-i*other.i, r*other.i + i*other.r}; }
@@ -227,8 +233,9 @@ template<typename T> struct cmplx {
                : Tres(r*other.r-i*other.i, r*other.i+i*other.r);
     }
 };
-template<typename T> void PMC(cmplx<T> &a, cmplx<T> &b,
-  const cmplx<T> &c, const cmplx<T> &d)
+template<typename T> inline void PMINPLACE(T &a, T &b)
+  { T t = a; a+=b; b=t-b; }
+template<typename T> void PMC(T &a, T &b, const T &c, const T &d)
   { a = c+d; b = c-d; }
 template<typename T> cmplx<T> conj(const cmplx<T> &a)
   { return {a.r, -a.i}; }
@@ -845,8 +852,6 @@ template <bool fwd, typename T> void ROTX135(T &a)
   else
     { auto tmp_=a.r; a.r=hsqt2*(-a.r-a.i); a.i=hsqt2*(tmp_-a.i); }
   }
-template<typename T> inline void PMINPLACE(T &a, T &b)
-  { T t = a; a.r+=b.r; a.i+=b.i; b.r=t.r-b.r; b.i=t.i-b.i; }
 
 template<bool fwd, typename T> void pass8 (size_t ido, size_t l1,
   const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
@@ -2359,12 +2364,15 @@ template<typename T0> class T_dcst23
     pocketfft_r<T0> fftplan;
     vector<T0> twiddle;
 
-    template<typename T> POCKETFFT_NOINLINE void exec2c(T c[], T0 fct,
-      bool ortho) const
+    template<typename T> POCKETFFT_NOINLINE void exec2(T c[], T0 fct,
+      bool ortho, bool cosine) const
       {
       constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
       size_t N=length();
       size_t NS2 = (N+1)/2;
+      if (!cosine)
+        for (size_t k=1; k<N; k+=2)
+          c[k] = -c[k];
       for (size_t i=2; i<N; i+=2)
         {
         T xim1 = T0(0.5)*(c[i-1]+c[i]);
@@ -2381,28 +2389,26 @@ template<typename T0> class T_dcst23
       if ((N&1)==0)
         c[NS2] = twiddle[NS2-1]*(c[NS2]+c[NS2]);
       for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
-        {
-        T tmp = c[k]+c[kc];
-        c[kc] = c[k]-c[kc];
-        c[k] = tmp;
-        }
+        PMINPLACE(c[k], c[kc]);
       c[0] *= 2;
+      if (!cosine)
+        for (size_t k=0, kc=N-1; k<kc; ++k, --kc)
+          swap(c[k], c[kc]);
       if (ortho) c[0]/=sqrt2;
       }
 
-    template<typename T> POCKETFFT_NOINLINE void exec3c(T c[], T0 fct,
-      bool ortho) const
+    template<typename T> POCKETFFT_NOINLINE void exec3(T c[], T0 fct,
+      bool ortho, bool cosine) const
       {
       constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
       size_t N=length();
       if (ortho) c[0]*=sqrt2;
       size_t NS2 = (N+1)/2;
+      if (!cosine)
+        for (size_t k=0, kc=N-1; k<NS2; ++k, --kc)
+          swap(c[k], c[kc]);
       for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
-        {
-        T tmp = c[k]-c[kc];
-        c[k] = c[k]+c[kc];
-        c[kc] = tmp;
-        }
+        PMINPLACE(c[k], c[kc]);
       if ((N&1)==0)
         c[NS2] = c[NS2]+c[NS2];
       for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
@@ -2420,33 +2426,9 @@ template<typename T0> class T_dcst23
         c[i] += c[i-1];
         c[i-1] = xim1;
         }
-      }
-
-    template<typename T> POCKETFFT_NOINLINE void exec2s(T c[], T0 fct,
-      bool ortho) const
-      {
-      constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
-      size_t N=length();
-      for (size_t k=1; k<N; k+=2)
-        c[k] = -c[k];
-      exec2c(c, fct, false);
-      for (size_t k=0, kc=N-1; k<kc; ++k, --kc)
-        swap(c[k], c[kc]);
-      if (ortho) c[0]/=sqrt2;
-      }
-
-    template<typename T> POCKETFFT_NOINLINE void exec3s(T c[], T0 fct,
-      bool ortho) const
-      {
-      constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
-      size_t N=length();
-      if (ortho) c[0]*=sqrt2;
-      size_t NS2 = N/2;
-      for (size_t k=0, kc=N-1; k<NS2; ++k, --kc)
-        swap(c[k], c[kc]);
-      exec3c(c, fct, false);
-      for (size_t k=1; k<N; k+=2)
-        c[k] = -c[k];
+      if (!cosine)
+        for (size_t k=1; k<N; k+=2)
+          c[k] = -c[k];
       }
 
   public:
@@ -2460,12 +2442,7 @@ template<typename T0> class T_dcst23
 
     template<typename T> POCKETFFT_NOINLINE void exec(T c[], T0 fct, bool ortho,
       int type, bool cosine) const
-      {
-      if (type==2)
-        cosine ? exec2c(c, fct, ortho) : exec2s(c, fct, ortho);
-      else
-        cosine ? exec3c(c, fct, ortho) : exec3s(c, fct, ortho);
-      }
+      { type==2 ? exec2(c, fct, ortho, cosine) : exec3(c, fct, ortho, cosine); }
 
     size_t length() const { return fftplan.length(); }
   };
@@ -2480,9 +2457,30 @@ template<typename T0> class T_dcst4
     unique_ptr<pocketfft_r<T0>> rfft;
     arr<cmplx<T0>> C2;
 
-    template<typename T> POCKETFFT_NOINLINE void exec4c(T c[], T0 fct) const
+  public:
+    POCKETFFT_NOINLINE T_dcst4(size_t length)
+      : N(length),
+        fft((N&1) ? nullptr : new pocketfft_c<T0>(N/2)),
+        rfft((N&1)? new pocketfft_r<T0>(N) : nullptr),
+        C2((N&1) ? 0 : N/2)
+      {
+      constexpr T0 pi = T0(3.141592653589793238462643383279502884197L);
+      if ((N&1)==0)
+        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,
+      bool /*ortho*/, int /*type*/, bool cosine) const
       {
       constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
+      size_t n2 = N/2;
+      if (!cosine)
+        for (size_t k=0, kc=N-1; k<n2; ++k, --kc)
+          swap(c[k], c[kc]);
       if (N&1)
         {
         // The following code is derived from the FFTW3 function apply_re11()
@@ -2491,7 +2489,6 @@ template<typename T0> class T_dcst4
 
         auto SGN_SET = [](T x, size_t i) {return (i%2) ? -x : x;};
         arr<T> y(N);
-        size_t n2 = N/2;
         {
         size_t i=0, m=n2;
         for (; m<N; ++i, m+=4)
@@ -2533,54 +2530,24 @@ template<typename T0> class T_dcst4
         {
         // even length algorithm from
         // https://www.appletonaudio.com/blog/2013/derivation-of-fast-dct-4-algorithm-based-on-dft/
-        arr<cmplx<T>> y(N/2);
-        for(size_t i=0; i<N/2; ++i)
+        arr<cmplx<T>> y(n2);
+        for(size_t i=0; i<n2; ++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)
+        for(size_t i=0, ic=n2-1; i<n2; ++i, --ic)
           {
-          c[2*i] = 2*y[i].r;
-          c[2*i+1] = -2*y[N/2-1-i].i;
+          c[2*i  ] =  2*(y[i ].r*C2[i ].r-y[i ].i*C2[i ].i);
+          c[2*i+1] = -2*(y[ic].i*C2[ic].r+y[ic].r*C2[ic].i);
           }
         }
+      if (!cosine)
+        for (size_t k=1; k<N; k+=2)
+          c[k] = -c[k];
       }
 
-    template<typename T> POCKETFFT_NOINLINE void exec4s(T c[], T0 fct) const
-      {
-      size_t N=length();
-      size_t NS2 = N/2;
-      for (size_t k=0, kc=N-1; k<NS2; ++k, --kc)
-        swap(c[k], c[kc]);
-      exec4c(c, fct);
-      for (size_t k=1; k<N; k+=2)
-        c[k] = -c[k];
-      }
-
-  public:
-    POCKETFFT_NOINLINE T_dcst4(size_t length)
-      : N(length),
-        fft((N&1) ? nullptr : new pocketfft_c<T0>(N/2)),
-        rfft((N&1)? new pocketfft_r<T0>(N) : nullptr),
-        C2((N&1) ? 0 : N/2)
-      {
-      constexpr T0 pi = T0(3.141592653589793238462643383279502884197L);
-      if ((N&1)==0)
-        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,
-      bool /*ortho*/, int /*type*/, bool cosine) const
-      { cosine ? exec4c(c, fct) : exec4s(c, fct); }
-
     size_t length() const { return N; }
   };
 
-- 
GitLab