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] 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