fork(2) download 
  1. // 2^(expn2)を有効数字disp10桁の10進数で指数表示する。
  2. // Usage: print_by_decimal <expn2> <disp10>
  3. #include <stdio.h>
  4. #include <stdlib.h>
  5. #include <stdint.h>
  6. #include <assert.h>
  7. #include <vector>
  8. #include <utility>
  9.  
  10. #include <time.h>
  11.  
  12. // 指数の型
  13. typedef uint64_t expn_t;
  14.  
  15. // 2乗したものを足し合わせてもuint64_tが桁あふれしない10進数
  16. #define BIG_DECIMAL		10000000ULL
  17.  
  18. static_assert(UINT64_MAX - (BIG_DECIMAL - 1) * (BIG_DECIMAL - 1) >= 2 * (BIG_DECIMAL - 1), "*** ERR ***");
  19. static_assert(2 * (BIG_DECIMAL - 1) + (BIG_DECIMAL - 1) * (BIG_DECIMAL - 1) <= UINT64_MAX, "*** ERR ***");
  20.  
  21. namespace
  22. {
  23. 	struct BigFloat
  24. 	{
  25. 		std::vector<uint32_t> frac;
  26. 		expn_t expn;
  27.  
  28. 		BigFloat() : expn(0) { }
  29. 		BigFloat(const size_t sz, const uint32_t inival = 0)
  30. 			: frac(sz, inival), expn(0) { }
  31.  
  32. 		size_t size() const { return frac.size(); }
  33. 		uint32_t operator[](const size_t k) const { return frac[k]; }
  34. 	};
  35.  
  36.  
  37. 	/// 先頭の0を切り詰めたサイズを返す。
  38. 	size_t truncated_size(const std::vector<uint32_t>& mbuf)
  39. 	{
  40. 		const size_t sz = mbuf.size();
  41. 		for (size_t i = sz; i > 0; i--) {
  42. 			if (mbuf[i - 1] != 0UL) {
  43. 				return i;
  44. 			}
  45. 		}
  46. 		return sz;
  47. 	}
  48.  
  49. 	/// 計算結果の記憶桁数に制限を設けた乗算
  50. 	/// メモリ消費を半分にするため、キャリー伝搬を単純に桁の積和の都度行う。
  51. 	/// (最上位の桁が0になることがあるため、記憶桁数bufszは厳密には精度を意味しない。)
  52. 	std::vector<uint32_t> mlen_mul(const std::vector<uint32_t>& lhs, const std::vector<uint32_t>& rhs,
  53. 		const expn_t bufsz_, const uint32_t a_priori_cy,
  54. 		expn_t& expn)
  55. 	{
  56. 		const size_t sz_lhs = truncated_size(lhs);
  57. 		const size_t sz_rhs = truncated_size(rhs);
  58.  
  59. 		const size_t dstsz = sz_lhs + sz_rhs;
  60. 		const size_t bufsz = std::min(dstsz, bufsz_);
  61. 		const size_t capable_min_idx = dstsz - bufsz;
  62.  
  63. 		// 指数出力
  64. 		expn = capable_min_idx;
  65.  
  66. 		std::vector<uint32_t> mbuf(bufsz, 0UL);
  67. 		for (expn_t i = 0; i < sz_lhs; i++) {
  68. 			uint32_t cy = 0UL;
  69. 			for (expn_t j = 0; j < sz_rhs; j++) {
  70. 				// 計算結果を格納すべき本来の位置
  71. 				const size_t dst_idx = i + j;
  72.  
  73. 				if (dst_idx >= capable_min_idx) {
  74. 					if (dst_idx == capable_min_idx && j > 0) {
  75. 						// アプリオリなキャリー設定
  76. 						// [capable_min_idx]より小さい桁の計算を行わないため、
  77. 						// mbuf[]の最下位桁に繰り上がる真のキャリーはわからない。
  78. 						cy = a_priori_cy;
  79. 					}
  80. 					const size_t k = dst_idx - capable_min_idx;
  81. 					const uint64_t val64 = (uint64_t)mbuf[k] + (uint64_t)lhs[i] * (uint64_t)rhs[j] + (uint64_t)cy;
  82. 					cy = (uint32_t)(val64 / BIG_DECIMAL);
  83. 					mbuf[k] = (uint32_t)(val64 - cy * BIG_DECIMAL);
  84. 				}
  85. 			}
  86. 			{
  87. 				const expn_t j = sz_rhs;
  88.  
  89. 				// 計算結果を格納すべき本来の位置
  90. 				const size_t dst_idx = i + j;
  91.  
  92. 				if (dst_idx >= capable_min_idx) {
  93. 					if (dst_idx == capable_min_idx) {
  94. 						// アプリオリなキャリー設定
  95. 						// [capable_min_idx]より小さい桁の計算を行わないため、
  96. 						// mbuf[]の最下位桁に繰り上がる真のキャリーはわからない。
  97. 						cy = a_priori_cy;
  98. 					}
  99. 					const size_t k = dst_idx - capable_min_idx;
  100. 					const uint64_t val64 = (uint64_t)mbuf[k] + (uint64_t)cy;
  101. 					cy = (uint32_t)(val64 / BIG_DECIMAL);
  102. 					mbuf[k] = (uint32_t)(val64 - cy * BIG_DECIMAL);
  103. 					assert(cy == 0);
  104. 				}
  105. 			}
  106. 		}
  107.  
  108. 		return mbuf;
  109. 	}
  110.  
  111. 	/// 浮動小数点演算による乗算
  112. 	BigFloat mlen_mul(const BigFloat& lhs, const BigFloat& rhs,
  113. 		const expn_t bufsz, const uint32_t a_priori_cy)
  114. 	{
  115. 		BigFloat result;
  116. 		expn_t res_expn;
  117. 		result.frac = mlen_mul(lhs.frac, rhs.frac, bufsz, a_priori_cy, res_expn);
  118. 		result.expn = lhs.expn + rhs.expn + res_expn;
  119. 		return result;
  120. 	}
  121.  
  122. 	/// 2^nをBIG_DECIMAL進数で計算する。
  123. 	BigFloat pow2(const expn_t n_, const size_t bufsz, const uint32_t a_priori_cy)
  124. 	{
  125. 		const time_t startTmg = time(NULL);
  126.  
  127. 		BigFloat mbuf((size_t)1, 1UL);
  128. 		BigFloat mbuf_work((size_t)1, 2UL);
  129. 		expn_t n = n_;
  130. 		for (;;) {
  131. 			if ((n & 0x1) != 0) {
  132. 				mbuf = mlen_mul(mbuf, mbuf_work, bufsz, a_priori_cy);
  133. 			}
  134. 			const time_t now = time(NULL);
  135. 			printf("n=%llu, mbuf sz=%zu --- %lld sec\n", n, mbuf.size(), now - startTmg);
  136. 			n >>= 1;
  137. 			if (n == 0ULL) {
  138. 				break;
  139. 			}
  140. 			mbuf_work = mlen_mul(mbuf_work, mbuf_work, bufsz, a_priori_cy);
  141. 		}
  142. 		return mbuf;
  143. 	}
  144.  
  145. 	/// 計算結果をファイルに出力する。
  146. 	void BigFloat_op_dump(const BigFloat& mbuf, FILE* const fp)
  147. 	{
  148. 		const size_t sz = mbuf.size();
  149. 		const size_t sqsz = ((sz + 9) / 10) * 10;
  150.  
  151. 		size_t cnt = 0;
  152. 		for (size_t i = sqsz; i > 0; i--) {
  153. 			if (i > sz) {
  154. 				fprintf(fp, " %7s", " ");
  155. 			}
  156. 			else {
  157. 				fprintf(fp, " %07lu", mbuf[i - 1]);
  158. 			}
  159. 			cnt++;
  160. 			if (cnt % 10 == 0) {
  161. 				if (i != 1) {
  162. 					fprintf(fp, "\n");
  163. 				}
  164. 				else {
  165. 					fprintf(fp, ". x 10^%llu\n", 7ULL * mbuf.expn);
  166. 				}
  167. 			}
  168. 		}
  169. 	}
  170.  
  171. #ifdef _WIN32
  172. 	/// 計算結果をファイルに出力する。
  173. 	bool write_to_file(const BigFloat& mbuf, const char* const fpath)
  174. 	{
  175. 		FILE* fp;
  176. 		errno_t errt;
  177.  
  178. 		errt = fopen_s(&fp, fpath, "w");
  179. 		if (errt != 0) {
  180. 			return false;
  181. 		}
  182.  
  183. 		BigFloat_op_dump(mbuf, fp);
  184.  
  185. 		fclose(fp);
  186. 		return true;
  187. 	}
  188. #endif
  189.  
  190. }	// namespace
  191.  
  192. int main(const int argc, char* argv[])
  193. {
  194. 	// 2のべきの指数
  195. 	expn_t expn2 = 17179869183LL;
  196. 	//expn_t expn2 = 171798;
  197. 	//expn_t expn2 = 1000;
  198. 	//expn_t expn2 = 64;
  199. 	//expn_t expn2 = 32;
  200. 	//expn_t expn2 = 16;
  201. 	//expn_t expn2 = 8;
  202. 	//expn_t expn2 = 0;
  203.  
  204. 	// 表示桁数(有効数字)
  205. 	expn_t disp10 = 64;
  206.  
  207. #ifdef _WIN32
  208. 	if (argc >= 2) {
  209. 		expn2 = _atoi64(argv[1]);
  210. 	}
  211. 	if (argc >= 3) {
  212. 		disp10 = _atoi64(argv[2]);
  213. 	}
  214. #else
  215. 	(void)(argc, argv);
  216. #endif
  217. 	printf("Calculating 2^%llu ...\n", expn2);
  218.  
  219. 	// テスト
  220. 	//BigFloat lhs((size_t)1, 9000000UL);
  221. 	//BigFloat rhs((size_t)1, 9000000UL);
  222. 	//BigFloat result = mlen_mul(lhs, rhs, 10, 0);
  223. 	//printf("result=");
  224. 	//BigFloat_op_dump(result, stdout);
  225.  
  226. 	const size_t bufsz = (disp10 + 6) / 7;
  227.  
  228. 	// 2^(expn2)の下限
  229. 	const BigFloat minv = pow2(expn2, bufsz, 0UL);
  230. 	BigFloat_op_dump(minv, stdout);
  231.  
  232. 	// 2^(expn2)の上限
  233. 	const BigFloat maxv = pow2(expn2, bufsz, 1UL);
  234. 	BigFloat_op_dump(maxv, stdout);
  235.  
  236. 	// 最終結果表示
  237. 	printf("\n\n\n");
  238. 	printf("2^%llu is between\n", expn2);
  239. 	BigFloat_op_dump(minv, stdout);
  240. 	BigFloat_op_dump(maxv, stdout);
  241.  
  242. #ifdef _WIN32
  243. 	// ファイル出力
  244. 	{
  245. 		char buf[_MAX_PATH];
  246. 		sprintf_s(buf, _countof(buf), "min_of_pow2(%llu)_in_decimal.txt", expn2);
  247. 		write_to_file(minv, buf);
  248. 		sprintf_s(buf, _countof(buf), "max_of_pow2(%llu)_in_decimal.txt", expn2);
  249. 		write_to_file(maxv, buf);
  250. 	}
  251. #endif
  252. }
  253.  
