阶乘和圆周率

无心人

发现50#用的程序是每10^7计数一次 所以以前说的每4秒 10^8错误 是40秒

无心人

当前界限下的阶乘77773146302!, 对应的对数0.497149872686531 开始数字3141592653534 检验： (15:09) gp > n = 77773146302 %26 = 77773146302 (07:55) gp > (n+1/2)*log(n) + 1/2*log(2*pi) - n + 1/(12*n) - 1/(360*n^3) %27 = 1872548868987.616634238479898278998546965437108456887088346 (07:55) gp > %27 / log(10) %28 = 813237640895.4971498726865435859265735865657796251777623316 (07:55) gp > %28 - floor(%28) %29 = 0.497149872686543585926573586565779625177762331512 (07:56) gp > 10^%29 %30 = 3.14159265353488687304651737046701462136950721154

gxqcn

原帖由 无心人 于 2008-12-17 14:08 发表 另外一个问题，似乎目前的库，对long double都不支持 我修改程序用long double，取不到对数值 而目前的double精度恐怕支持不了多久了

真正的 long double 有 80 bits，至少需要 10 字节存储， 所以一般编译器不支持。 但可以通过内嵌汇编得到80 bits精度的浮点数， 可惜FPU设计成栈模式，手工编写调试汇编代码非常累人（我之前曾写过一些）！ 以下是 FPU 中一些可能对本主题讨论话题有帮助的指令： FLDL2E Load the log base 2 of e (Napierian constant) FLDL2T Load the log base 2 of Ten FLDLG2 Load the log base 10 of 2 (common log of 2) FLDLN2 Load the log base e of 2 (natural log of 2) FLDPI Load the value of PI FSQRT Square root of ST(0) FYL2X Y*Log2X FYL2XP1 Y*Log2(X+1)

无心人

写汇编代码在 gcc里并不方便 在windows下编译GMP和MPFR也不方便 呵呵 两难

无心人

下一个选择就是5000亿左右了， 估计是极限了 仅仅在两个界限的对数值最后四位有差别 [2780, 4162] 每次加1000个，再调整一次 估计差别值在范围上很小 有错过的可能 ============================= 刚才加大了输出的精度 发现是6位的差别 我再去查下double的具体数值

无心人

今天有时间测试了下程序的时间

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

复制代码

=============================================== E:\math>factpi1 请输入起始搜索数字：100000000 请输入上界: 3.1415927 请输入下界: 3.1415926 起始位置：100000000 当前界限[0.49714986528586863, 0.49714987910989134] 高精度stirling函数执行时间(毫秒): 40.058000, 测试值：493.74586181 起始对数: 0.20876475177585671 每10^6次查找执行时间(秒): 0.460662 当前界限下的阶乘149832529!, 对应的对数0.497149873593380 ================================================= 对应机器P4 2.0/1G DDR400 软件XP 编译器MinGW + GCC 4.3.2 使用GMP 4.2.4 + MPFR 2.3.2， 编译为静态包 =================================================== 可以看到MPFR下，每1000次stirling函数计算仅0.04秒 而double下每1000*1000次叠加对数值是需要0.46秒 可计算出高精度浮点的运算时间仅是double叠加对数的100倍 所以在目前情况下，普通函数计算称为瓶颈 另外，根据执行时间判断 每次计算需要800-1000时钟周期

无心人

使用-ffast-math选项编译后的输出 请输入起始搜索数字：100000000 请输入上界: 3.1415927 请输入下界: 3.1415926 起始位置：100000000 当前界限[0.49714986528586863, 0.49714987910989134] 高精度stirling函数执行时间(毫秒): 30.043000, 测试值：493.74586181 起始对数: 0.20876475177585671 每10^6次查找执行时间(秒): 0.150216 当前界限下的阶乘143984543!, 对应的对数0.497149871271711

无心人

57#同样条件下 但修改程序如下

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;
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. hi = log10(hi);
82. low = log10(low);
84. printf("起始位置：%.0f\n当前界限[%.17f, %.17f]\n", start, low, hi);
85. a = start;
86. b = 1.0;
87. while (a >= 10.0)
88. {
89. a /= 10.0;
90. b *= 10.0; ////////////////////////////////////修改
91. }
93. lf = 2.0;
94. gettimeofday(&start_time, NULL);
95. for (i = 0; i < 1000; i ++)
96. lf += stirling(lf + 100000000.0);
97. gettimeofday(&end_time, NULL);
98. used_time = ((end_time.tv_sec - start_time.tv_sec)*1000000.0 + (end_time.tv_usec - start_time.tv_usec)) / 1000.0;
99. printf("高精度stirling函数执行时间(毫秒): %8.6lf, 测试值：%.8lf\n", used_time, lf);
101. lf = stirling(start);
102. printf("起始对数: %.17f\n", lf);
104. gettimeofday(&start_time, NULL);
106. while (1)
107. {
108. start += 1.0;
109. a = log10(start/b); //////////////////////////////修改
111. if (a >= 1.0)
112. {
113. b *= 10.0; //////////////////////////////修改
114. a -= 1.0;
115. }
117. i ++;
118. lf += a;
120. if (lf >= 1.0)
121. lf -= 1.0;
123. // printf("[%.0f, %.15f] ", start, lf);
125. if ((lf >= low) && (lf < hi))
126. {
127. find = 1;
128. break;
129. }
131. if (i >= 1000)
132. {
133. // lf = stirling(start);
134. i = 0;
135. j ++;
136. if (j >= 1000)
137. {
138. if (test1 == 0)
139. {
140. gettimeofday(&end_time, NULL);
141. used_time = ((end_time.tv_sec - start_time.tv_sec)*1000000.0 + (end_time.tv_usec - start_time.tv_usec))/1000000.0;
142. printf("每10^6次查找执行时间(秒): %.6lf\n", used_time);
143. test1 = 1;
144. }
145. }
146. }
147. }
149. if (find)
150. {
151. printf("\n当前界限下的阶乘%.0f!, 对应的对数%.15f", start, lf);
152. }
154. mpfr_clear(ln10);
155. mpfr_clear(cA);
156. mpfr_clear(cT);
157. return 1;
158. }

复制代码

请输入起始搜索数字：100000000 请输入上界: 3.1415927 请输入下界: 3.1415926 起始位置：100000000 当前界限[0.49714986528586863, 0.49714987910989134] 高精度stirling函数执行时间(毫秒): 30.043000, 测试值：493.74586181 起始对数: 0.20876475177585671 每10^6次查找执行时间(秒): 0.230332 当前界限下的阶乘149832529!, 对应的对数0.497149873593496

无心人

可以看到 乘法和除法在浮点运算下 对时间的影响是很大的 还有一个可以修改优化的地方是对取10的对数 改为取2的对数再乘以log10(2) 但不知道标准函数有取2的对数的函数么？ 还是需要自己写汇编

无心人

去掉-ffast-math，同时增加了对求10的对数的优化

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. }

复制代码

========================================= 请输入起始搜索数字：100000000 请输入上界: 3.1415927 请输入下界: 3.1415926 起始位置：100000000 当前界限[0.49714986528586863, 0.49714987910989134] 高精度stirling函数执行时间(毫秒): 30.044000, 测试值：493.74586181 起始对数: 0.20876475177585671 每10^6次查找执行时间(秒): 0.240346 当前界限下的阶乘149832529!, 对应的对数0.497149877671004