#include <complex>
#include <limits>
#include <cstdint>
#include <cstring>
#include <stdexcept>
#include <byteswap.h>

static const unsigned power_of_2 (size_t n)
{
	const unsigned table[37] = {
		32, 0, 1, 26, 2, 23, 27, 0, 3, 16, 24, 30, 28, 11, 0, 13,
		4, 7, 17, 0, 25, 22, 31, 15, 29, 10, 12, 6, 0, 21, 14, 9,
		5, 20, 8, 19, 18
	};

	if (n == 0 || n & (n - 1)) {
		return 0;
	} else {
		return table[n % 37];
	}
}

static unsigned revbin (unsigned x, unsigned ln)
{
#define R2(x) ((((x) >> 1) & 1) | (((x) & 1) << 1))
#define R4(x) (R2(((x) >> 2) & 3) | (R2((x) & 3) << 2))
#define R8(x) (R4(((x) >> 4) & 15) | (R4((x) & 15) << 4))
#define R1(x) R8(x), R8(x+1), R8(x+2), R8(x+3)
#define R(x) R1(16*x), R1(16*x+4), R1(16*x+8), R1(16*x+12)
	const unsigned char u8rev[256] = {
		R(0), R(1), R(2), R(3), R(4), R(5), R(6), R(7),
		R(8), R(9), R(10), R(11), R(12), R(13), R(14), R(15),
	};

	unsigned char * p = static_cast<unsigned char *> (static_cast<void *> (&x));

	for (unsigned i = 0; i < sizeof (x); ++i) {
		p[i] = u8rev[p[i]];
	}

	if (sizeof (x) == 4) {
		x = __bswap_32 (x);
	} else if (sizeof (x) == 8) {
		x = __bswap_64 (x);
	} else {
		for (unsigned i = 0; i < sizeof (x) / 2; ++i) {
			p[i] = p[sizeof (x) - i - 1];
		}
	}

	return (x >> (sizeof (x) * 8 - ln));
}

template<typename T>
class fft
{
	unsigned n;
	unsigned * permutecache;
	std::complex<T> * power;

	inline void permute (std::complex<T> *xs) const
	{
		for (unsigned i = 0; i < n; ++i) {
			if (i < permutecache[i])
				std::swap (xs[i], xs[permutecache[i]]);
		}
	}

	void ward (std::complex<T> *xs, bool backward) const
	{
		int back = backward ? -1 : 1;

		for (size_t b = 1; b < n; b *= 2) {
			for (size_t k = 0; k < b; ++k) {
				auto p = power [(back * n * k / (2 * b)) % n];
				for (size_t j = 0; j < n; j += 2 * b) {
					auto t = p * xs[j + k + b];
					xs[j + k + b] = xs[j + k] - t;
					xs[j + k] += t;
				}
			}
		}
	}

public:
	fft (unsigned n) : n(n)
	{
		unsigned p = power_of_2 (n);

		if (!p) {
			throw std::invalid_argument ("argument must be a power of 2");
		}

		permutecache = new unsigned[n];
		power = new std::complex<T>[n];

		std::complex<T> x (0, -2.0*M_PI/n);

		for (unsigned i = 0, j = 0; i < n; ++i) {
			permutecache[i] = revbin (i, p);
			power[i] = std::exp (x*static_cast<T> (i));
		}
	}

	~fft () {
		delete[] power;
		delete[] permutecache;
	}

	void forward (std::complex<T> *xs) const
	{
		permute (xs);
		ward (xs, false);
	}

	void backward (std::complex<T> *xs) const
	{
		permute (xs);
		ward (xs, true);

		for (size_t i = 0; i < n; ++i) {
			xs[i] /= n;
		}
	}
};