Success #stdin #stdout 0.01s 5368KB
comments (?)
stdin 
Standard input is empty
stdout 
Calculating 2^17179869183 ...
n=17179869183, mbuf sz=2 --- 0 sec
n=8589934591, mbuf sz=2 --- 0 sec
n=4294967295, mbuf sz=2 --- 0 sec
n=2147483647, mbuf sz=2 --- 0 sec
n=1073741823, mbuf sz=2 --- 0 sec
n=536870911, mbuf sz=4 --- 0 sec
n=268435455, mbuf sz=6 --- 0 sec
n=134217727, mbuf sz=10 --- 0 sec
n=67108863, mbuf sz=10 --- 0 sec
n=33554431, mbuf sz=10 --- 0 sec
n=16777215, mbuf sz=10 --- 0 sec
n=8388607, mbuf sz=10 --- 0 sec
n=4194303, mbuf sz=10 --- 0 sec
n=2097151, mbuf sz=10 --- 0 sec
n=1048575, mbuf sz=10 --- 0 sec
n=524287, mbuf sz=10 --- 0 sec
n=262143, mbuf sz=10 --- 0 sec
n=131071, mbuf sz=10 --- 0 sec
n=65535, mbuf sz=10 --- 0 sec
n=32767, mbuf sz=10 --- 0 sec
n=16383, mbuf sz=10 --- 0 sec
n=8191, mbuf sz=10 --- 0 sec
n=4095, mbuf sz=10 --- 0 sec
n=2047, mbuf sz=10 --- 0 sec
n=1023, mbuf sz=10 --- 0 sec
n=511, mbuf sz=10 --- 0 sec
n=255, mbuf sz=10 --- 0 sec
n=127, mbuf sz=10 --- 0 sec
n=63, mbuf sz=10 --- 0 sec
n=31, mbuf sz=10 --- 0 sec
n=15, mbuf sz=10 --- 0 sec
n=7, mbuf sz=10 --- 0 sec
n=3, mbuf sz=10 --- 0 sec
n=1, mbuf sz=10 --- 0 sec
 0000000 0000046 3718218 5095045 8438771 4531277 5985050 3531670 5893980 9439277. x 10^5171655888
