鱼C论坛

 找回密码
 立即注册
查看: 1338|回复: 7

c语言精确度问题

[复制链接]
发表于 2022-5-18 21:51:27 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能^_^

您需要 登录 才可以下载或查看,没有账号?立即注册

x
  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. void main( )
  19. { double  x;void NONO ();
  20.   printf("Input eps:") ;
  21.   scanf("%lf",&x); printf("\neps = %lf, PI=%lf\n", x, fun(x));
  22.   NONO();
  23. }

  24. void NONO ()
  25. {/* 本函数用于打开文件,输入数据,调用函数,输出数据,关闭文件。 */
  26.   FILE *fp, *wf ;
  27.   int i ;
  28.   double x ;

  29.   fp = fopen("in.dat","r") ;
  30.   wf = fopen("out.dat","w") ;
  31.   for(i = 0 ; i < 10 ; i++) {
  32.     fscanf(fp, "%lf", &x) ;
  33.     fprintf(wf, "%lf\n", fun(x)) ;
  34.   }
  35.   fclose(fp) ;
  36.   fclose(wf) ;
  37. }
复制代码

fun 函数里面哪里不对啊跟结果不大一样
[img]%5Bimg%5Dhttps://s3.bmp.ovh/imgs/2022/05/18/c4a7481e7128766d.jpg[/img][/img]
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复

使用道具 举报

发表于 2022-5-19 01:26:51 | 显示全部楼层
  1. #include <stdio.h>

  2. long double get_pi(long double eps) {
  3.     long double pi = 1;
  4.     for(size_t n = 1; ; ++n) {
  5.         long double a = 1;
  6.         for(size_t i = 1; i <= n; ++i) {
  7.             a *= i;
  8.         }
  9.         long double b = 1;
  10.         for(size_t i = 3; i <= 2 * n + 1; i += 2) {
  11.             b *= i;
  12.         }
  13.         long double temp = a / b;
  14.         if(temp < eps) break;
  15.         pi += temp;
  16.     }
  17.     return 2 * pi;
  18. }

  19. int main(void) {
  20.     printf("%.30Lf\n", get_pi(1e-4950L));
  21.     return 0;
  22. }
复制代码
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

发表于 2022-5-19 02:22:17 | 显示全部楼层
  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. long double get_pi(long double eps) {
  5.     mpfr_t pi, a, b, temp;
  6.     mpfr_inits2(256, pi, a, b, temp, (mpfr_ptr)0);
  7.     mpfr_set_ld(pi, 1, MPFR_RNDZ);
  8.     for(size_t n = 1; ; ++n) {
  9.         mpfr_set_ld(a, 1, MPFR_RNDZ);
  10.         for(size_t i = 1; i <= n; ++i) {
  11.             mpfr_mul_ui(a, a, i, MPFR_RNDZ);
  12.         }
  13.         mpfr_set_ld(b, 1, MPFR_RNDZ);
  14.         for(size_t i = 3; i <= 2 * n + 1; i += 2) {
  15.             mpfr_mul_ui(b, b, i, MPFR_RNDZ);
  16.         }
  17.         mpfr_div(temp, a, b, MPFR_RNDZ);
  18.         if(mpfr_cmp_ld(temp, eps) < 0) break;
  19.         mpfr_add(pi, pi, temp, MPFR_RNDZ);
  20.     }
  21.     long double result = mpfr_get_ld(pi, MPFR_RNDZ);
  22.     mpfr_clears(pi, a, b, temp, (mpfr_ptr)0);
  23.     mpfr_free_cache();
  24.     return 2 * result;
  25. }

  26. int main(void) {
  27.     printf("%.30Lf\n", get_pi(1e-4950L));
  28.     return 0;
  29. }
  30. $ gcc -Wall -O3 -o main main.c $(pkg-config --cflags --libs mpfr)
  31. $ time ./main
  32. 3.141592653589793238295968524909

  33. real        0m8.716s
  34. user        0m8.698s
  35. sys        0m0.001s
  36. $
复制代码


看起来这个算法到 3.141592653589793238 就是极限了
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

发表于 2022-5-19 02:31:24 | 显示全部楼层
  1. #include <stdio.h>
  2. #include <math.h>

  3. double fun(double eps) {
  4.     double pi, a = 1.0, b = 1.0, i = 1.0, t;
  5.     pi = 1.0;
  6.     t = 1.0;
  7.     //while(t >= eps) {
  8.     while(1) {
  9.         a *= i;
  10.         b *= (2 * i + 1);
  11.         t = a / b;
  12.         if(t < eps) break;
  13.         pi += t;
  14.         i++;
  15.     }
  16.     return (2 * pi);
  17. }

  18. // void main( )
  19. int main(void) {
  20.     double x;
  21.     void NONO();
  22.     //printf("Input eps:");
  23.     //scanf("%lf", &x);
  24.     x = 0.0005;
  25.     printf("\neps = %lf, PI=%lf\n", x, fun(x));
  26.     //NONO();
  27. }

  28. void NONO() { /* 本函数用于打开文件,输入数据,调用函数,输出数据,关闭文件。 */
  29.     FILE *fp, *wf;
  30.     int i;
  31.     double x;

  32.     fp = fopen("in.dat", "r");
  33.     wf = fopen("out.dat", "w");
  34.     for(i = 0; i < 10; i++) {
  35.         fscanf(fp, "%lf", &x);
  36.         fprintf(wf, "%lf\n", fun(x));
  37.     }
  38.     fclose(fp);
  39.     fclose(wf);
  40. }
