/*
This file is part of the ducc FFT library.

Copyright (C) 2010-2023 Max-Planck-Society
Copyright (C) 2019 Peter Bell

For the odd-sized DCT-IV transforms:
  Copyright (C) 2003, 2007-14 Matteo Frigo
  Copyright (C) 2003, 2007-14 Massachusetts Institute of Technology

Authors: Martin Reinecke, Peter Bell
*/

/* SPDX-License-Identifier: BSD-3-Clause OR GPL-2.0-or-later */

/*
All rights reserved.

Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:

* Redistributions of source code must retain the above copyright notice, this
  list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice, this
  list of conditions and the following disclaimer in the documentation and/or
  other materials provided with the distribution.
* Neither the name of the copyright holder nor the names of its contributors may
  be used to endorse or promote products derived from this software without
  specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

/*
 *  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
 */

#ifndef DUCC0_FFT_H
#define DUCC0_FFT_H

#include <cmath>
#include <cstddef>
#include <cstdlib>
#include <numeric>
#include <stdexcept>
#include <memory>
#include <vector>
#include <complex>
#include <algorithm>
#include "ducc0/infra/useful_macros.h"
#include "ducc0/infra/error_handling.h"
#include "ducc0/infra/threading.h"
#include "ducc0/infra/misc_utils.h"
#include "ducc0/infra/simd.h"
#include "ducc0/infra/mav.h"
#include "ducc0/infra/aligned_array.h"
#include "ducc0/math/cmplx.h"
#include "ducc0/math/unity_roots.h"
#include "ducc0/fft/fft1d.h"

/** \file fft.h
 *  Implementation of multi-dimensional Fast Fourier and related transforms
 *  \copyright Copyright (C) 2010-2021 Max-Planck-Society
 *  \copyright Copyright (C) 2019 Peter Bell
 *  \copyright
 *  \copyright For the odd-sized DCT-IV transforms:
 *  \copyright   Copyright (C) 2003, 2007-14 Matteo Frigo
 *  \copyright   Copyright (C) 2003, 2007-14 Massachusetts Institute of Technology
 *
 * \authors Martin Reinecke, Peter Bell
 */

namespace ducc0 {

namespace detail_fft {

// the next line is necessary to address some sloppy name choices in hipSYCL
using std::min, std::max;

template<typename T> constexpr inline size_t fft_simdlen
  = min<size_t>(8, native_simd<T>::size());
template<> constexpr inline size_t fft_simdlen<double>
  = min<size_t>(4, native_simd<double>::size());
template<> constexpr inline size_t fft_simdlen<float>
  = min<size_t>(8, native_simd<float>::size());
template<typename T> using fft_simd = typename simd_select<T,fft_simdlen<T>>::type;
template<typename T> constexpr inline bool fft_simd_exists = (fft_simdlen<T> > 1);

using shape_t=fmav_info::shape_t;
using stride_t=fmav_info::stride_t;

constexpr bool FORWARD  = true,
               BACKWARD = false;

struct util // hack to avoid duplicate symbols
  {
  static void sanity_check_axes(size_t ndim, const shape_t &axes)
    {
    if (ndim==1)
      {
      if ((axes.size()!=1) || (axes[0]!=0))
        throw std::invalid_argument("bad axes");
      return;
      }
    shape_t tmp(ndim,0);
    if (axes.empty()) throw std::invalid_argument("no axes specified");
    for (auto ax : axes)
      {
      if (ax>=ndim) throw std::invalid_argument("bad axis number");
      if (++tmp[ax]>1) throw std::invalid_argument("axis specified repeatedly");
      }
    }

  DUCC0_NOINLINE static void sanity_check_onetype(const fmav_info &a1,
    const fmav_info &a2, bool inplace, const shape_t &axes)
    {
    sanity_check_axes(a1.ndim(), axes);
    MR_assert(a1.conformable(a2), "array sizes are not conformable");
    if (inplace) MR_assert(a1.stride()==a2.stride(), "stride mismatch");
    }
  DUCC0_NOINLINE static void sanity_check_cr(const fmav_info &ac,
    const fmav_info &ar, const shape_t &axes)
    {
    sanity_check_axes(ac.ndim(), axes);
    MR_assert(ac.ndim()==ar.ndim(), "dimension mismatch");
    for (size_t i=0; i<ac.ndim(); ++i)
      MR_assert(ac.shape(i) == ((i==axes.back()) ? (ar.shape(i)/2+1) : ar.shape(i)),
        "axis length mismatch");
    }
  DUCC0_NOINLINE static void sanity_check_cr(const fmav_info &ac,
    const fmav_info &ar, const size_t axis)
    {
    if (axis>=ac.ndim()) throw std::invalid_argument("bad axis number");
    MR_assert(ac.ndim()==ar.ndim(), "dimension mismatch");
    for (size_t i=0; i<ac.ndim(); ++i)
      MR_assert(ac.shape(i) == ((i==axis) ? (ar.shape(i)/2+1) : ar.shape(i)),
        "axis length mismatch");
    }

  static size_t thread_count (size_t nthreads, const fmav_info &info,
    size_t axis, size_t vlen)
    {
    if (nthreads==1) return 1;
    size_t size = info.size();
    size_t parallel = size / (info.shape(axis) * vlen);
    if (info.shape(axis) < 1000)
      parallel /= 4;
    size_t max_threads = ducc0::adjust_nthreads(nthreads);
    return std::max(size_t(1), std::min(parallel, max_threads));
    }
  };


//
// sine/cosine transforms
//

template<typename T0> class T_dct1
  {
  private:
    pocketfft_r<T0> fftplan;

  public:
    DUCC0_NOINLINE T_dct1(size_t length, bool /*vectorize*/=false)
      : fftplan(2*(length-1)) {}

    template<typename T> DUCC0_NOINLINE T *exec(T c[], T buf[], T0 fct, bool ortho,
      int /*type*/, bool /*cosine*/, size_t nthreads=1) const
      {
      constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
      size_t N=fftplan.length(), n=N/2+1;
      if (ortho)
        { c[0]*=sqrt2; c[n-1]*=sqrt2; }
      auto tmp=&buf[0];
      tmp[0] = c[0];
      for (size_t i=1; i<n; ++i)
        tmp[i] = tmp[N-i] = c[i];
      auto res = fftplan.exec(tmp, &buf[N], fct, true, nthreads);
      c[0] = res[0];
      for (size_t i=1; i<n; ++i)
        c[i] = res[2*i-1];
      if (ortho)
        { c[0]*=sqrt2*T0(0.5); c[n-1]*=sqrt2*T0(0.5); }
      return c;
      }
    template<typename T> DUCC0_NOINLINE void exec_copyback(T c[], T buf[], T0 fct, bool ortho,
      int /*type*/, bool /*cosine*/, size_t nthreads=1) const
      {
      exec(c, buf, fct, ortho, 1, true, nthreads);
      }
    template<typename T> DUCC0_NOINLINE void exec(T c[], T0 fct, bool ortho,
      int /*type*/, bool /*cosine*/, size_t nthreads=1) const
      {
      quick_array<T> buf(bufsize());
      exec_copyback(c, buf.data(), fct, ortho, 1, true, nthreads);
      }

    size_t length() const { return fftplan.length()/2+1; }
    size_t bufsize() const { return fftplan.length()+fftplan.bufsize(); }
  };

template<typename T0> class T_dst1
  {
  private:
    pocketfft_r<T0> fftplan;

  public:
    DUCC0_NOINLINE T_dst1(size_t length, bool /*vectorize*/=false)
      : fftplan(2*(length+1)) {}

    template<typename T> DUCC0_NOINLINE T *exec(T c[], T buf[], T0 fct,
      bool /*ortho*/, int /*type*/, bool /*cosine*/, size_t nthreads=1) const
      {
      size_t N=fftplan.length(), n=N/2-1;
      auto tmp = &buf[0];
      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]; }
      auto res = fftplan.exec(tmp, buf+N, fct, true, nthreads);
      for (size_t i=0; i<n; ++i)
        c[i] = -res[2*i+2];
      return c;
      }
    template<typename T> DUCC0_NOINLINE void exec_copyback(T c[], T buf[], T0 fct,
      bool /*ortho*/, int /*type*/, bool /*cosine*/, size_t nthreads=1) const
      {
      exec(c, buf, fct, true, 1, false, nthreads);
      }
    template<typename T> DUCC0_NOINLINE void exec(T c[], T0 fct,
      bool /*ortho*/, int /*type*/, bool /*cosine*/, size_t nthreads) const
      {
      quick_array<T> buf(bufsize());
      exec_copyback(c, buf.data(), fct, true, 1, false, nthreads);
      }

    size_t length() const { return fftplan.length()/2-1; }
    size_t bufsize() const { return fftplan.length()+fftplan.bufsize(); }
  };

template<typename T0> class T_dcst23
  {
  private:
    pocketfft_r<T0> fftplan;
    std::vector<T0> twiddle;

  public:
    DUCC0_NOINLINE T_dcst23(size_t length, bool /*vectorize*/=false)
      : fftplan(length), twiddle(length)
      {
      UnityRoots<T0,Cmplx<T0>> tw(4*length);
      for (size_t i=0; i<length; ++i)
        twiddle[i] = tw[i+1].r;
      }

    template<typename T> DUCC0_NOINLINE T *exec(T c[], T buf[], T0 fct, bool ortho,
      int type, bool cosine, size_t nthreads=1) const
      {
      constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
      size_t N=length();
      size_t NS2 = (N+1)/2;
      if (type==2)
        {
        c[0] *= 2;
        if ((N&1)==0) c[N-1]*=2;
        if (cosine)
          for (size_t k=1; k<N-1; k+=2)
            MPINPLACE(c[k+1], c[k]);
        else
          for (size_t k=1; k<N-1; k+=2)
            PMINPLACE(c[k+1], c[k]);
        if ((!cosine) && ((N&1)==0))
          c[N-1] *= -1;
        auto res = fftplan.exec(c, buf, fct, false, nthreads);
        c[0] = res[0];
        for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
          {
          T t1 = twiddle[k-1]*res[kc]+twiddle[kc-1]*res[k];
          T t2 = twiddle[k-1]*res[k]-twiddle[kc-1]*res[kc];
          c[k] = T0(0.5)*(t1+t2); c[kc]=T0(0.5)*(t1-t2);
          }
        if ((N&1)==0)
          c[NS2] = res[NS2]*twiddle[NS2-1];
        if (!cosine)  // swap order completely
          for (size_t k=0, kc=N-1; k<kc; ++k, --kc)
            std::swap(c[k], c[kc]);
        if (ortho) c[0]*=sqrt2*T0(0.5);
        }
      else
        {
        if (ortho) c[0]*=sqrt2;
        if (!cosine)  // swap order completely
          for (size_t k=0, kc=N-1; k<NS2; ++k, --kc)
            std::swap(c[k], c[kc]);
        for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
          {
          T t1=c[k]+c[kc], t2=c[k]-c[kc];
          c[k] = twiddle[k-1]*t2+twiddle[kc-1]*t1;
          c[kc]= twiddle[k-1]*t1-twiddle[kc-1]*t2;
          }
        if ((N&1)==0)
          c[NS2] *= 2*twiddle[NS2-1];
        auto res = fftplan.exec(c, buf, fct, true, nthreads);
        if (res != c) // FIXME: not yet optimal
          copy_n(res, N, c);
        if ((!cosine) && ((N&1)==0))
          c[N-1] *= -1;
        if (cosine)
          for (size_t k=1; k<N-1; k+=2)
            MPINPLACE(c[k], c[k+1]);
        else
          for (size_t k=1; k<N-1; k+=2)
            PMINPLACE(c[k+1], c[k]);
        }
      return c;
      }
    template<typename T> DUCC0_NOINLINE void exec_copyback(T c[], T buf[], T0 fct,
      bool ortho, int type, bool cosine, size_t nthreads=1) const
      {
      exec(c, buf, fct, ortho, type, cosine, nthreads);
      }
    template<typename T> DUCC0_NOINLINE void exec(T c[], T0 fct, bool ortho,
      int type, bool cosine, size_t nthreads=1) const
      {
      quick_array<T> buf(bufsize());
      exec(c, &buf[0], fct, ortho, type, cosine, nthreads);
      }

    size_t length() const { return fftplan.length(); }
    size_t bufsize() const { return fftplan.bufsize(); }
  };

