c语言精确度问题

h-b-z-d-j-s-m

1. #include <stdio.h>
   
2. #include <math.h>
   
3. double  fun ( double  eps)
   
4. {
   
5.         double pi,a=1.0,b=1.0,i=1.0,t;
   
6.         pi=1.0;
   
7.         t=1.0;
   
8.         while(t>=eps)
   
9.         {
   
10.                 a*=i;
    
11.                 b*=(2*i+1);
    
12.                 t=a/b;
    
13.                 pi+=t;
    
14.                 i++;
    
15.         }
    
16.         return (2*pi);
    
17. 
    
18. 
    
19. }
    
20. 
    
21. void main( )
    
22. { double  x;void NONO ();
    
23.   printf("Input eps:") ;
    
24.   scanf("%lf",&x); printf("\neps = %lf, PI=%lf\n", x, fun(x));
    
25.   NONO();
    
26. }
    
27. 
    
28. void NONO ()
    
29. {/* 本函数用于打开文件，输入数据，调用函数，输出数据，关闭文件。 */
    
30.   FILE *fp, *wf ;
    
31.   int i ;
    
32.   double x ;
    
33. 
    
34.   fp = fopen("in.dat","r") ;
    
35.   wf = fopen("out.dat","w") ;
    
36.   for(i = 0 ; i < 10 ; i++) {
    
37.     fscanf(fp, "%lf", &x) ;
    
38.     fprintf(wf, "%lf\n", fun(x)) ;
    
39.   }
    
40.   fclose(fp) ;
    
41.   fclose(wf) ;
    
42. }

复制代码

fun 函数里面哪里不对啊跟结果不大一样
[img]%5Bimg%5Dhttps://s3.bmp.ovh/imgs/2022/05/18/c4a7481e7128766d.jpg[/img][/img]

人造人

1. #include <stdio.h>
   
2. 
   
3. long double get_pi(long double eps) {
   
4.     long double pi = 1;
   
5.     for(size_t n = 1; ; ++n) {
   
6.         long double a = 1;
   
7.         for(size_t i = 1; i <= n; ++i) {
   
8.             a *= i;
   
9.         }
   
10.         long double b = 1;
    
11.         for(size_t i = 3; i <= 2 * n + 1; i += 2) {
    
12.             b *= i;
    
13.         }
    
14.         long double temp = a / b;
    
15.         if(temp < eps) break;
    
16.         pi += temp;
    
17.     }
    
18.     return 2 * pi;
    
19. }
    
20. 
    
21. int main(void) {
    
22.     printf("%.30Lf\n", get_pi(1e-4950L));
    
23.     return 0;
    
24. }
    

复制代码

人造人

1. pi = 3.14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 58209 74944 59230 78164 06286 20899 86280 34825 34211 70679

复制代码

1. $ cat main.c
   
2. #include <stdio.h>
   
3. #include <mpfr.h>
   
4. 
   
5. long double get_pi(long double eps) {
   
6.     mpfr_t pi, a, b, temp;
   
7.     mpfr_inits2(256, pi, a, b, temp, (mpfr_ptr)0);
   
8.     mpfr_set_ld(pi, 1, MPFR_RNDZ);
   
9.     for(size_t n = 1; ; ++n) {
   
10.         mpfr_set_ld(a, 1, MPFR_RNDZ);
    
11.         for(size_t i = 1; i <= n; ++i) {
    
12.             mpfr_mul_ui(a, a, i, MPFR_RNDZ);
    
13.         }
    
14.         mpfr_set_ld(b, 1, MPFR_RNDZ);
    
15.         for(size_t i = 3; i <= 2 * n + 1; i += 2) {
    
16.             mpfr_mul_ui(b, b, i, MPFR_RNDZ);
    
17.         }
    
18.         mpfr_div(temp, a, b, MPFR_RNDZ);
    
19.         if(mpfr_cmp_ld(temp, eps) < 0) break;
    
20.         mpfr_add(pi, pi, temp, MPFR_RNDZ);
    
21.     }
    
22.     long double result = mpfr_get_ld(pi, MPFR_RNDZ);
    
23.     mpfr_clears(pi, a, b, temp, (mpfr_ptr)0);
    
24.     mpfr_free_cache();
    
25.     return 2 * result;
    
26. }
    
27. 
    
28. int main(void) {
    
29.     printf("%.30Lf\n", get_pi(1e-4950L));
    
30.     return 0;
    
31. }
    
32. $ gcc -Wall -O3 -o main main.c $(pkg-config --cflags --libs mpfr)
    
33. $ time ./main
    
34. 3.141592653589793238295968524909
    
35. 
    
36. real        0m8.716s
    
37. user        0m8.698s
    
38. sys        0m0.001s
    
39. $

复制代码

看起来这个算法到 3.141592653589793238 就是极限了

人造人

1. #include <stdio.h>
   
2. #include <math.h>
   
3. 
   
4. double fun(double eps) {
   
5.     double pi, a = 1.0, b = 1.0, i = 1.0, t;
   
6.     pi = 1.0;
   
7.     t = 1.0;
   
8.     //while(t >= eps) {
   
9.     while(1) {
   
10.         a *= i;
    
11.         b *= (2 * i + 1);
    
12.         t = a / b;
    
13.         if(t < eps) break;
    
14.         pi += t;
    
15.         i++;
    
16.     }
    
17.     return (2 * pi);
    
18. }
    
19. 
    
20. // void main( )
    
