fork download
copy 
  1. // 2^(2^34 - 1) = 2^17179869183 を計算する。
  2. #include <stdio.h>
  3. #include <stdlib.h>
  4. #include <stdint.h>
  5. #define _USE_MATH_DEFINES
  6. #include <math.h>
  7. #include <assert.h>
  8. #include <stdarg.h>
  9. #include <vector>
  10. #include <utility>
  11. #include <numeric>
  12.  
  13. #include <omp.h>
  14.  
  15. // [DBGSW] 中間結果表示
  16. //#define SHOW_ITM
  17.  
  18. // [DBGSW] I表示
  19. //#define SHOW_I
  20.  
  21. // [DBGSW] SIN_TABLE作成状況表示
  22. //#define SHOW_SIN_TABLE
  23.  
  24. // [DBGSW] 消費時間表示
  25. #define SHOW_CSM_TIME
  26.  
  27. // [DBGSW] 計算途中結果の保存
  28. //#define SAVE_ITM_RESULT
  29.  
  30. typedef uint32_t I_TYPE;
  31.  
  32. // [CONF] 基本SINテーブルの最大要素数
  33. #define BASIC_SINTBL_P_MAX		29
  34.  
  35. // [CONF] 多倍長バッファ1要素の10進数での桁数
  36. const int FIG = 5;
  37.  
  38. // [CONF] 10^FIG
  39. const I_TYPE RADIX = 100000;
  40.  
  41. // [CONF] 使用するCPUのL2キャッシュサイズ(2の冪 Byte/コア)の指数
  42. #ifdef NDEBUG
  43. #define EXPN_L2CHACHE_SZ			18			/* 256 KB */
  44. #else
  45. #define EXPN_L2CHACHE_SZ			8			/* 256 Byte (DEBUG) */
  46. #endif
  47.  
  48. // [CONF] 同一L2キャッシュを共用するスレッドの数(2の冪)の指数
  49. // Core i7 (HT有り)では1、Core i5 (HT無し)では0
  50. #define EXPN_L2_SHARING_THREADS		1
  51.  
  52. // [CONF] sizeof(複素数) (2の冪)の指数(double 2個で16 Byte)
  53. #define EXPN_SIZEOF_COMPLEX			4
  54.  
  55. static_assert((size_t)1 << EXPN_SIZEOF_COMPLEX == 2 * sizeof(double), "*** ERR ***");
  56.  
  57. // L2キャッシュにsrcとdstの両方が複素数が乗る最大要素数
  58. #define EXPN_CAPACITY	(EXPN_L2CHACHE_SZ - EXPN_SIZEOF_COMPLEX - EXPN_L2_SHARING_THREADS - 1)
  59.  
  60. #ifdef _WIN32
  61. int fopen_s_wrp(FILE** const p_fp, const char* const fpath, const char* const mode)
  62. {
  63. 	return (int)fopen_s(p_fp, fpath, mode);
  64. }
  65. int sprintf_s_wrp(char buf[], const size_t bufsz, const char* const fmt, ...)
  66. {
  67. 	va_list ap;
  68. 	int nRet;
  69.  
  70. 	va_start(ap, fmt);
  71. 	nRet = vsprintf_s(buf, bufsz, fmt, ap);
  72. 	va_end(ap);
  73.  
  74. 	return nRet;
  75. }
  76. #else
  77. int fopen_s_wrp(FILE** const p_fp, const char* const fpath, const char* const mode)
  78. {
  79. 	*p_fp = fopen(fpath, mode);
  80. 	return (*p_fp == NULL);
  81. }
  82. int sprintf_s_wrp(char buf[], const size_t bufsz, const char* const fmt, ...)
  83. {
  84. 	va_list ap;
  85. 	int nRet;
  86.  
  87. 	va_start(ap, fmt);
  88. 	nRet = vsprintf(buf, fmt, ap);
  89. 	va_end(ap);
  90.  
  91. 	return nRet;
  92. }
  93. #define _countof(arr)	(sizeof(arr) / sizeof(arr[0]))
  94. #endif
  95.  
  96. /// 実部と虚部が交代に並ぶ形式の複素数リストを表示する。
  97. void show_re_im_seq(const std::vector<double>& s, const int p)
  98. {
  99. 	const size_t hsz = (size_t)1 << (p - 1);
  100.  
  101. 	printf("Re:");
  102. 	for (size_t i = 0; i < hsz; i++) {
  103. 		printf(" % 12f", s[2 * i]);
  104. 	}
  105. 	printf("\n");
  106. 	printf("Im:");
  107. 	for (size_t i = 0; i < hsz; i++) {
  108. 		printf(" % 12f", s[2 * i + 1]);
  109. 	}
  110. 	printf("\n");
  111. }
  112.  
  113. /// 実部と虚部が交代に並ぶ形式の複素数リストを表示する。
  114. void show_re_im_seq(const std::vector<I_TYPE>& s, const int p)
  115. {
  116. 	char fmt[128];
  117. 	sprintf_s_wrp(fmt, _countof(fmt), "%%0%dlu", FIG);
  118.  
  119. 	const size_t hsz = (size_t)1 << (p - 1);
  120. 	printf("Re:");
  121. 	for (size_t i = 0; i < hsz; i++) {
  122. 		printf(fmt, s[2 * i]);
  123. 	}
  124. 	printf("\n");
  125. 	printf("Im:");
  126. 	for (size_t i = 0; i < hsz; i++) {
  127. 		printf(fmt, s[2 * i + 1]);
  128. 	}
  129. 	printf("\n");
  130. }
  131.  
  132. namespace priroot
  133. {
  134. 	// データの桁数の桁数を表す型
  135. 	typedef int expn_t;
  136.  
  137. 	// データの桁数を表す型
  138. 	typedef int64_t pow_t;
  139.  
  140. #ifdef SHOW_SIN_TABLE
  141. 	size_t g_sin_tblent_cnt = 0;
  142. #endif
  143.  
  144. 	/// 1の原始n乗根を表す。
  145. 	/// nの剰余計算を高速化するため、nは2^pに制限する。
  146. 	class basic_nth_root
  147. 	{
  148. 		/*const*/ expn_t sv_p;
  149. 		/*const*/ pow_t t360;
  150. 		/*const*/ pow_t m360;
  151. 		/*const*/ pow_t t090;
  152. 		std::vector<double> m_sin;
  153.  
  154. 	public:
  155. 		basic_nth_root() : sv_p(0), t360(0), m360(0), t090(0) { }
  156. 		basic_nth_root(const expn_t p)
  157. 			: sv_p(p), t360((pow_t)1 << p), m360(t360 - 1), t090(t360 / 4), m_sin((size_t)(t360 / 4 + 1))
  158. 		{
  159. 			if (p >= 2) {
  160. 				assert(t360 >= 4);
  161. 				t090 = t360 / 4;
  162. 				const double k = 0.5 * M_PI;
  163. 				for (pow_t i = 1; i < t090; i++) {
  164. 					m_sin[i] = sin((k * (double)i) / (double)t090);
  165. 				}
  166. #ifdef SHOW_SIN_TABLE
  167. 				g_sin_tblent_cnt += t090;
  168. 				printf("******* %s: p=%d table made. (+%zu, %zu)\n", __FUNCTION__, p, t090, g_sin_tblent_cnt);
  169. #endif
  170. 				m_sin[0] = 0.0;
  171. 				m_sin[t090] = 1.0;
  172. 			}
  173. 		}
  174.  
  175. 		expn_t p() const { return sv_p; }
  176.  
  177. 		pow_t n() const { return t360; }
  178.  
  179. 		/// 1の原始n乗根のk乗を取得する。
  180. 		void pow(const pow_t k_, double& re, double& im) const {
  181. 			assert(k_ >= 0);
  182. 			pow_t k = k_;
  183. 			k &= m360;		// nの剰余
  184. 			if (sv_p < 2) {
  185. 				assert(k == 0 || k == 1);
  186. 				im = 0.0;
  187. 				re = (k == 0) ? 1.0 : -1.0;
  188. 				return;
  189. 			}
  190. 			assert(0 <= k && k < 4 * t090);
  191. 			if (k <= t090) {
  192. 				re = m_sin[t090 - k];
  193. 				im = m_sin[k];
  194. 			}
  195. 			else if (k <= 2 * t090) {
  196. 				re = -m_sin[k - t090];			// = t090 - (2 * t090 - k)
  197. 				im = m_sin[2 * t090 - k];
  198. 			}
  199. 			else if (k <= 3 * t090) {
  200. 				re = -m_sin[3 * t090 - k];		// = t090 - (k - 2 * t090)
  201. 				im = -m_sin[k - 2 * t090];
  202. 			}
  203. 			else {
  204. 				re = m_sin[k - 3 * t090];		// = t090 - (4 * t090 - k)
  205. 				im = -m_sin[4 * t090 - k];
  206. 			}
  207. 		}
  208.  
  209. 	};
  210.  
  211. 	/// 1の原始n乗根を表す。
  212. 	class nth_root
  213. 	{
  214. 		expn_t max_p;
  215.  
  216. 		// 1の原始2^p乗根のr乗を収めたテーブルのベクタ
  217. 		// basic_tables[p]が1の原始2^p乗根のテーブルを表す。
  218. 		std::vector<basic_nth_root> basic_tables;
  219.  
  220. 		// 1の原始2^(p+m)乗根のみを収めたベクタ
  221. 		// high_order_root_tables[m]が1の原始2^(p+m)乗根の値。(こちらはベクタ要素がテーブルではなく値を表す。)
  222. 		std::vector<std::pair<double, double> > high_order_root_values;
  223.  
  224. 		/*
  225. 			(NOTE)
  226. 			1の原始n乗根をw_nと書くことにすると、{ b_n }を0か1からなる数列、rを適当な整数として
  227. 			　　k = b_0 + 2 * b_1 + (2^2) * b_2 + ... * (2^(m-1)) * b_(m-1) + { (2^m) * r }
  228. 			であるとき、
  229. 				(w_(2^(p+m)))^k = w_(2^(p+m))^b_0 * w_(2^(p+(m-1)))^b_1 * w_(2^(p+(m-2))^b_2 * ... * w_(2^(p+1))^b_(m-1) * { w_p^(rの2^pに対する剰余) }
  230. 			と書けるので、
  231. 				w_p^r (r=0..(2^p-1)) の値
  232. 				w_2^(p+q)  (q=1..m) の値
  233. 			があれば、(w^(p+m))^k を計算できる。
  234. 			前者は任意の指数rについて計算できねばならないため、テーブルで実現する場合はO(2^p)オーダーのメモリを要するが、
  235. 			後者はm個の値だけ記憶しておけば良いため、O(m)のオーダーで済む。そのかわりkが与えられる都度値の掛け算が発生する。
  236. 		*/
  237.  
  238. 	public:
  239. 		nth_root() : max_p(0) { }
  240.  
  241. 		/// 1の原始2^p乗根のテーブルを用意する。
  242. 		void alloc_table(const expn_t p) {
  243. 			if (p > max_p) {
  244. 				// (テーブル生成済のpより大きいp)
  245. 				// まずテーブル生成を検討する。
  246. 				if (p <= BASIC_SINTBL_P_MAX) {
  247. 					// (pが制限値以下)
  248. 					// [max_p+1..p]のテーブルを新規生成
  249. 					basic_tables.resize((size_t)(p + 1));
  250. 					for (expn_t i = max_p + 1; i <= p; i++) {
  251. 						basic_tables[i] = basic_nth_root(i);
  252. 					}
  253. 					max_p = p;
  254. 				}
  255. 			}
  256.  
  257. 			if (p > max_p) {
  258. 				// (任意指数に対するテーブル化が不能な大きさのp)
  259. 				// 指数k=1の値のみ記憶する。
  260. 				const expn_t m = p - max_p;
  261. 				if ((size_t)m >= high_order_root_values.size()) {
  262. 					high_order_root_values.resize((size_t)(m + 1));
  263. 					std::pair<double, double>& ent = high_order_root_values[m];		// Alias
  264. 					const pow_t big_n = (pow_t)1 << p;
  265. 					const double arg = (2 * M_PI) / big_n;
  266. 					ent.first = cos(arg);		// 実部
  267. 					ent.second = sin(arg);		// 虚部
  268. 				}
  269. 			}
  270. 		}
  271.  
  272. 		/// 1の原始2^p乗根のk乗を取得する。
  273. 		void pow(const expn_t p, const pow_t k, double& re, double& im) const {
  274. 			if (p <= max_p) {
  275. 				// (w_(2^p))^k算出
  276. 				basic_tables[p].pow(k, re, im);
  277. 			}
  278. 			else {
  279. 				const expn_t m = p - max_p;
  280. 				const pow_t r = k >> m;
  281. 				const pow_t b = k & (((pow_t)1 << m) - 1);
  282.  
  283. 				// (w_(2^max_p))^r算出
  284. 				basic_tables[max_p].pow(r, re, im);
  285.  
  286. 				// w_(2^(max_p + i))を掛け合わせる
  287. 				pow_t mask = (pow_t)1 << (m - 1);
  288. 				for (expn_t i = 1; i <= m; i++) {
  289. 					assert(mask > 0);
  290. 					if ((b & mask) != 0) {
  291. 						const std::pair<double, double>& ent = high_order_root_values[i];		// Alias
  292. 						const double tmp_re = re * ent.first - im * ent.second;
  293. 						const double tmp_im = re * ent.second + im * ent.first;
  294. 						re = tmp_re;
  295. 						im = tmp_im;
  296. 					}
  297. 					mask >>= 1;
  298. 				}
  299. 				assert(mask == 0);
  300. 			}
  301. 		}
  302.  
  303. 	};
  304.  
  305. }	// namespace priroot
  306.  
  307. namespace fftn
  308. {
  309. 	using namespace priroot;
  310.  
  311. 	/// ビット数bitwのBTR演算
  312. 	pow_t btr(const pow_t x, const expn_t bitw) {
  313. 		pow_t y = 0;
  314. 		for (expn_t i = 0; i <= (bitw - 1) / 2; i++) {
  315. 			const pow_t b = ((pow_t)1 << i) & x;
  316. 			y |= (b << ((bitw - 1) - 2 * i));
  317. 		}
  318. 		for (expn_t i = (bitw - 1) / 2 + 1; i < bitw; i++) {
  319. 			const pow_t b = ((pow_t)1 << i) & x;
  320. 			y |= (b >> (2 * i - (bitw - 1)));
  321. 		}
  322. 		return y;
  323. 	}
  324.  
  325. 	/// フーリエ変換
  326. 	/// DFTの原理式
  327. 	///   X(k) = Σ[j=0..t360-1]{ x(j) * w^(j*k) }  (k=0..t360-1)
  328. 	/// より
  329. 	///   X(k) = Σ[j=0..t180-1]{ x(2*j) * w^(2*j*k) } + Σ[j=0..t180-1]{ x(2*j+1) * w^((2*j+1)*k) }
  330. 	///        = Σ[j=0..t180-1]{ x(2*j) * w^(2*j*k) } + w^k * Σ[j=0..t180-1]{ x(2*j+1) * w^(2*j*k) }
  331. 	/// さらに、kをk+t180に置き換えると以下が成立する。
  332. 	///   X(k+t180) = Σ[j=0..t180-1]{ x(2*j) * w^(2*j*(k + t180) } + w^(k+t180) * Σ[j=0..t180-1]{ x(2*j+1) * w^(2*j*(k+t180) }
  333. 	///             = Σ[j=0..t180-1]{ x(2*j) * w^(2*j*k)         } - w^k        * Σ[j=0..t180-1]{ x(2*j+1) * w^(2*j*k) }
  334. 	/// よって、
  335. 	///   X_even(k) = Σ[j=0..t180-1]{ x(2*j  ) * w^(2*j*k) } (k=0..t180-1)
  336. 	///   X_odd(k)  = Σ[j=0..t180-1]{ x(2*j+1) * w^(2*j*k) } (k=0..t180-1)
  337. 	/// として、X(k) (k=0..t360)は以下のバタフライ演算で表せる。
  338. 	/// 　X(k     ) = X_even(k) + w^k * X_odd(k) (k=0..t180-1)
  339. 	///   X(k+t180) = X_even(k) - w^k * X_odd(k) (k=0..t180-1)
  340. 	/// これを、大元の時間領域データx(j)の実部がbuf[2*j]、虚部がbuf[2*j + 1]に格納されている前提で行う。
  341. 	void fouri(
  342. 		std::vector<double>::iterator buf,
  343. 		const expn_t p,
  344. 		nth_root& w)
  345. 	{
  346. 		// 再帰終端判定
  347. 		if (p == 0) {
  348. 			// (要素数2^0のフーリエ変換)
  349. 			// 単一要素x(j)のフーリエ変換はx(j)なので何もしない。
  350. 			return;
  351. 		}
  352.  
  353. 		const pow_t t360 = (pow_t)1 << p;
  354. 		const pow_t t180 = (pow_t)1 << (p - 1);
  355.  
  356. 		w.alloc_table(p);
  357.  
  358. 		// X_even(k)取得
  359. 		// [0..t180-1]に対し、w^2kを回転因子とするフーリエ変換をin-placeで実行する。
  360. 		fouri(buf, p - 1, w);
  361.  
  362. 		// X_odd(k)取得
  363. 		// [t180..t360-1]に対し、w^2kを回転因子とするフーリエ変換をin-placeで実行する。
  364. 		// t180に2を掛けているのは、実部と虚部の2要素セットでt001であることによる。
  365. 		fouri(buf + 2 * t180, p - 1, w);
  366.  
  367. 		// バタフライ演算
  368. 		// X(k     ) = X_even(k) + w^k * X_odd(k)
  369. 		// X(k+t180) = X_even(k) - w^k * X_odd(k)
  370. 		for (pow_t k = 0; k < t180; k++) {
  371. 			const pow_t ofs_even = k << 1;
  372. 			const pow_t ofs_odd = (k + t180) << 1;
  373.  
  374. 			// X_even(k)
  375. 			const double X_even_k_re = buf[ofs_even];
  376. 			const double X_even_k_im = buf[ofs_even + 1];
  377.  
  378. 			// X_odd(k)
  379. 			double X_odd_k_re = buf[ofs_odd];
  380. 			double X_odd_k_im = buf[ofs_odd + 1];
  381.  
  382. 			// X_odd(k) * w^k 算出、ここでwは1の原始t360乗根
  383. 			double wk_re, wk_im;
  384. 			w.pow(p, k, wk_re, wk_im);
  385. 			const double re = X_odd_k_re * wk_re - X_odd_k_im * wk_im;
  386. 			const double im = X_odd_k_re * wk_im + X_odd_k_im * wk_re;
  387.  
  388. 			// X(k     ) = X_even(k) + w^k * X_odd(k) = X_even(k) + (re + i * im)
  389. 			buf[ofs_even] = X_even_k_re + re;
  390. 			buf[ofs_even + 1] = X_even_k_im + im;
  391.  
  392. 			// X(k+t180) = X_even(k) - w^k * X_odd(k) = X_even(k) - (re + i * im)
  393. 			buf[ofs_odd] = X_even_k_re - re;
  394. 			buf[ofs_odd + 1] = X_even_k_im - im;
  395. 		}
  396. 	}
  397.  
  398. 	/// Bit reverse orderにin-placeで並べ替える。
  399. 	void sort_by_btr(
  400. 		std::vector<double>::iterator buf,
  401. 		const expn_t p)
  402. 	{
  403. 		const pow_t t360 = (pow_t)1 << p;
  404.  
  405. 		for (pow_t i = 0; i < t360; i++) {
  406. 			const pow_t btr_i = btr(i, p);
  407.  
  408. 			// (btr_i, i)と(i, btr_i)の両方についてswapすると元に戻ってしまうので、
  409. 			// 片方のみswap対象とする。
  410. 			if (btr_i < i) {
  411. 				std::swap(buf[2 * btr_i], buf[2 * i]);				// 実部
  412. 				std::swap(buf[2 * btr_i + 1], buf[2 * i + 1]);		// 虚部
  413. 			}
  414. 		}
  415. 	}
  416.  
  417. 	/// 逆フーリエ変換
  418. 	/// フーリエ変換とは逆の以下のバタフライ演算を行う。
  419. 	/// 　x(k     ) = (x_even(k) + w^(-k) * x_odd(k)) / t360 (k=0..t180-1)
  420. 	///   x(k+t180) = (x_even(k) - w^(-k) * x_odd(k)) / t360 (k=0..t180-1)
  421. 	/// これを、大元の周波数領域データX(k)の実部がbuf[2*k]、虚部がbuf[2*k+1]に格納されている前提で行う。
  422. 	void fouriinv(
  423. 		std::vector<double>::iterator buf,
  424. 		const expn_t p,
  425. 		nth_root& w)
  426. 	{
  427. 		// 再帰終端判定
  428. 		if (p == 0) {
  429. 			// (要素数2^0のフーリエ変換)
  430. 			// 単一要素x(j)のフーリエ変換はx(j)なので何もしない。
  431. 			return;
  432. 		}
  433.  
  434. 		const pow_t t360 = (pow_t)1 << p;
  435. 		const pow_t t180 = (pow_t)1 << (p - 1);
  436.  
  437. 		w.alloc_table(p);
  438.  
  439. 		// x_even(k)取得
  440. 		// [0..t180-1]に対し、w^2kを回転因子とするフーリエ変換をin-placeで実行する。
  441. 		fouriinv(buf, p - 1, w);
  442.  
  443. 		// x_odd(k)取得
  444. 		// [t180..t360-1]に対し、w^2kを回転因子とするフーリエ変換をin-placeで実行する。
  445. 		// t180に2を掛けているのは、実部と虚部の2要素セットでt001であることによる。
  446. 		fouriinv(buf + 2 * t180, p - 1, w);
  447.  
  448. 		// バタフライ演算
  449. 		// x(k     ) = x_even(k) + w^k * x_odd(k)
  450. 		// x(k+t180) = x_even(k) - w^k * x_odd(k)
  451. 		for (pow_t k = 0; k < t180; k++) {
  452. 			const pow_t ofs_even = k << 1;
  453. 			const pow_t ofs_odd = (k + t180) << 1;
  454.  
  455. 			// x_even(k)
  456. 			const double x_even_k_re = buf[ofs_even];
  457. 			const double x_even_k_im = buf[ofs_even + 1];
  458.  
  459. 			// x_odd(k)
  460. 			double x_odd_k_re = buf[ofs_odd];
  461. 			double x_odd_k_im = buf[ofs_odd + 1];
  462.  
  463. 			// x_odd(k) * w^(-k) 算出、ここでwは1の原始t360乗根
  464. 			double wk_re, wk_im;
  465. 			w.pow(p, k, wk_re, wk_im);
  466. 			const double re = x_odd_k_re * wk_re + x_odd_k_im * wk_im;
  467. 			const double im = x_odd_k_im * wk_re - x_odd_k_re * wk_im;
  468.  
  469. 			// X(k     ) = x_even(k) + w^k * x_odd(k) = x_even(k) + (re + i * im)
  470. 			buf[ofs_even] = (x_even_k_re + re) / 2.0;
  471. 			buf[ofs_even + 1] = (x_even_k_im + im) / 2.0;
  472.  
  473. 			// X(k+t180) = x_even(k) - w^k * x_odd(k) = x_even(k) - (re + i * im)
  474. 			buf[ofs_odd] = (x_even_k_re - re) / 2.0;
  475. 			buf[ofs_odd + 1] = (x_even_k_im - im) / 2.0;
  476. 		}
  477. 	}
  478.  
  479. }	// namespace fftn
  480.  
  481. // ブロックSix-Step FFTアルゴリズム
  482. // https://w...content-available-to-author-only...c.jp/public/VOL9/special/7.pdf
  483. // のC++実装
  484. // 実装の都合により、行数n1、列数n2、ブロックの行数nbは2の冪に制限する。
  485. // また、PCスタンドアロン実行用に行列転置処理を簡素化している。
  486. namespace fftb6
  487. {
  488. 	using namespace priroot;
  489.  
  490. 	/// 最適な行数(2の冪)の指数を返す。
  491. 	expn_t get_best_expn_n1(const expn_t p)
  492. 	{
  493. 		const expn_t exp_expn_n1 = (p + 1) / 2;
  494. 		return std::min(exp_expn_n1, EXPN_CAPACITY);
  495. 	}
  496.  
  497. 	/// 最適なブロックサイズ(2の冪)の指数を返す。
  498. 	expn_t get_best_expn_nb(const expn_t expn_n)
  499. 	{
  500. 		return std::min(expn_n, (EXPN_CAPACITY > expn_n) ? EXPN_CAPACITY - expn_n : 0);
  501. 	}
  502.  
  503. 	/// 転置処理
  504. 	/// 転置元配列src[]は (2^expn_n1)行(2^(p-expn_n1))列の2次元配列
  505. 	/// 転置先配列dst[]は (2^(p-expn_n1))行(2^expn_n1)列の2次元配列
  506. 	/// とそれぞれ解釈する。
  507. 	void transpose(
  508. 		std::vector<double>::const_iterator src,
  509. 		std::vector<double>::iterator dst,
  510. 		const expn_t p,
  511. 		const expn_t expn_n1)
  512. 	{
  513. 		const expn_t expn_n2 = p - expn_n1;
  514.  
  515. 		const pow_t t360 = (pow_t)1 << p;
  516. 		const pow_t n1 = (pow_t)1 << expn_n1;
  517. 		const pow_t n2 = (pow_t)1 << expn_n2;
  518.  
  519. 		// ストレートな転置アルゴリズム
  520. 		// 正方行列の転置への分解はPCスタンドアロン実行に向かないため行わない。
  521. #pragma omp parallel for
  522. 		for (pow_t i = 0; i < n1; i++) {
  523. 			for (pow_t j = 0; j < n2; j++) {
  524. 				const pow_t src_idx = n2 * i + j;
  525. 				const pow_t dst_idx = n1 * j + i;
  526.  
  527. 				dst[2 * dst_idx] = src[2 * src_idx];				// 実部
  528. 				dst[2 * dst_idx + 1] = src[2 * src_idx + 1];		// 虚部
  529. 			}
  530. 		}
  531. 	}
  532.  
  533. 	// 前方宣言
  534. 	void fouri_wrp(
  535. 		std::vector<double>::iterator buf,
  536. 		const expn_t p,
  537. 		nth_root& w);
  538. 	void fouriinv_wrp(
  539. 		std::vector<double>::iterator buf,
  540. 		const expn_t p,
  541. 		nth_root& w);
  542.  
  543. 	/// フーリエ変換
  544. 	void fourib6(
  545. 		std::vector<double>::iterator buf,
  546. 		const expn_t p,
  547. 		nth_root& w,
  548. 		const expn_t expn_n1)
  549. 	{
  550. 		// 行数が全体のサイズ以下であること
  551. 		if (!(0 <= expn_n1 && expn_n1 <= p)) {
  552. 			// ここに来たら呼び出し元のバグ
  553. 			abort();
  554. 		}
  555.  
  556. 		// 列数(2の冪)の指数
  557. 		const expn_t expn_n2 = p - expn_n1;
  558.  
  559. 		// ブロックサイズ
  560. 		const expn_t expn_nb = get_best_expn_nb(expn_n2);
  561.  
  562. 		// buf[0..t360-1]を複素数を要素とする2次元配列とみなす。
  563. 		// 入力に関しては j=n2*j1 + j2 (0≦j2＜n2) として、buf[j1, j2] (n1行n2列)
  564. 		// 出力に関しては k=n1*k2 + k1 (0≦k1＜n1) として、buf[k1, k2] (n1行n2列)
  565. 		// 出力の順序とメモリのアドレス順を合わせるには、buf[k1, k2]を転置してn2行n1列にする必要がある。
  566. 		const pow_t t360 = (pow_t)1 << p;
  567. 		const pow_t n1 = (pow_t)1 << expn_n1;
  568. 		const pow_t n2 = (pow_t)1 << expn_n2;
  569. 		const pow_t nb = (pow_t)1 << expn_nb;
  570.  
  571. 		// SINテーブル獲得
  572. 		// ここで指数pのテーブルを確保しておくことで、
  573. 		// マルチスレッド処理中にテーブルのstd::vector::realloc()が呼ばれなくして
  574. 		// 排他制御を省略できる。
  575. 		w.alloc_table(p);
  576.  
  577. 		// 6 step algorithm
  578. #pragma omp parallel for
  579. 		for (pow_t JJ = 0; JJ < n2; JJ += nb) {
  580. 			// 作業バッファ準備
  581. 			// 2を掛けているのは、実部と虚部の2要素セット1要素であることによる。
  582. 			// スレッド毎に確保する必要があるため、parallel for対象ブロック内で実施する。
  583. 			std::vector<double> work((size_t)(2 * nb * n1));
  584.  
  585. 			// ストレートな転置アルゴリズム
  586. 			// 正方行列の転置への分解はPCスタンドアロン実行に向かないため行わない。
  587. 			for (pow_t j = JJ; j < JJ + nb; j++) {
  588. 				for (pow_t i = 0; i < n1; i++) {
  589. 					// ここでのi、jはj1、j2にあたる。
  590. 					// 直下の処理は転置処理
  591. 					//   work[j2 - JJ, j1] = buf[j1, j2]
  592. 					// にあたる。
  593. 					const pow_t src_idx = n2 * i + j;
  594. 					const pow_t dst_idx = n1 * (j - JJ) + i;
  595. 					work[2 * dst_idx] = buf[2 * src_idx];				// 実部
  596. 					work[2 * dst_idx + 1] = buf[2 * src_idx + 1];		// 虚部
  597. 				}
  598. 			}
  599.  
  600. 			// Step 2: nb組のn1点multi row FFT
  601. 			for (pow_t i = 0; i < nb; i++) {
  602. 				// ここでのiはj2-JJに対応する。
  603. 				// 直下の処理はwork[j2 - JJ, j1]のj1に関するフーリエ変換
  604. 				//   work[j2 - JJ, k1] = Σ[j1=0..n1-1]{ work[j2 - JJ, j1] * (w_n1)^(j1 * k1) }
  605. 				//                     = Σ[j1=0..n1-1]{ buf[j1, j2] * (w_n1)^(j1 * k1) }
  606. 				// にあたる。ここで、w_n1は1の原始n1乗根。
  607. 				const pow_t row_head = 2 * n1 * i;
  608. 				fouri_wrp(work.begin() + row_head, expn_n1, w);
  609. 			}
  610.  
  611. 			// Step 3: ひねり係数をかける
  612. 			// Step 4: 転置
  613. 			for (pow_t j = JJ; j < JJ + nb; j++) {
  614. 				for (pow_t i = 0; i < n1; i++) {
  615. 					// ここでのi、jはk1、j2にあたる。
  616.  
  617. 					// 転置のsrcとdst
  618. 					const pow_t src_idx = n1 * (j - JJ) + i;	// n1 * (j2 - JJ) + k1
  619. 					const pow_t dst_idx = n2 * i + j;			// n2 * k1 + j2
  620.  
  621. 					// ひねり係数(w_n)^(j2*k1)
  622. 					// ここで、w_nは1の原始n乗根。
  623. 					double omg_re, omg_im;
  624. 					w.pow(p, j * i, omg_re, omg_im);
  625.  
  626. 					// ひねり係数を掛けて転置
  627. 					// 直下の処理は、ひねり係数を掛けて転置する処理
  628. 					//   buf[k1, j2] = work[j2 - JJ, k1] * (w_n)^(j2*k1)
  629. 					// にあたる。
  630. 					const double src_re = work[2 * src_idx];
  631. 					const double src_im = work[2 * src_idx + 1];
  632. 					buf[2 * dst_idx] = src_re * omg_re - src_im * omg_im;			// 実部
  633. 					buf[2 * dst_idx + 1] = src_im * omg_re + src_re * omg_im;		// 虚部
  634. 				}
  635. 			}
  636. 		}
  637.  
  638. 		// Setp 5: N1組のN2点multi row FFT
  639. 		// ブロックSix-step algorithmの引用元コードによれば
  640. 		// nb組のn2点multi row FFTに区分して行うべきところだが、
  641. 		// PCスタンドアロン実行に向かないため区分しない。
  642. #pragma omp parallel for
  643. 		for (pow_t i = 0; i < n1; i++) {
  644. 			// ここでのiはk1にあたる。
  645. 			// 直下の処理は、buf[k1, j2]のj2に関するフーリエ変換
  646. 			//   buf[k1, k2] = Σ[j2=0..n2-1]{ buf[k1, j2] * (w_n2)^(j2 * k2) }
  647. 			// にあたる。ここで、w_n2は1の原始n2乗根。
  648. 			const pow_t row_head = 2 * n2 * i;
  649. 			fouri_wrp(buf + row_head, expn_n2, w);
  650.  
  651. 			// Step 6: 転置
  652. 			// 本来はフーリエ変換結果X(k)がメモリ上で並ぶようにbuf[k1, k2]をn2行×n1列の2次元配列dst[,]に転置する工程だが、
  653. 			// buf[N2 * i + j] を X(N1 * j + i)と解釈することにして省略する。
  654. 		}
  655. 	}
  656.  
  657. 	/// フーリエ逆変換
  658. 	/// buf[]としては、fourib6()の結果を転置せずにそのまま与える前提。
  659. 	/// expn_n1もfourib6()に与えたものと同じ値を与える。
  660. 	void fourib6inv(
  661. 		std::vector<double>::iterator buf,
  662. 		const expn_t p,
  663. 		nth_root& w,
  664. 		const expn_t expn_n1)
  665. 	{
  666. 		// 行数が全体のサイズ以下であること
  667. 		if (!(0 <= expn_n1 && expn_n1 <= p)) {
  668. 			// ここに来たら呼び出し元のバグ
  669. 			abort();
  670. 		}
  671.  
  672. 		// 列数(2の冪)の指数
  673. 		const expn_t expn_n2 = p - expn_n1;
  674.  
  675. 		// ブロックサイズ
  676. 		const expn_t expn_nb = get_best_expn_nb(expn_n2);
  677.  
  678. 		// buf[0..t360-1]を複素数を要素とする2次元配列とみなす。
  679. 		// 入力に関しては k=n1*k2 + k1 (0≦k1＜n1) として、buf[k1, k2] (n1行n2列)
  680. 		// 出力に関しては j=n2*j1 + j2 (0≦j2＜n2) として、buf[j1, j2] (n1行n2列)
  681. 		// 入力の順序とメモリのアドレス順を合わせるには、buf[k1, k2]を転置してn2行n1列にする必要があるが
  682. 		// 当関数ではその転置を行わずに済ませる。
  683. 		// 出力の順序とメモリのアドレス順は一致する。
  684. 		const pow_t t360 = (pow_t)1 << p;
  685. 		const pow_t n1 = (pow_t)1 << expn_n1;
  686. 		const pow_t n2 = (pow_t)1 << expn_n2;
  687. 		const pow_t nb = (pow_t)1 << expn_nb;
  688.  
  689. 		// SINテーブル獲得
  690. 		// ここで指数pのテーブルを確保しておくことで、
  691. 		// マルチスレッド処理中にテーブルのstd::vector::realloc()が呼ばれなくして
  692. 		// 排他制御を省略できる。
  693. 		w.alloc_table(p);
  694.  
  695. 		// Modified 6 step algorithm
  696. 		// Setp 1: n1組のn2点multi row IFFT
  697. 		// ブロックSix-step algorithmの引用元コードによれば
  698. 		// nb組のn2点multi row FFTに区分して行うべきところだが、
  699. 		// PCスタンドアロン実行に向かないため区分しない。
  700. #pragma omp parallel for
  701. 		for (pow_t i = 0; i < n1; i++) {
  702. 			// ここでのiはk1にあたる。
  703. 			// 直下の処理はbuf[k1, k2]のk2に関する逆フーリエ変換
  704. 			//   buf[k1, j2] = (1/n2) * Σ[k2=0..n2-1]{ buf[k1, k2] * (w_n2)^(-(j2 * k2)) }
  705. 			// にあたる。ここで、w_n2は1の原始n2乗根。
  706. 			const pow_t row_head = 2 * n2 * i;
  707. 			fouriinv_wrp(buf + row_head, expn_n2, w);
  708. 		}
  709.  
  710. 		// Step 3: ひねり係数をかける
  711. 		// Step 4: 転置
  712. #pragma omp parallel for
  713. 		for (pow_t JJ = 0; JJ < n2; JJ += nb) {
  714. 			// 作業バッファ準備
  715. 			// 2を掛けているのは、実部と虚部の2要素セット1要素であることによる。
  716. 			std::vector<double> work((size_t)(2 * nb * n1));
  717.  
  718. 			// ストレートな転置アルゴリズム
  719. 			// 正方行列の転置への分解はPCスタンドアロン実行に向かないため行わない。
  720. 			for (pow_t j = JJ; j < JJ + nb; j++) {
  721. 				for (pow_t i = 0; i < n1; i++) {
  722. 					// ここでのi、jはk1、j2にあたる。
  723.  
  724. 					// 転置のsrcとdst
  725. 					const pow_t src_idx = n2 * i + j;
  726. 					const pow_t dst_idx = n1 * (j - JJ) + i;
  727.  
  728. 					// ひねり係数の逆数(w_n)^(j2*k1)
  729. 					// ここで、w_nは1の原始n乗根。
  730. 					double omg_re, omg_im;
  731. 					w.pow(p, j * i, omg_re, omg_im);
  732.  
  733. 					// ひねり係数を掛けて転置
  734. 					// 直下の処理は
  735. 					//   work[j2 - JJ, k1] = buf[k1, j2] * (w_n)^(-(j2*k1))
  736. 					// にあたる。
  737. 					const double src_re = buf[2 * src_idx];
  738. 					const double src_im = buf[2 * src_idx + 1];
  739. 					work[2 * dst_idx] = src_re * omg_re + src_im * omg_im;			// 実部
  740. 					work[2 * dst_idx + 1] = src_im * omg_re - src_re * omg_im;		// 虚部
  741. 				}
  742. 			}
  743.  
  744. 			// nb組のn1点multi row IFFT
  745. 			for (pow_t i = 0; i < nb; i++) {
  746. 				// ここでのiはj2-JJにあたる。
  747. 				// 直下の処理は、work[j2 - JJ, k1]のk1に関する逆フーリエ変換
  748. 				//   work[j2 - JJ, j2] = (1/n1) * Σ[k1=0..n1-1]{ work[j2 - JJ, k1] * (w_n1)^(j1 * k1) }
  749. 				//                     = (1/n1) * Σ[k1=0..n1-1]{ buf[k1, j2] * (w_n1)^(j1 * k1) }
  750. 				// にあたる。ここで、w_n1は1の原始n1乗根。
  751. 				const pow_t row_head = 2 * n1 * i;
  752. 				fouriinv_wrp(work.begin() + row_head, expn_n1, w);
  753. 			}
  754.  
  755. 			// IFFT結果をbuf[]の元の場所に書き戻す
  756. 			// ストレートな転置アルゴリズム
  757. 			// 正方行列の転置への分解はPCスタンドアロン実行に向かないため行わない。
  758. 			for (pow_t j = JJ; j < JJ + nb; j++) {
  759. 				for (pow_t i = 0; i < n1; i++) {
  760. 					// ここでのi、jはj1、j2にあたる。
  761. 					// 直下の処理は、転置処理
  762. 					//   buf[j1, j2] = work[j2 - JJ, j1]
  763. 					// にあたる。
  764. 					const pow_t src_idx = n1 * (j - JJ) + i;
  765. 					const pow_t dst_idx = n2 * i + j;
  766. 					buf[2 * dst_idx] = work[2 * src_idx];				// 実部
  767. 					buf[2 * dst_idx + 1] = work[2 * src_idx + 1];		// 虚部
  768. 				}
  769. 			}
  770. 		}
  771. 	}
  772.  
  773. 	/// フーリエ変換のwrapper
  774. 	void fouri_wrp(
  775. 		std::vector<double>::iterator buf,
  776. 		const expn_t p,
  777. 		nth_root& w)
  778. 	{
  779. 		if (p <= EXPN_CAPACITY) {
  780. 			// (FFT結果がL2キャッシュサイズの半分に収まる)
  781. 			// ストレートにfftする。
  782. 			fftn::sort_by_btr(buf, p);
  783. 			fftn::fouri(buf, p, w);
  784. 		}
  785. 		else {
  786. 			// (FFT結果がL2キャッシュサイズの半分に収まならい)
  787. 			// Blocked 6 step algorithmでFFTする
  788. 			// n1 = √(2^p)が最も効率が良いのでEXPN_CACACITYを超えない範囲で近い値にする。
  789. 			// (2^p)/n1がEXPN_CAPACITYを超える場合は、foourib6()から当関数にさらに再帰して細分化する。
  790. 			const int n1 = std::min((p + 1) / 2 + 1, EXPN_CAPACITY);
  791. 			// フーリエ変換
  792. 			fourib6(buf, p, w, n1);
  793. 			// 転置
  794. 			const pow_t t360 = (pow_t)1 << p;
  795. 			std::vector<double> tmp(buf, buf + 2 * t360);
  796. 			transpose(tmp.begin(), buf, p, n1);
  797. 		}
  798. 	}
  799.  
  800. 	/// 逆フーリエ変換のwrapper
  801. 	void fouriinv_wrp(
  802. 		std::vector<double>::iterator buf,
  803. 		const expn_t p,
  804. 		nth_root& w)
  805. 	{
  806. 		if (p <= EXPN_CAPACITY) {
  807. 			// (IFFT結果がL2キャッシュサイズの半分に収まる)
  808. 			// ストレートにifftする。
  809. 			fftn::sort_by_btr(buf, p);
  810. 			fftn::fouriinv(buf, p, w);
  811. 		}
  812. 		else {
  813. 			// (IFFT結果がL2キャッシュサイズの半分に収まならい)
  814. 			// Blocked 6 step algorithmでFFTする
  815. 			// n1 = √(2^p)が最も効率が良いのでEXPN_CACACITYを超えない範囲で近い値にする。
  816. 			// (2^p)/n1がEXPN_CAPACITYを超える場合は、foourib6()から当関数にさらに再帰して細分化する。
  817. 			const int n1 = std::min((p + 1) / 2 + 1, EXPN_CAPACITY);
  818. 			const pow_t t360 = (pow_t)1 << p;
  819. 			// 転置
  820. 			std::vector<double> tmp(buf, buf + (size_t)(2 * t360));
  821. 			transpose(tmp.begin(), buf, p, n1);
  822. 			// 逆フーリエ変換
  823. 			fourib6inv(buf, p, w, p - n1);
  824. 		}
  825. 	}
  826.  
  827. }	// namespace fftb6
  828.  
  829. namespace fftb6_4n
  830. {
  831. 	using namespace priroot;
  832.  
  833. 	/// 4n乗根によるフーリエ変換(Blocked 6-step algorithm版)
  834. 	/// x(j) (j=0..(2*t360)-1)から
  835. 	///   x'(j) = (x(j) + i * x(j+t360)) * u^(j/4) (j=0..t360-1)
  836. 	/// を構成し、x'(j)をDFTする。ここで、uは1の原始(4 * mag * t360)乗根 (mag ≧ 1)。
  837. 	void fourib6_4n(
  838. 		std::vector<double>::iterator buf,
  839. 		const expn_t p,
  840. 		nth_root& w,
  841. 		const expn_t expn_n1)
  842. 	{
  843. 		const pow_t t360 = (pow_t)1 << p;
  844.  
  845. 		w.alloc_table(p + 2);
  846.  
  847. 		// { x'(j) }構成
  848. 		// { x(j) + i * x(j+t360) } にu^(j/4)を掛ける。
  849. 		const expn_t p2 = p + 2;
  850. 		for (pow_t j = 0; j < t360; j++) {
  851. 			const double a = buf[2 * j];
  852. 			const double b = buf[2 * j + 1];
  853.  
  854. 			// u = w^4としてu^(j/4) = w^(j)
  855. 			double wj_re, wj_im;
  856. 			w.pow(p2, j, wj_re, wj_im);
  857. 			buf[2 * j] = a * wj_re - b * wj_im;
  858. 			buf[2 * j + 1] = a * wj_im + b * wj_re;
  859. 		}
  860.  
  861. 		// フーリエ変換
  862. 		fftb6::fourib6(buf, p, w, expn_n1);
  863. 	}
  864.  
  865. 	/// 4n乗根による逆フーリエ変換(Blocked 6-step algorithm版)
  866. 	void fourib6_4ninv(
  867. 		std::vector<double>::iterator buf,
  868. 		const expn_t p,
  869. 		nth_root& w,
  870. 		const expn_t expn_n1)
  871. 	{
  872. 		const pow_t t360 = (pow_t)1 << p;
  873.  
  874. 		w.alloc_table(p + 2);
  875.  
  876. 		// 逆フーリエ変換
  877. 		fftb6::fourib6inv(buf, p, w, expn_n1);
  878.  
  879. 		// { x(j) }復元
  880. 		const expn_t p2 = p + 2;
  881. 		for (pow_t j = 0; j < t360; j++) {
  882. 			const double a = buf[2 * j];
  883. 			const double b = buf[2 * j + 1];
  884.  
  885. 			// u = w^4としてu^(-j/4) = w^(-j)を掛ける
  886. 			double wj_re, wj_im;
  887. 			w.pow(p2, j, wj_re, wj_im);
  888. 			buf[2 * j] = a * wj_re + b * wj_im;
  889. 			buf[2 * j + 1] = b * wj_re - a * wj_im;
  890. 		}
  891. 	}
  892.  
  893. }	// namespace fft4n
  894.  
  895. namespace oper
  896. {
  897. 	using namespace fftb6_4n;
  898.  
  899. 	struct I
  900. 	{
  901. 		std::vector<I_TYPE> data;	// 基数RADIXの多倍長データ
  902. 		size_t ndigits;				// dataに入っている有効な桁数
  903.  
  904. 		I(const size_t cap)
  905. 			: data(cap), ndigits((size_t)0)
  906. 		{
  907. 			/*pass*/
  908. 		}
  909. 	};
  910.  
  911. 	struct F
  912. 	{
  913. 		std::vector<double> data;	// 基数[-(RADIX/2), RADIX/2)の多倍長データ
  914. 		int p;						// 2^p == (dataの要素数) となるp
  915.  
  916. 		F(const int p_)
  917. 			: data((size_t)1 << p_), p(0)
  918. 		{
  919. 			/*pass*/
  920. 		}
  921. 	};
  922.  
  923. 	/// 整数のフーリエ変換
  924. 	void i2f(F& a, const I& b, const int p, priroot::nth_root& w)
  925. 	{
  926. 		size_t i;
  927.  
  928. 		// 次数2^(p-1) = t360の4n乗根フーリエ変換を行う。
  929. 		// 結果は実部と虚部を合わせて2^p = (2 * t360)桁分となる。
  930. 		const pow_t t360 = (pow_t)1 << (p - 1);
  931.  
  932. 		assert(b.ndigits <= (size_t)t360);
  933. 		assert((size_t)(2 * t360) <= a.data.size());
  934.  
  935. #ifdef SHOW_I
  936. 		printf("%s: Before pre tr.:\n", __FUNCTION__);
  937. 		show_re_im_seq(b.data, p);
  938. #endif
  939.  
  940. 		// フーリエ変換の前処理
  941. 		// (bが表す整数) = Σ[i=0..a.ndigits-1]{ b.data[i] * RADIX^i } (0≦b.data[i]＜RADIX)
  942. 		// を、
  943. 		// (aが表す整数) = Σ[i=0..a.ndigits-1]{ a.data[i] * RADIX^i } (-RADIX / 2≦a.data[i]＜RADIX / 2)
  944. 		// に変換する。(a, bが同じ整数を表すように桁上げも行う。)
  945. 		a.p = p;
  946. 		int carry = 0;
  947. 		for (i = 0; i < b.ndigits - 1; i++) {
  948. 			// 実部虚部隣接化
  949. 			// b.data[t360..]を虚部に配置する。
  950. 			const size_t dst_i = (i < (size_t)t360) ? 2 * i : 2 * (i - (size_t)t360) + 1;
  951.  
  952. 			const int32_t t = b.data[i] + carry;	// 桁[i]の値
  953. 			if (t < (int32_t)RADIX / 2) {
  954. 				// (桁[i]の値∈[0, RADIX / 2) )
  955. 				// 何もしない
  956. 				a.data[dst_i] = t;
  957. 				carry = 0;
  958. 			}
  959. 			else {
  960. 				// (桁[i]の値∈[RADIX / 2, RADIX) )
  961. 				// 補正実施
  962. 				// まず桁[i] -= RADIX により桁[i]の値∈[0, RADIX / 2)とする。ただし
  963. 				// それだけでは(aが表す整数)が RADIX * RADIX^i = RADIX^(i+1)だけ減ってしまうので、
  964. 				// 減った分をcarryで一つ上の桁に回して補填する。
  965. 				a.data[dst_i] = (double)(t - (int32_t)RADIX);
  966. 				carry = 1;
  967. 			}
  968. 		}
  969. 		{
  970. 			// 実部虚部隣接化
  971. 			const size_t dst_i = (i < (size_t)t360) ? 2 * i : 2 * (i - (size_t)t360) + 1;
  972.  
  973. 			// 最上位桁
  974. 			a.data[dst_i] = (double)(b.data[i] + (I_TYPE)carry);
  975. 			i++;
  976. 		}
  977. 		for (; i < (size_t)(2 * t360); i++) {
  978. 			// 実部虚部隣接化
  979. 			const size_t dst_i = (i < (size_t)t360) ? 2 * i : 2 * (i - (size_t)t360) + 1;
  980.  
  981. 			a.data[dst_i] = 0.0;
  982. 		}
  983.  
  984. #if defined(SHOW_ITM) || defined(SHOW_I)
  985. 		printf("Before fouri4n():\n");
  986. 		show_re_im_seq(a.data, a.p);
  987. #endif
  988.  
  989. 		const expn_t expn_n2 = (p - 1) / 2;
  990. 		const expn_t expn_n1 = (p - 1) - expn_n2;
  991.  
  992. 		// フーリエ変換実施
  993. 		fftb6_4n::fourib6_4n(a.data.begin(), p - 1, w, expn_n1);
  994.  
  995. #ifdef SHOW_ITM
  996. 		printf("After fouri4n():\n");
  997. 		show_re_im_seq(a.data, a.p);
  998. #endif
  999. 	}
  1000.  
  1001. 	/// 整数のフーリエ変換結果の積(結果は整数)
  1002. 	void mul_f2i(I& a, const F& b, const F& c, F& d, double& err_max, priroot::nth_root& w) {
  1003. 		int64_t carry;
  1004. 		size_t i, t;
  1005.  
  1006. 		const int p = b.p;
  1007.  
  1008. 		// 次数2^(p-1) = t360の4n乗根フーリエ変換結果の積を逆フーリエ変換する。
  1009. 		// 結果は実部と虚部を合わせて2^p = (2 * t360)桁となる。
  1010. 		const pow_t t360 = (size_t)1 << (p - 1);
  1011.  
  1012. 		assert(p == c.p);
  1013. 		assert((size_t)(2 * t360) <= d.data.size());
  1014.  
  1015. 		// 積
  1016. 		for (i = 0; i < (size_t)t360; i++) {
  1017. 			const double b_i_re = b.data[2 * i];
  1018. 			const double b_i_im = b.data[2 * i + 1];
  1019. 			const double c_i_re = c.data[2 * i];
  1020. 			const double c_i_im = c.data[2 * i + 1];
  1021.  
  1022. 			// d[i] = b[i] * c[i]
  1023. 			const double d_i_re = b_i_re * c_i_re - b_i_im * c_i_im;
  1024. 			const double d_i_im = b_i_re * c_i_im + b_i_im * c_i_re;
  1025. 			d.data[2 * i] = d_i_re;
  1026. 			d.data[2 * i + 1] = d_i_im;
  1027. 		}
  1028.  
  1029. #ifdef SHOW_ITM
  1030. 		printf("Before fouri4ninv():\n");
  1031. 		show_re_im_seq(d.data, d.p);
  1032. #endif
  1033.  
  1034. 		const expn_t expn_n2 = (p - 1) / 2;
  1035. 		const expn_t expn_n1 = (p - 1) - expn_n2;
  1036.  
  1037. 		// 逆フーリエ変換
  1038. 		fftb6_4n::fourib6_4ninv(d.data.begin(), p - 1, w, expn_n1);
  1039.  
  1040. #ifdef SHOW_ITM
  1041. 		printf("After fouri4ninv():\n");
  1042. 		show_re_im_seq(d.data, d.p);
  1043. #endif
  1044.  
  1045. 		const size_t s = std::min(a.data.size(), (size_t)(2 * t360));
  1046.  
  1047. 		// Carry and fix処理
  1048. 		// ただし対象とする桁の並びd.data[]は、
  1049. 		// 桁の値が[0, RADIX)ではなく、[-RADIX / 2, RADIX / 2)である。
  1050. 		carry = 0;
  1051. 		t = 1;
  1052. 		//const double k = ldexp(1.0, -p + 1);		// k = 1.0 * 2^(-p + 1)
  1053. 		double k = 1.0;		// fouri4ninv()結果は2^pで割る必要無し。
  1054. 		for (i = 0; i < s; i++) {
  1055. 			// 実部虚部隣接化解除
  1056. 			// b.data[t360..]を虚部に配置する。
  1057. 			const size_t src_i = (i < (size_t)t360) ? 2 * i : 2 * (i - (size_t)t360) + 1;
  1058. 			//printf("t360=%lld, s=%zu, i=%zu, src_i=%zu\n", t360, s, i, src_i);
  1059.  
  1060. 			const double x = d.data[src_i] * k;
  1061. 			int64_t e = (int64_t)floor(x + 0.5);
  1062. 			err_max = std::max(err_max, fabs(x - e));
  1063. 			e += carry;
  1064. 			carry = e / RADIX;
  1065. 			int64_t f = e - carry * RADIX;
  1066. 			if (f < 0) {
  1067. 				f += RADIX;
  1068. 				carry--;
  1069. 			}
  1070. 			a.data[i] = (I_TYPE)f;
  1071. 			if (f != 0) {
  1072. 				t = i + 1;
  1073. 			}
  1074. 		}
  1075. 		a.ndigits = t;
  1076. 		if (carry != 0) {
  1077. 			if (a.ndigits < a.data.size()) {
  1078. 				assert(0 <= carry && carry < RADIX);
  1079. 				a.data[a.ndigits++] = (I_TYPE)carry;
  1080. 			}
  1081. 			else {
  1082. 				// ここに来たら桁あふれ
  1083. 				abort();
  1084. 			}
  1085. 		}
  1086. 		for (; i < s; i++) {
  1087. 			// 実部虚部隣接化解除
  1088. 			const size_t src_i = (i < (size_t)t360) ? 2 * i : 2 * i + 1;
  1089.  
  1090. 			const double x = d.data[src_i] * k;
  1091. 			const int64_t e = (int64_t)floor(x + 0.5);
  1092. 			err_max = std::max(err_max, fabs(x - e));
  1093. 			assert(e == 0);
  1094. 		}
  1095. 	}
  1096.  
  1097. 	/// 整数の和
  1098. 	void add_i(I& a, const I& b_, const I& c_) {
  1099. 		I_TYPE carry, d;
  1100. 		size_t i;
  1101.  
  1102. 		const I& b = (b_.ndigits > c_.ndigits) ? c_ : b_;	// Alias
  1103. 		const I& c = (b_.ndigits > c_.ndigits) ? b_ : c_;	// Alias
  1104. 		assert(b.ndigits <= c.ndigits);
  1105. 		assert(c.ndigits <= a.data.size());
  1106.  
  1107. 		carry = 0;
  1108. 		for (i = 0; i < b.ndigits; i++) {
  1109. 			d = b.data[i] + c.data[i] + carry;
  1110. 			carry = 0;
  1111. 			if (d >= RADIX) {
  1112. 				d -= RADIX;
  1113. 				carry++;
  1114. 			}
  1115. 			a.data[i] = d;
  1116. 		}
  1117. 		for (; i < c.ndigits; i++) {
  1118. 			d = c.data[i] + carry;
  1119. 			carry = 0;
  1120. 			if (d >= RADIX) {
  1121. 				d -= RADIX;
  1122. 				carry++;
  1123. 			}
  1124. 			a.data[i] = d;
  1125. 		}
  1126. 		a.ndigits = c.ndigits;
  1127. 		if (a.ndigits < a.data.size()) {
  1128. 			if (carry != 0) {
  1129. 				a.data[a.ndigits++] = carry;
  1130. 			}
  1131. 		}
  1132. 		else {
  1133. 			if (carry != 0) {
  1134. 				// ここに来たら桁あふれ
  1135. 				abort();
  1136. 			}
  1137. 		}
  1138. 	}
  1139.  
  1140. 	/// 整数のファイル出力
  1141. 	void print(const char* const filename, const I& a)
  1142. 	{
  1143. 		char format[10];
  1144. 		int nRet;
  1145.  
  1146. 		assert(a.data.size() > (size_t)0);
  1147. 		assert(a.ndigits <= a.data.size());
  1148.  
  1149. 		FILE* fp;
  1150. 		nRet = fopen_s_wrp(&fp, filename, "w");
  1151. 		if (nRet != 0) {
  1152. 			fflush(stdout);
  1153. 			fprintf(stderr, "ERROR: fopen_s_wrp() failed. (nRet=%d)\n", nRet);
  1154. 			abort();
  1155. 		}
  1156.  
  1157. 		sprintf_s_wrp(format, _countof(format), "%%0%du", FIG);
  1158.  
  1159. 		size_t s = a.ndigits;
  1160. 		fprintf(fp, "%lu", (unsigned long)a.data[--s]);
  1161. 		while (s != 0) {
  1162. 			fprintf(fp, format, a.data[--s]);
  1163. 		}
  1164. 		fclose(fp);
  1165. 	}
  1166.  
  1167. }	// namespace oper
  1168.  
  1169. namespace
  1170. {
  1171. 	// N ≧ log2(s1 * s2)を満たす最小のNを求める。
  1172. 	int s2p(const size_t s1, const size_t s2) {
  1173. 		int i;
  1174. 		for (i = 1; i < 64; i++) {
  1175. 			if (((size_t)1 << i) + 1 >= s1 + s2) {
  1176. 				// (2^i + 1 ≧ (s1の有効桁数) + (s2の有効桁数))
  1177. 				break;
  1178. 			}
  1179. 		}
  1180. 		return i;
  1181. 	}
  1182.  
  1183. }	// namespace
  1184.  
  1185. int main()
  1186. {
  1187. 	using namespace oper;
  1188.  
  1189. 	const int n = 34;
  1190.  
  1191. 	nth_root w;
  1192.  
  1193. 	const int max_p = s2p(((size_t)((1LLU << n) * log10(2)) / FIG + 1), 1);
  1194. 	const size_t cap = (size_t)1 << max_p;
  1195.  
  1196. 	I a(cap);
  1197. 	F x(max_p);
  1198.  
  1199. 	a.data[0] = 1;
  1200. 	a.ndigits = 1;
  1201.  
  1202. 	double err_max = 0.0;
  1203. #ifdef SHOW_CSM_TIME
  1204. 	std::vector<double> csm_time;
  1205. 	const double startTmg_0 = omp_get_wtime();
  1206. #endif
  1207. 	for (int i = 0; i < n; i++) {
  1208. #ifdef SHOW_CSM_TIME
  1209. 		const double startTmg = omp_get_wtime();
  1210. #endif
  1211.  
  1212. 		// 多倍長整数xのフーリエ変換
  1213. 		// p ≧ log2(a * a)を満たす最小のpを求め、
  1214. 		// aの次数2^pのフーリエ変換結果をxに格納する。
  1215. 		const int p = s2p(a.ndigits, a.ndigits);
  1216. 		i2f(x, a, p, w);
  1217.  
  1218. 		// x = x*x、および結果のxを整数に戻してaに格納する
  1219. 		mul_f2i(a, x, x, x, err_max, w);
  1220. 		add_i(a, a, a);
  1221.  
  1222. #ifdef SHOW_CSM_TIME
  1223. 		const double endTmg = omp_get_wtime();
  1224. 		csm_time.push_back(endTmg - startTmg);
  1225. #endif
  1226.  
  1227. 		// 途中経過出力
  1228. #ifdef SAVE_ITM_RESULT
  1229. 		/*pass*/
  1230. #else
  1231. 		if (i >= n - 1)
  1232. #endif
  1233. 		{
  1234. 			char filename[200];
  1235. 			sprintf_s_wrp(filename, _countof(filename), "%d.txt", i + 1);
  1236. 			print(filename, a);
  1237. 		}
  1238. 		printf("%d %f\n", i + 1, err_max);
  1239. 	}
  1240. #ifdef SHOW_CSM_TIME
  1241. 	const double endTmg_0 = omp_get_wtime();
  1242. 	printf("Total with file write = %10f sec\n", endTmg_0 - startTmg_0);
  1243. 	printf("Total only core time  = %10f sec\n", std::accumulate(csm_time.begin(), csm_time.end(), 0.0));
  1244. 	for (size_t i = 0; i < csm_time.size(); i++) {
  1245. 		printf("2^2^(%2zu) %10f sec\n", i + 1, csm_time[i]);
  1246. 	}
  1247. #endif
  1248.  
  1249. 	return 0;
  1250. }
  1251.  
Compilation error #stdin compilation error #stdout 0s 0KB
comments (?)
stdin
copy 
Standard input is empty
compilation info 
/usr/bin/ld: /home/JF5Sjr/ccrvN3nf.o: in function `main':
prog.cpp:(.text.startup+0x1c5): undefined reference to `omp_get_wtime'
/usr/bin/ld: prog.cpp:(.text.startup+0x1f1): undefined reference to `omp_get_wtime'
/usr/bin/ld: prog.cpp:(.text.startup+0x267): undefined reference to `omp_get_wtime'
/usr/bin/ld: prog.cpp:(.text.startup+0x2bc): undefined reference to `omp_get_wtime'
collect2: error: ld returned 1 exit status
stdout
copy 
Standard output is empty
https://ideone.com/ig7PuU
language:
C++14 (gcc 8.3)
created:
3 years ago
visibility:
secret

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!