鱼C论坛

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

c语言精确度问题

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

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

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

x
#include <stdio.h>
#include <math.h>
double  fun ( double  eps)
{
        double pi,a=1.0,b=1.0,i=1.0,t;
        pi=1.0;
        t=1.0;
        while(t>=eps)
        {
                a*=i;
                b*=(2*i+1);
                t=a/b;
                pi+=t;
                i++;
        }
        return (2*pi);


}

void main( )
{ double  x;void NONO ();
  printf("Input eps:") ;
  scanf("%lf",&x); printf("\neps = %lf, PI=%lf\n", x, fun(x));
  NONO();
}

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

  fp = fopen("in.dat","r") ;
  wf = fopen("out.dat","w") ;
  for(i = 0 ; i < 10 ; i++) {
    fscanf(fp, "%lf", &x) ;
    fprintf(wf, "%lf\n", fun(x)) ;
  }
  fclose(fp) ;
  fclose(wf) ;
}
fun 函数里面哪里不对啊跟结果不大一样
[img]%5Bimg%5Dhttps://s3.bmp.ovh/imgs/2022/05/18/c4a7481e7128766d.jpg[/img][/img]
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复

使用道具 举报

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

long double get_pi(long double eps) {
    long double pi = 1;
    for(size_t n = 1; ; ++n) {
        long double a = 1;
        for(size_t i = 1; i <= n; ++i) {
            a *= i;
        }
        long double b = 1;
        for(size_t i = 3; i <= 2 * n + 1; i += 2) {
            b *= i;
        }
        long double temp = a / b;
        if(temp < eps) break;
        pi += temp;
    }
    return 2 * pi;
}

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

使用道具 举报

发表于 2022-5-19 02:22:17 | 显示全部楼层
pi = 3.14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 58209 74944 59230 78164 06286 20899 86280 34825 34211 70679 
$ cat main.c
#include <stdio.h>
#include <mpfr.h>

long double get_pi(long double eps) {
    mpfr_t pi, a, b, temp;
    mpfr_inits2(256, pi, a, b, temp, (mpfr_ptr)0);
    mpfr_set_ld(pi, 1, MPFR_RNDZ);
    for(size_t n = 1; ; ++n) {
        mpfr_set_ld(a, 1, MPFR_RNDZ);
        for(size_t i = 1; i <= n; ++i) {
            mpfr_mul_ui(a, a, i, MPFR_RNDZ);
        }
        mpfr_set_ld(b, 1, MPFR_RNDZ);
        for(size_t i = 3; i <= 2 * n + 1; i += 2) {
            mpfr_mul_ui(b, b, i, MPFR_RNDZ);
        }
        mpfr_div(temp, a, b, MPFR_RNDZ);
        if(mpfr_cmp_ld(temp, eps) < 0) break;
        mpfr_add(pi, pi, temp, MPFR_RNDZ);
    }
    long double result = mpfr_get_ld(pi, MPFR_RNDZ);
    mpfr_clears(pi, a, b, temp, (mpfr_ptr)0);
    mpfr_free_cache();
    return 2 * result;
}

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

real        0m8.716s
user        0m8.698s
sys        0m0.001s
$

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

使用道具 举报

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

double fun(double eps) {
    double pi, a = 1.0, b = 1.0, i = 1.0, t;
    pi = 1.0;
    t = 1.0;
    //while(t >= eps) {
    while(1) {
        a *= i;
        b *= (2 * i + 1);
        t = a / b;
        if(t < eps) break;
        pi += t;
        i++;
    }
    return (2 * pi);
}

// void main( )
int main(void) {
    double x;
    void NONO();
    //printf("Input eps:");
    //scanf("%lf", &x);
    x = 0.0005;
    printf("\neps = %lf, PI=%lf\n", x, fun(x));
    //NONO();
}

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

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

使用道具 举报

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

void get_pi(mpfr_t *result, const mpfr_t eps) {
    mpfr_t pi, a, b, temp;
    mpfr_inits2(65536, pi, a, b, temp, (mpfr_ptr)0);
    mpfr_set_ld(pi, 1, MPFR_RNDZ);
    mpfr_set_ld(a, 1, MPFR_RNDZ);
    mpfr_set_ld(b, 1, MPFR_RNDZ);
    for(size_t n = 1; ; ++n) {
        mpfr_mul_ui(a, a, n, MPFR_RNDZ);
        mpfr_mul_ui(b, b, 2 * n + 1, MPFR_RNDZ);
        mpfr_div(temp, a, b, MPFR_RNDZ);
        if(mpfr_cmp(temp, eps) < 0) break;
        mpfr_add(pi, pi, temp, MPFR_RNDZ);
    }
    mpfr_mul_ui(pi, pi, 2, MPFR_RNDZ);
    mpfr_set(*result, pi, MPFR_RNDZ);
    mpfr_clears(pi, a, b, temp, (mpfr_ptr)0);
}

int main(void) {
    mpfr_t eps, result;
    mpfr_inits2(65536, eps, result, (mpfr_ptr)0);
    mpfr_set_str(eps, "1e-150", 10, MPFR_RNDZ);
    get_pi(&result, eps);
    mpfr_printf("%.100Rf\n", result);
    mpfr_clears(eps, result, (mpfr_ptr)0);
    mpfr_free_cache();
    return 0;
}
$ gcc-debug -o main main.c $(pkg-config --cflags --libs mpfr)
$ time ./main
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170680

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

使用道具 举报

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

void get_pi(mpfr_t result, double eps) {
    mpfr_t pi, a, b, temp;
    mpfr_inits2(65536, pi, a, b, temp, (mpfr_ptr)0);
    mpfr_set_ld(pi, 1, MPFR_RNDZ);
    mpfr_set_ld(a, 1, MPFR_RNDZ);
    mpfr_set_ld(b, 1, MPFR_RNDZ);
    for(size_t n = 1; ; ++n) {
        mpfr_mul_ui(a, a, n, MPFR_RNDZ);
        mpfr_mul_ui(b, b, 2 * n + 1, MPFR_RNDZ);
        mpfr_div(temp, a, b, MPFR_RNDZ);
        if(mpfr_cmp_d(temp, eps) < 0) break;
        mpfr_add(pi, pi, temp, MPFR_RNDZ);
    }
    mpfr_mul_ui(pi, pi, 2, MPFR_RNDZ);
    mpfr_set(result, pi, MPFR_RNDZ);
    mpfr_clears(pi, a, b, temp, (mpfr_ptr)0);
}

int main(void) {
    mpfr_t result;
    mpfr_init2(result, 65536);
    get_pi(result, 1e-150);
    mpfr_printf("%.100Rf\n", result);
    mpfr_clear(result);
    mpfr_free_cache();
    return 0;
}
$ gcc-debug -o main main.c $(pkg-config --cflags --libs mpfr)
$ time ./main
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170680

real        0m0.144s
user        0m0.130s
sys        0m0.014s
$
想知道小甲鱼最近在做啥?请访问 -> 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-11-17 16:46

Powered by Discuz! X3.4

© 2001-2023 Discuz! Team.

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