template<typename T0> class T_dcst4
  {
  private:
    size_t N;
    std::unique_ptr<pocketfft_c<T0>> fft;
    std::unique_ptr<pocketfft_r<T0>> rfft;
    quick_array<Cmplx<T0>> C2;
    size_t bufsz;

  public:
    DUCC0_NOINLINE T_dcst4(size_t length, bool /*vectorize*/=false)
      : N(length),
        fft((N&1) ? nullptr : make_unique<pocketfft_c<T0>>(N/2)),
        rfft((N&1)? make_unique<pocketfft_r<T0>>(N) : nullptr),
        C2((N&1) ? 0 : N/2),
        bufsz((N&1) ? (N+rfft->bufsize()) : (N+2*fft->bufsize()))
      {
      if ((N&1)==0)
        {
        UnityRoots<T0,Cmplx<T0>> tw(16*N);
        for (size_t i=0; i<N/2; ++i)
          C2[i] = tw[8*i+1].conj();
        }
      }

    template<typename T> DUCC0_NOINLINE T *exec(T c[], T buf[], T0 fct,
      bool /*ortho*/, int /*type*/, bool cosine, size_t nthreads) const
      {
      size_t n2 = N/2;
      if (!cosine)
        for (size_t k=0, kc=N-1; k<n2; ++k, --kc)
          std::swap(c[k], c[kc]);
      if (N&1)
        {
        // The following code is derived from the FFTW3 function apply_re11()
        // and is released under the 3-clause BSD license with friendly
        // permission of Matteo Frigo and Steven G. Johnson.

        auto y = buf;
        {
        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];
        for (; m<3*N; ++i, m+=4)
          y[i] = -c[m-2*N];
        for (; m<4*N; ++i, m+=4)
          y[i] = c[4*N-m-1];
        for (; i<N; ++i, m+=4)
          y[i] = c[m-4*N];
        }

        auto res = rfft->exec(y, y+N, fct, true, nthreads);
        {
        auto SGN = [](size_t i)
           {
           constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
           return (i&2) ? -sqrt2 : sqrt2;
           };
        c[n2] = res[0]*SGN(n2+1);
        size_t i=0, i1=1, k=1;
        for (; k<n2; ++i, ++i1, k+=2)
          {
          c[i    ] = res[2*k-1]*SGN(i1)     + res[2*k  ]*SGN(i);
          c[N -i1] = res[2*k-1]*SGN(N -i)   - res[2*k  ]*SGN(N -i1);
          c[n2-i1] = res[2*k+1]*SGN(n2-i)   - res[2*k+2]*SGN(n2-i1);
          c[n2+i1] = res[2*k+1]*SGN(n2+i+2) + res[2*k+2]*SGN(n2+i1);
          }
        if (k == n2)
          {
          c[i   ] = res[2*k-1]*SGN(i+1) + res[2*k]*SGN(i);
          c[N-i1] = res[2*k-1]*SGN(i+2) + res[2*k]*SGN(i1);
          }
        }

        // 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/
        auto y2 = reinterpret_cast<Cmplx<T> *>(buf);
        for(size_t i=0; i<n2; ++i)
          {
          y2[i].Set(c[2*i],c[N-1-2*i]);
          y2[i] *= C2[i];
          }

        auto res = fft->exec(y2, y2+N/2, fct, true, nthreads);
        for(size_t i=0, ic=n2-1; i<n2; ++i, --ic)
          {
          c[2*i  ] = T0( 2)*(res[i ].r*C2[i ].r-res[i ].i*C2[i ].i);
          c[2*i+1] = T0(-2)*(res[ic].i*C2[ic].r+res[ic].r*C2[ic].i);
          }
        }
      if (!cosine)
        for (size_t k=1; k<N; k+=2)
          c[k] = -c[k];
      return c;
      }
    template<typename T> DUCC0_NOINLINE void exec_copyback(T c[], T buf[], T0 fct,
      bool /*ortho*/, int /*type*/, bool cosine, size_t nthreads=1) const
      {
      exec(c, buf, fct, true, 4, cosine, nthreads);
      }
    template<typename T> DUCC0_NOINLINE void exec(T c[], T0 fct,
      bool /*ortho*/, int /*type*/, bool cosine, size_t nthreads=1) const
      {
      quick_array<T> buf(bufsize());
      exec(c, &buf[0], fct, true, 4, cosine, nthreads);
      }

    size_t length() const { return N; }
    size_t bufsize() const { return bufsz; }
  };


//
// multi-D infrastructure
//

template<typename T> std::shared_ptr<T> get_plan(size_t length, bool vectorize=false)
  {
#ifdef DUCC0_NO_FFT_CACHE
  return std::make_shared<T>(length, vectorize);
#else
  constexpr size_t nmax=10;
  struct entry { size_t n; bool vectorize; std::shared_ptr<T> ptr; };
  static std::array<entry, nmax> cache{{{0,0,nullptr}}};
  static std::array<size_t, nmax> last_access{{0}};
  static size_t access_counter = 0;
  static Mutex mut;

  auto find_in_cache = [&]() -> std::shared_ptr<T>
    {
    for (size_t i=0; i<nmax; ++i)
      if (cache[i].ptr && (cache[i].n==length) && (cache[i].vectorize==vectorize))
        {
        // no need to update if this is already the most recent entry
        if (last_access[i]!=access_counter)
          {
          last_access[i] = ++access_counter;
          // Guard against overflow
          if (access_counter == 0)
            last_access.fill(0);
          }
        return cache[i].ptr;
        }

    return nullptr;
    };

  {
  LockGuard lock(mut);

  auto p = find_in_cache();
  if (p) return p;
  }
  auto plan = std::make_shared<T>(length, vectorize);
  {
  LockGuard lock(mut);

  auto p = find_in_cache();
  if (p) return p;

  size_t lru = 0;
  for (size_t i=1; i<nmax; ++i)
    if (last_access[i] < last_access[lru])
      lru = i;

  cache[lru] = {length,vectorize, plan};
  last_access[lru] = ++access_counter;
  }
  return plan;
#endif
  }

template<size_t N> class multi_iter
  {
  private:
    shape_t shp, pos;
    stride_t str_i, str_o;
    size_t cshp_i, cshp_o, rem;
    ptrdiff_t cstr_i, cstr_o, sstr_i, sstr_o, p_ii, p_i[N], p_oi, p_o[N];
    bool uni_i, uni_o;

    void advance_i()
      {
      for (size_t i=0; i<pos.size(); ++i)
        {
        p_ii += str_i[i];
        p_oi += str_o[i];
        if (++pos[i] < shp[i])
          return;
        pos[i] = 0;
        p_ii -= ptrdiff_t(shp[i])*str_i[i];
        p_oi -= ptrdiff_t(shp[i])*str_o[i];
        }
      }

  public:
    multi_iter(const fmav_info &iarr, const fmav_info &oarr, size_t idim,
      size_t nshares, size_t myshare)
      : rem(iarr.size()/iarr.shape(idim)), sstr_i(0), sstr_o(0), p_ii(0), p_oi(0)
      {
      MR_assert(oarr.ndim()==iarr.ndim(), "dimension mismatch");
      MR_assert(iarr.ndim()>=1, "not enough dimensions");
      // Sort the extraneous dimensions in order of ascending output stride;
      // this should improve overall cache re-use and avoid clashes between
      // threads as much as possible.
      shape_t idx(iarr.ndim());
      std::iota(idx.begin(), idx.end(), 0);
      sort(idx.begin(), idx.end(),
        [&oarr](size_t i1, size_t i2) {return oarr.stride(i1) < oarr.stride(i2);});
      for (auto i: idx)
        if (i!=idim)
          {
          pos.push_back(0);
          MR_assert(iarr.shape(i)==oarr.shape(i), "shape mismatch");
          shp.push_back(iarr.shape(i));
          str_i.push_back(iarr.stride(i));
          str_o.push_back(oarr.stride(i));
          }
      MR_assert(idim<iarr.ndim(), "bad active dimension");
      cstr_i = iarr.stride(idim);
      cstr_o = oarr.stride(idim);
      cshp_i = iarr.shape(idim);
      cshp_o = oarr.shape(idim);

// collapse unneeded dimensions
      bool done = false;
      while(!done)
        {
        done=true;
        for (size_t i=1; i<shp.size(); ++i)
          if ((str_i[i] == str_i[i-1]*ptrdiff_t(shp[i-1]))
           && (str_o[i] == str_o[i-1]*ptrdiff_t(shp[i-1])))
            {
            shp[i-1] *= shp[i];
            str_i.erase(str_i.begin()+ptrdiff_t(i));
            str_o.erase(str_o.begin()+ptrdiff_t(i));
            shp.erase(shp.begin()+ptrdiff_t(i));
            pos.pop_back();
            done=false;
            }
        }
      if (pos.size()>0)
        {
        sstr_i = str_i[0];
        sstr_o = str_o[0];
        }

      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");
      auto [lo, hi] = calcShare(nshares, myshare, rem);
      size_t todo = hi-lo;

      size_t chunk = rem;
      for (size_t i2=0, i=pos.size()-1; i2<pos.size(); ++i2,--i)
        {
        chunk /= shp[i];
        size_t n_advance = lo/chunk;
        pos[i] += n_advance;
        p_ii += ptrdiff_t(n_advance)*str_i[i];
        p_oi += ptrdiff_t(n_advance)*str_o[i];
        lo -= n_advance*chunk;
        }
      MR_assert(lo==0, "must not happen");
      rem = todo;
      }
    void advance(size_t n)
      {
      if (rem<n) throw std::runtime_error("underrun");
      for (size_t i=0; i<n; ++i)
        {
        p_i[i] = p_ii;
        p_o[i] = p_oi;
        advance_i();
        }
      uni_i = uni_o = true;
      for (size_t i=1; i<n; ++i)
        {
        uni_i = uni_i && (p_i[i]-p_i[i-1] == sstr_i);
        uni_o = uni_o && (p_o[i]-p_o[i-1] == sstr_o);
        }
      rem -= n;
      }
    ptrdiff_t iofs(size_t i) const { return p_i[0] + ptrdiff_t(i)*cstr_i; }
    ptrdiff_t iofs(size_t j, size_t i) const { return p_i[j] + ptrdiff_t(i)*cstr_i; }
    ptrdiff_t iofs_uni(size_t j, size_t i) const { return p_i[0] + ptrdiff_t(j)*sstr_i + ptrdiff_t(i)*cstr_i; }
    ptrdiff_t oofs(size_t i) const { return p_o[0] + ptrdiff_t(i)*cstr_o; }
    ptrdiff_t oofs(size_t j, size_t i) const { return p_o[j] + ptrdiff_t(i)*cstr_o; }
    ptrdiff_t oofs_uni(size_t j, size_t i) const { return p_o[0] + ptrdiff_t(j)*sstr_o + ptrdiff_t(i)*cstr_o; }
    bool uniform_i() const { return uni_i; }
    ptrdiff_t unistride_i() const { return sstr_i; }
    bool uniform_o() const { return uni_o; }
    ptrdiff_t unistride_o() const { return sstr_o; }
    size_t length_in() const { return cshp_i; }
    size_t length_out() const { return cshp_o; }
    ptrdiff_t stride_in() const { return cstr_i; }
    ptrdiff_t stride_out() const { return cstr_o; }
    size_t remaining() const { return rem; }
    bool critical_stride_trans(size_t tsz) const
      {
      return ((abs<ptrdiff_t>(stride_in() *tsz)&4095)==0)
          || ((abs<ptrdiff_t>(stride_out()*tsz)&4095)==0);
      }
    bool critical_stride_other(size_t tsz) const
      {
      if (unistride_i()==0) return false;  // it's just one transform
      return ((abs<ptrdiff_t>(unistride_i()*tsz)&4095)==0)
          || ((abs<ptrdiff_t>(unistride_o()*tsz)&4095)==0);
      }
  };

