需要帮助修复近似 pi 的算法
Need help fixing an algorithm that approximates pi
我正在尝试为近似 pi 的算法编写 C 代码。它应该获得立方体的体积和该立方体内部球体的体积(球体的半径是立方体边的 1/2)。然后我应该用立方体的体积除以球体的体积再乘以 6 得到圆周率。
它正在工作,但它在应该获得卷的部分做了一些奇怪的事情。我认为这与我为近似值选择的增量有关。
使用边长为 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();
}
问题是 整数 的乘法 a * b * c
等同于加法 1 + 1 + 1 + 1 + ... + 1
a * b * c
次,对吧?
您正在添加 delta + delta + ...
(x / delta) * (y / delta) * (z / delta)
次。或者,换句话说,(x * y * z) / (delta ** 3)
次。
现在,delta
s 的总和与此相同:
delta * (1 + 1 + 1 + 1 + ...)
^^^^^^^^^^^^^^^^^^^^ (x * y * z) / (delta**3) times
因此,如果 delta
是 10 的幂,(x * y * z) / (delta**3)
将是一个整数,并且等于括号中 1 的总和(因为它与 product x * y * (z / (delta**3))
,其中最后一项是整数 - 请参阅此答案的第一句话)。因此,您的结果将如下所示:
delta * ( (x * y * z) / (delta ** 3) ) == (x * y * z) / (delta**2)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ the sum of ones
这就是您最终计算乘积 除以 delta
平方 的方法。
要解决这个问题,请将所有体积乘以 delta * delta
。
但是,我认为不可能将此逻辑用于不是 10 的幂的 delta
。事实上,对于 [=25=,代码会出现各种混乱] 和 l_cubo == 2
,例如:您将得到 9.261000000000061
而不是 8.
total+=delta;
if(in_esfera(x,y,z,r_esfera))
esfera+=delta;
total
和 esfera
是三维体积,而 delta
是一维长度。如果您要跟踪单位,则左侧有 m3,右侧有 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
请注意,此算法 "works" 用于任意 delta
值,但存在严重的准确性问题。它非常容易出现舍入问题。当delta
是2的幂时效果最佳:1/64.0
优于1/100.0
,例如:
v_sphere = 33.50365448; v_cube = 64.00000000
3.14096761
此外,如果您希望您的程序 运行 更快地删除所有这些打印输出!或者至少是内部循环中的那些...
我正在尝试为近似 pi 的算法编写 C 代码。它应该获得立方体的体积和该立方体内部球体的体积(球体的半径是立方体边的 1/2)。然后我应该用立方体的体积除以球体的体积再乘以 6 得到圆周率。
它正在工作,但它在应该获得卷的部分做了一些奇怪的事情。我认为这与我为近似值选择的增量有关。 使用边长为 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();
}
问题是 整数 的乘法 a * b * c
等同于加法 1 + 1 + 1 + 1 + ... + 1
a * b * c
次,对吧?
您正在添加 delta + delta + ...
(x / delta) * (y / delta) * (z / delta)
次。或者,换句话说,(x * y * z) / (delta ** 3)
次。
现在,delta
s 的总和与此相同:
delta * (1 + 1 + 1 + 1 + ...)
^^^^^^^^^^^^^^^^^^^^ (x * y * z) / (delta**3) times
因此,如果 delta
是 10 的幂,(x * y * z) / (delta**3)
将是一个整数,并且等于括号中 1 的总和(因为它与 product x * y * (z / (delta**3))
,其中最后一项是整数 - 请参阅此答案的第一句话)。因此,您的结果将如下所示:
delta * ( (x * y * z) / (delta ** 3) ) == (x * y * z) / (delta**2)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ the sum of ones
这就是您最终计算乘积 除以 delta
平方 的方法。
要解决这个问题,请将所有体积乘以 delta * delta
。
但是,我认为不可能将此逻辑用于不是 10 的幂的 delta
。事实上,对于 [=25=,代码会出现各种混乱] 和 l_cubo == 2
,例如:您将得到 9.261000000000061
而不是 8.
total+=delta; if(in_esfera(x,y,z,r_esfera)) esfera+=delta;
total
和 esfera
是三维体积,而 delta
是一维长度。如果您要跟踪单位,则左侧有 m3,右侧有 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
请注意,此算法 "works" 用于任意 delta
值,但存在严重的准确性问题。它非常容易出现舍入问题。当delta
是2的幂时效果最佳:1/64.0
优于1/100.0
,例如:
v_sphere = 33.50365448; v_cube = 64.00000000
3.14096761
此外,如果您希望您的程序 运行 更快地删除所有这些打印输出!或者至少是内部循环中的那些...