fork 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. // [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), 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)を計算できる。
  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. 				// 指数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 fft4n
  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*btr_j]、虚部がbuf[2*btr_j + 1]に格納されている前提で行う。
  313. 	/// ここで、t360=2^pとして、btr_jはjのビット幅pのbit reverse order。
  314. 	void fouri(
  315. 		std::vector<double>::iterator buf,
  316. 		const expn_t p,
  317. 		nth_root& w)
  318. 	{
  319. 		// 再帰終端判定
  320. 		if (p == 0) {
  321. 			// (要素数2^0のフーリエ変換)
  322. 			// 単一要素x(j)のフーリエ変換はx(j)なので何もしない。
  323. 			return;
  324. 		}
  325.  
  326. 		const pow_t t360 = (pow_t)1 << p;
  327. 		const pow_t t180 = (pow_t)1 << (p - 1);
  328.  
  329. 		w.alloc_table(p);
  330.  
  331. 		// X_even(k)取得
  332. 		// [0..t180-1]に対し、w^2kを回転因子とするフーリエ変換をin-placeで実行する。
  333. 		fouri(buf, p - 1, w);
  334.  
  335. 		// X_odd(k)取得
  336. 		// [t180..t360-1]に対し、w^2kを回転因子とするフーリエ変換をin-placeで実行する。
  337. 		// t180に2を掛けているのは、実部と虚部の2要素セットでt001であることによる。
  338. 		fouri(buf + 2 * t180, p - 1, w);
  339.  
  340. 		// バタフライ演算
  341. 		// X(k     ) = X_even(k) + w^k * X_odd(k)
  342. 		// X(k+t180) = X_even(k) - w^k * X_odd(k)
  343. 		for (pow_t k = 0; k < t180; k++) {
  344. 			const pow_t ofs_even = k << 1;
  345. 			const pow_t ofs_odd = (k + t180) << 1;
  346.  
  347. 			// X_even(k)
  348. 			const double X_even_k_re = buf[ofs_even];
  349. 			const double X_even_k_im = buf[ofs_even + 1];
  350.  
  351. 			// X_odd(k)
  352. 			double X_odd_k_re = buf[ofs_odd];
  353. 			double X_odd_k_im = buf[ofs_odd + 1];
  354.  
  355. 			// X_odd(k) * w^k 算出、ここでwは1の原始t360乗根
  356. 			double wk_re, wk_im;
  357. 			w.pow(p, k, wk_re, wk_im);
  358. 			const double re = X_odd_k_re * wk_re - X_odd_k_im * wk_im;
  359. 			const double im = X_odd_k_re * wk_im + X_odd_k_im * wk_re;
  360.  
  361. 			// X(k     ) = X_even(k) + w^k * X_odd(k) = X_even(k) + (re + i * im)
  362. 			buf[ofs_even] = X_even_k_re + re;
  363. 			buf[ofs_even + 1] = X_even_k_im + im;
  364.  
  365. 			// X(k+t180) = X_even(k) - w^k * X_odd(k) = X_even(k) - (re + i * im)
  366. 			buf[ofs_odd] = X_even_k_re - re;
  367. 			buf[ofs_odd + 1] = X_even_k_im - im;
  368. 		}
  369. 	}
  370.  
  371. 	/// Bit reverse orderにin-placeで並べ替える。
  372. 	void sort_by_btr(
  373. 		std::vector<double>::iterator buf,
  374. 		const expn_t p)
  375. 	{
  376. 		const pow_t t360 = (pow_t)1 << p;
  377.  
  378. 		for (pow_t i = 0; i < t360; i++) {
  379. 			const pow_t btr_i = btr(i, p);
  380.  
  381. 			// (btr_i, i)と(i, btr_i)の両方についてswapすると元に戻ってしまうので、
  382. 			// 片方のみswap対象とする。
  383. 			if (btr_i < i) {
  384. 				std::swap(buf[2 * btr_i], buf[2 * i]);				// 実部
  385. 				std::swap(buf[2 * btr_i + 1], buf[2 * i + 1]);		// 虚部
  386. 			}
  387. 		}
  388. 	}
  389.  
  390. 	/// 逆フーリエ変換
  391. 	/// フーリエ変換とは逆の以下のバタフライ演算を行う。
  392. 	/// 　x(k     ) = (x_even(k) + w^(-k) * x_odd(k)) / t360 (k=0..t180-1)
  393. 	///   x(k+t180) = (x_even(k) - w^(-k) * x_odd(k)) / t360 (k=0..t180-1)
  394. 	/// これを、大元の周波数領域データX(k)の実部がbuf[2*k]、虚部がbuf[2*k+1]に格納されている前提で行う。
  395. 	void fouriinv(
  396. 		std::vector<double>::iterator buf,
  397. 		const expn_t p,
  398. 		nth_root& w)
  399. 	{
  400. 		// 再帰終端判定
  401. 		if (p == 0) {
  402. 			// (要素数2^0のフーリエ変換)
  403. 			// 単一要素x(j)のフーリエ変換はx(j)なので何もしない。
  404. 			return;
  405. 		}
  406.  
  407. 		const pow_t t360 = (pow_t)1 << p;
  408. 		const pow_t t180 = (pow_t)1 << (p - 1);
  409.  
  410. 		w.alloc_table(p);
  411.  
  412. 		// x_even(k)取得
  413. 		// [0..t180-1]に対し、w^2kを回転因子とするフーリエ変換をin-placeで実行する。
  414. 		fouriinv(buf, p - 1, w);
  415.  
  416. 		// x_odd(k)取得
  417. 		// [t180..t360-1]に対し、w^2kを回転因子とするフーリエ変換をin-placeで実行する。
  418. 		// t180に2を掛けているのは、実部と虚部の2要素セットでt001であることによる。
  419. 		fouriinv(buf + 2 * t180, p - 1, w);
  420.  
  421. 		// バタフライ演算
  422. 		// x(k     ) = x_even(k) + w^k * x_odd(k)
  423. 		// x(k+t180) = x_even(k) - w^k * x_odd(k)
  424. 		for (pow_t k = 0; k < t180; k++) {
  425. 			const pow_t ofs_even = k << 1;
  426. 			const pow_t ofs_odd = (k + t180) << 1;
  427.  
  428. 			// x_even(k)
  429. 			const double x_even_k_re = buf[ofs_even];
  430. 			const double x_even_k_im = buf[ofs_even + 1];
  431.  
  432. 			// x_odd(k)
  433. 			double x_odd_k_re = buf[ofs_odd];
  434. 			double x_odd_k_im = buf[ofs_odd + 1];
  435.  
  436. 			// x_odd(k) * w^(-k) 算出、ここでwは1の原始t360乗根
  437. 			double wk_re, wk_im;
  438. 			w.pow(p, k, wk_re, wk_im);
  439. 			const double re = x_odd_k_re * wk_re + x_odd_k_im * wk_im;
  440. 			const double im = x_odd_k_im * wk_re - x_odd_k_re * wk_im;
  441.  
  442. 			// X(k     ) = x_even(k) + w^k * x_odd(k) = x_even(k) + (re + i * im)
  443. 			buf[ofs_even] = (x_even_k_re + re) / 2.0;
  444. 			buf[ofs_even + 1] = (x_even_k_im + im) / 2.0;
  445.  
  446. 			// X(k+t180) = x_even(k) - w^k * x_odd(k) = x_even(k) - (re + i * im)
  447. 			buf[ofs_odd] = (x_even_k_re - re) / 2.0;
  448. 			buf[ofs_odd + 1] = (x_even_k_im - im) / 2.0;
  449. 		}
  450. 	}
  451.  
  452. 	/// 4n乗根によるフーリエ変換
  453. 	/// x(j) (j=0..(2*t360)-1)から
  454. 	///   x'(j) = (x(j) + i * x(j+t360)) * u^(j/4) (j=0..t360-1)
  455. 	/// を構成し、x'(j)をDFTする。ここで、uは1の原始(4 * mag * t360)乗根 (mag ≧ 1)。
  456. 	void fouri4n(
  457. 		std::vector<double>::iterator buf,
  458. 		const expn_t p,
  459. 		nth_root& w)
  460. 	{
  461. 		const pow_t t360 = (pow_t)1 << p;
  462.  
  463. 		w.alloc_table(p + 2);
  464.  
  465. 		// { x'(j) }構成
  466. 		// { x(j) + i * x(j+t360) } にu^(j/4)を掛ける。
  467. 		const expn_t p2 = p + 2;
  468. 		for (pow_t j = 0; j < t360; j++) {
  469. 			const double a = buf[2 * j];
  470. 			const double b = buf[2 * j + 1];
  471.  
  472. 			// u = w^4としてu^(j/4) = w^(j)
  473. 			double wj_re, wj_im;
  474. 			w.pow(p2, j, wj_re, wj_im);
  475. 			buf[2 * j] = a * wj_re - b * wj_im;
  476. 			buf[2 * j + 1] = a * wj_im + b * wj_re;
  477. 		}
  478.  
  479. 		// { x'(j) }をbit reverse orderに変換
  480. 		sort_by_btr(buf, p);
  481.  
  482. 		// フーリエ変換
  483. 		fouri(buf, p, w);
  484. 	}
  485.  
  486. 	/// 4n乗根による逆フーリエ変換
  487. 	void fouri4ninv(
  488. 		std::vector<double>::iterator buf,
  489. 		const expn_t p,
  490. 		nth_root& w)
  491. 	{
  492. 		const pow_t t360 = (pow_t)1 << p;
  493.  
  494. 		w.alloc_table(p + 2);
  495.  
  496. 		// { X'(k) }をbit reverse orderに変換
  497. 		sort_by_btr(buf, p);
  498.  
  499. 		// 逆フーリエ変換
  500. 		fouriinv(buf, p, w);
  501.  
  502. 		// { x(j) }復元
  503. 		const expn_t p2 = p + 2;
  504. 		for (pow_t j = 0; j < t360; j++) {
  505. 			const double a = buf[2 * j];
  506. 			const double b = buf[2 * j + 1];
  507.  
  508. 			// u = w^4としてu^(-j/4) = w^(-j)を掛ける
  509. 			double wj_re, wj_im;
  510. 			w.pow(p2, j, wj_re, wj_im);
  511. 			buf[2 * j] = a * wj_re + b * wj_im;
  512. 			buf[2 * j + 1] = b * wj_re - a * wj_im;
  513. 		}
  514. 	}
  515.  
  516. }	// namespace fft4n
  517.  
  518. namespace oper
  519. {
  520. 	using namespace fft4n;
  521.  
  522. 	struct I
  523. 	{
  524. 		std::vector<I_TYPE> data;	// 基数RADIXの多倍長データ
  525. 		size_t ndigits;				// dataに入っている有効な桁数
  526.  
  527. 		I(const size_t cap)
  528. 			: data(cap), ndigits((size_t)0)
  529. 		{
  530. 			/*pass*/
  531. 		}
  532. 	};
  533.  
  534. 	struct F
  535. 	{
  536. 		std::vector<double> data;	// 基数[-(RADIX/2), RADIX/2)の多倍長データ
  537. 		int p;						// 2^p == (dataの要素数) となるp
  538.  
  539. 		F(const int p_)
  540. 			: data((size_t)1 << p_), p(0)
  541. 		{
  542. 			/*pass*/
  543. 		}
  544. 	};
  545.  
  546. 	/// 整数のフーリエ変換
  547. 	void i2f(F& a, const I& b, const int p, priroot::nth_root& w)
  548. 	{
  549. 		size_t i;
  550.  
  551. 		// 次数2^(p-1) = t360の4n乗根フーリエ変換を行う。
  552. 		// 結果は実部と虚部を合わせて2^p = (2 * t360)桁分となる。
  553. 		const pow_t t360 = (pow_t)1 << (p - 1);
  554.  
  555. 		assert(b.ndigits <= (size_t)t360);
  556. 		assert((size_t)(2 * t360) <= a.data.size());
  557.  
  558. #ifdef SHOW_I
  559. 		printf("%s: Before pre tr.:\n", __FUNCTION__);
  560. 		show_re_im_seq(b.data, p);
  561. #endif
  562.  
  563. 		// フーリエ変換の前処理
  564. 		// (bが表す整数) = Σ[i=0..a.ndigits-1]{ b.data[i] * RADIX^i } (0≦b.data[i]＜RADIX)
  565. 		// を、
  566. 		// (aが表す整数) = Σ[i=0..a.ndigits-1]{ a.data[i] * RADIX^i } (-RADIX / 2≦a.data[i]＜RADIX / 2)
  567. 		// に変換する。(a, bが同じ整数を表すように桁上げも行う。)
  568. 		a.p = p;
  569. 		int carry = 0;
  570. 		for (i = 0; i < b.ndigits - 1; i++) {
  571. 			// 実部虚部隣接化
  572. 			// b.data[t360..]を虚部に配置する。
  573. 			const size_t dst_i = (i < (size_t)t360) ? 2 * i : 2 * (i - (size_t)t360) + 1;
  574.  
  575. 			const int32_t t = b.data[i] + carry;	// 桁[i]の値
  576. 			if (t < (int32_t)RADIX / 2) {
  577. 				// (桁[i]の値∈[0, RADIX / 2) )
  578. 				// 何もしない
  579. 				a.data[dst_i] = t;
  580. 				carry = 0;
  581. 			}
  582. 			else {
  583. 				// (桁[i]の値∈[RADIX / 2, RADIX) )
  584. 				// 補正実施
  585. 				// まず桁[i] -= RADIX により桁[i]の値∈[0, RADIX / 2)とする。ただし
  586. 				// それだけでは(aが表す整数)が RADIX * RADIX^i = RADIX^(i+1)だけ減ってしまうので、
  587. 				// 減った分をcarryで一つ上の桁に回して補填する。
  588. 				a.data[dst_i] = (double)(t - (int32_t)RADIX);
  589. 				carry = 1;
  590. 			}
  591. 		}
  592. 		{
  593. 			// 実部虚部隣接化
  594. 			const size_t dst_i = (i < (size_t)t360) ? 2 * i : 2 * (i - (size_t)t360) + 1;
  595.  
  596. 			// 最上位桁
  597. 			a.data[dst_i] = (double)(b.data[i] + (I_TYPE)carry);
  598. 			i++;
  599. 		}
  600. 		for (; i < (size_t)(2 * t360); i++) {
  601. 			// 実部虚部隣接化
  602. 			const size_t dst_i = (i < (size_t)t360) ? 2 * i : 2 * (i - (size_t)t360) + 1;
  603.  
  604. 			a.data[dst_i] = 0.0;
  605. 		}
  606.  
  607. #if defined(SHOW_ITM) || defined(SHOW_I)
  608. 		printf("Before fouri4n():\n");
  609. 		show_re_im_seq(a.data, a.p);
  610. #endif
  611.  
  612. 		// フーリエ変換実施
  613. 		fft4n::fouri4n(a.data.begin(), p - 1, w);
  614.  
  615. #ifdef SHOW_ITM
  616. 		printf("After fouri4n():\n");
  617. 		show_re_im_seq(a.data, a.p);
  618. #endif
  619. 	}
  620.  
  621. 	/// 整数のフーリエ変換結果の積(結果は整数)
  622. 	void mul_f2i(I& a, const F& b, const F& c, F& d, double& err_max, priroot::nth_root& w) {
  623. 		int64_t carry;
  624. 		size_t i, t;
  625.  
  626. 		const int p = b.p;
  627.  
  628. 		// 次数2^(p-1) = t360の4n乗根フーリエ変換結果の積を逆フーリエ変換する。
  629. 		// 結果は実部と虚部を合わせて2^p = (2 * t360)桁となる。
  630. 		const pow_t t360 = (size_t)1 << (p - 1);
  631.  
  632. 		assert(p == c.p);
  633. 		assert((size_t)(2 * t360) <= d.data.size());
  634.  
  635. 		// 積
  636. 		for (i = 0; i < (size_t)t360; i++) {
  637. 			const double b_i_re = b.data[2 * i];
  638. 			const double b_i_im = b.data[2 * i + 1];
  639. 			const double c_i_re = c.data[2 * i];
  640. 			const double c_i_im = c.data[2 * i + 1];
  641.  
  642. 			// d[i] = b[i] * c[i]
  643. 			const double d_i_re = b_i_re * c_i_re - b_i_im * c_i_im;
  644. 			const double d_i_im = b_i_re * c_i_im + b_i_im * c_i_re;
  645. 			d.data[2 * i] = d_i_re;
  646. 			d.data[2 * i + 1] = d_i_im;
  647. 		}
  648.  
  649. #ifdef SHOW_ITM
  650. 		printf("Before fouri4ninv():\n");
  651. 		show_re_im_seq(d.data, d.p);
  652. #endif
  653.  
  654. 		// 逆フーリエ変換
  655. 		fft4n::fouri4ninv(d.data.begin(), p - 1, w);
  656.  
  657. #ifdef SHOW_ITM
  658. 		printf("After fouri4ninv():\n");
  659. 		show_re_im_seq(d.data, d.p);
  660. #endif
  661.  
  662. 		const size_t s = std::min(a.data.size(), (size_t)(2 * t360));
  663.  
  664. 		// Carry and fix処理
  665. 		// ただし対象とする桁の並びd.data[]は、
  666. 		// 桁の値が[0, RADIX)ではなく、[-RADIX / 2, RADIX / 2)である。
  667. 		carry = 0;
  668. 		t = 1;
  669. 		//const double k = ldexp(1.0, -p + 1);		// k = 1.0 * 2^(-p + 1)
  670. 		double k = 1.0;		// fouri4ninv()結果は2^pで割る必要無し。
  671. 		for (i = 0; i < s; i++) {
  672. 			// 実部虚部隣接化解除
  673. 			// b.data[t360..]を虚部に配置する。
  674. 			const size_t src_i = (i < (size_t)t360) ? 2 * i : 2 * (i - (size_t)t360) + 1;
  675. 			//printf("t360=%lld, s=%zu, i=%zu, src_i=%zu\n", t360, s, i, src_i);
  676.  
  677. 			const double x = d.data[src_i] * k;
  678. 			int64_t e = (int64_t)floor(x + 0.5);
  679. 			err_max = std::max(err_max, fabs(x - e));
  680. 			e += carry;
  681. 			carry = e / RADIX;
  682. 			int64_t f = e - carry * RADIX;
  683. 			if (f < 0) {
  684. 				f += RADIX;
  685. 				carry--;
  686. 			}
  687. 			a.data[i] = (I_TYPE)f;
  688. 			if (f != 0) {
  689. 				t = i + 1;
  690. 			}
  691. 		}
  692. 		a.ndigits = t;
  693. 		if (carry != 0) {
  694. 			if (a.ndigits < a.data.size()) {
  695. 				assert(0 <= carry && carry < RADIX);
  696. 				a.data[a.ndigits++] = (I_TYPE)carry;
  697. 			}
  698. 			else {
  699. 				// ここに来たら桁あふれ
  700. 				abort();
  701. 			}
  702. 		}
  703. 		for (; i < s; i++) {
  704. 			// 実部虚部隣接化解除
  705. 			const size_t src_i = (i < (size_t)t360) ? 2 * i : 2 * i + 1;
  706.  
  707. 			const double x = d.data[src_i] * k;
  708. 			const int64_t e = (int64_t)floor(x + 0.5);
  709. 			err_max = std::max(err_max, fabs(x - e));
  710. 			assert(e == 0);
  711. 		}
  712. 	}
  713.  
  714. 	/// 整数の和
  715. 	void add_i(I& a, const I& b_, const I& c_) {
  716. 		I_TYPE carry, d;
  717. 		size_t i;
  718.  
  719. 		const I& b = (b_.ndigits > c_.ndigits) ? c_ : b_;	// Alias
  720. 		const I& c = (b_.ndigits > c_.ndigits) ? b_ : c_;	// Alias
  721. 		assert(b.ndigits <= c.ndigits);
  722. 		assert(c.ndigits <= a.data.size());
  723.  
  724. 		carry = 0;
  725. 		for (i = 0; i < b.ndigits; i++) {
  726. 			d = b.data[i] + c.data[i] + carry;
  727. 			carry = 0;
  728. 			if (d >= RADIX) {
  729. 				d -= RADIX;
  730. 				carry++;
  731. 			}
  732. 			a.data[i] = d;
  733. 		}
  734. 		for (; i < c.ndigits; i++) {
  735. 			d = c.data[i] + carry;
  736. 			carry = 0;
  737. 			if (d >= RADIX) {
  738. 				d -= RADIX;
  739. 				carry++;
  740. 			}
  741. 			a.data[i] = d;
  742. 		}
  743. 		a.ndigits = c.ndigits;
  744. 		if (a.ndigits < a.data.size()) {
  745. 			if (carry != 0) {
  746. 				a.data[a.ndigits++] = carry;
  747. 			}
  748. 		}
  749. 		else {
  750. 			if (carry != 0) {
  751. 				// ここに来たら桁あふれ
  752. 				abort();
  753. 			}
  754. 		}
  755. 	}
  756.  
  757. 	/// 整数のファイル出力
  758. 	void print(const char* const filename, const I& a)
  759. 	{
  760. 		char format[10];
  761. 		int nRet;
  762.  
  763. 		assert(a.data.size() > (size_t)0);
  764. 		assert(a.ndigits <= a.data.size());
  765.  
  766. 		FILE* fp;
  767. 		nRet = fopen_s_wrp(&fp, filename, "w");
  768. 		if (nRet != 0) {
  769. 			fflush(stdout);
  770. 			fprintf(stderr, "ERROR: fopen_s_wrp() failed. (nRet=%d)\n", nRet);
  771. 			abort();
  772. 		}
  773.  
  774. 		sprintf_s_wrp(format, _countof(format), "%%0%du", FIG);
  775.  
  776. 		size_t s = a.ndigits;
  777. 		fprintf(fp, "%lu", a.data[--s]);
  778. 		while (s != 0) {
  779. 			fprintf(fp, format, a.data[--s]);
  780. 		}
  781. 		fclose(fp);
  782. 	}
  783.  
  784. }	// namespace oper
  785.  
  786. namespace
  787. {
  788. 	// N ≧ log2(s1 * s2)を満たす最小のNを求める。
  789. 	int s2p(const size_t s1, const size_t s2) {
  790. 		int i;
  791. 		for (i = 1; i < 64; i++) {
  792. 			if (((size_t)1 << i) + 1 >= s1 + s2) {
  793. 				// (2^i + 1 ≧ (s1の有効桁数) + (s2の有効桁数))
  794. 				break;
  795. 			}
  796. 		}
  797. 		return i;
  798. 	}
  799.  
  800. }	// namespace
  801.  
  802. int main()
  803. {
  804. 	using namespace oper;
  805.  
  806. 	const int n = 34;
  807.  
  808. 	nth_root w;
  809.  
  810. 	const int max_p = s2p(((size_t)((1LLU << n) * log10(2)) / FIG + 1), 1);
  811. 	const size_t cap = (size_t)1 << max_p;
  812.  
  813. 	I a(cap);
  814. 	F x(max_p);
  815.  
  816. 	a.data[0] = 1;
  817. 	a.ndigits = 1;
  818.  
  819. 	double err_max = 0.0;
  820. 	for (int i = 0; i < n; i++) {
  821.  
  822. 		// 多倍長整数xのフーリエ変換
  823. 		// p ≧ log2(a * a)を満たす最小のpを求め、
  824. 		// aの次数2^pのフーリエ変換結果をxに格納する。
  825. 		const int p = s2p(a.ndigits, a.ndigits);
  826. 		i2f(x, a, p, w);
  827.  
  828. 		// x = x*x、および結果のxを整数に戻してaに格納する
  829. 		mul_f2i(a, x, x, x, err_max, w);
  830. 		add_i(a, a, a);
  831.  
  832. 		// 途中経過出力
  833. 		{
  834. 			char filename[200];
  835. 			sprintf_s_wrp(filename, _countof(filename), "%d.txt", i + 1);
  836. 			print(filename, a);
  837. 		}
  838. 		printf("%d %f\n", i + 1, err_max);
  839. 	}
  840.  
  841. 	return 0;
  842. }
  843.  
Runtime error #stdin #stdout 0.69s 2095356KB
comments (?)
stdin 
Standard input is empty
stdout 
Standard output is empty
https://ideone.com/q0O74K
language:
C++14 (clang 8.0)
created:
3 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

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