template<typename T, typename T0> class TmpStorage
  {
  private:
    aligned_array<T> d;
    size_t dofs, dstride;

  public:
    TmpStorage(size_t n_trafo, size_t bufsize_data, size_t bufsize_trafo,
               size_t n_simultaneous, bool inplace)
      {
      if (inplace)
        {
        d.resize(bufsize_trafo);
        return;
        }
      constexpr auto vlen = fft_simdlen<T0>;
      // FIXME: when switching to C++20, use bit_floor(othersize)
      size_t buffct = std::min(vlen, n_trafo);
      size_t datafct = std::min(vlen, n_trafo);
      if (n_trafo>=n_simultaneous*vlen) datafct = n_simultaneous*vlen;
      dstride = bufsize_data;
      // critical stride avoidance
      if ((dstride&256)==0) dstride+=3;
      d.resize(buffct*(bufsize_trafo+17) + datafct*dstride);
      dofs = bufsize_trafo + 17;
      }

    template<typename T2> T2 *transformBuf()
      { return reinterpret_cast<T2 *>(d.data()); }
    template<typename T2> T2 *dataBuf()
      { return reinterpret_cast<T2 *>(d.data()) + dofs; }
    size_t data_stride() const
      { return dstride; }
  };

template<typename T2, typename T, typename T0> class TmpStorage2
  {
  private:
    TmpStorage<T, T0> &stg;

  public:
    using datatype = T2;
    TmpStorage2(TmpStorage<T,T0> &stg_): stg(stg_) {}

    T2 *transformBuf() { return stg.template transformBuf<T2>(); }
    T2 *dataBuf() { return stg.template dataBuf<T2>(); }
    size_t data_stride() const { return stg.data_stride(); }
  };

// Yes, this looks strange. But this is currently the only way I found to
// stop compilers from vectorizing the copying loops and messing up the ordering
// of the memory accesses, which is really important here.
template <typename Titer, typename Ts> DUCC0_NOINLINE void copy_inputx2(const Titer &it,
  const cfmav<Cmplx<Ts>> &src, Ts *DUCC0_RESTRICT dst, size_t vlen)
  {
  for (size_t i=0; i<it.length_in(); ++i)
    for (size_t j=0; j<vlen; ++j)
      {
      dst[2*i*vlen+j     ] = src.raw(it.iofs(j,i)).r;
      dst[2*i*vlen+j+vlen] = src.raw(it.iofs(j,i)).i;
      }
  }
template <typename Titer, typename Ts> DUCC0_NOINLINE void copy_inputx(const Titer &it,
  const cfmav<Cmplx<Ts>> &src, Ts *DUCC0_RESTRICT dst, size_t vlen)
  {
  if (it.stride_in()==1)
    return copy_inputx2(it, src, dst, vlen);
  for (size_t i=0; i<it.length_in(); ++i)
    for (size_t j=0; j<vlen; ++j)
      {
      dst[2*i*vlen+j     ] = src.raw(it.iofs(j,i)).r;
      dst[2*i*vlen+j+vlen] = src.raw(it.iofs(j,i)).i;
      }
  }
template <typename Tsimd, typename Titer> DUCC0_NOINLINE void copy_input(const Titer &it,
  const cfmav<Cmplx<typename Tsimd::value_type>> &src, Cmplx<Tsimd> *DUCC0_RESTRICT dst)
  {
  constexpr auto vlen=Tsimd::size();
  copy_inputx(it, src, reinterpret_cast<typename Tsimd::value_type *>(dst),vlen);
  }

template <typename Tsimd, typename Titer> DUCC0_NOINLINE void copy_input(const Titer &it,
  const cfmav<typename Tsimd::value_type> &src, Tsimd *DUCC0_RESTRICT dst)
  {
  constexpr auto vlen=Tsimd::size();
  for (size_t i=0; i<it.length_in(); ++i)
    for (size_t j=0; j<vlen; ++j)
      dst[i][j] = src.raw(it.iofs(j,i));
  }

template <typename Titer, typename T> DUCC0_NOINLINE void copy_input(const Titer &it,
  const cfmav<T> &src, T *DUCC0_RESTRICT dst)
  {
  if (dst == &src.raw(it.iofs(0))) return;  // in-place
  for (size_t i=0; i<it.length_in(); ++i)
    dst[i] = src.raw(it.iofs(i));
  }

template<typename Titer,typename Ts> DUCC0_NOINLINE void copy_outputx2(const Titer &it,
  const Ts *DUCC0_RESTRICT src, vfmav<Cmplx<Ts>> &dst, size_t vlen)
  {
  Cmplx<Ts> * DUCC0_RESTRICT ptr = dst.data();
  for (size_t i=0; i<it.length_out(); ++i)
    for (size_t j=0; j<vlen; ++j)
        ptr[it.oofs(j,i)].Set(src[i*2*vlen+j],src[i*2*vlen+j+vlen]);
  }
template<typename Titer,typename Ts> DUCC0_NOINLINE void copy_outputx(const Titer &it,
  const Ts *DUCC0_RESTRICT src, vfmav<Cmplx<Ts>> &dst, size_t vlen)
  {
  if (it.stride_out()==1)
    return copy_outputx2(it,src,dst,vlen);
  Cmplx<Ts> * DUCC0_RESTRICT ptr = dst.data();
  for (size_t i=0; i<it.length_out(); ++i)
    for (size_t j=0; j<vlen; ++j)
        ptr[it.oofs(j,i)].Set(src[i*2*vlen+j],src[i*2*vlen+j+vlen]);
  }
template<typename Tsimd, typename Titer> DUCC0_NOINLINE void copy_output(const Titer &it,
  const Cmplx<Tsimd> *DUCC0_RESTRICT src, vfmav<Cmplx<typename Tsimd::value_type>> &dst)
  {
  constexpr auto vlen=Tsimd::size();
  copy_outputx(it, reinterpret_cast<const typename Tsimd::value_type *>(src), dst, vlen);
  }

template<typename Tsimd, typename Titer> DUCC0_NOINLINE void copy_output(const Titer &it,
  const Tsimd *DUCC0_RESTRICT src, vfmav<typename Tsimd::value_type> &dst)
  {
  constexpr auto vlen=Tsimd::size();
  auto ptr=dst.data();
  for (size_t i=0; i<it.length_out(); ++i)
    for (size_t j=0; j<vlen; ++j)
      ptr[it.oofs(j,i)] = src[i][j];
  }

template<typename T, typename Titer> DUCC0_NOINLINE void copy_output(const Titer &it,
  const T *DUCC0_RESTRICT src, vfmav<T> &dst)
  {
  auto ptr=dst.data();
  if (src == &dst.raw(it.oofs(0))) return;  // in-place
  for (size_t i=0; i<it.length_out(); ++i)
    ptr[it.oofs(i)] = src[i];
  }
template <typename Tsimd, typename Titer> DUCC0_NOINLINE void copy_input(const Titer &it,
  const cfmav<Cmplx<typename Tsimd::value_type>> &src, Cmplx<Tsimd> *dst, size_t nvec, size_t vstr)
  {
  constexpr auto vlen=Tsimd::size();
  for (size_t i=0; i<it.length_in(); ++i)
    for (size_t j0=0; j0<nvec; ++j0)
      for (size_t j1=0; j1<vlen; ++j1)
        {
        dst[j0*vstr+i].r[j1] = src.raw(it.iofs(j0*vlen+j1,i)).r;
        dst[j0*vstr+i].i[j1] = src.raw(it.iofs(j0*vlen+j1,i)).i;
        }
  }
template <typename T, typename Titer> DUCC0_NOINLINE void copy_input(const Titer &it,
  const cfmav<Cmplx<T>> &src, Cmplx<T> *dst, size_t nvec, size_t vstr)
  {
  for (size_t i=0; i<it.length_in(); ++i)
    for (size_t j0=0; j0<nvec; ++j0)
      dst[j0*vstr+i] = src.raw(it.iofs(j0,i));
  }

template <typename Tsimd, typename Titer> DUCC0_NOINLINE void copy_input(const Titer &it,
  const cfmav<typename Tsimd::value_type> &src, Tsimd *dst, size_t nvec, size_t vstr)
  {
  constexpr auto vlen=Tsimd::size();
  for (size_t i=0; i<it.length_in(); ++i)
    for (size_t j0=0; j0<nvec; ++j0)
      for (size_t j1=0; j1<vlen; ++j1)
        dst[j0*vstr+i][j1] = src.raw(it.iofs(j0*vlen+j1,i));
  }

template <typename T, typename Titer> DUCC0_NOINLINE void copy_input(const Titer &it,
  const cfmav<T> &src, T *dst, size_t nvec, size_t vstr)
  {
  for (size_t i=0; i<it.length_in(); ++i)
    for (size_t j0=0; j0<nvec; ++j0)
      dst[j0*vstr+i] = src.raw(it.iofs(j0,i));
  }

template<typename Tsimd, typename Titer> DUCC0_NOINLINE void copy_output(const Titer &it,
  const Cmplx<Tsimd> *src, vfmav<Cmplx<typename Tsimd::value_type>> &dst, size_t nvec, size_t vstr)
  {
  constexpr auto vlen=Tsimd::size();
  Cmplx<typename Tsimd::value_type> * DUCC0_RESTRICT ptr = dst.data();
  for (size_t i=0; i<it.length_out(); ++i)
    for (size_t j0=0; j0<nvec; ++j0)
      for (size_t j1=0; j1<vlen; ++j1)
        ptr[it.oofs(j0*vlen+j1,i)].Set(src[j0*vstr+i].r[j1],src[j0*vstr+i].i[j1]);
  }
template<typename T, typename Titer> DUCC0_NOINLINE void copy_output(const Titer &it,
  const Cmplx<T> *src, vfmav<Cmplx<T>> &dst, size_t nvec, size_t vstr)
  {
  Cmplx<T> * DUCC0_RESTRICT ptr = dst.data();
  for (size_t i=0; i<it.length_out(); ++i)
    for (size_t j0=0; j0<nvec; ++j0)
      ptr[it.oofs(j0,i)] = src[j0*vstr+i];
  }
