Avoid some type-punning and clean up pffft a bit

1 files changed, 148 insertions, 161 deletions

diff --git a/common/pffft.cpp b/common/pffft.cpp
index 0c2dd940..06ae66ec 100644
--- a/common/pffft.cpp
+++ b/common/pffft.cpp
@@ -57,14 +57,18 @@
 
 #include "pffft.h"
 
+#include <array>
 #include <assert.h>
 #include <cmath>
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
+#include <vector>
 
+#include "albit.h"
 #include "almalloc.h"
 #include "alnumbers.h"
+#include "vector.h"
 
 #if defined(__GNUC__)
 #define ALWAYS_INLINE(return_type) inline return_type __attribute__ ((always_inline))
@@ -1328,10 +1332,9 @@ struct PFFFT_Setup {
     int     Ncvec; // nb of complex simd vectors (N/4 if PFFFT_COMPLEX, N/8 if PFFFT_REAL)
     int ifac[15];
     pffft_transform_t transform;
-    float *e;    // points into 'data' , N/4*3 elements
-    float *twiddle; // points into 'data', N/4 elements
 
-    alignas(MALLOC_V4SF_ALIGNMENT) v4sf data[1];
+    float *twiddle; // N/4 elements
+    alignas(MALLOC_V4SF_ALIGNMENT) v4sf e[1]; // N/4*3 elements
 };
 
 PFFFT_Setup *pffft_new_setup(int N, pffft_transform_t transform)
@@ -1348,7 +1351,7 @@ PFFFT_Setup *pffft_new_setup(int N, pffft_transform_t transform)
         assert((N%(SIMD_SZ*SIMD_SZ)) == 0);
 
     const unsigned int Ncvec = static_cast<unsigned>(transform == PFFFT_REAL ? N/2 : N)/SIMD_SZ;
-    size_t storelen{offsetof(PFFFT_Setup, data[0]) + (2u*Ncvec * sizeof(v4sf))};
+    size_t storelen{offsetof(PFFFT_Setup, e[0]) + (2u*Ncvec * sizeof(v4sf))};
 
     void *store{al_calloc(MALLOC_V4SF_ALIGNMENT, storelen)};
     if(!store) return nullptr;
@@ -1358,39 +1361,28 @@ PFFFT_Setup *pffft_new_setup(int N, pffft_transform_t transform)
     s->transform = transform;
     /* nb of complex simd vectors */
     s->Ncvec = static_cast<int>(Ncvec);
-    s->e = reinterpret_cast<float*>(&s->data[0]);
-    s->twiddle = reinterpret_cast<float*>(&s->data[2u*Ncvec*(SIMD_SZ-1)/SIMD_SZ]);
+    s->twiddle = reinterpret_cast<float*>(&s->e[2u*Ncvec*(SIMD_SZ-1)/SIMD_SZ]);
 
-    if(transform == PFFFT_REAL)
+    if constexpr(SIMD_SZ > 1)
     {
+        al::vector<float,16> e(2u*Ncvec*(SIMD_SZ-1));
         for(int k=0; k < s->Ncvec; ++k)
         {
-            int i = k/SIMD_SZ;
-            int j = k%SIMD_SZ;
-            for(int m=0; m < SIMD_SZ-1; ++m)
+            size_t i{static_cast<size_t>(k) / SIMD_SZ};
+            size_t j{static_cast<size_t>(k) % SIMD_SZ};
+            for(size_t m{0};m < SIMD_SZ-1;++m)
             {
-                const double A = -2.0*al::numbers::pi*(m+1)*k / N;
-                s->e[(2*(i*3 + m) + 0) * SIMD_SZ + j] = static_cast<float>(std::cos(A));
-                s->e[(2*(i*3 + m) + 1) * SIMD_SZ + j] = static_cast<float>(std::sin(A));
+                const double A = -2.0*al::numbers::pi*static_cast<double>(m+1)*k / N;
+                e[(2*(i*3 + m) + 0)*SIMD_SZ + j] = static_cast<float>(std::cos(A));
+                e[(2*(i*3 + m) + 1)*SIMD_SZ + j] = static_cast<float>(std::sin(A));
             }
         }
-        rffti1_ps(N/SIMD_SZ, s->twiddle, s->ifac);
+        std::memcpy(s->e, e.data(), e.size()*sizeof(float));
     }
+    if(transform == PFFFT_REAL)
+        rffti1_ps(N/SIMD_SZ, s->twiddle, s->ifac);
     else