21. int main(void) {
    
22.     double x;
    
23.     void NONO();
    
24.     //printf("Input eps:");
    
25.     //scanf("%lf", &x);
    
26.     x = 0.0005;
    
27.     printf("\neps = %lf, PI=%lf\n", x, fun(x));
    
28.     //NONO();
    
29. }
    
30. 
    
31. void NONO() { /* 本函数用于打开文件，输入数据，调用函数，输出数据，关闭文件。 */
    
32.     FILE *fp, *wf;
    
33.     int i;
    
34.     double x;
    
35. 
    
36.     fp = fopen("in.dat", "r");
    
37.     wf = fopen("out.dat", "w");
    
38.     for(i = 0; i < 10; i++) {
    
39.         fscanf(fp, "%lf", &x);
    
40.         fprintf(wf, "%lf\n", fun(x));
    
41.     }
    
42.     fclose(fp);
    
43.     fclose(wf);
    
44. }
    

复制代码

人造人

原来是mpfr_t转long double出现了问题
其实根本不需要1e-4950，1e-150就可以精确到小数点后100位

1. $ cat main.c
   
2. #include <stdio.h>
   
3. #include <mpfr.h>
   
4. 
   
5. void get_pi(mpfr_t *result, const mpfr_t eps) {
   
6.     mpfr_t pi, a, b, temp;
   
7.     mpfr_inits2(65536, pi, a, b, temp, (mpfr_ptr)0);
   
8.     mpfr_set_ld(pi, 1, MPFR_RNDZ);
   
9.     mpfr_set_ld(a, 1, MPFR_RNDZ);
   
10.     mpfr_set_ld(b, 1, MPFR_RNDZ);
    
11.     for(size_t n = 1; ; ++n) {
    
12.         mpfr_mul_ui(a, a, n, MPFR_RNDZ);
    
13.         mpfr_mul_ui(b, b, 2 * n + 1, MPFR_RNDZ);
    
14.         mpfr_div(temp, a, b, MPFR_RNDZ);
    
15.         if(mpfr_cmp(temp, eps) < 0) break;
    
16.         mpfr_add(pi, pi, temp, MPFR_RNDZ);
    
17.     }
    
18.     mpfr_mul_ui(pi, pi, 2, MPFR_RNDZ);
    
19.     mpfr_set(*result, pi, MPFR_RNDZ);
    
20.     mpfr_clears(pi, a, b, temp, (mpfr_ptr)0);
    
21. }
    
22. 
    
23. int main(void) {
    
24.     mpfr_t eps, result;
    
25.     mpfr_inits2(65536, eps, result, (mpfr_ptr)0);
    
26.     mpfr_set_str(eps, "1e-150", 10, MPFR_RNDZ);
    
27.     get_pi(&result, eps);
    
28.     mpfr_printf("%.100Rf\n", result);
    
29.     mpfr_clears(eps, result, (mpfr_ptr)0);
    
30.     mpfr_free_cache();
    
31.     return 0;
    
32. }
    
33. $ gcc-debug -o main main.c $(pkg-config --cflags --libs mpfr)
    
34. $ time ./main
    
35. 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170680
    
36. 
    
37. real        0m0.164s
    
38. user        0m0.150s
    
39. sys        0m0.014s
    
40. $

复制代码

人造人

1. $ cat main.c
   
2. #include <stdio.h>
   
3. #include <mpfr.h>
   
4. 
   
5. void get_pi(mpfr_t result, double eps) {
   
6.     mpfr_t pi, a, b, temp;
   
7.     mpfr_inits2(65536, pi, a, b, temp, (mpfr_ptr)0);
   
8.     mpfr_set_ld(pi, 1, MPFR_RNDZ);
   
9.     mpfr_set_ld(a, 1, MPFR_RNDZ);
   
10.     mpfr_set_ld(b, 1, MPFR_RNDZ);
    
11.     for(size_t n = 1; ; ++n) {
    
12.         mpfr_mul_ui(a, a, n, MPFR_RNDZ);
    
13.         mpfr_mul_ui(b, b, 2 * n + 1, MPFR_RNDZ);
    
14.         mpfr_div(temp, a, b, MPFR_RNDZ);
    
15.         if(mpfr_cmp_d(temp, eps) < 0) break;
    
16.         mpfr_add(pi, pi, temp, MPFR_RNDZ);
    
17.     }
    
18.     mpfr_mul_ui(pi, pi, 2, MPFR_RNDZ);
    
19.     mpfr_set(result, pi, MPFR_RNDZ);
    
20.     mpfr_clears(pi, a, b, temp, (mpfr_ptr)0);
    
21. }
    
22. 
    
23. int main(void) {
    
24.     mpfr_t result;
    
25.     mpfr_init2(result, 65536);
    
26.     get_pi(result, 1e-150);
    
27.     mpfr_printf("%.100Rf\n", result);
    
28.     mpfr_clear(result);
    
29.     mpfr_free_cache();
    
30.     return 0;
    
31. }
    
32. $ gcc-debug -o main main.c $(pkg-config --cflags --libs mpfr)
    
33. $ time ./main
    
34. 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170680
    
35. 
    
36. real        0m0.144s
    
37. user        0m0.130s
    
38. sys        0m0.014s
    
39. $

复制代码

顶级太阳

新人求助，你们是怎么把这个代码写到帖子里的？尝试了发帖和回复里面的各个功能都没能实现。求各位指点。

wp231957

顶级太阳 发表于 2022-5-19 09:19
新人求助，你们是怎么把这个代码写到帖子里的？尝试了发帖和回复里面的各个功能都没能实现。求各位指点。

有一对尖括号，你点它就对了