template<typename Tsimd, typename Titer> DUCC0_NOINLINE void copy_output(const Titer &it,
  const Tsimd *src, vfmav<typename Tsimd::value_type> &dst, size_t nvec, size_t vstr)
  {
  constexpr auto vlen=Tsimd::size();
  typename Tsimd::value_type * DUCC0_RESTRICT ptr = dst.data();
  for (size_t i=0; i<it.length_out(); ++i)
    for (size_t j0=0; j0<nvec; ++j0)
      for (size_t j1=0; j1<vlen; ++j1)
        ptr[it.oofs(j0*vlen+j1,i)] = src[j0*vstr+i][j1];
  }
template<typename T, typename Titer> DUCC0_NOINLINE void copy_output(const Titer &it,
  const T *src, vfmav<T> &dst, size_t nvec, size_t vstr)
  {
  T * DUCC0_RESTRICT ptr = dst.data();
  for (size_t i=0; i<it.length_out(); ++i)
    for (size_t j0=0; j0<nvec; ++j0)
      ptr[it.oofs(j0,i)] = src[j0*vstr+i];
  }


template <typename T, size_t vlen> struct add_vec
  { using type = typename simd_select<T, vlen>::type; };
template <typename T, size_t vlen> struct add_vec<Cmplx<T>, vlen>
  { using type = Cmplx<typename simd_select<T, vlen>::type>; };
template <typename T, size_t vlen> using add_vec_t = typename add_vec<T, vlen>::type;

template<typename Tplan, typename T, typename T0, typename Exec>
DUCC0_NOINLINE void general_nd(const cfmav<T> &in, vfmav<T> &out,
  const shape_t &axes, T0 fct, size_t nthreads, const Exec &exec,
  const bool /*allow_inplace*/=true)
  {
  if ((in.ndim()==1)&&(in.stride(0)==1)&&(out.stride(0)==1))
    {
    auto plan = get_plan<Tplan>(in.shape(0), true);
    exec.exec_simple(in.data(), out.data(), *plan, fct, nthreads);
    return;
    }
  std::shared_ptr<Tplan> plan;
  size_t nth1d = (in.ndim()==1) ? nthreads : 1;

  for (size_t iax=0; iax<axes.size(); ++iax)
    {
    size_t len=in.shape(axes[iax]);
    if ((!plan) || (len!=plan->length()))
      plan = get_plan<Tplan>(len, in.ndim()==1);

    execParallel(util::thread_count(nthreads, in, axes[iax], fft_simdlen<T0>),
      [&](Scheduler &sched)
      {
      constexpr auto vlen = fft_simdlen<T0>;
      constexpr size_t nmax = 16;
      const auto &tin(iax==0? in : out);
      multi_iter<nmax> it(tin, out, axes[iax], sched.num_threads(), sched.thread_num());

      //size_t working_set_size = (2-(tin.data()==out.data()))*len*sizeof(T)+plan->bufsize()*sizeof(T);
      size_t working_set_size = sizeof(T)*(len+plan->bufsize());
      size_t n_simul=1;
      while(n_simul*2*working_set_size<=262144) n_simul*=2;
      n_simul = min(n_simul, vlen);
      size_t n_bunch = n_simul;
      if ((in.stride(axes[iax])!=1) || (out.stride(axes[iax])!=1))  // might make sense to bunch
        //while ((n_bunch<nmax) && ((n_bunch*sizeof(T)<=32) || (n_bunch*2*working_set_size<=262144)))
        while ((n_bunch<nmax) && (n_bunch*sizeof(T)<=32))
          n_bunch*=2;
      bool inplace = (in.stride(axes[iax])==1) && (out.stride(axes[iax])==1) && (n_bunch==1);
      MR_assert(n_bunch<=nmax, "must not happen");
      TmpStorage<T,T0> storage(in.size()/len, len, plan->bufsize(), (n_bunch+vlen-1)/vlen, inplace);

      // first, do all possible steps of size n_bunch, then n_simul
      if (n_bunch>1)
        {
#ifndef DUCC0_NO_SIMD
        if constexpr (vlen>1)
          {
          constexpr size_t lvlen = vlen;
          if (n_simul>=lvlen)
            {
            if (it.remaining()>=n_bunch)
              {
              TmpStorage2<add_vec_t<T, lvlen>,T,T0> storage2(storage);
              while (it.remaining()>=n_bunch)
                {
                it.advance(n_bunch);
                exec.exec_n(it, tin, out, storage2, *plan, fct, n_bunch/lvlen, nth1d);
                }
              }
            }
          if (n_simul==lvlen)
            {
            if (it.remaining()>=lvlen)
              {
              TmpStorage2<add_vec_t<T, lvlen>,T,T0> storage2(storage);
              while (it.remaining()>=lvlen)
                {
                it.advance(lvlen);
                exec(it, tin, out, storage2, *plan, fct, nth1d);
                }
              }
            }
          }
        if constexpr ((vlen>2) && (simd_exists<T0,vlen/2>))
          {
          constexpr size_t lvlen = vlen/2;
          if (n_simul>=lvlen)
            {
            if (it.remaining()>=n_bunch)
              {
              TmpStorage2<add_vec_t<T, lvlen>,T,T0> storage2(storage);
              while (it.remaining()>=n_bunch)
                {
                it.advance(n_bunch);
                exec.exec_n(it, tin, out, storage2, *plan, fct, n_bunch/lvlen, nth1d);
                }
              }
            }
          if (n_simul==lvlen)
            {
            if (it.remaining()>=lvlen)
              {
              TmpStorage2<add_vec_t<T, lvlen>,T,T0> storage2(storage);
              while (it.remaining()>=lvlen)
                {
                it.advance(lvlen);
                exec(it, tin, out, storage2, *plan, fct, nth1d);
                }
              }
            }
          }
        if constexpr ((vlen>4) && (simd_exists<T0,vlen/4>))
          {
          constexpr size_t lvlen = vlen/4;
          if (n_simul>=lvlen)
            {
            if (it.remaining()>=n_bunch)
              {
              TmpStorage2<add_vec_t<T, lvlen>,T,T0> storage2(storage);
              while (it.remaining()>=n_bunch)
                {
                it.advance(n_bunch);
                exec.exec_n(it, tin, out, storage2, *plan, fct, n_bunch/lvlen, nth1d);
                }
              }
            }
          if (n_simul==lvlen)
            {
            if (it.remaining()>=lvlen)
              {
              TmpStorage2<add_vec_t<T, lvlen>,T,T0> storage2(storage);
              while (it.remaining()>=lvlen)
                {
                it.advance(lvlen);
                exec(it, tin, out, storage2, *plan, fct, nth1d);
                }
              }
            }
          }
#endif
        {
        TmpStorage2<T,T,T0> storage2(storage);
        while (it.remaining()>=n_bunch)
          {
          it.advance(n_bunch);
          exec.exec_n(it, tin, out, storage2, *plan, fct, n_bunch, nth1d);
          }
        }
        }
        {
        TmpStorage2<T,T,T0> storage2(storage);
        while (it.remaining()>0)
          {
          it.advance(1);
          exec(it, tin, out, storage2, *plan, fct, nth1d, inplace);
          }
        }
      });  // end of parallel region
    fct = T0(1); // factor has been applied, use 1 for remaining axes
    }
  }

struct ExecC2C
  {
  bool forward;

  template <typename T0, typename Tstorage, typename Titer> DUCC0_NOINLINE void operator() (
    const Titer &it, const cfmav<Cmplx<T0>> &in,
    vfmav<Cmplx<T0>> &out, Tstorage &storage, const pocketfft_c<T0> &plan, T0 fct,
    size_t nthreads, bool inplace=false) const
    {
    using T = typename Tstorage::datatype;
    if constexpr(is_same<Cmplx<T0>, T>::value)
      if (inplace)
        {
        if (in.data()!=out.data())
          copy_input(it, in, out.data()+it.oofs(0));
        plan.exec_copyback(out.data()+it.oofs(0), storage.transformBuf(), fct, forward, nthreads);
        return;
        }
    T *buf1=storage.transformBuf(), *buf2=storage.dataBuf();
    copy_input(it, in, buf2);
    auto res = plan.exec(buf2, buf1, fct, forward, nthreads);
    copy_output(it, res, out);
    }
  template <typename T0, typename Tstorage, typename Titer> DUCC0_NOINLINE void exec_n (
    const Titer &it, const cfmav<Cmplx<T0>> &in,
    vfmav<Cmplx<T0>> &out, Tstorage &storage, const pocketfft_c<T0> &plan, T0 fct, size_t nvec,
    size_t nthreads) const
    {
    using T = typename Tstorage::datatype;
    size_t dstr = storage.data_stride();
    T *buf1=storage.transformBuf(), *buf2=storage.dataBuf();
    copy_input(it, in, buf2, nvec, dstr);
    for (size_t i=0; i<nvec; ++i)
      plan.exec_copyback(buf2+i*dstr, buf1, fct, forward, nthreads);
    copy_output(it, buf2, out, nvec, dstr);
    }
  template <typename T0> DUCC0_NOINLINE void exec_simple (
    const Cmplx<T0> *in, Cmplx<T0> *out, const pocketfft_c<T0> &plan, T0 fct,
    size_t nthreads) const
    {
    if (in!=out) copy_n(in, plan.length(), out);
    plan.exec(out, fct, forward, nthreads);
    }
  };

struct ExecHartley
  {
  template <typename T0, typename Tstorage, typename Titer> DUCC0_NOINLINE void operator() (
    const Titer &it, const cfmav<T0> &in, vfmav<T0> &out,
    Tstorage &storage, const pocketfft_hartley<T0> &plan, T0 fct, size_t nthreads,
    bool inplace=false) const
    {
    using T = typename Tstorage::datatype;
    if constexpr(is_same<T0, T>::value)
      if (inplace)
        {
        if (in.data()!=out.data())
          copy_input(it, in, out.data()+it.oofs(0));
        plan.exec_copyback(out.data()+it.oofs(0), storage.transformBuf(), fct, nthreads);
        return;
        }
    T *buf1=storage.transformBuf(), *buf2=storage.dataBuf();
    copy_input(it, in, buf2);
    auto res = plan.exec(buf2, buf1, fct, nthreads);
    copy_output(it, res, out);
    }
  template <typename T0, typename Tstorage, typename Titer> DUCC0_NOINLINE void exec_n (
    const Titer &it, const cfmav<T0> &in,
    vfmav<T0> &out, Tstorage &storage, const pocketfft_hartley<T0> &plan, T0 fct, size_t nvec,
    size_t nthreads) const
    {
    using T = typename Tstorage::datatype;
    size_t dstr = storage.data_stride();
    T *buf1=storage.transformBuf(), *buf2=storage.dataBuf();
    copy_input(it, in, buf2, nvec, dstr);
    for (size_t i=0; i<nvec; ++i)
      plan.exec_copyback(buf2+i*dstr, buf1, fct, nthreads);
    copy_output(it, buf2, out, nvec, dstr);
    }
  template <typename T0> DUCC0_NOINLINE void exec_simple (
    const T0 *in, T0 *out, const pocketfft_hartley<T0> &plan, T0 fct,
    size_t nthreads) const
    {
    if (in!=out) copy_n(in, plan.length(), out);
    plan.exec(out, fct, nthreads);
    }
  };