-    {
-        for(int k=0; k < s->Ncvec; ++k)
-        {
-            int i = k/SIMD_SZ;
-            int j = k%SIMD_SZ;
-            for(int m=0; m < SIMD_SZ-1; ++m)
-            {
-                const double A = -2.0*al::numbers::pi*(m+1)*k / N;
-                s->e[(2*(i*3 + m) + 0)*SIMD_SZ + j] = static_cast<float>(std::cos(A));
-                s->e[(2*(i*3 + m) + 1)*SIMD_SZ + j] = static_cast<float>(std::sin(A));
-            }
-        }
         cffti1_ps(N/SIMD_SZ, s->twiddle, s->ifac);
-    }
 
     /* check that N is decomposable with allowed prime factors */
     int m = 1;
@@ -1455,13 +1447,14 @@ static void unreversed_copy(int N, const v4sf *in, v4sf *out, int out_stride)
 
 void pffft_zreorder(PFFFT_Setup *setup, const float *in, float *out, pffft_direction_t direction)
 {
-    const int N = setup->N, Ncvec = setup->Ncvec;
-    const v4sf *vin = reinterpret_cast<const v4sf*>(in);
-    v4sf *vout = reinterpret_cast<v4sf*>(out);
     assert(in != out);
+
+    const int N{setup->N}, Ncvec{setup->Ncvec};
+    const v4sf *vin{reinterpret_cast<const v4sf*>(in)};
+    v4sf *vout{reinterpret_cast<v4sf*>(out)};
     if(setup->transform == PFFFT_REAL)
     {
-        const int dk = N/32;
+        const int dk{N/32};
         if(direction == PFFFT_FORWARD)
         {
             for(int k=0; k < dk; ++k)
@@ -1469,8 +1462,8 @@ void pffft_zreorder(PFFFT_Setup *setup, const float *in, float *out, pffft_direc
                 INTERLEAVE2(vin[k*8 + 0], vin[k*8 + 1], vout[2*(0*dk + k) + 0], vout[2*(0*dk + k) + 1]);
                 INTERLEAVE2(vin[k*8 + 4], vin[k*8 + 5], vout[2*(2*dk + k) + 0], vout[2*(2*dk + k) + 1]);
             }
-            reversed_copy(dk, vin+2, 8, reinterpret_cast<v4sf*>(out + N/2));
-            reversed_copy(dk, vin+6, 8, reinterpret_cast<v4sf*>(out + N));
+            reversed_copy(dk, vin+2, 8, vout + N/SIMD_SZ/2);
+            reversed_copy(dk, vin+6, 8, vout + N/SIMD_SZ);
         }
         else
         {
@@ -1479,8 +1472,8 @@ void pffft_zreorder(PFFFT_Setup *setup, const float *in, float *out, pffft_direc
                 UNINTERLEAVE2(vin[2*(0*dk + k) + 0], vin[2*(0*dk + k) + 1], vout[k*8 + 0], vout[k*8 + 1]);
                 UNINTERLEAVE2(vin[2*(2*dk + k) + 0], vin[2*(2*dk + k) + 1], vout[k*8 + 4], vout[k*8 + 5]);
             }
-            unreversed_copy(dk, reinterpret_cast<const v4sf*>(in + N/4), reinterpret_cast<v4sf*>(out + N - 6*SIMD_SZ), -8);
-            unreversed_copy(dk, reinterpret_cast<const v4sf*>(in + 3*N/4), reinterpret_cast<v4sf*>(out + N - 2*SIMD_SZ), -8);
+            unreversed_copy(dk, vin + N/SIMD_SZ/4, vout + N/SIMD_SZ - 6, -8);
+            unreversed_copy(dk, vin + 3*N/SIMD_SZ/4, vout + N/SIMD_SZ - 2, -8);
         }
     }
     else
@@ -1504,31 +1497,29 @@ void pffft_zreorder(PFFFT_Setup *setup, const float *in, float *out, pffft_direc
     }
 }
 
