阶乘和圆周率

无心人

可以发现上面几个例子的结果是两个不同的输出 测试发现正确的输出是 149832529!!!! 但计算器计算发现，似乎此时的对数值精度并不高

无心人

重新修改成

2. #include <gmp.h>
3. #include <mpfr.h>
4. #include <stdio.h>
5. #include <math.h>
6. #include <sys/time.h>
8. #define DefaultBitLength 128
10. mpfr_t ln10, cA, cT;
11. //=1/2ln(2π)+(n+1/2)ln(n)-n+1/12n
12. //n > 10^8 误差<1/360Z^3 = 10^(-26)
13. long double stirling(double n)
14. {
15. mpfr_t a, b, c, s;
16. long double r;
18. mpfr_init2(a, DefaultBitLength);
19. mpfr_init2(b, DefaultBitLength);
20. mpfr_set_d(a, n, GMP_RNDD); //a=n
21. mpfr_init2(c, DefaultBitLength);
22. mpfr_init2(s, DefaultBitLength);
24. mpfr_div(s, cT, a, GMP_RNDD); //s = 1/12a = 1/12n
25. mpfr_sub(s, s, a, GMP_RNDD); //s = s - a = - n + 1/12n
27. mpfr_log(c, a, GMP_RNDD); //c = ln(n)
29. mpfr_set_d(b, 0.5, GMP_RNDD);
30. mpfr_add(b, a, b, GMP_RNDD);
31. mpfr_mul(b, b, c, GMP_RNDD); //b = (n + 1/2)ln(n)
33. mpfr_add(s, s, b, GMP_RNDD); //s = (n + 1/2)ln(n) - n + 1/12n
35. mpfr_add(s, s, cA, GMP_RNDD); //最终结果s=1/2ln(2*pi) + (n + 1/2)ln(n) - n + 1/12n
37. mpfr_mul(s, s, ln10, GMP_RNDD); //转换为log10
38. mpfr_floor(b, s);
39. mpfr_sub(s, s, b, GMP_RNDD); //只要小数
41. r = mpfr_get_ld(s, GMP_RNDD);
43. mpfr_clear(a);
44. mpfr_clear(b);
45. mpfr_clear(c);
46. mpfr_clear(s);
47. return r;
48. }
50. int main(void)
51. {
52. double a, b, low, hi, lf, t;
53. double start, used_time, LT2;
54. int find = 0, test1 = 0;
55. long i = 0, j = 0;
56. struct timeval start_time, end_time;
58. mpfr_init2(ln10, DefaultBitLength);
59. mpfr_init2(cA, DefaultBitLength);
60. mpfr_init2(cT, DefaultBitLength);
61. mpfr_set_ui(ln10, 10, GMP_RNDD);
62. mpfr_set_ui(cT, 1, GMP_RNDD);
64. mpfr_log(ln10, ln10, GMP_RNDD);
65. mpfr_div(ln10, cT, ln10, GMP_RNDD);
67. mpfr_const_pi(cA, GMP_RNDD);
68. mpfr_mul_ui(cA, cA, 2, GMP_RNDD);
69. mpfr_log(cA, cA, GMP_RNDD);
70. mpfr_div_ui(cA, cA, 2, GMP_RNDD);
72. mpfr_div_ui(cT, cT, 12, GMP_RNDD);
74. printf("请输入起始搜索数字：");
75. scanf("%lf", &start);
76. printf("请输入上界: ");
77. scanf("%lf", &hi);
78. printf("请输入下界: ");
79. scanf("%lf", &low);
81. LT2 = log10(2);
82. hi = log10(hi);
83. low = log10(low);
85. printf("起始位置：%.0f\n当前界限[%.17f, %.17f]\n", start, low, hi);
86. a = start;
87. b = 1.0;
88. while (a >= 10.0)
89. {
90. a /= 10.0;
91. b *= 10.0; ///////////////////////////////////
92. }
94. lf = 2.0;
95. gettimeofday(&start_time, NULL);
96. for (i = 0; i < 1000; i ++)
97. lf += stirling(lf + 100000000.0);
98. gettimeofday(&end_time, NULL);
99. used_time = ((end_time.tv_sec - start_time.tv_sec)*1000000.0 + (end_time.tv_usec - start_time.tv_usec)) / 1000.0;
100. printf("高精度stirling函数执行时间(毫秒): %8.6lf, 测试值：%.8lf\n", used_time, lf);
102. lf = stirling(start);
103. printf("起始对数: %.17f\n", lf);
105. gettimeofday(&start_time, NULL);
107. while (1)
108. {
109. start += 1.0;
110. a = log2(start/b) * LT2; ///////////////////////////
112. if (a >= 1.0)
113. {
114. b *= 10.0; ////////////////////////////////
115. a -= 1.0;
116. }
118. i ++;
119. lf += a;
121. if (lf >= 1.0)
122. lf -= 1.0;
124. // printf("[%.0f, %.15f] ", start, lf);
126. if ((lf >= low) && (lf < hi))
127. {
128. find = 1;
129. break;
130. }
132. if (i >= 1000)
133. {
134. // lf = stirling(start);
135. i = 0;
136. j ++;
137. if (j >= 1000)
138. {
139. if (test1 == 0)
140. {
141. gettimeofday(&end_time, NULL);
142. used_time = ((end_time.tv_sec - start_time.tv_sec)*1000000.0 + (end_time.tv_usec - start_time.tv_usec))/1000000.0;
143. printf("每10^6次查找执行时间(秒): %.6lf\n", used_time);
144. test1 = 1;
145. }
146. }
147. }
148. }
150. if (find)
151. {
152. printf("\n当前界限下的阶乘%.0f!, 对应的对数%.15f", start, lf);
153. }
155. mpfr_clear(ln10);
156. mpfr_clear(cA);
157. mpfr_clear(cT);
158. return 1;
159. }

