# 
# FFT and convolution test (Python)
# 
# Copyright (c) 2026 Project Nayuki. (MIT License)
# https://www.nayuki.io/page/free-small-fft-in-multiple-languages
# 
# Permission is hereby granted, free of charge, to any person obtaining a copy of
# this software and associated documentation files (the "Software"), to deal in
# the Software without restriction, including without limitation the rights to
# use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
# the Software, and to permit persons to whom the Software is furnished to do so,
# subject to the following conditions:
# - The above copyright notice and this permission notice shall be included in
#   all copies or substantial portions of the Software.
# - The Software is provided "as is", without warranty of any kind, express or
#   implied, including but not limited to the warranties of merchantability,
#   fitness for a particular purpose and noninfringement. In no event shall the
#   authors or copyright holders be liable for any claim, damages or other
#   liability, whether in an action of contract, tort or otherwise, arising from,
#   out of or in connection with the Software or the use or other dealings in the
#   Software.
# 

import cmath, math, random, unittest
import fft


class FftTest(unittest.TestCase):
	
	# ---- Test functions ----
	
	def test_fourier_power_of_2(self) -> None:
		for i in range(0, 12 + 1):
			self._test_fft(2 ** i)
	
	
	def test_fourier_small(self) -> None:
		for i in range(0, 30):
			self._test_fft(i)
	
	
	def test_fourier_diverse(self) -> None:
		prev: int = 0
		for i in range(100 + 1):
			n: int = int(round(1500 ** (i / 100)))
			if n > prev:
				self._test_fft(n)
				prev = n
	
	
	def test_convolution_power_of_2(self) -> None:
		for i in range(0, 12 + 1):
			self._test_convolution(2 ** i)
	
	
	def test_convolution_diverse(self) -> None:
		prev: int = 0
		for i in range(100 + 1):
			n: int = int(round(1500 ** (i / 100)))
			if n > prev:
				self._test_convolution(n)
				prev = n
		
	
	def _test_fft(self, size: int) -> None:
		input: list[complex] = FftTest._random_vector(size)
		expect: list[complex] = FftTest._naive_dft(input, False)
		actual: list[complex] = fft.transform(input, False)
		err: float = FftTest._log10_rms_err(expect, actual)
		
		actual = [(x / size) for x in expect]
		actual = fft.transform(actual, True)
		err = max(FftTest._log10_rms_err(input, actual), err)
		print(f"fftsize={size:4d}  logerr={err:5.1f}")
	
	
	def _test_convolution(self, size: int) -> None:
		input0: list[complex] = FftTest._random_vector(size)
		input1: list[complex] = FftTest._random_vector(size)
		expect: list[complex] = FftTest._naive_convolution(input0, input1)
		actual: list[complex] = fft.convolve(input0, input1)
		print(f"convsize={size:4d}  logerr={FftTest._log10_rms_err(expect, actual):5.1f}")
	
	
	# ---- Naive reference computation functions ----
	
	@staticmethod
	def _naive_dft(input: list[complex], inverse: bool) -> list[complex]:
		n: int = len(input)
		output: list[complex] = []
		if n == 0:
			return output
		coef: float = (2 if inverse else -2) * math.pi / n
		for k in range(n):  # For each output element
			s: complex = 0
			for t in range(n):  # For each input element
				s += input[t] * cmath.rect(1, (t * k % n) * coef)
			output.append(s)
		return output
	
	
	@staticmethod
	def _naive_convolution(xvec: list[complex], yvec: list[complex]) -> list[complex]:
		assert len(xvec) == len(yvec)
		n: int = len(xvec)
		result: list[complex] = [0] * n
		for i in range(n):
			for j in range(n):
				result[(i + j) % n] += xvec[i] * yvec[j]
		return result
	
	
	# ---- Utility functions ----
	
	@staticmethod
	def _log10_rms_err(xvec: list[complex], yvec: list[complex]) -> float:
		global _max_log_err
		assert len(xvec) == len(yvec)
		err: float = 10.0**(-99 * 2)
		for (x, y) in zip(xvec, yvec):
			err += abs(x - y) ** 2
		err = math.sqrt(err / max(len(xvec), 1))  # Now this is a root mean square (RMS) error
		err = math.log10(err)
		_max_log_err = max(err, _max_log_err)
		return err
	
	
	@staticmethod
	def _random_vector(n: int) -> list[complex]:
		return [complex(random.uniform(-1.0, 1.0), random.uniform(-1.0, 1.0)) for _ in range(n)]



if __name__ == "__main__":
	_max_log_err: float = -math.inf
	unittest.main(exit=False)
	print()
	print(f"Max log err = {_max_log_err:.1f}")
	print(f"Test {'passed' if (_max_log_err < -10) else 'failed'}")