struct ExecFHT
  {
  template <typename T0, typename Tstorage, typename Titer> DUCC0_NOINLINE void operator() (
    const Titer &it, const cfmav<T0> &in, vfmav<T0> &out,
    Tstorage &storage, const pocketfft_fht<T0> &plan, T0 fct, size_t nthreads,
    bool inplace=false) const
    {
    using T = typename Tstorage::datatype;
    if constexpr(is_same<T0, T>::value)
      if (inplace)
        {
        if (in.data()!=out.data())
          copy_input(it, in, out.data()+it.oofs(0));
        plan.exec_copyback(out.data()+it.oofs(0), storage.transformBuf(), fct, nthreads);
        return;
        }
    T *buf1=storage.transformBuf(), *buf2=storage.dataBuf();
    copy_input(it, in, buf2);
    auto res = plan.exec(buf2, buf1, fct, nthreads);
    copy_output(it, res, out);
    }
  template <typename T0, typename Tstorage, typename Titer> DUCC0_NOINLINE void exec_n (
    const Titer &it, const cfmav<T0> &in,
    vfmav<T0> &out, Tstorage &storage, const pocketfft_fht<T0> &plan, T0 fct, size_t nvec,
    size_t nthreads) const
    {
    using T = typename Tstorage::datatype;
    size_t dstr = storage.data_stride();
    T *buf1=storage.transformBuf(), *buf2=storage.dataBuf();
    copy_input(it, in, buf2, nvec, dstr);
    for (size_t i=0; i<nvec; ++i)
      plan.exec_copyback(buf2+i*dstr, buf1, fct, nthreads);
    copy_output(it, buf2, out, nvec, dstr);
    }
  template <typename T0> DUCC0_NOINLINE void exec_simple (
    const T0 *in, T0 *out, const pocketfft_fht<T0> &plan, T0 fct,
    size_t nthreads) const
    {
    if (in!=out) copy_n(in, plan.length(), out);
    plan.exec(out, fct, nthreads);
    }
  };

struct ExecFFTW
  {
  bool forward;

  template <typename T0, typename Tstorage, typename Titer> DUCC0_NOINLINE void operator() (
    const Titer &it, const cfmav<T0> &in, vfmav<T0> &out,
    Tstorage &storage, const pocketfft_fftw<T0> &plan, T0 fct, size_t nthreads,
    bool inplace=false) const
    {
    using T = typename Tstorage::datatype;
    if constexpr(is_same<T0, T>::value)
      if (inplace)
        {
        if (in.data()!=out.data())
          copy_input(it, in, out.data()+it.oofs(0));
        plan.exec_copyback(out.data()+it.oofs(0), storage.transformBuf(), fct, forward, nthreads);
        return;
        }
    T *buf1=storage.transformBuf(), *buf2=storage.dataBuf();
    copy_input(it, in, buf2);
    auto res = plan.exec(buf2, buf1, fct, forward, nthreads);
    copy_output(it, res, out);
    }
  template <typename T0, typename Tstorage, typename Titer> DUCC0_NOINLINE void exec_n (
    const Titer &it, const cfmav<T0> &in,
    vfmav<T0> &out, Tstorage &storage, const pocketfft_fftw<T0> &plan, T0 fct, size_t nvec,
    size_t nthreads) const
    {
    using T = typename Tstorage::datatype;
    size_t dstr = storage.data_stride();
    T *buf1=storage.transformBuf(), *buf2=storage.dataBuf();
    copy_input(it, in, buf2, nvec, dstr);
    for (size_t i=0; i<nvec; ++i)
      plan.exec_copyback(buf2+i*dstr, buf1, fct, forward, nthreads);
    copy_output(it, buf2, out, nvec, dstr);
    }
  template <typename T0> DUCC0_NOINLINE void exec_simple (
    const T0 *in, T0 *out, const pocketfft_fftw<T0> &plan, T0 fct,
    size_t nthreads) const
    {
    if (in!=out) copy_n(in, plan.length(), out);
    plan.exec(out, fct, forward, nthreads);
    }
  };

struct ExecDcst
  {
  bool ortho;
  int type;
  bool cosine;

  template <typename T0, typename Tstorage, typename Tplan, typename Titer>
  DUCC0_NOINLINE void operator() (const Titer &it, const cfmav<T0> &in,
    vfmav <T0> &out, Tstorage &storage, const Tplan &plan, T0 fct, size_t nthreads,
    bool inplace=false) const
    {
    using T = typename Tstorage::datatype;
    if constexpr(is_same<T0, T>::value)
      if (inplace)
        {
        if (in.data()!=out.data())
          copy_input(it, in, out.data()+it.oofs(0));
        plan.exec_copyback(out.data()+it.oofs(0), storage.transformBuf(), fct, ortho, type, cosine, nthreads);
        return;
        }
    T *buf1=storage.transformBuf(), *buf2=storage.dataBuf();
    copy_input(it, in, buf2);
    auto res = plan.exec(buf2, buf1, fct, ortho, type, cosine, nthreads);
    copy_output(it, res, out);
    }
  template <typename T0, typename Tstorage, typename Tplan, typename Titer> DUCC0_NOINLINE void exec_n (
    const Titer &it, const cfmav<T0> &in,
    vfmav<T0> &out, Tstorage &storage, const Tplan &plan, T0 fct, size_t nvec,
    size_t nthreads) const
    {
    using T = typename Tstorage::datatype;
    size_t dstr = storage.data_stride();
    T *buf1=storage.transformBuf(), *buf2=storage.dataBuf();
    copy_input(it, in, buf2, nvec, dstr);
    for (size_t i=0; i<nvec; ++i)
      plan.exec_copyback(buf2+i*dstr, buf1, fct, ortho, type, cosine, nthreads);
    copy_output(it, buf2, out, nvec, dstr);
    }
  template <typename T0, typename Tplan> DUCC0_NOINLINE void exec_simple (
    const T0 *in, T0 *out, const Tplan &plan, T0 fct,
    size_t nthreads) const
    {
    if (in!=out) copy_n(in, plan.length(), out);
    plan.exec(out, fct, ortho, type, cosine, nthreads);
    }
  };

template<typename T> DUCC0_NOINLINE void general_r2c(
  const cfmav<T> &in, vfmav<Cmplx<T>> &out, size_t axis, bool forward, T fct,
  size_t nthreads)
  {
  size_t nth1d = (in.ndim()==1) ? nthreads : 1;
  auto plan = std::make_unique<pocketfft_r<T>>(in.shape(axis));
  size_t len=in.shape(axis);
  execParallel(
    util::thread_count(nthreads, in, axis, fft_simdlen<T>),
    [&](Scheduler &sched) {
    constexpr auto vlen = fft_simdlen<T>;
    TmpStorage<T,T> storage(in.size()/len, len, plan->bufsize(), 1, false);
    multi_iter<vlen> it(in, out, axis, sched.num_threads(), sched.thread_num());
#ifndef DUCC0_NO_SIMD
    if constexpr (vlen>1)
      {
      TmpStorage2<add_vec_t<T, vlen>,T,T> storage2(storage);
      auto dbuf = storage2.dataBuf();
      auto tbuf = storage2.transformBuf();
      while (it.remaining()>=vlen)
        {
        it.advance(vlen);
        copy_input(it, in, dbuf);
        auto res = plan->exec(dbuf, tbuf, fct, true, nth1d);
        auto vout = out.data();
        for (size_t j=0; j<vlen; ++j)
          vout[it.oofs(j,0)].Set(res[0][j]);
        size_t i=1, ii=1;
        if (forward)
          for (; i<len-1; i+=2, ++ii)
            for (size_t j=0; j<vlen; ++j)
              vout[it.oofs(j,ii)].Set(res[i][j], res[i+1][j]);
        else
          for (; i<len-1; i+=2, ++ii)
            for (size_t j=0; j<vlen; ++j)
              vout[it.oofs(j,ii)].Set(res[i][j], -res[i+1][j]);
        if (i<len)
          for (size_t j=0; j<vlen; ++j)
            vout[it.oofs(j,ii)].Set(res[i][j]);
        }
      }
    if constexpr (vlen>2)
      if constexpr (simd_exists<T,vlen/2>)
        if (it.remaining()>=vlen/2)
          {
          TmpStorage2<add_vec_t<T, vlen/2>,T,T> storage2(storage);
          auto dbuf = storage2.dataBuf();
          auto tbuf = storage2.transformBuf();
          it.advance(vlen/2);
          copy_input(it, in, dbuf);
          auto res = plan->exec(dbuf, tbuf, fct, true, nth1d);
          auto vout = out.data();
          for (size_t j=0; j<vlen/2; ++j)
            vout[it.oofs(j,0)].Set(res[0][j]);
          size_t i=1, ii=1;
          if (forward)
            for (; i<len-1; i+=2, ++ii)
              for (size_t j=0; j<vlen/2; ++j)
                vout[it.oofs(j,ii)].Set(res[i][j], res[i+1][j]);
          else
            for (; i<len-1; i+=2, ++ii)
              for (size_t j=0; j<vlen/2; ++j)
                vout[it.oofs(j,ii)].Set(res[i][j], -res[i+1][j]);
          if (i<len)
            for (size_t j=0; j<vlen/2; ++j)
              vout[it.oofs(j,ii)].Set(res[i][j]);
          }
    if constexpr (vlen>4)
      if constexpr( simd_exists<T,vlen/4>)
        if (it.remaining()>=vlen/4)
          {
          TmpStorage2<add_vec_t<T, vlen/4>,T,T> storage2(storage);
          auto dbuf = storage2.dataBuf();
          auto tbuf = storage2.transformBuf();
          it.advance(vlen/4);
          copy_input(it, in, dbuf);
          auto res = plan->exec(dbuf, tbuf, fct, true, nth1d);
          auto vout = out.data();
          for (size_t j=0; j<vlen/4; ++j)
            vout[it.oofs(j,0)].Set(res[0][j]);
          size_t i=1, ii=1;
          if (forward)
            for (; i<len-1; i+=2, ++ii)
              for (size_t j=0; j<vlen/4; ++j)
                vout[it.oofs(j,ii)].Set(res[i][j], res[i+1][j]);
          else
            for (; i<len-1; i+=2, ++ii)
              for (size_t j=0; j<vlen/4; ++j)
                vout[it.oofs(j,ii)].Set(res[i][j], -res[i+1][j]);
          if (i<len)
            for (size_t j=0; j<vlen/4; ++j)
              vout[it.oofs(j,ii)].Set(res[i][j]);
          }
#endif
    {
    TmpStorage2<T,T,T> storage2(storage);
    auto dbuf = storage2.dataBuf();
    auto tbuf = storage2.transformBuf();
    while (it.remaining()>0)
      {
      it.advance(1);
      copy_input(it, in, dbuf);
      auto res = plan->exec(dbuf, tbuf, fct, true, nth1d);
      auto vout = out.data();
      vout[it.oofs(0)].Set(res[0]);
      size_t i=1, ii=1;
      if (forward)
        for (; i<len-1; i+=2, ++ii)
          vout[it.oofs(ii)].Set(res[i], res[i+1]);
      else
        for (; i<len-1; i+=2, ++ii)
          vout[it.oofs(ii)].Set(res[i], -res[i+1]);
      if (i<len)
        vout[it.oofs(ii)].Set(res[i]);
      }
    }
    });  // end of parallel region
  }