n=17179869183, mbuf sz=2 --- 0 sec
n=8589934591, mbuf sz=2 --- 0 sec
n=4294967295, mbuf sz=2 --- 0 sec
n=2147483647, mbuf sz=2 --- 0 sec
n=1073741823, mbuf sz=2 --- 0 sec
n=536870911, mbuf sz=4 --- 0 sec
n=268435455, mbuf sz=6 --- 0 sec
n=134217727, mbuf sz=10 --- 0 sec
n=67108863, mbuf sz=10 --- 0 sec
n=33554431, mbuf sz=10 --- 0 sec
n=16777215, mbuf sz=10 --- 0 sec
n=8388607, mbuf sz=10 --- 0 sec
n=4194303, mbuf sz=10 --- 0 sec
n=2097151, mbuf sz=10 --- 0 sec
n=1048575, mbuf sz=10 --- 0 sec
n=524287, mbuf sz=10 --- 0 sec
n=262143, mbuf sz=10 --- 0 sec
n=131071, mbuf sz=10 --- 0 sec
n=65535, mbuf sz=10 --- 0 sec
n=32767, mbuf sz=10 --- 0 sec
n=16383, mbuf sz=10 --- 0 sec
n=8191, mbuf sz=10 --- 0 sec
n=4095, mbuf sz=10 --- 0 sec
n=2047, mbuf sz=10 --- 0 sec
n=1023, mbuf sz=10 --- 0 sec
n=511, mbuf sz=10 --- 0 sec
n=255, mbuf sz=10 --- 0 sec
n=127, mbuf sz=10 --- 0 sec
n=63, mbuf sz=10 --- 0 sec
n=31, mbuf sz=10 --- 0 sec
n=15, mbuf sz=10 --- 0 sec
n=7, mbuf sz=10 --- 0 sec
n=3, mbuf sz=10 --- 0 sec
n=1, mbuf sz=10 --- 0 sec
 0000000 0000046 3718218 5095045 8438771 4531277 5985050 3531670 5894606 3902528. x 10^5171655888



2^17179869183 is between
 0000000 0000046 3718218 5095045 8438771 4531277 5985050 3531670 5893980 9439277. x 10^5171655888
 0000000 0000046 3718218 5095045 8438771 4531277 5985050 3531670 5894606 3902528. x 10^5171655888
https://ideone.com/R3szJK
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!