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