template<typename T> DUCC0_NOINLINE void general_c2r(
  const cfmav<Cmplx<T>> &in, vfmav<T> &out, size_t axis, bool forward, T fct,
  size_t nthreads)
  {
  size_t nth1d = (in.ndim()==1) ? nthreads : 1;
  auto plan = std::make_unique<pocketfft_r<T>>(out.shape(axis));
  size_t len=out.shape(axis);
  execParallel(
    util::thread_count(nthreads, in, axis, fft_simdlen<T>),
    [&](Scheduler &sched) {
      constexpr auto vlen = fft_simdlen<T>;
      TmpStorage<T,T> storage(out.size()/len, len, plan->bufsize(), 1, false);
      multi_iter<vlen> it(in, out, axis, sched.num_threads(), sched.thread_num());
#ifndef DUCC0_NO_SIMD
      if constexpr (vlen>1)
        {
        TmpStorage2<add_vec_t<T, vlen>,T,T> storage2(storage);
        auto dbuf = storage2.dataBuf();
        auto tbuf = storage2.transformBuf();
        while (it.remaining()>=vlen)
          {
          it.advance(vlen);
          for (size_t j=0; j<vlen; ++j)
            dbuf[0][j]=in.raw(it.iofs(j,0)).r;
          {
          size_t i=1, ii=1;
          if (forward)
            for (; i<len-1; i+=2, ++ii)
              for (size_t j=0; j<vlen; ++j)
                {
                dbuf[i  ][j] =  in.raw(it.iofs(j,ii)).r;
                dbuf[i+1][j] = -in.raw(it.iofs(j,ii)).i;
                }
          else
            for (; i<len-1; i+=2, ++ii)
              for (size_t j=0; j<vlen; ++j)
                {
                dbuf[i  ][j] = in.raw(it.iofs(j,ii)).r;
                dbuf[i+1][j] = in.raw(it.iofs(j,ii)).i;
                }
          if (i<len)
            for (size_t j=0; j<vlen; ++j)
              dbuf[i][j] = in.raw(it.iofs(j,ii)).r;
          }
          auto res = plan->exec(dbuf, tbuf, fct, false, nth1d);
          copy_output(it, res, out);
          }
        }
      if constexpr (vlen>2)
        if constexpr (simd_exists<T,vlen/2>)
          if (it.remaining()>=vlen/2)
            {
            TmpStorage2<add_vec_t<T, vlen/2>,T,T> storage2(storage);
            auto dbuf = storage2.dataBuf();
            auto tbuf = storage2.transformBuf();
            it.advance(vlen/2);
            for (size_t j=0; j<vlen/2; ++j)
              dbuf[0][j]=in.raw(it.iofs(j,0)).r;
            {
            size_t i=1, ii=1;
            if (forward)
              for (; i<len-1; i+=2, ++ii)
                for (size_t j=0; j<vlen/2; ++j)
                  {
                  dbuf[i  ][j] =  in.raw(it.iofs(j,ii)).r;
                  dbuf[i+1][j] = -in.raw(it.iofs(j,ii)).i;
                  }
            else
              for (; i<len-1; i+=2, ++ii)
                for (size_t j=0; j<vlen/2; ++j)
                  {
                  dbuf[i  ][j] = in.raw(it.iofs(j,ii)).r;
                  dbuf[i+1][j] = in.raw(it.iofs(j,ii)).i;
                  }
            if (i<len)
              for (size_t j=0; j<vlen/2; ++j)
                dbuf[i][j] = in.raw(it.iofs(j,ii)).r;
            }
            auto res = plan->exec(dbuf, tbuf, fct, false, nth1d);
            copy_output(it, res, out);
            }
      if constexpr (vlen>4)
        if constexpr(simd_exists<T,vlen/4>)
          if (it.remaining()>=vlen/4)
            {
            TmpStorage2<add_vec_t<T, vlen/4>,T,T> storage2(storage);
            auto dbuf = storage2.dataBuf();
            auto tbuf = storage2.transformBuf();
            it.advance(vlen/4);
            for (size_t j=0; j<vlen/4; ++j)
              dbuf[0][j]=in.raw(it.iofs(j,0)).r;
            {
            size_t i=1, ii=1;
            if (forward)
              for (; i<len-1; i+=2, ++ii)
                for (size_t j=0; j<vlen/4; ++j)
                  {
                  dbuf[i  ][j] =  in.raw(it.iofs(j,ii)).r;
                  dbuf[i+1][j] = -in.raw(it.iofs(j,ii)).i;
                  }
            else
              for (; i<len-1; i+=2, ++ii)
                for (size_t j=0; j<vlen/4; ++j)
                  {
                  dbuf[i  ][j] = in.raw(it.iofs(j,ii)).r;
                  dbuf[i+1][j] = in.raw(it.iofs(j,ii)).i;
                  }
            if (i<len)
              for (size_t j=0; j<vlen/4; ++j)
                dbuf[i][j] = in.raw(it.iofs(j,ii)).r;
            }
            auto res = plan->exec(dbuf, tbuf, fct, false, nth1d);
            copy_output(it, res, out);
            }
#endif
      {
      TmpStorage2<T,T,T> storage2(storage);
      auto dbuf = storage2.dataBuf();
      auto tbuf = storage2.transformBuf();
      while (it.remaining()>0)
        {
        it.advance(1);
        dbuf[0]=in.raw(it.iofs(0)).r;
        {
        size_t i=1, ii=1;
        if (forward)
          for (; i<len-1; i+=2, ++ii)
            {
            dbuf[i  ] =  in.raw(it.iofs(ii)).r;
            dbuf[i+1] = -in.raw(it.iofs(ii)).i;
            }
        else
          for (; i<len-1; i+=2, ++ii)
            {
            dbuf[i  ] = in.raw(it.iofs(ii)).r;
            dbuf[i+1] = in.raw(it.iofs(ii)).i;
            }
        if (i<len)
          dbuf[i] = in.raw(it.iofs(ii)).r;
        }
        auto res = plan->exec(dbuf, tbuf, fct, false, nth1d);
        copy_output(it, res, out);
        }
      }
    });  // end of parallel region
  }

struct ExecR2R
  {
  bool r2c, forward;

  template <typename T0, typename Tstorage, typename Titer> DUCC0_NOINLINE void operator() (
    const Titer &it, const cfmav<T0> &in, vfmav<T0> &out, Tstorage &storage,
    const pocketfft_r<T0> &plan, T0 fct, size_t nthreads,
    bool inplace=false) const
    {
    using T = typename Tstorage::datatype;
    if constexpr(is_same<T0, T>::value)
      if (inplace)
        {
        T *buf1=storage.transformBuf(), *buf2=out.data();
        if (in.data()!=buf2)
          copy_input(it, in, buf2);
        if ((!r2c) && forward)
          for (size_t i=2; i<it.length_out(); i+=2)
            buf2[i] = -buf2[i];
        plan.exec_copyback(buf2, buf1, fct, r2c, nthreads);
        if (r2c && (!forward))
          for (size_t i=2; i<it.length_out(); i+=2)
            buf2[i] = -buf2[i];
        return;
        }

    T *buf1=storage.transformBuf(), *buf2=storage.dataBuf();
    copy_input(it, in, buf2);
    if ((!r2c) && forward)
      for (size_t i=2; i<it.length_out(); i+=2)
        buf2[i] = -buf2[i];
    auto res = plan.exec(buf2, buf1, fct, r2c, nthreads);
    if (r2c && (!forward))
      for (size_t i=2; i<it.length_out(); i+=2)
        res[i] = -res[i];
    copy_output(it, res, out);
    }
  template <typename T0, typename Tstorage, typename Titer> DUCC0_NOINLINE void exec_n (
    const Titer &it, const cfmav<T0> &in,
    vfmav<T0> &out, Tstorage &storage, const pocketfft_r<T0> &plan, T0 fct, size_t nvec,
    size_t nthreads) const
    {
    using T = typename Tstorage::datatype;
    size_t dstr = storage.data_stride();
    T *buf1=storage.transformBuf(), *buf2=storage.dataBuf();
    copy_input(it, in, buf2, nvec, dstr);
    if ((!r2c) && forward)
      for (size_t k=0; k<nvec; ++k)
        for (size_t i=2; i<it.length_out(); i+=2)
          buf2[i+k*dstr] = -buf2[i+k*dstr];
    for (size_t i=0; i<nvec; ++i)
      plan.exec_copyback(buf2+i*dstr, buf1, fct, r2c, nthreads);
    if (r2c && (!forward))
      for (size_t k=0; k<nvec; ++k)
        for (size_t i=2; i<it.length_out(); i+=2)
          buf2[i+k*dstr] = -buf2[i+k*dstr];
    copy_output(it, buf2, out, nvec, dstr);
    }
  template <typename T0> DUCC0_NOINLINE void exec_simple (
    const T0 *in, T0 *out, const pocketfft_r<T0> &plan, T0 fct,
    size_t nthreads) const
    {
    if (in!=out) copy_n(in, plan.length(), out);
    if ((!r2c) && forward)
      for (size_t i=2; i<plan.length(); i+=2)
        out[i] = -out[i];
    plan.exec(out, fct, r2c, nthreads);
    if (r2c && (!forward))
      for (size_t i=2; i<plan.length(); i+=2)
        out[i] = -out[i];
    }
  };

/// Complex-to-complex Fast Fourier Transform
/** This executes a Fast Fourier Transform on \a in and stores the result in
 *  \a out.
 *
 *  \a in and \a out must have identical shapes; they may point to the same
 *  memory; in this case their strides must also be identical.
 *
 *  \a axes specifies the axes over which the transform is carried out.
 *
 *  If \a forward is true, a minus sign will be used in the exponent.
 *
 *  No normalization factors will be applied by default; if multiplication by
 *  a constant is desired, it can be supplied in \a fct.
 *
 *  If the underlying array has more than one dimension, the computation will
 *  be distributed over \a nthreads threads.
 */
template<typename T> DUCC0_NOINLINE void c2c(const cfmav<std::complex<T>> &in,
  vfmav<std::complex<T>> &out, const shape_t &axes, bool forward,
  T fct, size_t nthreads=1)
  {
  util::sanity_check_onetype(in, out, in.data()==out.data(), axes);
  if (in.size()==0) return;
  const auto &in2(reinterpret_cast<const cfmav<Cmplx<T> >&>(in));
  auto &out2(reinterpret_cast<vfmav<Cmplx<T> >&>(out));
  if ((axes.size()>1) && (in.data()!=out.data())) // optimize axis order
    for (size_t i=1; i<axes.size(); ++i)
      if ((in.stride(i)==1)&&(out.stride(i)==1))
        {
        shape_t axes2(axes);
        swap(axes2[0],axes2[i]);
        general_nd<pocketfft_c<T>>(in2, out2, axes2, fct, nthreads, ExecC2C{forward});
        return;
        }
  general_nd<pocketfft_c<T>>(in2, out2, axes, fct, nthreads, ExecC2C{forward});
  }

/// Fast Discrete Cosine Transform
/** This executes a DCT on \a in and stores the result in \a out.
 *
 *  \a in and \a out must have identical shapes; they may point to the same
 *  memory; in this case their strides must also be identical.
 *
 *  \a axes specifies the axes over which the transform is carried out.
 *
 *  If \a forward is true, a DCT is computed, otherwise an inverse DCT.
 *
 *  \a type specifies the desired type (1-4) of the transform.
 *
 *  No normalization factors will be applied by default; if multiplication by
 *  a constant is desired, it can be supplied in \a fct.
 *
 *  If \a ortho is true, the first and last array entries are corrected (if
 *  necessary) to allow an orthonormalized transform.
 *
 *  If the underlying array has more than one dimension, the computation will
 *  be distributed over \a nthreads threads.
 */