复制代码

编译命令 gcc factPi1.c -o factpi1.exe -lmpfr -lgmp -lm ====================================================== 请输入起始搜索数字：100000000 请输入上界: 3.1415927 请输入下界: 3.1415926 起始位置：100000000 当前界限[0.49714986528586863, 0.49714987910989134] 高精度stirling函数执行时间(毫秒): 30.043000, 测试值：493.74586181 起始对数: 0.20876475177585671 每10^6次查找执行时间(秒): 0.300432 当前界限下的阶乘149832529!, 对应的对数0.497149873638380 ============================================ 计算器得到的值 0.4971498735933375434005 ============================================ 直接使用10的对数 其他条件不变的输出 即使用 a = log10(start/b);// * LT2;取代a = log2(start / b) * LT2; ============================================ 请输入起始搜索数字：100000000 请输入上界: 3.1415927 请输入下界: 3.1415926 起始位置：100000000 当前界限[0.49714986528586863, 0.49714987910989134] 高精度stirling函数执行时间(毫秒): 30.044000, 测试值：493.74586181 起始对数: 0.20876475177585671 每10^6次查找执行时间(秒): 0.450648 当前界限下的阶乘149832529!, 对应的对数0.497149873595150 =========================================== 精度大些，但消耗时间多了 即使使用log10(2)的高精度值直接乘而不是用对数函数 也不能增大精度，所以似乎库函数存在某些优化？？？？

无心人

程序可以确定为62#程序但 使用a = log10(start/b);// * LT2;取代a = log2(start / b) * LT2; 此时 不同的优化结果如下 ========================================= 不使用任何优化 每10^6次查找执行时间(秒): 0.470677 当前界限下的阶乘149832529!, 对应的对数0.4971 4987 3595 150 ========================================= 带-O2 每10^6次查找执行时间(秒): 0.450648 当前界限下的阶乘149832529!, 对应的对数0.4971 4987 3593 083 =========================================== 带-O2 -ffast-math 每10^6次查找执行时间(秒): 0.160230 当前界限下的阶乘149832529!, 对应的对数0.4971 4987 3593 496 =========================================== 对比计算器得到的数值0.4971 4987 3593 3375 4340 05 ========================================== 反而是最后一种情况又快， 精度又高 呵呵

无心人

修改了代码 使用gcc -O2 -ffast-math -march=pentium4 factpi1.c -o factpi1.exe -lmpfr -lgmp

