h-b-z-d-j-s-m 发表于 2022-5-18 21:51:27

c语言精确度问题

#include <stdio.h>
#include <math.h>
doublefun ( doubleeps)
{
        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( )
{ doublex;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 函数里面哪里不对啊跟结果不大一样
%5Bimg%5Dhttps://s3.bmp.ovh/imgs/2022/05/18/c4a7481e7128766d.jpg

人造人 发表于 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;
}

人造人 发表于 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 就是极限了

人造人 发表于 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);
}

人造人 发表于 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
$

人造人 发表于 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
$

顶级太阳 发表于 2022-5-19 09:19:19

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

wp231957 发表于 2022-5-19 09:51:50

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

有一对尖括号，你点它就对了
页: [1]
查看完整版本: c语言精确度问题