c - C 中的乘法结果错误

Question

我正在编写 C 代码，使用阿基米德方法计算 pi，精度为 80 位小数。 This video可以阐明该方法。

我基本上是通过从直径为 2、边长为 1 的正六边形开始计算 Pi 并执行 Pi = 周长/直径，并将尺寸数量加倍，这样我就可以得到一个多边形每次我的程序计算它时都更接近一个圆。

这是代码:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int main()
{
    long double pi, circunferencia, s, s1, a, b, diametro;
    int lado = 6;

    s=1;
    diametro = 2;

    while (lado<=(6*pow(2,27))){

        printf("numero de lados =%d\n", lado);

        a= sqrt(1-(s/2)*(s/2));
      //  printf("a= %.80Lf\n", a);
        b= 1 - a;
      //  printf("b= %.80Lf\n", b);
        s1= sqrt ((b*b)+((s/2)*(s/2)));
      //  printf("c= %.80Lf\n", s1);


        lado = lado*2;
        s = s1;
        circunferencia = lado * s1;
        pi = circunferencia/diametro;
        printf("\npi= %.80Lf\n\n\n", pi);

    }

    return 0;
}

使用 2²⁷，我得到了大约 15 位小数的精度，并且 lado = 805306368。如果我使用 2²⁸ 或更高，乘法从这部分开始出错，我开始得到负数。一切都变得一团糟，lado突然被设置为0并且变成了无限循环。

我是 C 和编程的初学者，如果你们能给我推荐一些关于这个主题的读物，以了解它到底是如何工作的以及为什么我会遇到这个错误，那就太棒了。

我正在使用 Code::Blocks 版本 17.12、GNU GCC 编译器(我启用了“让 gcc 遵循 2011 ISO C 语言标准 [-std=c11]”选项，因为我无法请勿使用“%Lf”来引用long double变量)和Windows 10 x64操作系统。

Answer 1

If i do 2^28 or higher, the multiplication goes wrong from this part,

lado<=(6*pow(2,27)) 没问题，但 lado = lado*2; 大约在那时开始溢出 int 数学。

一旦lado > INT_MAX/2(可能是1,073,741,823)，int就会发生溢出，任何事情都可能发生。在OP的例子中，lado*2最终变成0。

代码可以使用 lado

// int lado = 6;
long long lado = 6;
...
  // printf("numero de lados =%d\n", lado);
  printf("numero de lados =%lld\n", lado);

这将修补立即循环问题大约另外 20 或 所以循环。

---

我怀疑代码能否用 long long 实现“以 80 位小数精度计算 pi”(可能是 17 - 34 位小数精度)。

---

OP 报告“我的老师说我可以使用 long double 存储最多 80 位数字”确定性意味着 long double 为 80 位 x86 extended precision format 或大约 20 位十进制数字的精度。

我建议使用 long double 函数而不是 long double 函数。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int main() {
  long double pi, circunferencia, s, s1, a, b, diametro;
  long long lado = 6;

  s = 1;
  diametro = 2;

  while (lado <= (6 * pow(2, 27))) {
    printf("numero de lados =%lld\t", lado);

    // a = sqrt(1 - (s / 2) * (s / 2));
    a = sqrtl(1 - (s / 2) * (s / 2));
    b = 1 - a;
    // s1 = sqrt((b * b) + ((s / 2) * (s / 2)));
    s1 = hypotl(b, s / 2);

    lado = lado * 2;
    s = s1;
    circunferencia = lado * s1;
    pi = circunferencia / diametro;
    printf("pi= %.40Lf\n", pi);
  }
  printf("best                      \tpi= %.40Lf\n", acosl(-1.0));
  return 0;
}

示例输出

numero de lados =6  pi= 3.1058285412302491482368360653509853364085
numero de lados =12 pi= 3.1326286132812381970275117382129792531487
numero de lados =24 pi= 3.1393502030468672071242958176995330177306
....
numero de lados =201326592  pi= 3.1415926535897932057699033503439522974077
numero de lados =402653184  pi= 3.1415926535897932296223511450250498455716
numero de lados =805306368  pi= 3.1415926535897932356938833109438746760134
best                        pi= 3.1415926535897932385128089594061862044327
// in my case good to about     1 23456789012345678 (18 digits)