fork(1) download
copy 
  1. #include <iostream>
  2. #include <vector>
  3.  
  4. // Возведение a в степень b по модулю p
  5. int powmod (int a, int b, int p)
  6. {
  7. 	int res = 1;
  8. 	while (b)
  9. 		if (b & 1)
  10. 			res = int ((long long) res * a % p),  --b;
  11. 		else
  12. 			a = int ((long long) a * a % p),  b >>= 1;
  13. 	return res;
  14. }
  15.  
  16. // Вычисление обратного к a элемента по модулю n
  17. int inverse(int a, int n)
  18. {
  19. 	int b0 = n, t, q;
  20. 	int x0 = 0, x1 = 1;
  21. 	if (n == 1) return 1;
  22. 	while (a > 1) {
  23. 		q = a / n;
  24. 		t = n, n = a % n, a = t;
  25. 		t = x0, x0 = x1 - q * x0, x1 = t;
  26. 	}
  27. 	if (x1 < 0) x1 += b0;
  28. 	return x1;
  29. }
  30.  
  31. // Вычисление примитивного корня числового поля
  32. // (Стандартный DFT использует комплексные числа, поэтому изменим алгоритм для работы в Zp)
  33. // Ограничение: p - простое
  34. int generator (int p)
  35. {
  36. 	std::vector<int> fact;
  37. 	int phi = p - 1,  n = phi;
  38. 	for (int i = 2; i * i <= n; ++i)
  39. 		if (n % i == 0)
  40. 		{
  41. 			fact.push_back (i);
  42. 			while (n % i == 0)
  43. 				n /= i;
  44. 		}
  45. 	if (n > 1)
  46. 		fact.push_back (n);
  47.  
  48. 	for (int res = 2; res <= p; ++res) {
  49. 		bool ok = true;
  50. 		for (size_t i = 0; i < fact.size() && ok; ++i)
  51. 			ok &= powmod (res, phi / fact[i], p) != 1;
  52. 		if (ok)  return res;
  53. 	}
  54. 	return -1;
  55. }
  56.  
  57. // Реализация прямого DFT над Zp
  58. // Принимает вектор длины n, эта длина используется при вычислении
  59. // примитивного корня поля
  60. // Вычисления приводятся по mod (n + 1), т.к. длина вектора на 1 меньше p,
  61. // где p: Zp
  62. std::vector<int> forward_dft (std::vector<int>& a)
  63. {
  64. 	std::vector<int> result;
  65. 	int prim_root = generator (a.size());	// Корень степени n из единицы
  66. 	for (int i = 0; i < a.size(); i++)
  67. 	{
  68. 		int sum   = 0;
  69. 		int power = 0;
  70. 		for (int j = 0; j < a.size(); j++)
  71. 		{
  72. 			int power_of_root = powmod (prim_root, power, a.size() + 1);
  73. 			sum	  += power_of_root * a[j];	// Накапливаем сумму произведений элементов строки матрицы на вектор
  74. 			power += i;						// Увеличиваем степень, в которую возводится примитивный корень из единицы
  75. 		}
  76. 		result.push_back(sum % (a.size() + 1)); // Набрали сумму, приводим по mod p, p = длина вектора + 1
  77. 	}
  78. 	return result;
  79. }
  80.  
  81.  
  82. std::vector<int> inversed_dft (std::vector<int>& a)
  83. {
  84. 	std::vector<int> result;
  85. 	int prim_root = generator (a.size());	// Корень степени n из единицы
  86. 	int inv_prim_root = inverse (prim_root, a.size() + 1);	// Обратный к найденному примитивному корню
  87. 	int inv_n = inverse(a.size(), a.size() + 1);
  88. 	for (int i = 0; i < a.size(); i++)
  89. 	{
  90. 		int sum   = 0;
  91. 		int power = 0;
  92. 		for (int j = 0; j < a.size(); j++)
  93. 		{
  94. 			int power_of_inv_root = powmod (inv_prim_root, power, a.size() + 1);
  95. 			sum	  += inv_n * power_of_inv_root * a[j];	// Накапливаем сумму произведений элементов строки матрицы на вектор
  96. 			power += i;							// Увеличиваем степень, в которую возводится примитивный корень из единицы
  97. 		}
  98. 		result.push_back(sum % (a.size() + 1)); // Набрали сумму, приводим по mod p, p = длина вектора + 1
  99. 	}
  100. 	return result;
  101. }
  102.  
  103. // Вывод сообщений
  104. void debug(const std::string& str, std::vector<int>& v)
  105. {
  106. 	std::cout << "[DEBUG]:" << str << std::endl;
  107. 	for(int i = 0; i < v.size(); i++)
  108. 		std::cout << v[i] << " ";
  109. 	std::cout << std::endl;
  110. }
  111.  
  112. int main()
  113. {
  114. 	int a[] = {4, 3, 2, 1};
  115. 	std::vector<int> v(std::begin(a), std::end(a));
  116.  
  117. 	debug("Source vector", v);
  118.  
  119. 	std::vector<int> result = forward_dft(v);
  120.  
  121. 	debug("Forward DFT", result);
  122.  
  123. 	result = inversed_dft(result);
  124.  
  125. 	debug("Inversed DFT", result);
  126. }
Success #stdin #stdout 0s 3464KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
[DEBUG]:Source vector
4 3 2 1 
[DEBUG]:Forward DFT
0 1 2 3 
[DEBUG]:Inversed DFT
4 3 2 1
https://ideone.com/jhd6r7
language:
C++14 (gcc 8.3)
created:
8 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!