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