复制代码
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

发表于 2022-5-19 03:01:44 | 显示全部楼层
原来是mpfr_t转long double出现了问题
其实根本不需要1e-4950,1e-150就可以精确到小数点后100位

  1. $ cat main.c
  2. #include <stdio.h>
  3. #include <mpfr.h>

  4. void get_pi(mpfr_t *result, const mpfr_t eps) {
  5.     mpfr_t pi, a, b, temp;
  6.     mpfr_inits2(65536, pi, a, b, temp, (mpfr_ptr)0);
  7.     mpfr_set_ld(pi, 1, MPFR_RNDZ);
  8.     mpfr_set_ld(a, 1, MPFR_RNDZ);
  9.     mpfr_set_ld(b, 1, MPFR_RNDZ);
  10.     for(size_t n = 1; ; ++n) {
  11.         mpfr_mul_ui(a, a, n, MPFR_RNDZ);
  12.         mpfr_mul_ui(b, b, 2 * n + 1, MPFR_RNDZ);
  13.         mpfr_div(temp, a, b, MPFR_RNDZ);
  14.         if(mpfr_cmp(temp, eps) < 0) break;
  15.         mpfr_add(pi, pi, temp, MPFR_RNDZ);
  16.     }
  17.     mpfr_mul_ui(pi, pi, 2, MPFR_RNDZ);
  18.     mpfr_set(*result, pi, MPFR_RNDZ);
  19.     mpfr_clears(pi, a, b, temp, (mpfr_ptr)0);
  20. }

  21. int main(void) {
  22.     mpfr_t eps, result;
  23.     mpfr_inits2(65536, eps, result, (mpfr_ptr)0);
  24.     mpfr_set_str(eps, "1e-150", 10, MPFR_RNDZ);
  25.     get_pi(&result, eps);
  26.     mpfr_printf("%.100Rf\n", result);
  27.     mpfr_clears(eps, result, (mpfr_ptr)0);
  28.     mpfr_free_cache();
  29.     return 0;
  30. }
  31. $ gcc-debug -o main main.c $(pkg-config --cflags --libs mpfr)
  32. $ time ./main
  33. 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170680

  34. real        0m0.164s
  35. user        0m0.150s
  36. sys        0m0.014s
  37. $
复制代码
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

发表于 2022-5-19 03:06:55 | 显示全部楼层
  1. $ cat main.c
  2. #include <stdio.h>
  3. #include <mpfr.h>

  4. void get_pi(mpfr_t result, double eps) {
  5.     mpfr_t pi, a, b, temp;
  6.     mpfr_inits2(65536, pi, a, b, temp, (mpfr_ptr)0);
  7.     mpfr_set_ld(pi, 1, MPFR_RNDZ);
  8.     mpfr_set_ld(a, 1, MPFR_RNDZ);
  9.     mpfr_set_ld(b, 1, MPFR_RNDZ);
  10.     for(size_t n = 1; ; ++n) {
  11.         mpfr_mul_ui(a, a, n, MPFR_RNDZ);
  12.         mpfr_mul_ui(b, b, 2 * n + 1, MPFR_RNDZ);
  13.         mpfr_div(temp, a, b, MPFR_RNDZ);
  14.         if(mpfr_cmp_d(temp, eps) < 0) break;
  15.         mpfr_add(pi, pi, temp, MPFR_RNDZ);
  16.     }
  17.     mpfr_mul_ui(pi, pi, 2, MPFR_RNDZ);
  18.     mpfr_set(result, pi, MPFR_RNDZ);
  19.     mpfr_clears(pi, a, b, temp, (mpfr_ptr)0);
  20. }

  21. int main(void) {
  22.     mpfr_t result;
  23.     mpfr_init2(result, 65536);
  24.     get_pi(result, 1e-150);
  25.     mpfr_printf("%.100Rf\n", result);
  26.     mpfr_clear(result);
  27.     mpfr_free_cache();
  28.     return 0;
  29. }
  30. $ gcc-debug -o main main.c $(pkg-config --cflags --libs mpfr)
  31. $ time ./main
  32. 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170680

  33. real        0m0.144s
  34. user        0m0.130s
  35. sys        0m0.014s
  36. $
复制代码
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

发表于 2022-5-19 09:19:19 | 显示全部楼层
新人求助,你们是怎么把这个代码写到帖子里的?尝试了发帖和回复里面的各个功能都没能实现。求各位指点。
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

发表于 2022-5-19 09:51:50 From FishC Mobile | 显示全部楼层
顶级太阳 发表于 2022-5-19 09:19
新人求助,你们是怎么把这个代码写到帖子里的?尝试了发帖和回复里面的各个功能都没能实现。求各位指点。

有一对尖括号,你点它就对了
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

小黑屋|手机版|Archiver|鱼C工作室 ( 粤ICP备18085999号-1 | 粤公网安备 44051102000585号)

GMT+8, 2024-3-28 19:44

Powered by Discuz! X3.4

© 2001-2023 Discuz! Team.

快速回复 返回顶部 返回列表