2. #include <gmp.h>
3. #include <mpfr.h>
4. #include <stdio.h>
5. #include <stdlib.h>
6. #include <math.h>
7. #include <sys/time.h>
9. #define DefaultBitLength 128
10. #define TEST
12. mpfr_t ln10, cA, cT;
14. FILE * fp;
16. //=1/2ln(2π)+(n+1/2)ln(n)-n+1/12n
17. //n > 10^8 误差<1/360Z^3 = 10^(-26)
18. long double stirling(double n)
19. {
20. mpfr_t a, b, c, s;
21. long double r;
23. mpfr_init2(a, DefaultBitLength);
24. mpfr_init2(b, DefaultBitLength);
25. mpfr_set_d(a, n, GMP_RNDD); //a=n
26. mpfr_init2(c, DefaultBitLength);
27. mpfr_init2(s, DefaultBitLength);
29. mpfr_div(s, cT, a, GMP_RNDD); //s = 1/12a = 1/12n
30. mpfr_sub(s, s, a, GMP_RNDD); //s = s - a = - n + 1/12n
32. mpfr_log(c, a, GMP_RNDD); //c = ln(n)
34. mpfr_set_d(b, 0.5, GMP_RNDD);
35. mpfr_add(b, a, b, GMP_RNDD);
36. mpfr_mul(b, b, c, GMP_RNDD); //b = (n + 1/2)ln(n)
38. mpfr_add(s, s, b, GMP_RNDD); //s = (n + 1/2)ln(n) - n + 1/12n
40. mpfr_add(s, s, cA, GMP_RNDD); //最终结果s=1/2ln(2*pi) + (n + 1/2)ln(n) - n + 1/12n
42. mpfr_mul(s, s, ln10, GMP_RNDD); //转换为log10
43. mpfr_floor(b, s);
44. mpfr_sub(s, s, b, GMP_RNDD); //只要小数
46. r = mpfr_get_ld(s, GMP_RNDD);
48. mpfr_clear(a);
49. mpfr_clear(b);
50. mpfr_clear(c);
51. mpfr_clear(s);
52. return r;
53. }
55. int main(void)
56. {
57. double a, b, low, hi, lf, t;
58. double start, used_time;
59. int find = 0, test1 = 0;
60. long i, j = 0;
61. struct timeval start_time, end_time;
63. mpfr_init2(ln10, DefaultBitLength);
64. mpfr_init2(cA, DefaultBitLength);
65. mpfr_init2(cT, DefaultBitLength);
66. mpfr_set_ui(ln10, 10, GMP_RNDD);
67. mpfr_set_ui(cT, 1, GMP_RNDD);
69. mpfr_log(ln10, ln10, GMP_RNDD);
70. mpfr_div(ln10, cT, ln10, GMP_RNDD);
72. mpfr_const_pi(cA, GMP_RNDD);
73. mpfr_mul_ui(cA, cA, 2, GMP_RNDD);
74. mpfr_log(cA, cA, GMP_RNDD);
75. mpfr_div_ui(cA, cA, 2, GMP_RNDD);
77. mpfr_div_ui(cT, cT, 12, GMP_RNDD);
79. printf("请输入起始搜索数字：");
80. scanf("%lf", &start);
81. printf("请输入上界: ");
82. scanf("%lf", &hi);
83. printf("请输入下界: ");
84. scanf("%lf", &low);
85. fp = fopen("factpi.log", "a+t");
86. printf("起始位置：%.0f\n当前界限[%.17f, %.17f]\n", start, low, hi);
87. fprintf(fp, "\n///////////////////////////////////////////////////////////////\n");
88. fprintf(fp, "起始位置：%.0f\n当前界限[%.17f, %.17f]\n", start, low, hi);
90. hi = log10(hi);
91. low = log10(low);
93. printf("当前对数界限[%.17f, %.17f]\n", low, hi);
94. fprintf(fp, "当前对数界限[%.17f, %.17f]\n", low, hi);
95. a = start;
96. b = 1.0;
97. while (a >= 10.0)
98. {
99. a /= 10.0;
100. b *= 10.0;
101. }
103. lf = 2.0;
104. gettimeofday(&start_time, NULL);
105. for (i = 0; i < 1000; i ++)
106. lf += stirling(lf + 100000000.0);
107. gettimeofday(&end_time, NULL);
108. used_time = ((end_time.tv_sec - start_time.tv_sec)*1000000.0 + (end_time.tv_usec - start_time.tv_usec)) / 1000.0;
109. printf("高精度stirling函数执行时间(毫秒): %8.6lf, 测试值：%.8lf\n", used_time, lf);
110. fprintf(fp, "高精度stirling函数执行时间(毫秒): %8.6lf, 测试值：%.8lf\n", used_time, lf);
112. lf = stirling(start);
113. printf("起始对数: %.17f\n", lf);
114. fprintf(fp, "起始对数: %.17f\n", lf);
115. fflush(fp);
116. gettimeofday(&start_time, NULL);
118. i = 0;
119. while (1)
120. {
121. if ((a = log10((start += 1.0)/b)) >= 1.0)
122. {
123. b *= 10.0;
124. a -= 1.0;
125. }
127. if ((lf += a) >= 1.0)
128. lf -= 1.0;
130. // printf("[%.0f, %.15f] ", start, lf);
132. if ((lf >= low) && (lf < hi))
133. {
134. find = 1;
135. break;
136. }
138. if (( ++i) >= 1000)
139. {
140. lf = stirling(start);
141. i = 0;
142. j ++;
143. #ifdef TEST
144. if (j >= 1000)
145. {
146. if (test1 == 0)
147. {
148. gettimeofday(&end_time, NULL);
149. used_time = ((end_time.tv_sec - start_time.tv_sec)*1000000.0 + (end_time.tv_usec - start_time.tv_usec))/1000000.0;
150. printf("每10^6次查找执行时间(秒): %.6lf\n", used_time);
151. test1 = 1;
152. }
153. }
154. #else
155. if (j >= 100000000)
156. {
157. j = 0;
158. fprintf(fp, "当前搜索数值: %.0f\n", start);
159. fflush(fp);
160. }
161. #endif
162. }
163. }
165. if (find)
166. {
167. printf("\n当前界限下的阶乘%.0f!, 对应的对数%.15f", start, lf);
168. fprintf(fp, "\n当前界限下的阶乘%.0f!, 对应的对数%.15f", start, lf);
169. }
171. mpfr_clear(ln10);
172. mpfr_clear(cA);
173. mpfr_clear(cT);
174. fclose(fp);
175. return 1;
176. }