-void pffft_cplx_finalize(int Ncvec, const v4sf *in, v4sf *out, const v4sf *e)
+void pffft_cplx_finalize(const int Ncvec, const v4sf *in, v4sf *out, const v4sf *e)
 {
-    const int dk = Ncvec/SIMD_SZ; // number of 4x4 matrix blocks
-    v4sf r0, i0, r1, i1, r2, i2, r3, i3;
-    v4sf sr0, dr0, sr1, dr1, si0, di0, si1, di1;
     assert(in != out);
+
+    const int dk{Ncvec/SIMD_SZ}; // number of 4x4 matrix blocks
     for(int k=0; k < dk; ++k)
     {
-        r0 = in[8*k+0]; i0 = in[8*k+1];
-        r1 = in[8*k+2]; i1 = in[8*k+3];
-        r2 = in[8*k+4]; i2 = in[8*k+5];
-        r3 = in[8*k+6]; i3 = in[8*k+7];
+        v4sf r0{in[8*k+0]}, i0{in[8*k+1]};
+        v4sf r1{in[8*k+2]}, i1{in[8*k+3]};
+        v4sf r2{in[8*k+4]}, i2{in[8*k+5]};
+        v4sf r3{in[8*k+6]}, i3{in[8*k+7]};
         VTRANSPOSE4(r0,r1,r2,r3);
         VTRANSPOSE4(i0,i1,i2,i3);
         VCPLXMUL(r1,i1,e[k*6+0],e[k*6+1]);
         VCPLXMUL(r2,i2,e[k*6+2],e[k*6+3]);
         VCPLXMUL(r3,i3,e[k*6+4],e[k*6+5]);
 
-        sr0 = VADD(r0,r2); dr0 = VSUB(r0, r2);
-        sr1 = VADD(r1,r3); dr1 = VSUB(r1, r3);
-        si0 = VADD(i0,i2); di0 = VSUB(i0, i2);
-        si1 = VADD(i1,i3); di1 = VSUB(i1, i3);
+        v4sf sr0{VADD(r0,r2)}, dr0{VSUB(r0, r2)};
+        v4sf sr1{VADD(r1,r3)}, dr1{VSUB(r1, r3)};
+        v4sf si0{VADD(i0,i2)}, di0{VSUB(i0, i2)};
+        v4sf si1{VADD(i1,i3)}, di1{VSUB(i1, i3)};
 
-        /*
-         * transformation for each column is:
+        /* transformation for each column is:
          *
          * [1   1   1   1   0   0   0   0]   [r0]
          * [1   0  -1   0   0  -1   0   1]   [r1]
@@ -1550,23 +1541,22 @@ void pffft_cplx_finalize(int Ncvec, const v4sf *in, v4sf *out, const v4sf *e)
     }
 }
 
-void pffft_cplx_preprocess(int Ncvec, const v4sf *in, v4sf *out, const v4sf *e)
+void pffft_cplx_preprocess(const int Ncvec, const v4sf *in, v4sf *out, const v4sf *e)
 {
-    const int dk = Ncvec/SIMD_SZ; // number of 4x4 matrix blocks
-    v4sf r0, i0, r1, i1, r2, i2, r3, i3;
-    v4sf sr0, dr0, sr1, dr1, si0, di0, si1, di1;
     assert(in != out);
+
+    const int dk{Ncvec/SIMD_SZ}; // number of 4x4 matrix blocks
     for(int k=0; k < dk; ++k)
     {
-        r0 = in[8*k+0]; i0 = in[8*k+1];
-        r1 = in[8*k+2]; i1 = in[8*k+3];
-        r2 = in[8*k+4]; i2 = in[8*k+5];
-        r3 = in[8*k+6]; i3 = in[8*k+7];
+        v4sf r0{in[8*k+0]}, i0{in[8*k+1]};
+        v4sf r1{in[8*k+2]}, i1{in[8*k+3]};
+        v4sf r2{in[8*k+4]}, i2{in[8*k+5]};
+        v4sf r3{in[8*k+6]}, i3{in[8*k+7]};
 
-        sr0 = VADD(r0,r2); dr0 = VSUB(r0, r2);
-        sr1 = VADD(r1,r3); dr1 = VSUB(r1, r3);
-        si0 = VADD(i0,i2); di0 = VSUB(i0, i2);
-        si1 = VADD(i1,i3); di1 = VSUB(i1, i3);
+        v4sf sr0{VADD(r0,r2)}, dr0{VSUB(r0, r2)};
+        v4sf sr1{VADD(r1,r3)}, dr1{VSUB(r1, r3)};
+        v4sf si0{VADD(i0,i2)}, di0{VSUB(i0, i2)};
+        v4sf si1{VADD(i1,i3)}, di1{VSUB(i1, i3)};
 
         r0 = VADD(sr0, sr1); i0 = VADD(si0, si1);
         r1 = VSUB(dr0, di1); i1 = VADD(di0, dr1);
@@ -1645,20 +1635,18 @@ static ALWAYS_INLINE(void) pffft_real_finalize_4x4(const v4sf *in0, const v4sf *
 
 static NEVER_INLINE(void) pffft_real_finalize(int Ncvec, const v4sf *in, v4sf *out, const v4sf *e)
 {
-    static constexpr float s = al::numbers::sqrt2_v<float>/2.0f;
-    const int dk = Ncvec/SIMD_SZ; // number of 4x4 matrix blocks
-    /* fftpack order is f0r f1r f1i f2r f2i ... f(n-1)r f(n-1)i f(n)r */
-
-    v4sf_union cr, ci, *uout = reinterpret_cast<v4sf_union*>(out);
-    v4sf save = in[7], zero=VZERO();
-    float xr0, xi0, xr1, xi1, xr2, xi2, xr3, xi3;
+    static constexpr float s{al::numbers::sqrt2_v<float>/2.0f};
 
-    cr.v = in[0]; ci.v = in[Ncvec*2-1];
     assert(in != out);
+    const int dk{Ncvec/SIMD_SZ}; // number of 4x4 matrix blocks
+    /* fftpack order is f0r f1r f1i f2r f2i ... f(n-1)r f(n-1)i f(n)r */
+
+    const v4sf zero{VZERO()};
+    const auto cr = al::bit_cast<std::array<float,SIMD_SZ>>(in[0]);
+    const auto ci = al::bit_cast<std::array<float,SIMD_SZ>>(in[Ncvec*2-1]);
     pffft_real_finalize_4x4(&zero, &zero, in+1, e, out);
 
-    /*
-     * [cr0 cr1 cr2 cr3 ci0 ci1 ci2 ci3]
+    /* [cr0 cr1 cr2 cr3 ci0 ci1 ci2 ci3]
      *
      * [Xr(1)]  ] [1   1   1   1   0   0   0   0]
      * [Xr(N/4) ] [0   0   0   0   1   s   0  -s]
@@ -1670,29 +1658,26 @@ static NEVER_INLINE(void) pffft_real_finalize(int Ncvec, const v4sf *in, v4sf *o
      * [Xi(3N/4)] [0   0   0   0   0  -s   1  -s]
      */
 
-    xr0=(cr.f[0]+cr.f[2]) + (cr.f[1]+cr.f[3]); uout[0].f[0] = xr0;
-    xi0=(cr.f[0]+cr.f[2]) - (cr.f[1]+cr.f[3]); uout[1].f[0] = xi0;
-    xr2=(cr.f[0]-cr.f[2]);                     uout[4].f[0] = xr2;
-    xi2=(cr.f[3]-cr.f[1]);                     uout[5].f[0] = xi2;
-    xr1= ci.f[0] + s*(ci.f[1]-ci.f[3]);        uout[2].f[0] = xr1;
-    xi1=-ci.f[2] - s*(ci.f[1]+ci.f[3]);        uout[3].f[0] = xi1;
-    xr3= ci.f[0] - s*(ci.f[1]-ci.f[3]);        uout[6].f[0] = xr3;
-    xi3= ci.f[2] - s*(ci.f[1]+ci.f[3]);        uout[7].f[0] = xi3;
-
-    for(int k=1; k < dk; ++k)
-    {
-        v4sf save_next = in[8*k+7];
-        pffft_real_finalize_4x4(&save, &in[8*k+0], in + 8*k+1, e + k*6, out + k*8);
-        save = save_next;
-    }
+    auto *uout = reinterpret_cast<v4sf_union*>(out);
+    const float xr0{(cr[0]+cr[2]) + (cr[1]+cr[3])}; uout[0].f[0] = xr0;
+    const float xi0{(cr[0]+cr[2]) - (cr[1]+cr[3])}; uout[1].f[0] = xi0;
+    const float xr2{(cr[0]-cr[2])};                 uout[4].f[0] = xr2;
+    const float xi2{(cr[3]-cr[1])};                 uout[5].f[0] = xi2;
+    const float xr1{ ci[0] + s*(ci[1]-ci[3])};      uout[2].f[0] = xr1;
+    const float xi1{-ci[2] - s*(ci[1]+ci[3])};      uout[3].f[0] = xi1;
+    const float xr3{ ci[0] - s*(ci[1]-ci[3])};      uout[6].f[0] = xr3;
+    const float xi3{ ci[2] - s*(ci[1]+ci[3])};      uout[7].f[0] = xi3;
+
+    for(int k{1};k < dk;++k)
+        pffft_real_finalize_4x4(&in[8*k-1], &in[8*k+0], in + 8*k+1, e + k*6, out + k*8);
 }
 
 static ALWAYS_INLINE(void) pffft_real_preprocess_4x4(const v4sf *in, const v4sf *e, v4sf *out,
     int first)
 {
     v4sf r0=in[0], i0=in[1], r1=in[2], i1=in[3], r2=in[4], i2=in[5], r3=in[6], i3=in[7];
-    /*
-     * transformation for each column is:
+
+    /* transformation for each column is:
      *
      * [1   1   1   1   0   0   0   0]   [r0]
      * [1   0   0  -1   0  -1  -1   0]   [r1]
@@ -1741,22 +1726,22 @@ static ALWAYS_INLINE(void) pffft_real_preprocess_4x4(const v4sf *in, const v4sf
 static NEVER_INLINE(void) pffft_real_preprocess(int Ncvec, const v4sf *in, v4sf *out, const v4sf *e)
 {
     static constexpr float s = al::numbers::sqrt2_v<float>;
-    const int dk = Ncvec/SIMD_SZ; // number of 4x4 matrix blocks
-    /* fftpack order is f0r f1r f1i f2r f2i ... f(n-1)r f(n-1)i f(n)r */
 
-    v4sf_union Xr, Xi, *uout = reinterpret_cast<v4sf_union*>(out);
-    float cr0, ci0, cr1, ci1, cr2, ci2, cr3, ci3;
     assert(in != out);
-    for(int k=0; k < 4; ++k)
+    const int dk{Ncvec/SIMD_SZ}; // number of 4x4 matrix blocks
+    /* fftpack order is f0r f1r f1i f2r f2i ... f(n-1)r f(n-1)i f(n)r */
+
+    std::array<float,SIMD_SZ> Xr, Xi;
+    for(size_t k{0};k < 4;++k)
     {
-        Xr.f[k] = reinterpret_cast<const float*>(in)[8*k];
-        Xi.f[k] = reinterpret_cast<const float*>(in)[8*k+4];
+        /* TODO: Use _mm_cvtss_f32 or equivalent. */
+        Xr[k] = al::bit_cast<std::array<float,SIMD_SZ>>(in[4*k])[0];
+        Xi[k] = al::bit_cast<std::array<float,SIMD_SZ>>(in[4*k + 1])[0];
     }
 
     pffft_real_preprocess_4x4(in, e, out+1, 1); // will write only 6 values
 
-    /*
-     * [Xr0 Xr1 Xr2 Xr3 Xi0 Xi1 Xi2 Xi3]
+    /* [Xr0 Xr1 Xr2 Xr3 Xi0 Xi1 Xi2 Xi3]
      *
      * [cr0] [1   0   2   0   1   0   0   0]
      * [cr1] [1   0   0   0  -1   0  -2   0]
@@ -1767,57 +1752,60 @@ static NEVER_INLINE(void) pffft_real_preprocess(int Ncvec, const v4sf *in, v4sf
      * [ci2] [0   0   0   0   0  -2   0   2]
      * [ci3] [0  -s   0   s   0  -s   0  -s]
      */
-    for(int k=1; k < dk; ++k)
+    for(int k{1};k < dk;++k)
         pffft_real_preprocess_4x4(in+8*k, e + k*6, out-1+k*8, 0);
 
-    cr0=(Xr.f[0]+Xi.f[0]) + 2*Xr.f[2]; uout[0].f[0] = cr0;
-    cr1=(Xr.f[0]-Xi.f[0]) - 2*Xi.f[2]; uout[0].f[1] = cr1;
-    cr2=(Xr.f[0]+Xi.f[0]) - 2*Xr.f[2]; uout[0].f[2] = cr2;
-    cr3=(Xr.f[0]-Xi.f[0]) + 2*Xi.f[2]; uout[0].f[3] = cr3;
-    ci0= 2*(Xr.f[1]+Xr.f[3]);                       uout[2*Ncvec-1].f[0] = ci0;
-    ci1= s*(Xr.f[1]-Xr.f[3]) - s*(Xi.f[1]+Xi.f[3]); uout[2*Ncvec-1].f[1] = ci1;
-    ci2= 2*(Xi.f[3]-Xi.f[1]);                       uout[2*Ncvec-1].f[2] = ci2;
-    ci3=-s*(Xr.f[1]-Xr.f[3]) - s*(Xi.f[1]+Xi.f[3]); uout[2*Ncvec-1].f[3] = ci3;
+    /* TODO: Use _mm_set_ps or equivalent. */
+    auto *uout = reinterpret_cast<v4sf_union*>(out);
+    const float cr0{(Xr[0]+Xi[0]) + 2*Xr[2]}; uout[0].f[0] = cr0;
+    const float cr1{(Xr[0]-Xi[0]) - 2*Xi[2]}; uout[0].f[1] = cr1;
+    const float cr2{(Xr[0]+Xi[0]) - 2*Xr[2]}; uout[0].f[2] = cr2;
+    const float cr3{(Xr[0]-Xi[0]) + 2*Xi[2]}; uout[0].f[3] = cr3;
+    const float ci0{ 2*(Xr[1]+Xr[3])};                   uout[2*Ncvec-1].f[0] = ci0;
+    const float ci1{ s*(Xr[1]-Xr[3]) - s*(Xi[1]+Xi[3])}; uout[2*Ncvec-1].f[1] = ci1;
+    const float ci2{ 2*(Xi[3]-Xi[1])};                   uout[2*Ncvec-1].f[2] = ci2;
+    const float ci3{-s*(Xr[1]-Xr[3]) - s*(Xi[1]+Xi[3])}; uout[2*Ncvec-1].f[3] = ci3;
 }
 
 
 void pffft_transform_internal(PFFFT_Setup *setup, const float *finput, float *foutput,
     v4sf *scratch, pffft_direction_t direction, int ordered)
 {
-    const int Ncvec  = setup->Ncvec;
-    const int nf_odd = (setup->ifac[1] & 1);
-
-    // temporary buffer is allocated on the stack if the scratch pointer is NULL
     assert(scratch != nullptr);
+    assert(VALIGNED(finput) && VALIGNED(foutput) && VALIGNED(scratch));
 
-    const v4sf *vinput = reinterpret_cast<const v4sf*>(finput);
-    v4sf *voutput      = reinterpret_cast<v4sf*>(foutput);
-    v4sf *buff[2]      = { voutput, scratch };
-    int ib = (nf_odd ^ ordered ? 1 : 0);
+    const int Ncvec{setup->Ncvec};
+    const int nf_odd{setup->ifac[1] & 1};
 
-    assert(VALIGNED(finput) && VALIGNED(foutput));
+    auto *vinput = reinterpret_cast<const v4sf*>(finput);
+    auto *voutput = reinterpret_cast<v4sf*>(foutput);
+    assert(voutput != scratch);
 
-    //assert(finput != foutput);
+    v4sf *buff[2]{voutput, scratch};
+    int ib{(nf_odd ^ ordered) ? 1 : 0};
     if(direction == PFFFT_FORWARD)
     {
+        /* Swap the initial work buffer for forward FFTs, which helps avoid an
+         * extra copy for output.
+         */
         ib = !ib;
         if(setup->transform == PFFFT_REAL)
         {
-            ib = (rfftf1_ps(Ncvec*2, vinput, buff[ib], buff[!ib], setup->twiddle, &setup->ifac[0]) == buff[0] ? 0 : 1);
-            pffft_real_finalize(Ncvec, buff[ib], buff[!ib], reinterpret_cast<v4sf*>(setup->e));
+            ib = (rfftf1_ps(Ncvec*2, vinput, buff[ib], buff[!ib], setup->twiddle, setup->ifac) == buff[1]);
+            pffft_real_finalize(Ncvec, buff[ib], buff[!ib], setup->e);
         }
         else
         {
-            v4sf *tmp = buff[ib];
+            v4sf *tmp{buff[ib]};
             for(int k=0; k < Ncvec; ++k)
-            {
                 UNINTERLEAVE2(vinput[k*2], vinput[k*2+1], tmp[k*2], tmp[k*2+1]);
-            }
-            ib = (cfftf1_ps(Ncvec, buff[ib], buff[!ib], buff[ib], setup->twiddle, &setup->ifac[0], -1.0f) == buff[0] ? 0 : 1);
-            pffft_cplx_finalize(Ncvec, buff[ib], buff[!ib], reinterpret_cast<v4sf*>(setup->e));
+
+            ib = (cfftf1_ps(Ncvec, buff[ib], buff[!ib], buff[ib], setup->twiddle, setup->ifac, -1.0f) == buff[1]);
+            pffft_cplx_finalize(Ncvec, buff[ib], buff[!ib], setup->e);
         }
         if(ordered)
-            pffft_zreorder(setup, reinterpret_cast<float*>(buff[!ib]), reinterpret_cast<float*>(buff[ib]), PFFFT_FORWARD);
+            pffft_zreorder(setup, reinterpret_cast<float*>(buff[!ib]),
+                reinterpret_cast<float*>(buff[ib]), PFFFT_FORWARD);
         else
             ib = !ib;
     }
@@ -1828,23 +1816,22 @@ void pffft_transform_internal(PFFFT_Setup *setup, const float *finput, float *fo
 
         if(ordered)
         {
-            pffft_zreorder(setup, reinterpret_cast<const float*>(vinput), reinterpret_cast<float*>(buff[ib]), PFFFT_BACKWARD);
+            pffft_zreorder(setup, reinterpret_cast<const float*>(vinput),
+                reinterpret_cast<float*>(buff[ib]), PFFFT_BACKWARD);
             vinput = buff[ib];
             ib = !ib;
         }
         if(setup->transform == PFFFT_REAL)
         {
-            pffft_real_preprocess(Ncvec, vinput, buff[ib], reinterpret_cast<v4sf*>(setup->e));
-            ib = (rfftb1_ps(Ncvec*2, buff[ib], buff[0], buff[1], setup->twiddle, &setup->ifac[0]) == buff[0] ? 0 : 1);
+            pffft_real_preprocess(Ncvec, vinput, buff[ib], setup->e);
+            ib = (rfftb1_ps(Ncvec*2, buff[ib], buff[0], buff[1], setup->twiddle, setup->ifac) == buff[1]);
         }
         else
         {
-            pffft_cplx_preprocess(Ncvec, vinput, buff[ib], reinterpret_cast<v4sf*>(setup->e));
-            ib = (cfftf1_ps(Ncvec, buff[ib], buff[0], buff[1],  setup->twiddle, &setup->ifac[0], +1.0f) == buff[0] ? 0 : 1);
-            for(int k=0; k < Ncvec; ++k)
-            {
+            pffft_cplx_preprocess(Ncvec, vinput, buff[ib], setup->e);
+            ib = (cfftf1_ps(Ncvec, buff[ib], buff[0], buff[1],  setup->twiddle, setup->ifac, +1.0f) == buff[1]);
+            for(int k{0};k < Ncvec;++k)
                 INTERLEAVE2(buff[ib][k*2], buff[ib][k*2+1], buff[ib][k*2], buff[ib][k*2+1]);
-            }
         }
     }
 
@@ -1852,14 +1839,13 @@ void pffft_transform_internal(PFFFT_Setup *setup, const float *finput, float *fo
     {
         /* extra copy required -- this situation should only happen when finput == foutput */
         assert(finput==foutput);
-        for(int k=0; k < Ncvec; ++k)
+        for(int k{0};k < Ncvec;++k)
         {
-            v4sf a = buff[ib][2*k], b = buff[ib][2*k+1];
+            v4sf a{buff[ib][2*k]}, b{buff[ib][2*k+1]};
             voutput[2*k] = a; voutput[2*k+1] = b;
         }
         ib = !ib;
     }
-    assert(buff[ib] == voutput);
 }
 
 void pffft_zconvolve_accumulate(PFFFT_Setup *s, const float *a, const float *b, float *ab,
@@ -1867,10 +1853,10 @@ void pffft_zconvolve_accumulate(PFFFT_Setup *s, const float *a, const float *b,
 {
     assert(VALIGNED(a) && VALIGNED(b) && VALIGNED(ab));
 
-    const int Ncvec = s->Ncvec;
-    const v4sf *RESTRICT va = reinterpret_cast<const v4sf*>(a);
-    const v4sf *RESTRICT vb = reinterpret_cast<const v4sf*>(b);
-    v4sf *RESTRICT vab = reinterpret_cast<v4sf*>(ab);
+    const int Ncvec{s->Ncvec};
+    const v4sf *RESTRICT va{reinterpret_cast<const v4sf*>(a)};
+    const v4sf *RESTRICT vb{reinterpret_cast<const v4sf*>(b)};
+    v4sf *RESTRICT vab{reinterpret_cast<v4sf*>(ab)};
 
 #ifdef __arm__
     __builtin_prefetch(va);
@@ -1891,18 +1877,19 @@ void pffft_zconvolve_accumulate(PFFFT_Setup *s, const float *a, const float *b,
 #endif
 
 #ifndef ZCONVOLVE_USING_INLINE_ASM
-    const v4sf vscal = LD_PS1(scaling);
+    const v4sf vscal{LD_PS1(scaling)};
 #endif
-    float ar1 = reinterpret_cast<const v4sf_union*>(va)[0].f[0];
-    float ai1 = reinterpret_cast<const v4sf_union*>(va)[1].f[0];
-    float br1 = reinterpret_cast<const v4sf_union*>(vb)[0].f[0];
-    float bi1 = reinterpret_cast<const v4sf_union*>(vb)[1].f[0];
-    float abr1 = reinterpret_cast<v4sf_union*>(vab)[0].f[0];
-    float abi1 = reinterpret_cast<v4sf_union*>(vab)[1].f[0];
+    /* TODO: Use _mm_cvtss_f32 or equivalent. */
+    float ar1{reinterpret_cast<const v4sf_union*>(va)[0].f[0]};
+    float ai1{reinterpret_cast<const v4sf_union*>(va)[1].f[0]};
+    float br1{reinterpret_cast<const v4sf_union*>(vb)[0].f[0]};
+    float bi1{reinterpret_cast<const v4sf_union*>(vb)[1].f[0]};
+    float abr1{reinterpret_cast<v4sf_union*>(vab)[0].f[0]};
+    float abi1{reinterpret_cast<v4sf_union*>(vab)[1].f[0]};
 
 #ifdef ZCONVOLVE_USING_INLINE_ASM // inline asm version, unfortunately miscompiled by clang 3.2, at least on ubuntu.. so this will be restricted to gcc
-    const float *a_ = a, *b_ = b; float *ab_ = ab;
-    int N = Ncvec;
+    const float *a_{a}, *b_{b}; float *ab_{ab};
+    int N{Ncvec};
     asm volatile("mov         r8, %2                  \n"
                 "vdup.f32    q15, %4                 \n"
                 "1:                                  \n"
@@ -1941,9 +1928,8 @@ void pffft_zconvolve_accumulate(PFFFT_Setup *s, const float *a, const float *b,
 
     for(int i=0; i < Ncvec; i += 2)
     {
-        v4sf ar4, ai4, br4, bi4;
-        ar4 = va[2*i+0]; ai4 = va[2*i+1];
-        br4 = vb[2*i+0]; bi4 = vb[2*i+1];
+        v4sf ar4{va[2*i+0]}, ai4{va[2*i+1]};
+        v4sf br4{vb[2*i+0]}, bi4{vb[2*i+1]};
         VCPLXMUL(ar4, ai4, br4, bi4);
         vab[2*i+0] = VMADD(ar4, vscal, vab[2*i+0]);
         vab[2*i+1] = VMADD(ai4, vscal, vab[2*i+1]);
@@ -2000,22 +1986,24 @@ void pffft_zreorder_nosimd(PFFFT_Setup *setup, const float *in, float *out,
 void pffft_transform_internal_nosimd(PFFFT_Setup *setup, const float *input, float *output,
     float *scratch, pffft_direction_t direction, int ordered)
 {
-    const int Ncvec  = setup->Ncvec;
-    const int nf_odd = (setup->ifac[1] & 1);
+    const int Ncvec{setup->Ncvec};
+    const int nf_odd{setup->ifac[1] & 1};
 
     assert(scratch != nullptr);
 
+    /* z-domain data for complex transforms is already ordered without SIMD. */
     if(setup->transform == PFFFT_COMPLEX)
-        ordered = 0; // it is always ordered.
-    int ib = (nf_odd ^ ordered ? 1 : 0);
-    float *buff[2] = { output, scratch };
+        ordered = 0;
+
+    float *buff[2]{output, scratch};
+    int ib{(nf_odd ^ ordered) ? 1 : 0};
 
     if(direction == PFFFT_FORWARD)
     {
         if(setup->transform == PFFFT_REAL)
-            ib = (rfftf1_ps(Ncvec*2, input, buff[ib], buff[!ib], setup->twiddle, &setup->ifac[0]) == buff[0] ? 0 : 1);
+            ib = (rfftf1_ps(Ncvec*2, input, buff[ib], buff[!ib], setup->twiddle, setup->ifac) == buff[1]);
         else
-            ib = (cfftf1_ps(Ncvec, input, buff[ib], buff[!ib], setup->twiddle, &setup->ifac[0], -1.0f) == buff[0] ? 0 : 1);
+            ib = (cfftf1_ps(Ncvec, input, buff[ib], buff[!ib], setup->twiddle, setup->ifac, -1.0f) == buff[1]);
         if(ordered)
         {
             pffft_zreorder(setup, buff[ib], buff[!ib], PFFFT_FORWARD);
@@ -2034,9 +2022,9 @@ void pffft_transform_internal_nosimd(PFFFT_Setup *setup, const float *input, flo
             ib = !ib;
         }
         if(setup->transform == PFFFT_REAL)
-            ib = (rfftb1_ps(Ncvec*2, input, buff[ib], buff[!ib],  setup->twiddle, &setup->ifac[0]) == buff[0] ? 0 : 1);
+            ib = (rfftb1_ps(Ncvec*2, input, buff[ib], buff[!ib],  setup->twiddle, setup->ifac) == buff[1]);
         else
-            ib = (cfftf1_ps(Ncvec, input, buff[ib], buff[!ib], setup->twiddle, &setup->ifac[0], +1.0f) == buff[0] ? 0 : 1);
+            ib = (cfftf1_ps(Ncvec, input, buff[ib], buff[!ib], setup->twiddle, setup->ifac, +1.0f) == buff[1]);
     }
     if(buff[ib] != output)
     {
@@ -2049,7 +2037,6 @@ void pffft_transform_internal_nosimd(PFFFT_Setup *setup, const float *input, flo
         }
         ib = !ib;
     }
-    assert(buff[ib] == output);
 }
 
 #define pffft_zconvolve_accumulate_nosimd pffft_zconvolve_accumulate