template<typename T> DUCC0_NOINLINE void dct(const cfmav<T> &in, vfmav<T> &out,
  const shape_t &axes, int type, T fct, bool ortho, size_t nthreads=1)
  {
  if ((type<1) || (type>4)) throw std::invalid_argument("invalid DCT type");
  util::sanity_check_onetype(in, out, in.data()==out.data(), axes);
  if (in.size()==0) return;
  const ExecDcst exec{ortho, type, true};
  if (type==1)
    general_nd<T_dct1<T>>(in, out, axes, fct, nthreads, exec);
  else if (type==4)
    general_nd<T_dcst4<T>>(in, out, axes, fct, nthreads, exec);
  else
    general_nd<T_dcst23<T>>(in, out, axes, fct, nthreads, exec);
  }

/// Fast Discrete Sine Transform
/** This executes a DST on \a in and stores the result in \a out.
 *
 *  \a in and \a out must have identical shapes; they may point to the same
 *  memory; in this case their strides must also be identical.
 *
 *  \a axes specifies the axes over which the transform is carried out.
 *
 *  If \a forward is true, a DST is computed, otherwise an inverse DST.
 *
 *  \a type specifies the desired type (1-4) of the transform.
 *
 *  No normalization factors will be applied by default; if multiplication by
 *  a constant is desired, it can be supplied in \a fct.
 *
 *  If \a ortho is true, the first and last array entries are corrected (if
 *  necessary) to allow an orthonormalized transform.
 *
 *  If the underlying array has more than one dimension, the computation will
 *  be distributed over \a nthreads threads.
 */
template<typename T> DUCC0_NOINLINE void dst(const cfmav<T> &in, vfmav<T> &out,
  const shape_t &axes, int type, T fct, bool ortho, size_t nthreads=1)
  {
  if ((type<1) || (type>4)) throw std::invalid_argument("invalid DST type");
  util::sanity_check_onetype(in, out, in.data()==out.data(), axes);
  if (in.size()==0) return;
  const ExecDcst exec{ortho, type, false};
  if (type==1)
    general_nd<T_dst1<T>>(in, out, axes, fct, nthreads, exec);
  else if (type==4)
    general_nd<T_dcst4<T>>(in, out, axes, fct, nthreads, exec);
  else
    general_nd<T_dcst23<T>>(in, out, axes, fct, nthreads, exec);
  }

template<typename T> DUCC0_NOINLINE void r2c(const cfmav<T> &in,
  vfmav<std::complex<T>> &out, size_t axis, bool forward, T fct,
  size_t nthreads=1)
  {
  util::sanity_check_cr(out, in, axis);
  if (in.size()==0) return;
  auto &out2(reinterpret_cast<vfmav<Cmplx<T>>&>(out));
  general_r2c(in, out2, axis, forward, fct, nthreads);
  }

template<typename T> DUCC0_NOINLINE void r2c(const cfmav<T> &in,
  vfmav<std::complex<T>> &out, const shape_t &axes,
  bool forward, T fct, size_t nthreads=1)
  {
  util::sanity_check_cr(out, in, axes);
  if (in.size()==0) return;
  r2c(in, out, axes.back(), forward, fct, nthreads);
  if (axes.size()==1) return;

  auto newaxes = shape_t{axes.begin(), --axes.end()};
  c2c(out, out, newaxes, forward, T(1), nthreads);
  }

template<typename T> DUCC0_NOINLINE void c2r(const cfmav<std::complex<T>> &in,
  vfmav<T> &out,  size_t axis, bool forward, T fct, size_t nthreads=1)
  {
  util::sanity_check_cr(in, out, axis);
  if (in.size()==0) return;
  const auto &in2(reinterpret_cast<const cfmav<Cmplx<T>>&>(in));
  general_c2r(in2, out, axis, forward, fct, nthreads);
  }

template<typename T> DUCC0_NOINLINE void c2r(const cfmav<std::complex<T>> &in,
  vfmav<T> &out, const shape_t &axes, bool forward, T fct,
  size_t nthreads=1)
  {
  if (axes.size()==1)
    return c2r(in, out, axes[0], forward, fct, nthreads);
  util::sanity_check_cr(in, out, axes);
  if (in.size()==0) return;
  auto atmp(vfmav<std::complex<T>>::build_noncritical(in.shape(), UNINITIALIZED));
  auto newaxes = shape_t{axes.begin(), --axes.end()};
  c2c(in, atmp, newaxes, forward, T(1), nthreads);
  c2r(atmp, out, axes.back(), forward, fct, nthreads);
  }

template<typename T> DUCC0_NOINLINE void c2r_mut(vfmav<std::complex<T>> &in,
  vfmav<T> &out, const shape_t &axes, bool forward, T fct,
  size_t nthreads=1)
  {
  if (axes.size()==1)
    return c2r(in, out, axes[0], forward, fct, nthreads);
  util::sanity_check_cr(in, out, axes);
  if (in.size()==0) return;
  auto newaxes = shape_t{axes.begin(), --axes.end()};
  c2c(in, in, newaxes, forward, T(1), nthreads);
  c2r(in, out, axes.back(), forward, fct, nthreads);
  }

template<typename T> DUCC0_NOINLINE void r2r_fftpack(const cfmav<T> &in,
  vfmav<T> &out, const shape_t &axes, bool real2hermitian, bool forward,
  T fct, size_t nthreads=1)
  {
  util::sanity_check_onetype(in, out, in.data()==out.data(), axes);
  if (in.size()==0) return;
  general_nd<pocketfft_r<T>>(in, out, axes, fct, nthreads,
    ExecR2R{real2hermitian, forward});
  }

template<typename T> DUCC0_NOINLINE void r2r_fftw(const cfmav<T> &in,
  vfmav<T> &out, const shape_t &axes, bool forward,
  T fct, size_t nthreads=1)
  {
  util::sanity_check_onetype(in, out, in.data()==out.data(), axes);
  if (in.size()==0) return;
  general_nd<pocketfft_fftw<T>>(in, out, axes, fct, nthreads,
    ExecFFTW{forward});
  }

template<typename T> DUCC0_NOINLINE void r2r_separable_hartley(const cfmav<T> &in,
  vfmav<T> &out, const shape_t &axes, T fct, size_t nthreads=1)
  {
  util::sanity_check_onetype(in, out, in.data()==out.data(), axes);
  if (in.size()==0) return;
  general_nd<pocketfft_hartley<T>>(in, out, axes, fct, nthreads,
    ExecHartley{}, false);
  }

template<typename T> DUCC0_NOINLINE void r2r_separable_fht(const cfmav<T> &in,
  vfmav<T> &out, const shape_t &axes, T fct, size_t nthreads=1)
  {
  util::sanity_check_onetype(in, out, in.data()==out.data(), axes);
  if (in.size()==0) return;
  general_nd<pocketfft_fht<T>>(in, out, axes, fct, nthreads,
    ExecFHT{}, false);
  }

template<typename T0, typename T1, typename Func> void hermiteHelper(size_t idim, ptrdiff_t iin,
  ptrdiff_t iout0, ptrdiff_t iout1, const cfmav<T0> &c,
  vfmav<T1> &r, const shape_t &axes, Func func, size_t nthreads)
  {
  auto cstr=c.stride(idim), str=r.stride(idim);
  auto len=r.shape(idim);

  if (idim+1==c.ndim())  // last dimension, not much gain in parallelizing
    {
    if (idim==axes.back())  // halfcomplex axis
      for (size_t i=0,ic=0; i<len/2+1; ++i,ic=len-i)
        func (c.raw(iin+i*cstr), r.raw(iout0+i*str), r.raw(iout1+ic*str));
    else if (find(axes.begin(), axes.end(), idim) != axes.end())  // FFT axis
      for (size_t i=0,ic=0; i<len; ++i,ic=len-i)
        func (c.raw(iin+i*cstr), r.raw(iout0+i*str), r.raw(iout1+ic*str));
    else  // non-FFT axis
      for (size_t i=0; i<len; ++i)
        func (c.raw(iin+i*cstr), r.raw(iout0+i*str), r.raw(iout1+i*str));
    }
  else
    {
    if (idim==axes.back())
      {
      if (nthreads==1)
        for (size_t i=0,ic=0; i<len/2+1; ++i,ic=len-i)
          hermiteHelper(idim+1, iin+i*cstr, iout0+i*str, iout1+ic*str, c, r, axes, func, 1);
      else
        execParallel(0, len/2+1, nthreads, [&](size_t lo, size_t hi)
          {
          for (size_t i=lo,ic=(i==0?0:len-i); i<hi; ++i,ic=len-i)
            hermiteHelper(idim+1, iin+i*cstr, iout0+i*str, iout1+ic*str, c, r, axes, func, 1);
          });
      }
    else if (find(axes.begin(), axes.end(), idim) != axes.end())
      {
      if (nthreads==1)
        for (size_t i=0,ic=0; i<len; ++i,ic=len-i)
          hermiteHelper(idim+1, iin+i*cstr, iout0+i*str, iout1+ic*str, c, r, axes, func, 1);
      else
        execParallel(0, len/2+1, nthreads, [&](size_t lo, size_t hi)
          {
          for (size_t i=lo,ic=(i==0?0:len-i); i<hi; ++i,ic=len-i)
            {
            size_t io0=iout0+i*str, io1=iout1+ic*str;
            hermiteHelper(idim+1, iin+i*cstr, io0, io1, c, r, axes, func, 1);
            if (i!=ic)
              hermiteHelper(idim+1, iin+ic*cstr, io1, io0, c, r, axes, func, 1);
            }
          });
      }
    else
      {
      if (nthreads==1)
        for (size_t i=0; i<len; ++i)
          hermiteHelper(idim+1, iin+i*cstr, iout0+i*str, iout1+i*str, c, r, axes, func, 1);
      else
        execParallel(0, len, nthreads, [&](size_t lo, size_t hi)
          {
          for (size_t i=lo; i<hi; ++i)
            hermiteHelper(idim+1, iin+i*cstr, iout0+i*str, iout1+i*str, c, r, axes, func, 1);
          });
      }
    }
  }

template<typename T> void oscarize(vfmav<T> &data, size_t ax0, size_t ax1,
  size_t nthreads)
  {
  auto nu=data.shape(ax0), nv=data.shape(ax1);
  if ((nu<3)||(nv<3)) return;
  vector<slice> slc(data.ndim());
  slc[ax0] = slice(1,(nu+1)/2);
  slc[ax1] = slice(1,(nv+1)/2);
  auto all = subarray(data, slc);
  slc[ax0] = slice(nu-1,nu/2,-1);
  auto ahl = subarray(data, slc);
  slc[ax1] = slice(nv-1,nv/2,-1);
  auto ahh = subarray(data, slc);
  slc[ax0] = slice(1,(nu+1)/2);
  auto alh = subarray(data, slc);
  mav_apply([](T &ll, T &hl, T &hh, T &lh)
    {
    T tll=ll, thl=hl, tlh=lh, thh=hh;
    T v = T(0.5)*(tll+tlh+thl+thh);
    ll = v-thh;
    hl = v-tlh;
    lh = v-thl;
    hh = v-tll;
    }, nthreads, all, ahl, ahh, alh);
  }

