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.  
  12. // [DBGSW] 中間結果表示
  13. //#define SHOW_ITM
  14.  
  15. // [DBGSW] I表示
  16. //#define SHOW_I
  17.  
  18. // [DBGSW] SIN_TABLE作成状況表示
  19. //#define SHOW_SIN_TABLE
  20.  
  21. typedef uint32_t I_TYPE;
  22.  
  23. // [CONF] 基本SINテーブルの最大要素数
  24. #define BASIC_SINTBL_P_MAX		29
  25.  
  26. // [CONF] 多倍長バッファ1要素の10進数での桁数
  27. const int FIG = 5;
  28.  
  29. // [CONF] 10^FIG
  30. const I_TYPE RADIX = 100000;
  31.  
  32. #ifdef _WIN32
  33. int fopen_s_wrp(FILE** const p_fp, const char* const fpath, const char* const mode)
  34. {
  35. 	return (int)fopen_s(p_fp, fpath, mode);
  36. }
  37. int sprintf_s_wrp(char buf[], const size_t bufsz, const char* const fmt, ...)
  38. {
  39. 	va_list ap;
  40. 	int nRet;
  41.  
  42. 	va_start(ap, fmt);
  43. 	nRet = vsprintf_s(buf, bufsz, fmt, ap);
  44. 	va_end(ap);
  45.  
  46. 	return nRet;
  47. }
  48. #else
  49. int fopen_s_wrp(FILE** const p_fp, const char* const fpath, const char* const mode)
  50. {
  51. 	*p_fp = fopen(fpath, mode);
  52. 	return (*p_fp == NULL);
  53. }
  54. int sprintf_s_wrp(char buf[], const size_t bufsz, const char* const fmt, ...)
  55. {
  56. 	va_list ap;
  57. 	int nRet;
  58.  
  59. 	va_start(ap, fmt);
  60. 	nRet = vsprintf(buf, fmt, ap);
  61. 	va_end(ap);
  62.  
  63. 	return nRet;
  64. }
  65. #define _countof(arr)	(sizeof(arr) / sizeof(arr[0]))
  66. #endif
  67.  
  68. /// 実部と虚部が交代に並ぶ形式の複素数リストを表示する。
  69. void show_re_im_seq(const std::vector<double>& s, const int p)
  70. {
  71. 	const size_t hsz = (size_t)1 << (p - 1);
  72.  
  73. 	printf("Re:");
  74. 	for (size_t i = 0; i < hsz; i++) {
  75. 		printf(" % 12f", s[2 * i]);
  76. 	}
  77. 	printf("\n");
  78. 	printf("Im:");
  79. 	for (size_t i = 0; i < hsz; i++) {
  80. 		printf(" % 12f", s[2 * i + 1]);
  81. 	}
  82. 	printf("\n");
  83. }
  84.  
  85. /// 実部と虚部が交代に並ぶ形式の複素数リストを表示する。
  86. void show_re_im_seq(const std::vector<I_TYPE>& s, const int p)
  87. {
  88. 	char fmt[128];
  89. 	sprintf_s_wrp(fmt, _countof(fmt), "%%0%dlu", FIG);
  90.  
  91. 	const size_t hsz = (size_t)1 << (p - 1);
  92. 	printf("Re:");
  93. 	for (size_t i = 0; i < hsz; i++) {
  94. 		printf(fmt, s[2 * i]);
  95. 	}
  96. 	printf("\n");
  97. 	printf("Im:");
  98. 	for (size_t i = 0; i < hsz; i++) {
  99. 		printf(fmt, s[2 * i + 1]);
  100. 	}
  101. 	printf("\n");
  102. }
  103.  
  104. namespace priroot
  105. {
  106. 	// データの桁数の桁数を表す型
  107. 	typedef int expn_t;
  108.  
  109. 	// データの桁数を表す型
  110. 	typedef int64_t pow_t;
  111.  
  112. #ifdef SHOW_SIN_TABLE
  113. 	size_t g_sin_tblent_cnt = 0;
  114. #endif
  115.  
  116. 	/// 1の原始n乗根を表す。
  117. 	/// nの剰余計算を高速化するため、nは2^pに制限する。
  118. 	class basic_nth_root
  119. 	{
  120. 		/*const*/ expn_t sv_p;
  121. 		/*const*/ pow_t t360;
  122. 		/*const*/ pow_t m360;
  123. 		/*const*/ pow_t t090;
  124. 		std::vector<double> m_sin;
  125.  
  126. 	public:
  127. 		basic_nth_root() : sv_p(0), t360(0), t090(0), m360(0) { }
  128. 		basic_nth_root(const expn_t p)
  129. 			: sv_p(p), t360((pow_t)1 << p), m360(t360 - 1), t090(t360 / 4), m_sin((size_t)(t360 / 4 + 1))
  130. 		{
  131. 			if (p >= 2) {
  132. 				assert(t360 >= 4);
  133. 				t090 = t360 / 4;
  134. 				const double k = 0.5 * M_PI;
  135. 				for (pow_t i = 1; i < t090; i++) {
  136. 					m_sin[i] = sin((k * (double)i) / (double)t090);
  137. 				}
  138. #ifdef SHOW_SIN_TABLE
  139. 				g_sin_tblent_cnt += t090;
  140. 				printf("******* %s: p=%d table made. (+%zu, %zu)\n", __FUNCTION__, p, t090, g_sin_tblent_cnt);
  141. #endif
  142. 				m_sin[0] = 0.0;
  143. 				m_sin[t090] = 1.0;
  144. 			}
  145. 		}
  146.  
  147. 		expn_t p() const { return sv_p; }
  148.  
  149. 		pow_t n() const { return t360; }
  150.  
  151. 		/// 1の原始n乗根のk乗を取得する。
  152. 		void pow(const pow_t k_, double& re, double& im) const {
  153. 			assert(k_ >= 0);
  154. 			pow_t k = k_;
  155. 			k &= m360;		// nの剰余
  156. 			if (sv_p < 2) {
  157. 				assert(k == 0 || k == 1);
  158. 				im = 0.0;
  159. 				re = (k == 0) ? 1.0 : -1.0;
  160. 				return;
  161. 			}
  162. 			assert(0 <= k && k < 4 * t090);
  163. 			if (k <= t090) {
  164. 				re = m_sin[t090 - k];
  165. 				im = m_sin[k];
  166. 			}
  167. 			else if (k <= 2 * t090) {
  168. 				re = -m_sin[k - t090];			// = t090 - (2 * t090 - k)
  169. 				im = m_sin[2 * t090 - k];
  170. 			}
  171. 			else if (k <= 3 * t090) {
  172. 				re = -m_sin[3 * t090 - k];		// = t090 - (k - 2 * t090)
  173. 				im = -m_sin[k - 2 * t090];
  174. 			}
  175. 			else {
  176. 				re = m_sin[k - 3 * t090];		// = t090 - (4 * t090 - k)
  177. 				im = -m_sin[4 * t090 - k];
  178. 			}
  179. 		}
  180.  
  181. 	};
  182.  
  183. 	/// 1の原始n乗根を表す。
  184. 	class nth_root
  185. 	{
  186. 		expn_t max_p;
  187.  
  188. 		// 1の原始2^p乗根のr乗を収めたテーブルのベクタ
  189. 		// basic_tables[p]が1の原始2^p乗根のテーブルを表す。
  190. 		std::vector<basic_nth_root> basic_tables;
  191.  
  192. 		// 1の原始2^(p+m)乗根のみを収めたベクタ
  193. 		// high_order_root_tables[m]が1の原始2^(p+m)乗根の値。(こちらはベクタ要素がテーブルではなく値を表す。)
  194. 		std::vector<std::pair<double, double> > high_order_root_values;
  195.  
  196. 		/*
  197. 			(NOTE)
  198. 			1の原始n乗根をw_nと書くことにすると、{ b_n }を0か1からなる数列、rを適当な整数として
  199. 			　　k = b_0 + 2 * b_1 + (2^2) * b_2 + ... * (2^(m-1)) * b_(m-1) + { (2^m) * r }
  200. 			であるとき、
  201. 				(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に対する剰余) }
  202. 			と書けるので、
  203. 				w_p^r (r=0..(2^p-1)) の値
  204. 				w_2^(p+q)  (q=1..m) の値
  205. 			があれば、(w^(p+m))^k を計算できる。
  206. 			前者は任意の指数rについて計算できねばならないため、テーブルで実現する場合はO(2^p)オーダーのメモリを要するが、
  207. 			後者はm個の値だけ記憶しておけば良いため、O(m)のオーダーで済む。そのかわりkが与えられる都度値の掛け算が発生する。
  208. 		*/
  209.  
  210. 	public:
  211. 		nth_root() : max_p(0) { }
  212.  
  213. 		/// 1の原始2^p乗根のテーブルを用意する。
  214. 		void alloc_table(const expn_t p) {
  215. 			if (p > max_p) {
  216. 				// (テーブル生成済のpより大きいp)
  217. 				// まずテーブル生成を検討する。
  218. 				if (p <= BASIC_SINTBL_P_MAX) {
  219. 					// (pが制限値以下)
  220. 					// [max_p+1..p]のテーブルを新規生成
  221. 					basic_tables.resize((size_t)(p + 1));
  222. 					for (expn_t i = max_p + 1; i <= p; i++) {
  223. 						basic_tables[i] = basic_nth_root(i);
  224. 					}
  225. 					max_p = p;
  226. 				}
  227. 			}
  228.  
  229. 			if (p > max_p) {
  230. 				// (任意指数に対するテーブル化が不能な大きさのp)
  231. 				// 指数k=1の値のみ記憶する。
  232. 				const expn_t m = p - max_p;
  233. 				if ((size_t)m >= high_order_root_values.size()) {
  234. 					high_order_root_values.resize((size_t)(m + 1));
  235. 					std::pair<double, double>& ent = high_order_root_values[m];		// Alias
  236. 					const pow_t big_n = (pow_t)1 << p;
  237. 					const double arg = (2 * M_PI) / big_n;
  238. 					ent.first = cos(arg);		// 実部
  239. 					ent.second = sin(arg);		// 虚部
  240. 				}
  241. 			}
  242. 		}
  243.  
  244. 		/// 1の原始2^p乗根のk乗を取得する。
  245. 		void pow(const expn_t p, const pow_t k, double& re, double& im) const {
  246. 			if (p <= max_p) {
  247. 				// (w_(2^p))^k算出
  248. 				basic_tables[p].pow(k, re, im);
  249. 			}
  250. 			else {
  251. 				const expn_t m = p - max_p;
  252. 				const pow_t r = k >> m;
  253. 				const pow_t b = k & (((pow_t)1 << m) - 1);
  254.  
  255. 				// (w_(2^max_p))^r算出
  256. 				basic_tables[max_p].pow(r, re, im);
  257.  
  258. 				// w_(2^(max_p + i))を掛け合わせる
  259. 				pow_t mask = (pow_t)1 << (m - 1);
  260. 				for (expn_t i = 1; i <= m; i++) {
  261. 					assert(mask > 0);
  262. 					if ((b & mask) != 0) {
  263. 						const std::pair<double, double>& ent = high_order_root_values[i];		// Alias
  264. 						const double tmp_re = re * ent.first - im * ent.second;
  265. 						const double tmp_im = re * ent.second + im * ent.first;
  266. 						re = tmp_re;
  267. 						im = tmp_im;
  268. 					}
  269. 					mask >>= 1;
  270. 				}
  271. 				assert(mask == 0);
  272. 			}
  273. 		}
  274.  
  275. 	};
  276.  
  277. }	// namespace priroot
  278.  
  279. namespace fftn
  280. {
  281. 	using namespace priroot;
  282.  
  283. 	/// ビット数bitwのBTR演算
  284. 	pow_t btr(const pow_t x, const expn_t bitw) {
  285. 		pow_t y = 0;
  286. 		for (expn_t i = 0; i <= (bitw - 1) / 2; i++) {
  287. 			const pow_t b = ((pow_t)1 << i) & x;
  288. 			y |= (b << ((bitw - 1) - 2 * i));
  289. 		}
  290. 		for (expn_t i = (bitw - 1) / 2 + 1; i < bitw; i++) {
  291. 			const pow_t b = ((pow_t)1 << i) & x;
  292. 			y |= (b >> (2 * i - (bitw - 1)));
  293. 		}
  294. 		return y;
  295. 	}
  296.  
  297. 	/// フーリエ変換
  298. 	/// DFTの原理式
  299. 	///   X(k) = Σ[j=0..t360-1]{ x(j) * w^(j*k) }  (k=0..t360-1)
  300. 	/// より
  301. 	///   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) }
  302. 	///        = Σ[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) }
  303. 	/// さらに、kをk+t180に置き換えると以下が成立する。
  304. 	///   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) }
  305. 	///             = Σ[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) }
  306. 	/// よって、
  307. 	///   X_even(k) = Σ[j=0..t180-1]{ x(2*j  ) * w^(2*j*k) } (k=0..t180-1)
  308. 	///   X_odd(k)  = Σ[j=0..t180-1]{ x(2*j+1) * w^(2*j*k) } (k=0..t180-1)
  309. 	/// として、X(k) (k=0..t360)は以下のバタフライ演算で表せる。
  310. 	/// 　X(k     ) = X_even(k) + w^k * X_odd(k) (k=0..t180-1)
  311. 	///   X(k+t180) = X_even(k) - w^k * X_odd(k) (k=0..t180-1)
  312. 	/// これを、大元の時間領域データx(j)の実部がbuf[2*j]、虚部がbuf[2*j + 1]に格納されている前提で行う。
  313. 	void fouri(
  314. 		std::vector<double>::iterator buf,
  315. 		const expn_t p,
  316. 		nth_root& w)
  317. 	{
  318. 		// 再帰終端判定
  319. 		if (p == 0) {
  320. 			// (要素数2^0のフーリエ変換)
  321. 			// 単一要素x(j)のフーリエ変換はx(j)なので何もしない。
  322. 			return;
  323. 		}
  324.  
  325. 		const pow_t t360 = (pow_t)1 << p;
  326. 		const pow_t t180 = (pow_t)1 << (p - 1);
  327.  
  328. 		w.alloc_table(p);
  329.  
  330. 		// X_even(k)取得
  331. 		// [0..t180-1]に対し、w^2kを回転因子とするフーリエ変換をin-placeで実行する。
  332. 		fouri(buf, p - 1, w);
  333.  
  334. 		// X_odd(k)取得
  335. 		// [t180..t360-1]に対し、w^2kを回転因子とするフーリエ変換をin-placeで実行する。
  336. 		// t180に2を掛けているのは、実部と虚部の2要素セットでt001であることによる。
  337. 		fouri(buf + 2 * t180, p - 1, w);
  338.  
  339. 		// バタフライ演算
  340. 		// X(k     ) = X_even(k) + w^k * X_odd(k)
  341. 		// X(k+t180) = X_even(k) - w^k * X_odd(k)
  342. 		for (pow_t k = 0; k < t180; k++) {
  343. 			const pow_t ofs_even = k << 1;
  344. 			const pow_t ofs_odd = (k + t180) << 1;
  345.  
  346. 			// X_even(k)
  347. 			const double X_even_k_re = buf[ofs_even];
  348. 			const double X_even_k_im = buf[ofs_even + 1];
  349.  
  350. 			// X_odd(k)
  351. 			double X_odd_k_re = buf[ofs_odd];
  352. 			double X_odd_k_im = buf[ofs_odd + 1];
  353.  
  354. 			// X_odd(k) * w^k 算出、ここでwは1の原始t360乗根
  355. 			double wk_re, wk_im;
  356. 			w.pow(p, k, wk_re, wk_im);
  357. 			const double re = X_odd_k_re * wk_re - X_odd_k_im * wk_im;
  358. 			const double im = X_odd_k_re * wk_im + X_odd_k_im * wk_re;
  359.  
  360. 			// X(k     ) = X_even(k) + w^k * X_odd(k) = X_even(k) + (re + i * im)
  361. 			buf[ofs_even] = X_even_k_re + re;
  362. 			buf[ofs_even + 1] = X_even_k_im + im;
  363.  
  364. 			// X(k+t180) = X_even(k) - w^k * X_odd(k) = X_even(k) - (re + i * im)
  365. 			buf[ofs_odd] = X_even_k_re - re;
  366. 			buf[ofs_odd + 1] = X_even_k_im - im;
  367. 		}
  368. 	}
  369.  
  370. 	/// Bit reverse orderにin-placeで並べ替える。
  371. 	void sort_by_btr(
  372. 		std::vector<double>::iterator buf,
  373. 		const expn_t p)
  374. 	{
  375. 		const pow_t t360 = (pow_t)1 << p;
  376.  
  377. 		for (pow_t i = 0; i < t360; i++) {
  378. 			const pow_t btr_i = btr(i, p);
  379.  
  380. 			// (btr_i, i)と(i, btr_i)の両方についてswapすると元に戻ってしまうので、
  381. 			// 片方のみswap対象とする。
  382. 			if (btr_i < i) {
  383. 				std::swap(buf[2 * btr_i], buf[2 * i]);				// 実部
  384. 				std::swap(buf[2 * btr_i + 1], buf[2 * i + 1]);		// 虚部
  385. 			}
  386. 		}
  387. 	}
  388.  
  389. 	/// 逆フーリエ変換
  390. 	/// フーリエ変換とは逆の以下のバタフライ演算を行う。
  391. 	/// 　x(k     ) = (x_even(k) + w^(-k) * x_odd(k)) / t360 (k=0..t180-1)
  392. 	///   x(k+t180) = (x_even(k) - w^(-k) * x_odd(k)) / t360 (k=0..t180-1)
  393. 	/// これを、大元の周波数領域データX(k)の実部がbuf[2*k]、虚部がbuf[2*k+1]に格納されている前提で行う。
  394. 	void fouriinv(
  395. 		std::vector<double>::iterator buf,
  396. 		const expn_t p,
  397. 		nth_root& w)
  398. 	{
  399. 		// 再帰終端判定
  400. 		if (p == 0) {
  401. 			// (要素数2^0のフーリエ変換)
  402. 			// 単一要素x(j)のフーリエ変換はx(j)なので何もしない。
  403. 			return;
  404. 		}
  405.  
  406. 		const pow_t t360 = (pow_t)1 << p;
  407. 		const pow_t t180 = (pow_t)1 << (p - 1);
  408.  
  409. 		w.alloc_table(p);
  410.  
  411. 		// x_even(k)取得
  412. 		// [0..t180-1]に対し、w^2kを回転因子とするフーリエ変換をin-placeで実行する。
  413. 		fouriinv(buf, p - 1, w);
  414.  
  415. 		// x_odd(k)取得
  416. 		// [t180..t360-1]に対し、w^2kを回転因子とするフーリエ変換をin-placeで実行する。
  417. 		// t180に2を掛けているのは、実部と虚部の2要素セットでt001であることによる。
  418. 		fouriinv(buf + 2 * t180, p - 1, w);
  419.  
  420. 		// バタフライ演算
  421. 		// x(k     ) = x_even(k) + w^k * x_odd(k)
  422. 		// x(k+t180) = x_even(k) - w^k * x_odd(k)
  423. 		for (pow_t k = 0; k < t180; k++) {
  424. 			const pow_t ofs_even = k << 1;
  425. 			const pow_t ofs_odd = (k + t180) << 1;
  426.  
  427. 			// x_even(k)
  428. 			const double x_even_k_re = buf[ofs_even];
  429. 			const double x_even_k_im = buf[ofs_even + 1];
  430.  
  431. 			// x_odd(k)
  432. 			double x_odd_k_re = buf[ofs_odd];
  433. 			double x_odd_k_im = buf[ofs_odd + 1];
  434.  
  435. 			// x_odd(k) * w^(-k) 算出、ここでwは1の原始t360乗根
  436. 			double wk_re, wk_im;
  437. 			w.pow(p, k, wk_re, wk_im);
  438. 			const double re = x_odd_k_re * wk_re + x_odd_k_im * wk_im;
  439. 			const double im = x_odd_k_im * wk_re - x_odd_k_re * wk_im;
  440.  
  441. 			// X(k     ) = x_even(k) + w^k * x_odd(k) = x_even(k) + (re + i * im)
  442. 			buf[ofs_even] = (x_even_k_re + re) / 2.0;
  443. 			buf[ofs_even + 1] = (x_even_k_im + im) / 2.0;
  444.  
  445. 			// X(k+t180) = x_even(k) - w^k * x_odd(k) = x_even(k) - (re + i * im)
  446. 			buf[ofs_odd] = (x_even_k_re - re) / 2.0;
  447. 			buf[ofs_odd + 1] = (x_even_k_im - im) / 2.0;
  448. 		}
  449. 	}
  450.  
  451. }	// namespace fftn
  452.  
  453. namespace fft4n
  454. {
  455. 	using namespace priroot;
  456.  
  457. 	/// 4n乗根によるフーリエ変換
  458. 	/// x(j) (j=0..(2*t360)-1)から
  459. 	///   x'(j) = (x(j) + i * x(j+t360)) * u^(j/4) (j=0..t360-1)
  460. 	/// を構成し、x'(j)をDFTする。ここで、uは1の原始(4 * mag * t360)乗根 (mag ≧ 1)。
  461. 	void fouri4n(
  462. 		std::vector<double>::iterator buf,
  463. 		const expn_t p,
  464. 		nth_root& w)
  465. 	{
  466. 		const pow_t t360 = (pow_t)1 << p;
  467.  
  468. 		w.alloc_table(p + 2);
  469.  
  470. 		// { x'(j) }構成
  471. 		// { x(j) + i * x(j+t360) } にu^(j/4)を掛ける。
  472. 		const expn_t p2 = p + 2;
  473. 		for (pow_t j = 0; j < t360; j++) {
  474. 			const double a = buf[2 * j];
  475. 			const double b = buf[2 * j + 1];
  476.  
  477. 			// u = w^4としてu^(j/4) = w^(j)
  478. 			double wj_re, wj_im;
  479. 			w.pow(p2, j, wj_re, wj_im);
  480. 			buf[2 * j] = a * wj_re - b * wj_im;
  481. 			buf[2 * j + 1] = a * wj_im + b * wj_re;
  482. 		}
  483.  
  484. 		// { x'(j) }をbit reverse orderに変換
  485. 		fftn::sort_by_btr(buf, p);
  486.  
  487. 		// フーリエ変換
  488. 		fftn::fouri(buf, p, w);
  489. 	}
  490.  
  491. 	/// 4n乗根による逆フーリエ変換
  492. 	void fouri4ninv(
  493. 		std::vector<double>::iterator buf,
  494. 		const expn_t p,
  495. 		nth_root& w)
  496. 	{
  497. 		const pow_t t360 = (pow_t)1 << p;
  498.  
  499. 		w.alloc_table(p + 2);
  500.  
  501. 		// { X'(k) }をbit reverse orderに変換
  502. 		fftn::sort_by_btr(buf, p);
  503.  
  504. 		// 逆フーリエ変換
  505. 		fftn::fouriinv(buf, p, w);
  506.  
  507. 		// { x(j) }復元
  508. 		const expn_t p2 = p + 2;
  509. 		for (pow_t j = 0; j < t360; j++) {
  510. 			const double a = buf[2 * j];
  511. 			const double b = buf[2 * j + 1];
  512.  
  513. 			// u = w^4としてu^(-j/4) = w^(-j)を掛ける
  514. 			double wj_re, wj_im;
  515. 			w.pow(p2, j, wj_re, wj_im);
  516. 			buf[2 * j] = a * wj_re + b * wj_im;
  517. 			buf[2 * j + 1] = b * wj_re - a * wj_im;
  518. 		}
  519. 	}
  520.  
  521. }	// namespace fft4n
  522.  
  523. namespace oper
  524. {
  525. 	using namespace fft4n;
  526.  
  527. 	struct I
  528. 	{
  529. 		std::vector<I_TYPE> data;	// 基数RADIXの多倍長データ
  530. 		size_t ndigits;				// dataに入っている有効な桁数
  531.  
  532. 		I(const size_t cap)
  533. 			: data(cap), ndigits((size_t)0)
  534. 		{
  535. 			/*pass*/
  536. 		}
  537. 	};
  538.  
  539. 	struct F
  540. 	{
  541. 		std::vector<double> data;	// 基数[-(RADIX/2), RADIX/2)の多倍長データ
  542. 		int p;						// 2^p == (dataの要素数) となるp
  543.  
  544. 		F(const int p_)
  545. 			: data((size_t)1 << p_), p(0)
  546. 		{
  547. 			/*pass*/
  548. 		}
  549. 	};
  550.  
  551. 	/// 整数のフーリエ変換
  552. 	void i2f(F& a, const I& b, const int p, priroot::nth_root& w)
  553. 	{
  554. 		size_t i;
  555.  
  556. 		// 次数2^(p-1) = t360の4n乗根フーリエ変換を行う。
  557. 		// 結果は実部と虚部を合わせて2^p = (2 * t360)桁分となる。
  558. 		const pow_t t360 = (pow_t)1 << (p - 1);
  559.  
  560. 		assert(b.ndigits <= (size_t)t360);
  561. 		assert((size_t)(2 * t360) <= a.data.size());
  562.  
  563. #ifdef SHOW_I
  564. 		printf("%s: Before pre tr.:\n", __FUNCTION__);
  565. 		show_re_im_seq(b.data, p);
  566. #endif
  567.  
  568. 		// フーリエ変換の前処理
  569. 		// (bが表す整数) = Σ[i=0..a.ndigits-1]{ b.data[i] * RADIX^i } (0≦b.data[i]＜RADIX)
  570. 		// を、
  571. 		// (aが表す整数) = Σ[i=0..a.ndigits-1]{ a.data[i] * RADIX^i } (-RADIX / 2≦a.data[i]＜RADIX / 2)
  572. 		// に変換する。(a, bが同じ整数を表すように桁上げも行う。)
  573. 		a.p = p;
  574. 		int carry = 0;
  575. 		for (i = 0; i < b.ndigits - 1; i++) {
  576. 			// 実部虚部隣接化
  577. 			// b.data[t360..]を虚部に配置する。
  578. 			const size_t dst_i = (i < (size_t)t360) ? 2 * i : 2 * (i - (size_t)t360) + 1;
  579.  
  580. 			const int32_t t = b.data[i] + carry;	// 桁[i]の値
  581. 			if (t < (int32_t)RADIX / 2) {
  582. 				// (桁[i]の値∈[0, RADIX / 2) )
  583. 				// 何もしない
  584. 				a.data[dst_i] = t;
  585. 				carry = 0;
  586. 			}
  587. 			else {
  588. 				// (桁[i]の値∈[RADIX / 2, RADIX) )
  589. 				// 補正実施
  590. 				// まず桁[i] -= RADIX により桁[i]の値∈[0, RADIX / 2)とする。ただし
  591. 				// それだけでは(aが表す整数)が RADIX * RADIX^i = RADIX^(i+1)だけ減ってしまうので、
  592. 				// 減った分をcarryで一つ上の桁に回して補填する。
  593. 				a.data[dst_i] = (double)(t - (int32_t)RADIX);
  594. 				carry = 1;
  595. 			}
  596. 		}
  597. 		{
  598. 			// 実部虚部隣接化
  599. 			const size_t dst_i = (i < (size_t)t360) ? 2 * i : 2 * (i - (size_t)t360) + 1;
  600.  
  601. 			// 最上位桁
  602. 			a.data[dst_i] = (double)(b.data[i] + (I_TYPE)carry);
  603. 			i++;
  604. 		}
  605. 		for (; i < (size_t)(2 * t360); i++) {
  606. 			// 実部虚部隣接化
  607. 			const size_t dst_i = (i < (size_t)t360) ? 2 * i : 2 * (i - (size_t)t360) + 1;
  608.  
  609. 			a.data[dst_i] = 0.0;
  610. 		}
  611.  
  612. #if defined(SHOW_ITM) || defined(SHOW_I)
  613. 		printf("Before fouri4n():\n");
  614. 		show_re_im_seq(a.data, a.p);
  615. #endif
  616.  
  617. 		// フーリエ変換実施
  618. 		fft4n::fouri4n(a.data.begin(), p - 1, w);
  619.  
  620. #ifdef SHOW_ITM
  621. 		printf("After fouri4n():\n");
  622. 		show_re_im_seq(a.data, a.p);
  623. #endif
  624. 	}
  625.  
  626. 	/// 整数のフーリエ変換結果の積(結果は整数)
  627. 	void mul_f2i(I& a, const F& b, const F& c, F& d, double& err_max, priroot::nth_root& w) {
  628. 		int64_t carry;
  629. 		size_t i, t;
  630.  
  631. 		const int p = b.p;
  632.  
  633. 		// 次数2^(p-1) = t360の4n乗根フーリエ変換結果の積を逆フーリエ変換する。
  634. 		// 結果は実部と虚部を合わせて2^p = (2 * t360)桁となる。
  635. 		const pow_t t360 = (size_t)1 << (p - 1);
  636.  
  637. 		assert(p == c.p);
  638. 		assert((size_t)(2 * t360) <= d.data.size());
  639.  
  640. 		// 積
  641. 		for (i = 0; i < (size_t)t360; i++) {
  642. 			const double b_i_re = b.data[2 * i];
  643. 			const double b_i_im = b.data[2 * i + 1];
  644. 			const double c_i_re = c.data[2 * i];
  645. 			const double c_i_im = c.data[2 * i + 1];
  646.  
  647. 			// d[i] = b[i] * c[i]
  648. 			const double d_i_re = b_i_re * c_i_re - b_i_im * c_i_im;
  649. 			const double d_i_im = b_i_re * c_i_im + b_i_im * c_i_re;
  650. 			d.data[2 * i] = d_i_re;
  651. 			d.data[2 * i + 1] = d_i_im;
  652. 		}
  653.  
  654. #ifdef SHOW_ITM
  655. 		printf("Before fouri4ninv():\n");
  656. 		show_re_im_seq(d.data, d.p);
  657. #endif
  658.  
  659. 		// 逆フーリエ変換
  660. 		fft4n::fouri4ninv(d.data.begin(), p - 1, w);
  661.  
  662. #ifdef SHOW_ITM
  663. 		printf("After fouri4ninv():\n");
  664. 		show_re_im_seq(d.data, d.p);
  665. #endif
  666.  
  667. 		const size_t s = std::min(a.data.size(), (size_t)(2 * t360));
  668.  
  669. 		// Carry and fix処理
  670. 		// ただし対象とする桁の並びd.data[]は、
  671. 		// 桁の値が[0, RADIX)ではなく、[-RADIX / 2, RADIX / 2)である。
  672. 		carry = 0;
  673. 		t = 1;
  674. 		//const double k = ldexp(1.0, -p + 1);		// k = 1.0 * 2^(-p + 1)
  675. 		double k = 1.0;		// fouri4ninv()結果は2^pで割る必要無し。
  676. 		for (i = 0; i < s; i++) {
  677. 			// 実部虚部隣接化解除
  678. 			// b.data[t360..]を虚部に配置する。
  679. 			const size_t src_i = (i < (size_t)t360) ? 2 * i : 2 * (i - (size_t)t360) + 1;
  680. 			//printf("t360=%lld, s=%zu, i=%zu, src_i=%zu\n", t360, s, i, src_i);
  681.  
  682. 			const double x = d.data[src_i] * k;
  683. 			int64_t e = (int64_t)floor(x + 0.5);
  684. 			err_max = std::max(err_max, fabs(x - e));
  685. 			e += carry;
  686. 			carry = e / RADIX;
  687. 			int64_t f = e - carry * RADIX;
  688. 			if (f < 0) {
  689. 				f += RADIX;
  690. 				carry--;
  691. 			}
  692. 			a.data[i] = (I_TYPE)f;
  693. 			if (f != 0) {
  694. 				t = i + 1;
  695. 			}
  696. 		}
  697. 		a.ndigits = t;
  698. 		if (carry != 0) {
  699. 			if (a.ndigits < a.data.size()) {
  700. 				assert(0 <= carry && carry < RADIX);
  701. 				a.data[a.ndigits++] = (I_TYPE)carry;
  702. 			}
  703. 			else {
  704. 				// ここに来たら桁あふれ
  705. 				abort();
  706. 			}
  707. 		}
  708. 		for (; i < s; i++) {
  709. 			// 実部虚部隣接化解除
  710. 			const size_t src_i = (i < (size_t)t360) ? 2 * i : 2 * i + 1;
  711.  
  712. 			const double x = d.data[src_i] * k;
  713. 			const int64_t e = (int64_t)floor(x + 0.5);
  714. 			err_max = std::max(err_max, fabs(x - e));
  715. 			assert(e == 0);
  716. 		}
  717. 	}
  718.  
  719. 	/// 整数の和
  720. 	void add_i(I& a, const I& b_, const I& c_) {
  721. 		I_TYPE carry, d;
  722. 		size_t i;
  723.  
  724. 		const I& b = (b_.ndigits > c_.ndigits) ? c_ : b_;	// Alias
  725. 		const I& c = (b_.ndigits > c_.ndigits) ? b_ : c_;	// Alias
  726. 		assert(b.ndigits <= c.ndigits);
  727. 		assert(c.ndigits <= a.data.size());
  728.  
  729. 		carry = 0;
  730. 		for (i = 0; i < b.ndigits; i++) {
  731. 			d = b.data[i] + c.data[i] + carry;
  732. 			carry = 0;
  733. 			if (d >= RADIX) {
  734. 				d -= RADIX;
  735. 				carry++;
  736. 			}
  737. 			a.data[i] = d;
  738. 		}
  739. 		for (; i < c.ndigits; i++) {
  740. 			d = c.data[i] + carry;
  741. 			carry = 0;
  742. 			if (d >= RADIX) {
  743. 				d -= RADIX;
  744. 				carry++;
  745. 			}
  746. 			a.data[i] = d;
  747. 		}
  748. 		a.ndigits = c.ndigits;
  749. 		if (a.ndigits < a.data.size()) {
  750. 			if (carry != 0) {
  751. 				a.data[a.ndigits++] = carry;
  752. 			}
  753. 		}
  754. 		else {
  755. 			if (carry != 0) {
  756. 				// ここに来たら桁あふれ
  757. 				abort();
  758. 			}
  759. 		}
  760. 	}
  761.  
  762. 	/// 整数のファイル出力
  763. 	void print(const char* const filename, const I& a)
  764. 	{
  765. 		char format[10];
  766. 		int nRet;
  767.  
  768. 		assert(a.data.size() > (size_t)0);
  769. 		assert(a.ndigits <= a.data.size());
  770.  
  771. 		FILE* fp;
  772. 		nRet = fopen_s_wrp(&fp, filename, "w");
  773. 		if (nRet != 0) {
  774. 			fflush(stdout);
  775. 			fprintf(stderr, "ERROR: fopen_s_wrp() failed. (nRet=%d)\n", nRet);
  776. 			abort();
  777. 		}
  778.  
  779. 		sprintf_s_wrp(format, _countof(format), "%%0%du", FIG);
  780.  
  781. 		size_t s = a.ndigits;
  782. 		fprintf(fp, "%lu", a.data[--s]);
  783. 		while (s != 0) {
  784. 			fprintf(fp, format, a.data[--s]);
  785. 		}
  786. 		fclose(fp);
  787. 	}
  788.  
  789. }	// namespace oper
  790.  
  791. namespace
  792. {
  793. 	// N ≧ log2(s1 * s2)を満たす最小のNを求める。
  794. 	int s2p(const size_t s1, const size_t s2) {
  795. 		int i;
  796. 		for (i = 1; i < 64; i++) {
  797. 			if (((size_t)1 << i) + 1 >= s1 + s2) {
  798. 				// (2^i + 1 ≧ (s1の有効桁数) + (s2の有効桁数))
  799. 				break;
  800. 			}
  801. 		}
  802. 		return i;
  803. 	}
  804.  
  805. }	// namespace
  806.  
  807. int main()
  808. {
  809. 	using namespace oper;
  810.  
  811. 	const int n = 34;
  812.  
  813. 	nth_root w;
  814.  
  815. 	const int max_p = s2p(((size_t)((1LLU << n) * log10(2)) / FIG + 1), 1);
  816. 	const size_t cap = (size_t)1 << max_p;
  817.  
  818. 	I a(cap);
  819. 	F x(max_p);
  820.  
  821. 	a.data[0] = 1;
  822. 	a.ndigits = 1;
  823.  
  824. 	double err_max = 0.0;
  825. 	for (int i = 0; i < n; i++) {
  826.  
  827. 		// 多倍長整数xのフーリエ変換
  828. 		// p ≧ log2(a * a)を満たす最小のpを求め、
  829. 		// aの次数2^pのフーリエ変換結果をxに格納する。
  830. 		const int p = s2p(a.ndigits, a.ndigits);
  831. 		i2f(x, a, p, w);
  832.  
  833. 		// x = x*x、および結果のxを整数に戻してaに格納する
  834. 		mul_f2i(a, x, x, x, err_max, w);
  835. 		add_i(a, a, a);
  836.  
  837. 		// 途中経過出力
  838. 		{
  839. 			char filename[200];
  840. 			sprintf_s_wrp(filename, _countof(filename), "%d.txt", i + 1);
  841. 			print(filename, a);
  842. 		}
  843. 		printf("%d %f\n", i + 1, err_max);
  844. 	}
  845.  
  846. 	return 0;
  847. }
  848.  
Runtime error #stdin #stdout 0.78s 2095224KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/xHSzFD
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!