需要帮助修复一个近似于 pi 的算法答案

【问题标题】：Need help fixing an algorithm that approximates pi需要帮助修复一个近似于 pi 的算法
【发布时间】：2020-04-10 17:33:48
【问题描述】：

我正在尝试为近似 pi 的算法编写 C 代码。它应该得到立方体的体积和立方体内球体的体积（球体的半径是立方体边的 1/2）。然后我应该将立方体的体积除以球体的体积并乘以 6 得到 pi。

它正在工作，但它在应该获得卷的部分做了一些奇怪的事情。我认为这与我为近似值选择的增量有关。用 4 边的立方体代替 64 的体积，它给了我 6400。用球体代替 33，它给了我 3334。一些东西。

有人能弄清楚吗？这是代码（我注释了相关部分）：

#include <stdio.h>      

int in_esfera(double x, double y, double z, double r_esfera){
    double dist = (x-r_esfera)*(x-r_esfera) + (y-r_esfera)*(y-r_esfera) + (z-r_esfera)*(z-r_esfera);

    return  dist <= (r_esfera)*(r_esfera) ? 1 : 0;   
}   

double get_pi(double l_cubo){   
    double r_esfera = l_cubo/2;   
    double total = 0;
    double esfera = 0;    
//this is delta, for the precision. If I set it to 1E anything less than -1 the program continues endlessly. Is this normal?
    double delta = (1E-1);   

    for(double x = 0; x < l_cubo; x+=delta){
        printf("x => %f; delta => %.6f\n",x,delta);
        for(double y = 0; y <l_cubo; y+=delta){
            printf("y => %f; delta => %.6f\n",y,delta);
            for(double z = 0; z < l_cubo; z+=delta){
                printf("z => %f; delta => %.6f\n",z,delta);
                total+=delta;
                if(in_esfera(x,y,z,r_esfera))
                    esfera+=delta;
            }
        }
    }

    //attempt at fixing this
        //esfera/=delta;
        //total/=delta;
    //

//This printf displays the volumes. Notice how the place of the point is off. If delta isn't a power of 10 the values are completely wrong.   
    printf("v_sphere = %.8f; v_cube = %.8f\n",esfera,total);   

    return (esfera)/(total)*6;
}   

void teste_pi(){        
    double l_cubo = 4;    
    double pi = get_pi(l_cubo);

    printf("%.8f\n",pi);
}   

int main(){   
    teste_pi();
}

【问题讨论】：

您的get_pi 函数类似于O((l_cubo / delta)**3)（立方复杂性！难怪它“无休止地继续”）。它应该计算立方体和球体的体积吗？
@ForceBru 部分。首先是它的作用（检查该函数的最后一个 printf），然后返回比率 * 6（或 pi）。
事实证明，您需要将两个卷乘以(delta * delta)，如esfera *= delta * delta; 来解决此问题。这是简单的数学，但我不太确定如何用语言表达lol
@ForceBru 即使 delta 不是 10 的幂，它也有效。谢谢！我很好奇它为什么会起作用......
@ForceBru 如果他们打算降低时间复杂度，他们可以选择方形和圆形。此外，似乎稍微降低了错误。

标签： c optimization pi approximation

【解决方案1】：

total+=delta;
if(in_esfera(x,y,z,r_esfera))
    esfera+=delta;

total 和 esfera 是三维体积，而 delta 是一维长度。如果您要跟踪单位，则左侧有 m³，右侧有 m。单位不兼容。

要修复它，请使用立方体 delta，这样您就可以在概念上累积微小的立方体而不是微小的线条。

total+=delta*delta*delta;
if(in_esfera(x,y,z,r_esfera))
    esfera+=delta*delta*delta;

这样做可以修复输出，也适用于delta 的任何值：

v_sphere = 33.37400000; v_cube = 64.00000000
3.12881250

请注意，此算法对任意 delta 值“有效”，但存在严重的准确性问题。它非常容易出现舍入问题。当delta 是 2 的幂时效果最佳：1/64.0 优于 1/100.0，例如：

v_sphere = 33.50365448; v_cube = 64.00000000
3.14096761

此外，如果您希望程序运行得更快，请摆脱所有这些打印输出！或者至少是内部循环中的那些......

【讨论】：

【解决方案2】：

问题是像a * b * c 这样的整数乘法与将1 + 1 + 1 + 1 + ... + 1 a * b * c 相加是一样的，对吧？

您正在添加delta + delta + ... (x / delta) * (y / delta) * (z / delta) 次。或者，换句话说，(x * y * z) / (delta ** 3) 次。

现在，deltas 的总和与此相同：

delta * (1 + 1 + 1 + 1 + ...)
         ^^^^^^^^^^^^^^^^^^^^ (x * y * z) / (delta**3) times

所以，如果delta 是 10 的幂，(x * y * z) / (delta**3) 将是一个整数，它等于括号中 1 的总和（因为它与积 x * y * (z / (delta**3))，其中最后一项是整数 - 请参阅此答案的第一句话）。因此，您的结果将如下所示：

delta * ( (x * y * z) / (delta ** 3) ) == (x * y * z) / (delta**2)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ the sum of ones

这就是您最终计算乘积的方式除以delta 平方。

要解决这个问题，请将所有卷乘以 delta * delta。

但是，我认为对于不是 10 的幂的deltas 使用此逻辑是不可能的。事实上，对于 delta == 0.21 和 l_cubo == 2，代码会出现各种混乱，例如：你会得到 9.261000000000061 而不是 8。

【讨论】：