template<typename T> void r2r_genuine_hartley(const cfmav<T> &in,
  vfmav<T> &out, const shape_t &axes, T fct, size_t nthreads=1)
  {
  if (axes.size()==1)
    return r2r_separable_hartley(in, out, axes, fct, nthreads);
  if (axes.size()==2)
    {
    r2r_separable_hartley(in, out, axes, fct, nthreads);
    oscarize(out, axes[0], axes[1], nthreads);
    return;
    }
  util::sanity_check_onetype(in, out, in.data()==out.data(), axes);
  if (in.size()==0) return;
  shape_t tshp(in.shape());
  tshp[axes.back()] = tshp[axes.back()]/2+1;
  auto atmp(vfmav<std::complex<T>>::build_noncritical(tshp, UNINITIALIZED));
  r2c(in, atmp, axes, true, fct, nthreads);
  hermiteHelper(0, 0, 0, 0, atmp, out, axes, [](const std::complex<T> &c, T &r0, T &r1)
    {
    r0 = c.real()+c.imag();
    r1 = c.real()-c.imag();
    }, nthreads);
  }

template<typename T> void r2r_genuine_fht(const cfmav<T> &in,
  vfmav<T> &out, const shape_t &axes, T fct, size_t nthreads=1)
  {
  if (axes.size()==1)
    return r2r_separable_fht(in, out, axes, fct, nthreads);
  if (axes.size()==2)
    {
    r2r_separable_fht(in, out, axes, fct, nthreads);
    oscarize(out, axes[0], axes[1], nthreads);
    return;
    }
  util::sanity_check_onetype(in, out, in.data()==out.data(), axes);
  if (in.size()==0) return;
  shape_t tshp(in.shape());
  tshp[axes.back()] = tshp[axes.back()]/2+1;
  auto atmp(vfmav<std::complex<T>>::build_noncritical(tshp, UNINITIALIZED));
  r2c(in, atmp, axes, true, fct, nthreads);
  hermiteHelper(0, 0, 0, 0, atmp, out, axes, [](const std::complex<T> &c, T &r0, T &r1)
    {
    r0 = c.real()-c.imag();
    r1 = c.real()+c.imag();
    }, nthreads);
  }

template<typename Tplan, typename T0, typename T, typename Exec>
DUCC0_NOINLINE void general_convolve_axis(const cfmav<T> &in, vfmav<T> &out,
  const size_t axis, const cmav<T,1> &kernel, size_t nthreads,
  const Exec &exec)
  {
  std::unique_ptr<Tplan> plan1, plan2;

  size_t l_in=in.shape(axis), l_out=out.shape(axis);
  MR_assert(kernel.size()==l_in, "bad kernel size");
  plan1 = std::make_unique<Tplan>(l_in);
  plan2 = std::make_unique<Tplan>(l_out);
  size_t bufsz = max(plan1->bufsize(), plan2->bufsize());

  vmav<T,1> fkernel({kernel.shape(0)}, UNINITIALIZED);
  for (size_t i=0; i<kernel.shape(0); ++i)
    fkernel(i) = kernel(i);
  plan1->exec(fkernel.data(), T0(1)/T0(l_in), true, nthreads);

  execParallel(
    util::thread_count(nthreads, in, axis, fft_simdlen<T0>),
    [&](Scheduler &sched) {
      constexpr auto vlen = fft_simdlen<T0>;
      TmpStorage<T,T0> storage(in.size()/l_in, l_in+l_out, bufsz, 1, false);
      multi_iter<vlen> it(in, out, axis, sched.num_threads(), sched.thread_num());
#ifndef DUCC0_NO_SIMD
      if constexpr (vlen>1)
        {
        TmpStorage2<add_vec_t<T, vlen>,T,T0> storage2(storage);
        while (it.remaining()>=vlen)
          {
          it.advance(vlen);
          exec(it, in, out, storage2, *plan1, *plan2, fkernel);
          }
        }
      if constexpr (vlen>2)
        if constexpr (simd_exists<T,vlen/2>)
          if (it.remaining()>=vlen/2)
            {
            TmpStorage2<add_vec_t<T, vlen/2>,T,T0> storage2(storage);
            it.advance(vlen/2);
            exec(it, in, out, storage2, *plan1, *plan2, fkernel);
            }
      if constexpr (vlen>4)
        if constexpr (simd_exists<T,vlen/4>)
          if (it.remaining()>=vlen/4)
            {
            TmpStorage2<add_vec_t<T, vlen/4>,T,T0> storage2(storage);
            it.advance(vlen/4);
            exec(it, in, out, storage2, *plan1, *plan2, fkernel);
            }
#endif
      {
      TmpStorage2<T,T,T0> storage2(storage);
      while (it.remaining()>0)
        {
        it.advance(1);
        exec(it, in, out, storage2, *plan1, *plan2, fkernel);
        }
      }
    });  // end of parallel region
  }

struct ExecConv1R
  {
  template <typename T0, typename Tstorage, typename Titer> void operator() (
    const Titer &it, const cfmav<T0> &in, vfmav<T0> &out,
    Tstorage &storage, const pocketfft_r<T0> &plan1, const pocketfft_r<T0> &plan2,
    const cmav<T0,1> &fkernel) const
    {
    using T = typename Tstorage::datatype;
    size_t l_in = plan1.length(),
           l_out = plan2.length(),
           l_min = std::min(l_in, l_out);
    T *buf1=storage.transformBuf(), *buf2=storage.dataBuf();
    copy_input(it, in, buf2);
    plan1.exec_copyback(buf2, buf1, T0(1), true);
    auto res = buf2;
    {
    res[0] *= fkernel(0);
    size_t i;
    for (i=1; 2*i<l_min; ++i)
      {
      Cmplx<T> t1(res[2*i-1], res[2*i]);
      Cmplx<T0> t2(fkernel(2*i-1), fkernel(2*i));
      auto t3 = t1*t2;
      res[2*i-1] = t3.r;
      res[2*i] = t3.i;
      }
    if (2*i==l_min)
      {
      if (l_min<l_out) // padding
        res[2*i-1] *= fkernel(2*i-1)*T0(0.5);
      else if (l_min<l_in) // truncation
        {
        Cmplx<T> t1(res[2*i-1], res[2*i]);
        Cmplx<T0> t2(fkernel(2*i-1), fkernel(2*i));
        res[2*i-1] = (t1*t2).r*T0(2);
        }
      else
        res[2*i-1] *= fkernel(2*i-1);
      }
    }
    for (size_t i=l_in; i<l_out; ++i) res[i] = T(0);
    res = plan2.exec(res, buf1, T0(1), false);
    copy_output(it, res, out);
    }
  };
struct ExecConv1C
  {
  template <typename T0, typename Tstorage, typename Titer> void operator() (
    const Titer &it, const cfmav<Cmplx<T0>> &in, vfmav<Cmplx<T0>> &out,
    Tstorage &storage, const pocketfft_c<T0> &plan1, const pocketfft_c<T0> &plan2,
    const cmav<Cmplx<T0>,1> &fkernel) const
    {
    using T = typename Tstorage::datatype;
    size_t l_in = plan1.length(),
           l_out = plan2.length(),
           l_min = std::min(l_in, l_out);
    T *buf1=storage.transformBuf(), *buf2=storage.dataBuf();
    copy_input(it, in, buf2);
    auto res = plan1.exec(buf2, buf1, T0(1), true);
    auto res2 = buf2+l_in;
    {
    res2[0] = res[0]*fkernel(0);
    size_t i;
    for (i=1; 2*i<l_min; ++i)
      {
      res2[i] = res[i]*fkernel(i);
      res2[l_out-i] = res[l_in-i]*fkernel(l_in-i);
      }
    if (2*i==l_min)
      {
      if (l_min<l_out) // padding
        res2[l_out-i] = res2[i] = res[i]*fkernel(i)*T0(.5);
      else if (l_min<l_in) // truncation
        res2[i] = res[i]*fkernel(i) + res[l_in-i]*fkernel(l_in-i);
      else
        res2[i] = res[i]*fkernel(i);
      ++i;
      }
    for (; 2*i<=l_out; ++i)
      res2[i] = res2[l_out-i] = T(0,0);
    }
    res = plan2.exec(res2, buf1, T0(1), false);
    copy_output(it, res, out);
    }
  };

/// Convolution and zero-padding/truncation along one axis
/** This performs a circular convolution with the kernel \a kernel on axis
 *  \a axis of \a in, applies the necessary zero-padding/truncation on this
 *  axis to give it the length \a out.shape(axis),and returns the result
 *  in \a out.
 *
 *  The main purpose of this routine is efficiency: the combination of the above
 *  operations can be carried out more quickly than running the individual
 *  operations in succession.
 *
 *  \a in and \a out must have identical shapes, with the possible exception
 *  of the axis \a axis; they may point to the same memory; in this case all
 *  of their strides must be identical.
 *
 *  \a axis specifies the axis over which the operation is carried out.
 *
 *  \a kernel must have the same length as \a in.shape(axis); it must be
 *  provided in the same domain as \a in (i.e. not pre-transformed).
 *
 *  If \a in has more than one dimension, the computation will
 *  be distributed over \a nthreads threads.
 */
template<typename T> DUCC0_NOINLINE void convolve_axis(const cfmav<T> &in,
  vfmav<T> &out, size_t axis, const cmav<T,1> &kernel, size_t nthreads=1)
  {
  MR_assert(axis<in.ndim(), "bad axis number");
  MR_assert(in.ndim()==out.ndim(), "dimensionality mismatch");
  if (in.data()==out.data())
    MR_assert(in.stride()==out.stride(), "strides mismatch");
  for (size_t i=0; i<in.ndim(); ++i)
    if (i!=axis)
      MR_assert(in.shape(i)==out.shape(i), "shape mismatch");
  if (in.size()==0) return;
  general_convolve_axis<pocketfft_r<T>, T>(in, out, axis, kernel, nthreads,
    ExecConv1R());
  }
template<typename T> DUCC0_NOINLINE void convolve_axis(const cfmav<complex<T>> &in,
  vfmav<complex<T>> &out, size_t axis, const cmav<complex<T>,1> &kernel,
  size_t nthreads=1)
  {
  MR_assert(axis<in.ndim(), "bad axis number");
  MR_assert(in.ndim()==out.ndim(), "dimensionality mismatch");
  if (in.data()==out.data())
    MR_assert(in.stride()==out.stride(), "strides mismatch");
  for (size_t i=0; i<in.ndim(); ++i)
    if (i!=axis)
      MR_assert(in.shape(i)==out.shape(i), "shape mismatch");
  if (in.size()==0) return;
  const auto &in2(reinterpret_cast<const cfmav<Cmplx<T>>&>(in));
  auto &out2(reinterpret_cast<vfmav<Cmplx<T>>&>(out));
  const auto &kernel2(reinterpret_cast<const cmav<Cmplx<T>,1>&>(kernel));
  general_convolve_axis<pocketfft_c<T>, T>(in2, out2, axis, kernel2, nthreads,
    ExecConv1C());
  }

} // namespace detail_fft

using detail_fft::FORWARD;
using detail_fft::BACKWARD;
using detail_fft::c2c;
using detail_fft::c2r;
using detail_fft::c2r_mut;
using detail_fft::r2c;
using detail_fft::r2r_fftpack;
using detail_fft::r2r_fftw;
using detail_fft::r2r_separable_hartley;
using detail_fft::r2r_genuine_hartley;
using detail_fft::r2r_separable_fht;
using detail_fft::r2r_genuine_fht;
using detail_fft::dct;
using detail_fft::dst;
using detail_fft::convolve_axis;

} // namespace ducc0

#endif // POCKETFFT_HDRONLY_H