libsimple_compress/source/simple/compress/fft.hpp

#ifndef SIMPLE_COMPRESS_FFT_HPP
#define SIMPLE_COMPRESS_FFT_HPP

#include <iterator> // std::begin std::end
#include <cassert> // assert
#include <cstddef> // std::size_t

#include "simple/support/algorithm/utils.hpp" // support::distance
#include "simple/support/algorithm/numeric.hpp" // support::midpoint
#include "simple/support/range.hpp" // support::range

namespace simple::compress
{

	template <typename I, typename O, typename Rotate, typename RotationStep>
	constexpr void fft(Rotate rotate, RotationStep r_step, I in, std::size_t in_step, O&& out)
	{
		using std::begin; using std::end;

		assert(support::distance(out) != 0);
		if(support::distance(out) == 1)
		{
			*begin(out) = *in;
			return;
		}

		auto even_out = support::range{begin(out), support::midpoint(out)};
		auto odd_out = support::range{support::midpoint(out), end(out)};

		fft(rotate, rotate(r_step, r_step), in, in_step*2, even_out);
		fft(rotate, rotate(r_step, r_step), in + in_step, in_step*2, odd_out);

		auto even = begin(even_out);
		auto odd = begin(odd_out);

		// NOTE: do first step of the loop manually so that we don't need to initialize identity, don't have powerful enough type system for that yet
		{
			auto odd_ = *odd;
			*odd = *even - odd_;
			*even = *even + odd_;
			++even; ++odd;
		}

		for(auto r = r_step; even != end(even_out); ++even, ++odd, r = rotate(r, r_step))
		{
			auto odd_ = rotate(*odd, r);
			*odd = *even - odd_;
			*even = *even + odd_;
		}
	}

	template <typename I, typename O, typename Rotate, typename RotationStep>
	constexpr void fft(Rotate rotate, RotationStep step, const I& in, O&& out)
	{
		using std::begin;
		// TODO
		// assert distance(in) = distance(out) and is power of 2
		// assert rstep is tau / distance(in) ??
		fft(rotate, step, begin(in), 1, out);
	}


} // namespace simple::compress

#endif /* end of include guard */