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. #ifdef _WIN32
  22. int fopen_s_wrp(FILE** const p_fp, const char* const fpath, const char* const mode)
  23. {
  24. 	return (int)fopen_s(p_fp, fpath, mode);
  25. }
  26. int sprintf_s_wrp(char buf[], const size_t bufsz, const char* const fmt, ...)
  27. {
  28. 	va_list ap;
  29. 	int nRet;
  30.  
  31. 	va_start(ap, fmt);
  32. 	nRet = vsprintf_s(buf, bufsz, fmt, ap);
  33. 	va_end(ap);
  34.  
  35. 	return nRet;
  36. }
  37. #else
  38. int fopen_s_wrp(FILE** const p_fp, const char* const fpath, const char* const mode)
  39. {
  40. 	*p_fp = fopen(fpath, mode);
  41. 	return (*p_fp == NULL);
  42. }
  43. int sprintf_s_wrp(char buf[], const size_t bufsz, const char* const fmt, ...)
  44. {
  45. 	va_list ap;
  46. 	int nRet;
  47.  
  48. 	va_start(ap, fmt);
  49. 	nRet = vsprintf(buf, fmt, ap);
  50. 	va_end(ap);
  51.  
  52. 	return nRet;
  53. }
  54. #define _countof(arr)	(sizeof(arr) / sizeof(arr[0]))
  55. #endif
  56.  
  57. namespace fft
  58. {
  59. 	void fr(
  60. 		std::vector<double>::iterator s1,
  61. 		std::vector<double>::iterator s2,
  62. 		std::vector<double>::iterator end,
  63. 		const size_t add,
  64. 		const size_t t90,
  65. 		const double sin1, const double cos1)
  66. 	{
  67. 		for (; s1 < end; (s1 += add, s2 += add)) {
  68. #ifdef SHOW_ITM
  69. 			printf("%s: end-s1=%lld, add=%zu\n", __FUNCTION__, end - s1, add);
  70. #endif
  71. 			const double r0 = s1[0];
  72. 			const double r1 = s2[0];
  73. 			const double r2 = s1[t90];
  74. 			const double r3 = s2[t90];
  75.  
  76. 			// 1のN乗根wを cos1 + i * sin1 として、次式を計算
  77. 			// (r2 + i * r3) * w
  78. 			//   = r2 * cos1 + r2 * i * sin1 + r3 * i * cos1 - r3 * sin1
  79. 			//   = (r2 * cos1 - r3 * sin1) + i * (r2 * sin1 + r3 * cos1)
  80. 			const double m0 = r2 * cos1 - r3 * sin1;		// 実部
  81. 			const double m1 = r2 * sin1 + r3 * cos1;		// 虚部
  82.  
  83. 			s1[0] = r0 + m0;
  84. 			s2[0] = r0 - m0;
  85. 			s1[t90] = -r1 + m1;
  86. 			s2[t90] = r1 + m1;
  87. 		}
  88. 	}
  89.  
  90. 	void fri(
  91. 		std::vector<double>::iterator s1,
  92. 		std::vector<double>::iterator s2,
  93. 		std::vector<double>::iterator end,
  94. 		const size_t add,
  95. 		const size_t t90,
  96. 		const double sin1, const double cos1)
  97. 	{
  98. 		for (; s1 < end; (s1 += add, s2 += add)) {
  99. 			const double r0 = s1[0];
  100. 			const double r1 = s2[0];
  101. 			const double r2 = s1[t90];
  102. 			const double r3 = s2[t90];
  103.  
  104. 			const double m0 = r0 + r1;
  105. 			const double m1 = -r2 + r3;
  106. 			const double m2 = r0 - r1;
  107. 			const double m3 = r2 + r3;
  108.  
  109. 			s1[0] = m0;
  110. 			s2[0] = m1;
  111. 			s1[t90] = m2 * cos1 + m3 * sin1;
  112. 			s2[t90] = -m2 * sin1 + m3 * cos1;
  113. 		}
  114. 	}
  115.  
  116. 	void fr0(
  117. 		std::vector<double>::iterator s1,
  118. 		std::vector<double>::iterator end,
  119. 		const size_t add,
  120. 		const size_t t90)
  121. 	{
  122. 		for (; s1 < end; s1 += add) {
  123. 			const double r0 = s1[0];
  124. 			const double r1 = s1[t90];
  125.  
  126. 			s1[0] = r0 + r1;
  127. 			s1[t90] = r0 - r1;
  128. 		}
  129. 	}
  130.  
  131. 	void bxchg(
  132. 		std::vector<double>& s,
  133. 		const int p)
  134. 	{
  135. 		// 2^pを8等分した境界位置
  136. 		const size_t t045 = (size_t)1 << (p - 3);	// 1 * 2^(p-3)
  137. 		const size_t t090 = t045 + t045;			// 2 * 2^(p-3)
  138. 		const size_t t180 = t090 + t090;			// 4 * 2^(p-3)
  139. 		const size_t t270 = t180 + t090;			// 6 * 2^(p-3)
  140. 		/*
  141. 			(計算例)
  142. 			p=8のとき
  143. 			t045 =  32 = B'00100000
  144. 			t090 =  64 = B'01000000
  145. 			t180 = 128 = B'10000000
  146. 			t270 = 192 = B'11000000
  147. 			t4 = t090 - 4 = 60
  148. 			t5 = t4       = 60
  149. 		*/
  150.  
  151. 		size_t t4 = t090 - 4;		// 2 * 2^(p-3) - 4	--- p=4のとき0、p=5のとき4、p=6のとき28
  152. 		size_t t5 = t4;				// 2 * 2^(p-3) - 4  --- (同上)
  153. 		for (;;) {
  154. 			std::swap(s[t4 + 1], s[t180 + t5]);		// 180 degを挟んだ2点を交換
  155. 			std::swap(s[t4 + 2], s[t090 + t5]);
  156. 			std::swap(s[t4 + 3], s[t270 + t5]);
  157. 			std::swap(s[t090 + t4 + 1], s[t180 + t5 + 2]);
  158. 			std::swap(s[t090 + t4 + 3], s[t270 + t5 + 2]);
  159. 			std::swap(s[t180 + t4 + 3], s[t270 + t5 + 1]);
  160.  
  161. 			if (t4 < t5) {
  162. 				std::swap(s[t4], s[t5]);
  163. 				std::swap(s[t090 + t4 + 2], s[t090 + t5 + 2]);
  164. 				std::swap(s[t180 + t4 + 1], s[t180 + t5 + 1]);
  165. 				std::swap(s[t270 + t4 + 3], s[t270 + t5 + 3]);
  166. 			}
  167. 			if (t4 == 0) {
  168. 				break;
  169. 			}
  170.  
  171. 			t4 -= 4;
  172. 			size_t t6 = t045;
  173. 			for (;;) {
  174. 				t5 ^= t6;
  175. 				if ((t5 & t6) == 0) {
  176. 					break;
  177. 				}
  178. 				t6 >>= 1;
  179. 			}
  180. 		}
  181. 	}
  182.  
  183. 	void bitexchg(
  184. 		std::vector<double>& s,
  185. 		const int p)
  186. 	{
  187. 		if (p <= 1) {
  188. 			/*pass*/
  189. 		}
  190. 		else if (p == 2) {
  191. 			std::swap(s[1], s[2]);
  192. 		}
  193. 		else if (p == 3) {
  194. 			std::swap(s[1], s[4]);
  195. 			std::swap(s[3], s[6]);
  196. 		}
  197. 		else {
  198. 			bxchg(s, p);
  199. 		}
  200. 	}
  201.  
  202. #ifdef SHOW_SIN_TABLE
  203. 	size_t g_sin_tblent_cnt = 0;
  204. #endif
  205.  
  206. 	/// sin関数のテーブル化
  207. 	/// SIN_TABLE[p]に以下のリストを設定する。
  208. 	/// SIN_TABLE[2] = [ sin(0) ]
  209. 	/// SIN_TABLE[3] = [ sin(0), sin(45 deg) ]
  210. 	/// SIN_TABLE[4] = [ sin(0), sin(90 / 4 deg), sin(45 deg), sin(90 * 3 / 4 deg) ]
  211. 	/// ...
  212. 	void sin_table(
  213. 		const int p,
  214. 		std::vector<std::vector<double> >& SIN_TABLE)
  215. 	{
  216. 		const size_t sv_sz = SIN_TABLE.size();
  217. 		if (p >= SIN_TABLE.size()) {
  218. 			// (pに対するテーブルが未作成)
  219. 			// テーブルを収める領域を確保
  220. 			SIN_TABLE.resize((size_t)(p + 1));
  221. 		}
  222. 		std::vector<double>& sin_tbl = SIN_TABLE[p];		// Alias
  223.  
  224. 		if (p >= 2 && sv_sz <= p) {
  225. 			// (pが2以上かつpに対するテーブルが未作成だった)
  226. 			// pに対するテーブルの領域を確保し、数値格納
  227. 			const size_t t090 = (size_t)1 << (p - 2);
  228. 			sin_tbl.resize(t090);
  229. 			const double k = 0.5 * M_PI;
  230. 			for (size_t i = 0; i < t090; i++) {
  231. 				sin_tbl[i] = sin((k * i) / t090);
  232. 			}
  233. #ifdef SHOW_SIN_TABLE
  234. 			g_sin_tblent_cnt += t090;
  235. 			printf("******* %s: p=%d table made. (+%zu, %zu; sv_sz=%zu)\n", __FUNCTION__, p, t090, g_sin_tblent_cnt, sv_sz);
  236. #endif
  237. 		}
  238. 	}
  239.  
  240. 	/// フーリエ変換
  241. 	void fouri(
  242. 		std::vector<double>& s,
  243. 		const int p,
  244. 		std::vector<std::vector<double> >& SIN_TABLE)
  245. 	{
  246. #ifdef SHOW_ITM
  247. 		printf("%s: p=%d: Called.\n", __FUNCTION__, p);
  248. #endif
  249. 		sin_table(p, SIN_TABLE);
  250. 		const std::vector<double>& table = SIN_TABLE[p];
  251. 		const size_t t360 = (size_t)1 << p;				// = 2^p
  252. 		const size_t t090 = t360 >> 2;					// = 2^(p-1)
  253.  
  254. 		bitexchg(s, p);
  255.  
  256. 		for (int i = 1; i <= p; i++) {					// i = 1, 2, 3, ..., p
  257. 			const size_t r360 = (size_t)1 << i;			// = 2^i
  258. 			const size_t r180 = r360 >> 1;				// = 2^(i-1)
  259. 			const size_t r090 = r180 >> 1;				// = 2^(i-2)
  260.  
  261. 			fr0(s.begin(), s.begin() + t360, r360, r180);
  262.  
  263. 			size_t k1 = 1;
  264. 			size_t k2 = r180 - 1;						// = 2^(i-1)-1
  265. 			const size_t t2 = (size_t)1 << (p - i);		// = 2^(p-i) = (2^p) / (2^i) = t360 / r360 (幾何学的意味)
  266. 			size_t t0 = t2;								// = t360 / r360 (幾何学的意味)
  267. 			size_t t1 = t090 - t2;						// = t090 - t360 / r360 (幾何学的意味)
  268. 			for (; k1 < r090; (k1++, k2--)) {			// k1 = 1, 2, 3, ..., 2^(i-2)-1; k2 = 2^(i-1)-1, 2^(i-1)-2, 2^(i-1)-3, ..., 2^(i-2)+1
  269. 				const double sin1 = table[t0];			// 1の(4*p)乗根の(2^(p-i))乗の虚部
  270. 				const double cos1 = table[t1];			// 1の(4*p)乗根の(2^(p-i))乗の実部
  271. #ifdef SHOW_ITM
  272. 				{
  273. 					const double degree = (2 * M_PI) / atan2(sin1, cos1);
  274. 					printf("%s: i=%d, k1=%zu, k2=%zu, r090=%zu, w^%f=1, t0=%zu, t1=%zu\n",
  275. 						__FUNCTION__,
  276. 						i, k1, k2, r090, degree,
  277. 						t0, t1);
  278. 					const double ideg = floor(degree + 0.5);
  279. 					if (std::abs(degree - ideg) > 0.1) {
  280. 						printf("*** fraction ***\n");
  281. 					}
  282. 				}
  283. #endif
  284. 				fr(s.begin() + k1, s.begin() + k2, s.begin() + t360, r360, r180, sin1, cos1);
  285. 				t0 += t2;	// t0 = t360/r360, (2*t360)/r360, (3*t360)/r360, ...,  ((2^(i-2)-1)*2^p)/2^i, ..., ((2^(i-2)-1)*t360)/r360
  286. 				t1 -= t2;
  287. 			}
  288. 		}
  289. #ifdef SHOW_ITM
  290. 		printf("%s: p=%d: Done.\n", __FUNCTION__, p);
  291. #endif
  292. 	}
  293.  
  294. 	/// 逆フーリエ変換
  295. 	void fouriinv(
  296. 		std::vector<double>& s,
  297. 		const int p,
  298. 		std::vector<std::vector<double> >& SIN_TABLE)
  299. 	{
  300. 		sin_table(p, SIN_TABLE);
  301. 		const std::vector<double>& table = SIN_TABLE[p];
  302. 		const size_t t360 = (size_t)1 << p;
  303. 		const size_t t090 = t360 >> 2;
  304.  
  305. 		for (int i = p; i >= 1; i--) {
  306. 			const size_t r360 = (size_t)1 << i;
  307. 			const size_t r180 = r360 >> 1;
  308. 			const size_t r090 = r180 >> 1;
  309.  
  310. 			fr0(s.begin(), s.begin() + t360, r360, r180);
  311.  
  312. 			size_t k1 = 1;
  313. 			size_t k2 = r180 - 1;
  314. 			size_t t2 = (size_t)1 << (p - i);
  315. 			size_t t0 = t2;
  316. 			size_t t1 = t090 - t2;
  317. 			for (; k1 < r090; k1++, k2--) {
  318. 				const double sin1 = table[t0];
  319. 				const double cos1 = table[t1];
  320. 				fri(s.begin() + k1, s.begin() + k2, s.begin() + t360, r360, r180, sin1, cos1);
  321. 				t0 += t2;
  322. 				t1 -= t2;
  323. 			}
  324. 		}
  325. 		bitexchg(s, p);
  326. 	}
  327.  
  328. }	// namespace fft
  329.  
  330. typedef uint32_t I_TYPE;
  331.  
  332. const int FIG = 3;
  333. const I_TYPE RADIX = 1000;
  334.  
  335. namespace oper
  336. {
  337. 	using namespace fft;
  338.  
  339. 	struct I
  340. 	{
  341. 		std::vector<I_TYPE> data;	// 基数RADIXの多倍長データ
  342. 		size_t ndigits;				// dataに入っている有効な桁数
  343.  
  344. 		I(const size_t cap)
  345. 			: data(cap), ndigits((size_t)0)
  346. 		{
  347. 			/*pass*/
  348. 		}
  349. 	};
  350.  
  351. 	struct F
  352. 	{
  353. 		std::vector<double> data;	// 基数[-(RADIX/2), RADIX/2)の多倍長データ
  354. 		int p;						// 2^p == (dataの要素数) となるp
  355.  
  356. 		F(const int p_)
  357. 			: data((size_t)1 << p_), p(0)
  358. 		{
  359. 			/*pass*/
  360. 		}
  361. 	};
  362.  
  363. 	/// 整数のフーリエ変換
  364. 	void i2f(F& a, const I& b, const int p, std::vector<std::vector<double> >& SIN_TABLE)
  365. 	{
  366. 		size_t i;
  367.  
  368. 		const size_t t360 = (size_t)1 << p;
  369.  
  370. 		assert(b.ndigits <= t360);
  371. 		assert(t360 <= a.data.size());
  372.  
  373. #ifdef SHOW_I
  374. 		printf("%s: Before pre tr.:\n", __FUNCTION__);
  375. 		for (size_t ii = 0; ii < b.ndigits; ii++) {
  376. 			printf(" %03ld", b.data[ii]);
  377. 		}
  378. 		printf("\n");
  379. 		for (size_t ii = b.ndigits; ii < b.data.size(); ii++) {
  380. 			if (b.data[ii] != 0) {
  381. 				abort();
  382. 			}
  383. 		}
  384. #endif
  385.  
  386. 		// フーリエ変換の前処理
  387. 		// (bが表す整数) = Σ[i=0..a.ndigits-1]{ b.data[i] * RADIX^i } (0≦b.data[i]＜RADIX)
  388. 		// を、
  389. 		// (aが表す整数) = Σ[i=0..a.ndigits-1]{ a.data[i] * RADIX^i } (-RADIX / 2≦a.data[i]＜RADIX / 2)
  390. 		// に変換する。(a, bが同じ整数を表すように桁上げも行う。)
  391. 		a.p = p;
  392. 		int carry = 0;
  393. 		for (i = 0; i < b.ndigits - 1; i++) {
  394. 			const int32_t t = b.data[i] + carry;	// 桁[i]の値
  395. 			if (t < (int32_t)RADIX / 2) {
  396. 				// (桁[i]の値∈[0, RADIX / 2) )
  397. 				// 何もしない
  398. 				a.data[i] = t;
  399. 				carry = 0;
  400. 			}
  401. 			else {
  402. 				// (桁[i]の値∈[RADIX / 2, RADIX) )
  403. 				// 補正実施
  404. 				// まず桁[i] -= RADIX により桁[i]の値∈[0, RADIX / 2)とする。ただし
  405. 				// それだけでは(aが表す整数)が RADIX * RADIX^i = RADIX^(i+1)だけ減ってしまうので、
  406. 				// 減った分をcarryで一つ上の桁に回して補填する。
  407. 				a.data[i] = (double)(t - (int32_t)RADIX);
  408. 				carry = 1;
  409. 			}
  410. 		}
  411. 		// 最上位桁
  412. 		a.data[i] = (double)(b.data[i] + (I_TYPE)carry);
  413. 		i++;
  414. 		for (; i < t360; i++) {
  415. 			a.data[i] = 0.0;
  416. 		}
  417.  
  418. #ifdef SHOW_I
  419. 		{
  420. 			size_t n;
  421. 			for (n = a.data.size(); n > 0; n--) {
  422. 				if (std::abs(a.data[n - 1]) >= 0.000001) {
  423. 					break;
  424. 				}
  425. 			}
  426. 			printf("%s: After pre tr.:\n", __FUNCTION__);
  427. 			for (size_t ii = 0; ii < n; ii++) {
  428. 				printf(" %03.0f", a.data[ii]);
  429. 			}
  430. 			printf("\n");
  431. 		}
  432. #endif
  433.  
  434. 		// フーリエ変換実施
  435. 		fouri(a.data, p, SIN_TABLE);
  436. 	}
  437.  
  438. 	/// 整数のフーリエ変換結果の積(結果は整数)
  439. 	void mul_f2i(I& a, const F& b, const F& c, F& d, double& err_max, std::vector<std::vector<double> >& SIN_TABLE)
  440. 	{
  441. 		size_t i, t, t180, t360, k1, k2;
  442. 		int64_t carry;
  443.  
  444. 		const int p = b.p;
  445. 		t180 = (size_t)1 << (p - 1);
  446. 		t360 = t180 + t180;
  447.  
  448. 		assert(p == c.p);
  449. 		assert(t360 <= b.data.size());
  450. 		assert(t360 <= c.data.size());
  451. 		assert(t360 <= d.data.size());
  452.  
  453. 		d.data[0] = b.data[0] * c.data[0] * 0.5;
  454. 		d.data[t180] = b.data[t180] * c.data[t180] * 0.5;
  455.  
  456. 		k1 = 1;
  457. 		k2 = t360 - 1;
  458. 		for (; k1 < t180; (k1++, k2--)) {
  459. 			const double re = b.data[k1] * c.data[k1] - b.data[k2] * c.data[k2];
  460. 			const double im = b.data[k1] * c.data[k2] + b.data[k2] * c.data[k1];
  461. 			d.data[k1] = re;
  462. 			d.data[k2] = im;
  463. 		}
  464.  
  465. 		// 逆フーリエ変換
  466. 		fouriinv(d.data, p, SIN_TABLE);
  467.  
  468. 		const size_t s = std::min(a.data.size(), t360);
  469.  
  470. 		// Carry and fix処理
  471. 		// ただし対象とする桁の並びd.data[]は、
  472. 		// 桁の値が[0, RADIX)ではなく、[-RADIX / 2, RADIX / 2)である。
  473. 		carry = 0;
  474. 		t = 1;
  475. 		const double k = ldexp(1.0, -p + 1);		// k = 1.0 * 2^(-p + 1)
  476. 		for (i = 0; i < s; i++) {
  477. 			const double x = d.data[i] * k;
  478. 			int64_t e = (int64_t)floor(x + 0.5);
  479. 			err_max = std::max(err_max, fabs(x - e));
  480. 			e += carry;
  481. 			carry = e / RADIX;
  482. 			int64_t f = e - carry * RADIX;
  483. 			if (f < 0) {
  484. 				f += RADIX;
  485. 				carry--;
  486. 			}
  487. 			a.data[i] = (I_TYPE)f;
  488. 			if (f != 0) {
  489. 				t = i + 1;
  490. 			}
  491. 		}
  492. 		a.ndigits = t;
  493. 		if (carry != 0) {
  494. 			if (a.ndigits < a.data.size()) {
  495. 				assert(0 <= carry && carry < RADIX);
  496. 				a.data[a.ndigits++] = (I_TYPE)carry;
  497. 			}
  498. 			else {
  499. 				// ここに来たら桁あふれ
  500. 				abort();
  501. 			}
  502. 		}
  503. 		for (; i < t360; i++) {
  504. 			const double x = d.data[i] * k;
  505. 			const int64_t e = (int64_t)floor(x + 0.5);
  506. 			err_max = std::max(err_max, fabs(x - e));
  507. 			assert(e == 0);
  508. 		}
  509. 	}
  510.  
  511. 	/// 整数の和
  512. 	void add_i(I& a, const I& b_, const I& c_) {
  513. 		I_TYPE carry, d;
  514. 		size_t i;
  515.  
  516. 		const I& b = (b_.ndigits > c_.ndigits) ? c_ : b_;	// Alias
  517. 		const I& c = (b_.ndigits > c_.ndigits) ? b_ : c_;	// Alias
  518. 		assert(b.ndigits <= c.ndigits);
  519. 		assert(c.ndigits <= a.data.size());
  520.  
  521. 		carry = 0;
  522. 		for (i = 0; i < b.ndigits; i++) {
  523. 			d = b.data[i] + c.data[i] + carry;
  524. 			carry = 0;
  525. 			if (d >= RADIX) {
  526. 				d -= RADIX;
  527. 				carry++;
  528. 			}
  529. 			a.data[i] = d;
  530. 		}
  531. 		for (; i < c.ndigits; i++) {
  532. 			d = c.data[i] + carry;
  533. 			carry = 0;
  534. 			if (d >= RADIX) {
  535. 				d -= RADIX;
  536. 				carry++;
  537. 			}
  538. 			a.data[i] = d;
  539. 		}
  540. 		a.ndigits = c.ndigits;
  541. 		if (a.ndigits < a.data.size()) {
  542. 			if (carry != 0) {
  543. 				a.data[a.ndigits++] = carry;
  544. 			}
  545. 		}
  546. 		else {
  547. 			// ここに来たら桁あふれ
  548. 			abort();
  549. 		}
  550. 	}
  551.  
  552. 	/// 整数のファイル出力
  553. 	void print(const char* const filename, const I& a)
  554. 	{
  555. 		char format[10];
  556. 		int nRet;
  557.  
  558. 		assert(a.data.size() > (size_t)0);
  559. 		assert(a.ndigits <= a.data.size());
  560.  
  561. 		FILE* fp;
  562. 		nRet = fopen_s_wrp(&fp, filename, "w");
  563. 		if (nRet != 0) {
  564. 			fflush(stdout);
  565. 			fprintf(stderr, "ERROR: fopen_s_wrp() failed. (nRet=%d)\n", nRet);
  566. 			abort();
  567. 		}
  568.  
  569. 		sprintf_s_wrp(format, _countof(format), "%%0%du", FIG);
  570.  
  571. 		size_t s = a.ndigits;
  572. 		fprintf(fp, "%lu", a.data[--s]);
  573. 		while (s != 0) {
  574. 			fprintf(fp, format, a.data[--s]);
  575. 		}
  576. 		fclose(fp);
  577. 	}
  578.  
  579. }	// namespace oper
  580.  
  581. namespace
  582. {
  583. 	// N ≧ log2(s1 * s2)を満たす最小のNを求める。
  584. 	int s2p(const size_t s1, const size_t s2) {
  585. 		int i;
  586. 		for (i = 1; i < 64; i++) {
  587. 			if (((size_t)1 << i) + 1 >= s1 + s2) {
  588. 				// (2^i + 1 ≧ (s1の有効桁数) + (s2の有効桁数))
  589. 				break;
  590. 			}
  591. 		}
  592. 		return i;
  593. 	}
  594.  
  595. }	// namespace
  596.  
  597. bool fft_test()
  598. {
  599. 	using namespace fft;
  600. 	using namespace oper;
  601.  
  602. 	const int max_p = 5;
  603. 	const size_t cap = (size_t)1 << max_p;
  604.  
  605. 	I a(cap);
  606. 	F x(max_p);
  607.  
  608. 	x.p = max_p;
  609. 	for (size_t i = 0; i < cap; i++) {
  610. 		x.data[i] = (double)i;
  611. 	}
  612.  
  613. 	// フーリエ変換実施
  614. 	std::vector<std::vector<double> > SIN_TABLE;	// working object for fouri()
  615. 	fouri(x.data, x.p, SIN_TABLE);
  616.  
  617. 	return true;
  618. }
  619.  
  620. int main()
  621. {
  622. 	using namespace oper;
  623.  
  624. #ifdef SHOW_ITM
  625. 	if (fft_test()) {
  626. 		return 0;
  627. 	}
  628. #endif
  629.  
  630. 	std::vector<std::vector<double> > SIN_TABLE;
  631.  
  632. 	const int n = 34;
  633.  
  634. 	const int max_p = s2p(((size_t)((1LLU << n) * log10(2)) / FIG + 1), 1);
  635. 	const size_t cap = (size_t)1 << max_p;
  636.  
  637. 	I a(cap);
  638. 	F x(max_p);
  639.  
  640. 	a.data[0] = 1;
  641. 	a.ndigits = 1;
  642.  
  643. 	double err_max = 0.0;
  644. 	for (int i = 0; i < n; i++) {
  645.  
  646. 		// 多倍長整数xのフーリエ変換
  647. 		// p ≧ log2(a * a)を満たす最小のpを求め、
  648. 		// aの次数2^pのフーリエ変換結果をxに格納する。
  649. 		const int p = s2p(a.ndigits, a.ndigits);
  650. 		i2f(x, a, p, SIN_TABLE);
  651.  
  652. 		// x = x*x、および結果のxを整数に戻してaに格納する
  653. 		mul_f2i(a, x, x, x, err_max, SIN_TABLE);
  654. 		add_i(a, a, a);
  655.  
  656. 		// 途中経過出力
  657. 		{
  658. 			char filename[200];
  659. 			sprintf_s_wrp(filename, _countof(filename), "%d.txt", i + 1);
  660. 			print(filename, a);
  661. 		}
  662. 		printf("%d %f\n", i + 1, err_max);
  663. 	}
  664.  
  665. 	return 0;
  666. }
  667.  
Runtime error #stdin #stdout 0.87s 2090640KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/oDgjXr
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!