复制代码

========================================================== 请输入起始搜索数字：100000000 请输入上界: 3.1415927 请输入下界: 3.1415926 起始位置：100000000 当前界限[3.14159260000000010, 3.14159269999999990] 当前对数界限[0.49714986528586863, 0.49714987910989134] 高精度stirling函数执行时间(毫秒): 30.044000, 测试值：493.74586181 起始对数: 0.20876475177585671 每10^6次查找执行时间(秒): 0.170244 当前界限下的阶乘149832529!, 对应的对数0.497149873593338 ======================================================== 经过查询汇编代码，程序已经被gcc优化到了很好的程度 相对以前，执行时间已经缩小了很多了 使用了P6以上的fcomi和直接的log指令， 另外，其最终的对数值精度提高了 ========================================================== P4 2.0 /1G DDR/XP SP3

无心人

d = 314159265358, n = 1154318938997!, 对应的对数0.497149872693353，首数字3141592653584138 ============================================================= 通过PARI验证 (09:52) gp > n = 1154318938997 %38 = 1154318938997 (07:44) gp > (n+1/2)*log(n) + 1/2*log(2*pi) - n + 1/(12*n) - 1/(360*n^3) %39 = 30906348934959.82738115133341720201016520726795684547468467 (07:44) gp > %39 / log(10) %40 = 13422456798229.49714987269335219103203453332675692373542468 (07:44) gp > %40 - floor(%40) %41 = 0.497149872693352191032034533326756923735424682255 (07:44) gp > 10^(%41) %42 = 3.14159265358413885452821840174232968895743730673

无心人

用来计算的机器是 Core 2 XEON 1.6G / 1G DDR2 下一个将在10万亿附近产生 目前增加到10000个数字校正一次 每10 ^8的时间是7.8125秒 如果还是每1000校正一次 则是9.375秒 ========================== 差距是两天 呵呵 ========================== 依据是对最新求得的结果1154318938997 使用10000间隔校正 从之前的11543000000重新搜索 结果对数是0.497149872693352 而GP得到的结果是0.497149872693353 =========================== 目前的对数范围 [0.49714987269402422 , 0.49714987269416244] 差距是13822 但因为所有代码均被gcc优化成汇编代码,而gcc默认是63位浮点 所以实际上的差距是1382200 我想足够10000间隔校正而不丢失分辨精度吧 ============================= 呵呵,求出下个结果,就不再计算了 再下个至少需要100天

无心人

在这十天的空隙里 先贴点以10^n开始的 d = 1, n = 5 d = 10, n = 27 d = 100, n = 197 d = 1000, n = 197 d = 10000, n = 7399 //////////////////////////////////////////// 呵呵,下面的似乎是错误的 看来需要校正结果

无心人

竟然昨天下午5:37就得到了数据 ============================= 请输入起始搜索数字：1154318938997 请输入上界: 3.141592653590 请输入下界: 3.141592653589 起始位置：1154318938997 当前界限[3.14159265358900000, 3.14159265359000010] 当前对数界限[0.49714987269402422, 0.49714987269416244] 高精度stirling函数执行时间(毫秒): 15.625000, 测试值：493.74586181 起始对数: 0.49714987269335220 每10^6次查找执行时间(秒): 0.078125 当前界限下的阶乘1544607046599!, 对应的对数0.497149872694156 ================================================= GP的验证 (11:04) gp > stirling1(n) = (n+1/2)*log(n) + 1/2*log(2*pi) - n + 1/(12*n) - 1/( 360*n^3) (11:04) gp > stirling(n) = stirling1(n) / log(10) - floor(stirling1(n) / log(10 )) (11:06) gp > stirling(1544607046599) %57 = 0.497149872694155502964270816144690329633515005612 (07:54) gp > 10^%57 %58 = 3.14159265358994983986193148792491264569132697520

gxqcn

看兄弟比较孤寂的，回复一个。 可惜没有这么好的软硬件条件，所以只能欣赏而已。

无心人

呵呵 这个情况 导致了我还有继续计算下去的必要 再算一个，直到结果超越10万亿