求定积分有什么好的算法吗？

pupilzeng

不过你提供的函数还是没法满足我所要的精度。
这时由c语言原始数据类型精度限制的。
我使用了gmp库，可以达到很高的精度了，起码小数点后几十位是没问题的。而且很容易扩充成更高精度的。

1. 
   
2. #include<stdio.h>
   
3. #include<gmp.h>
   
4. #include<stdlib.h>
   
5. 
   
6. #define PREC 328
   
7. 
   
8. void normal_cdf(mpf_t, const mpf_t, unsigned long int);
   
9. 
   
10. int
    
11. main(int argc, char * argv[])
    
12. {
    
13.         double var;
    
14.         unsigned long int eff;
    
15.         var = atof(argv[1]);
    
16.         eff = atof(argv[2]);
    
17.         mpf_t x, result;
    
18.         mpf_init(result);
    
19.         mpf_init_set_d(x, var);
    
20.         normal_cdf(result, x, eff);
    
21.         char format[50];
    
22.         sprintf(format, "norma_cdf(%%f) = %%.%dFe\n", eff);
    
23.         gmp_printf(format, var, result);
    
24.         mpf_clear(result);       
    
25.         mpf_clear(x);
    
26.         return 0;
    
27. }
    
28. void
    
29. normal_cdf(mpf_t result,const mpf_t x, unsigned long int eff)
    
30. {       
    
31.         mpz_t n;
    
32.         mpf_t  a, error;
    
33.         mpq_t q;
    
34.         mpz_init(n);
    
35.         mpf_init(a);
    
36.         mpq_init(q);
    
37.         mpf_init(error);
    
38.         mpf_set_prec(a, PREC);
    
39.         mpf_set_prec(result, PREC);
    
40.         mpf_set_prec(error, PREC);
    
41.        
    
42.         mpf_set_str(a, "0.3989422804014326779399460599343818684758586311649346576659258296706579258993018385012523339073069364", 10);
    
43. //                printf("error\n");
    
44. //        gmp_printf("a = %.100Ff\n", a);
    
45. //        gmp_printf("error = %40Ff\n", error);
    
46. //        gmp_printf("x =%20Ff \n", x);
    
47.         mpf_set_d(result, 0.5);
    
48. //        gmp_printf("result =%20Ff \n", result);
    
49.         mpf_mul(a, a, x);
    
50. //        gmp_printf("a = %.100Ff\n", a);
    
51.         mpf_t dist, temp;
    
52.         mpf_init(dist);
    
53.         mpf_init(temp);
    
54.         mpf_set_prec(dist, PREC);
    
55.         mpf_set_prec(temp, PREC);
    
56.         mpz_t num, den;
    
57.         mpz_init(num);
    
58.         mpz_init(den);
    
59.         mpf_set_d(error, 10);
    
60. //        gmp_printf("error =%20Ff \n", error);
    
61.         mpf_pow_ui(error, error, eff + 1);
    
62. //        gmp_printf("error =%20Ff \n", error);
    
63.         mpf_ui_div(error, 1, error);
    
64. //        gmp_printf("error =%20Ff \n", error);
    
65.         do
    
66.         {
    
67. //                gmp_printf("\nthe %Zd-th loop \n", n);
    
68.                 mpf_add(result, result, a);
    
69. //                gmp_printf("result =%20Ff \n", result);
    
70.                 mpz_add_ui(n, n, 1);
    
71.                 mpz_mul_ui(num, n, 2);
    
72.                 mpz_sub_ui(num, num, 1);
    
73. //                gmp_printf("num =%Zd \n", num);
    
74.                 mpz_add_ui(den, num, 2);
    
75.                 mpz_mul(den, den, n);
    
76.                 mpz_mul_ui(den, den, 2);
    
77. //                gmp_printf("den =%Zd \n", den);
    
78.                 mpq_set_num(q, num);
    
79.                 mpq_set_den(q, den);
    
80.                 mpf_set_q(temp, q);
    
81. //                gmp_printf("temp = %20Ff\n", temp);
    
82.                 mpf_mul(temp, temp, x);
    
83.                 mpf_mul(temp, temp, x);
    
84.                 mpf_neg(temp, temp);
    
85.                 mpf_mul(a, a, temp);
    
86. //                gmp_printf("a = %.20Ff\n", a);
    
87.                 mpf_set(dist, a);
    
88. //                gmp_printf("dist =%20Ff \n", dist);
    
89.                 mpf_div(dist, dist, result);
    
90.                 mpf_abs(dist, dist);
    
91. //                gmp_printf("dist =%.70Ff \n", dist);
    
92.         }while(mpf_cmp(dist, error) > 0);
    
93. //        gmp_printf("phi_cdf(%f) = %.4Fe \n", var, result);
    
94.         mpf_clear(a);
    
95.         mpq_clear(q);
    
96.         mpf_clear(dist);
    
97.         mpf_clear(temp);
    
98.         mpf_clear(error);
    
99.         mpz_clear(n);
    
100.         mpz_clear(num);
     
101.         mpz_clear(den);
     
102. }
     

复制代码

编译时需要libgmp3-dev库，使用gcc -lgmp -o normal_cdf normal_cdf.c编译。
运行方法：normal_cdf var prec
第一个参数：要求概率函数的自变量。
第二个参数：结果的精确度，即小数点后面的位数。

这个级数的收敛速度不错，是你自己推导出来的吗？
佩服ing

无双

不错

lordbyorn

Compare to your code, my code looks like pseudocode.

这个级数的收敛速度不错，是你自己推导出来的吗？

Not very hard. :p

pupilzeng

最初由 Sandy 发表
f:=dy/dx不够精确，I like some cool methods

其实要提高精度单从算法上来说并不难，主要是受到c语言本身的数据精度影响。
所以要特别高的精度的话，还是要利用任意精度的数学库。

soloforce

浮点数运算密集的话可以考虑用 mpfr 数学库，基于gmp的，不过精度控制和舍入误差比较合理。