多种解法计算圆周率π.docx
《多种解法计算圆周率π.docx》由会员分享，可在线阅读，更多相关《多种解法计算圆周率π.docx（14页珍藏版）》请在冰豆网上搜索。
多种解法计算圆周率π
课程设计报告
学院、系：
吉林大学珠海学院计算机科学与技术系
专业名称：
软件工程
课程设计科目
C语言程序课程设计
所在班级：
学生学号：
学生姓名：
指导教师：
曾志平
完成时间：
2012年3月-5月
题目：
多种解法计算圆周率π
一、设计任务与目标
此问题求圆周率在之前所学的C语言中已经接触过，但是算法单一，采用级数，而且收敛较慢，故运行时间较长，此程序设计要解决的问题是如何实现高精度的运算，如何对结果进行输出，并且尝试采用不同的方法进行求解圆周率。
本次上机实践所使用的平台和相关软件。
平台：
windows7
相关软件：
VC++6.0
二、方案设计与论证
随机数法求圆周率可以利用计算机中随机数函数模拟出两个0~1之间的浮点型点（x,y）,建立直角坐标系思想，利用边长为1的正方形内切半径为0.5圆的方程（x-0.5）*（x-0.5）+（y-0.5）*（y-0.5）<=0.25判断点是否在圆内，用计数器b保存在的点，如此模拟5000次。
用落入圆内的点数b的4倍除以总的扔的点数N用，可大概求出圆周率的值。
一般来说，根据概率思想，N值越大，模拟次数越多，其求出来圆周率的值越接近真实的圆周率。
祖冲之迭代法因为圆内接正六边形边长等于半径的思想，故可以从正六边形出发，不断迭代，当正多边形边数增加时，其周长也逐渐逼近圆的周长，反过来即可求的一定精度的圆周率，设圆内接正六边形的边长为b,边数为i，利用公式
进行迭代运算，为了提高精度，算法中对公式进行分开运算，求得边数为2i的圆内接正多边形后得出其周长，运用迭代后的正多边形周长减去迭代前的正多边形周长，获取其精度程度。
如果最终求出的圆内接正多边形的周长，即接近圆的周长，最后利用数学公式即可以求出圆周率。
输出其最后迭代后得圆内接正多边形的边数和圆周率即可。
用级数法求圆周率，定义一个a[400]的数组用于存储计算结果，从结果出发，因为其要输出一定精度的圆周率，若n值太大会造成计算冗余，利用数学中不等式确定其n项，用一个循环从n~1计算每一项的值并存储，1+n/（2n+1）用数组模拟手工乘除加法，除法1/（2n+1），相除后商为a[0]，然后将余数乘以10，作为被除数再除以除数取商为十分位，存于a[1],如此类推；乘法则每个数组乘以n，如果满十则向前面数组进一，再取其个位存储；加法因为加1，则直接可以加到a[0]上。
然后保存计算结果后用该值计算n-1项，如此重复一直到第1项后按照格式输出每个数组元素即可。
附加题圆外切正多边形与圆内接正多边形算法相似，但是圆半径的长度为1时，其外切正六边形的边长为2*sqrt（3）/3，并且其迭代公式也相应的进行了改变，b=2*（sqrt（b*b+4）-2）/b;由于圆外切正多边形迭代后的周长小于迭代前的周长，故控制精度时用迭代前的正多边形周长减去迭代后的正多边形边长。
三、程序框图或流程图，程序清单与调用关系
由主函数分别调用voidsuiji（）（即课程第一小问，随机数函数），voidneijie（）（第二小问，内接正多边形法），voidjishu（）（第三小问，级数法）,voidwaiqie（）（外切正多边形法）.
由于jishu函数的流程图过长，造成排版不便，故进行了分段和省略控制输出等部分流程图。
四、全部源程序清单
#include
#include
#include
#include
#defineN5000
voidmain（）
{
voidjishu（）;
voidneijie（）;
voidsuiji（）;
voidwaiqie（）;
suiji（）;
printf（"-------------------------------------------------------\n"）;
neijie（）;
printf（"-------------------------------------------------------\n"）;
jishu（）;
printf（"-------------------------------------------------------\n"）;
waiqie（）;
printf（"-------------------------------------------------------\n"）;
}
voidsuiji（）
{
doublex,y;
inta=0,b=0;
srand（time（0））;
while（a++<=N）//投5000次点
{
x=（double）rand（）/RAND_MAX;//产生0~1之间的浮点数
y=（double）rand（）/RAND_MAX;//产生0~1之间的浮点数
if（（x-0.5）*（x-0.5）+（y-0.5）*（y-0.5）<=0.25）//判断所产生的点是否在圆内
b++;//汇总落在圆内的点数
}
printf（"用随机数法求得π=%lf\n",4.0*b/N）;
}
voidneijie（）
{
doublee,b=1,d;//b：
为正多边形边长
longinti;//i：
正多边形边数
for（i=6;;i*=2）//正多边形边数加倍
{
d=1.0-sqrt（1.0-b*b/4）;//计算圆内接正多边形的边长
b=sqrt（b*b/4+d*d）;
if（2*i*b-i*e<1e-15）//2*i*b迭代后的周长，i*e原来的周长
break;//精度达1e-15则停止计算
e=b;//保存本次正多边形的边长作为下一次精度控制的依据
}
printf（"用圆内接正多边形计算π=%.15lf\n迭代后得正多边形的边数为：
%d\n",i*b,i）;
}
voidwaiqie（）
{
doublee,b=2*sqrt（3）/3,d;//b：
为正多边形边长
longinti;//i：
正多边形边数
for（i=6;;i*=2）//正多边形边数加倍
{
e=b;//e：
暂存边数为加倍前的边长
d=sqrt（1.0+b*b/4）-1;//计算圆外切正多边形的边长
b=2*（b*b/4-d*d）/b;
if（i*e-2*i*b<1e-15）
break;//精度达1e-15则停止计算
e=b;//保存本次正多边形的边长作为下一次精度控制的依据
}
printf（"附加题:
外切正多边形计算π=%.15lf\n迭代后得正多边形的边数为：
%d\n",i*b,i）;
}
voidjishu（）
{
floats;
intx,b,c,i,j,n,e,a[400];
for（s=0,n=1;n<=400;n++）//通过小数精确程度对项数n进行确定
{
s=s+log（（2*n+1）/n）;
if（s/log（10）>100）break;
}
for（i=0;i<=101;i++）a[i]=0;
for（x=1,b=n;b>=1;b--）
{
c=2*b+1;//首先用除法计算1/2n+1的高精度商
for（i=0;i<=100;i++）
{
a[i]=x/c;x=（x%c）*10+a[i+1];/*用数组模拟手工除法，相除后商为a[0]，然后将余数添0，再除以除数取商为十分位，存于a[1],如此类推*/
}
a[101]=x/c;
for（j=0,i=101;i>=0;i--）//计算n/（2n+1）的高精度值
{
a[i]=a[i]*b+j;//数组模拟手工乘法，每个数组乘以n
j=a[i]/10;//如果満十则对前面数组进1
a[i]=a[i]%10;//而该数组取个位
}
a[0]=a[0]+1;//计算1+n/（2n+1）的高精度值
x=a[0];
}
for（j=0,i=101;i>=0;i--）
{
a[i]=a[i]*2+j;j=a[i]/10;a[i]=a[i]%10;//因为公式左边还有个1/2，对结果整体乘以二，模拟手工乘法
}
printf（"用级数方法求得π=%d.",a[0]）;//控制格式并输出结果
for（e=0,i=1;i<=100;i++）
{
printf（"%d",a[i]）;
e++;
if（e%30==0）
printf（"\n"）;
}
printf（"\n"）;
}五、程序运行结果测试与分析
由输出结果可以看出，随机数法求得圆周率浮动较大，也许是N值设为5000还过于小造成结果误差大，然后把N值设为50000后继续测试如下
当N值为50000时，测试结果较为接近圆周率且比较数值浮动较小，其余小问测试结果基本接近真实圆周率的值，基本完成程序设计目标。
六、结论与心得
1．设计中遇到的问题及解决过程
在设计过程中遇到第三小问需要保存100位圆周率的时候，和同学讨论后才明白结果可以用数组保存并输出，
附加小问时圆外切正多边形时刚来采用的是d=sqrt（1.0+b*b/4）-1；b=2*（sqrt（b*b+4）-2）/b，后来发现误差很大，改用公式b=2*（b*b/4-d*d）/b才能够接近，思考后发现b=2*（b*b/4-d*d）/b（误差较小）将d代入，得b=2*（b*b/4-（sqrt（b*b/4+1）-1）*（sqrt（b*b/4+1）-1））/b;（就这样写误差还是比较小的，几乎跟内切法的精度一样），再化简得b=2*（b*b/4-（b*b/4+1）-1+2*sqrt（b*b/4+1））/b;（这个时候误差就大了，i只有49152），再化简b=2*（2*sqrt（b*b/4+1）-2）/b即b=2*（sqrt（b*b+4）-2）/b（误差依然比较大）,所以断定d*d为误差所在。
我的理解：
无论如何，b,d等这些值是一定有偏差的，b=2*（b*b/4-d*d）/b=2*（b/2+d）*（b/2-d）/b,由于d进行了开根处理，误差微乎其微，几乎可以忽略不计，即可以当他们是常数，这个式子的误差级别y=ax+b（a,b为常量），b=2*（sqrt（b*b+4）-2）/b的误差级别y=a/x+b（a,b为常量），x在（0，1]区间显然后者的误差比较大）；
调试程序的时候会遇到小错误，但经过反复修改后程序能够运行。
2．设计体会和收获。
此次程序设计明白了对数据类型及其精度的重要性，有些科学研究需要高精度的数据，同时明白通向解决问题的道路是多种多样的，我们应当用不同方向和方法去解决问题。
在我们解决不了问题的时候，要尝试着去分析问题回归问题的最基本。
七、参考资料
[1],XX百科（祖冲之迭代法）
[2]杨克昌编著，计算机程序设计典型例题精解，国防科技大学大学出版社，1999年3月，第297-299页）
八、致谢
致谢曾志平老师对程序提出的改进和思路；
课程设计成绩评定表
对课程设计工作过程的简短介绍和自我评价
学生签名：
2012年4月11日
（以下由评定小组教师填写）
质量评价指标（在相应栏目打√）
评价项目
评价质量
优秀
良好
中等
及格
不及格
工作量和态度
实验、计算可靠性
文字和图表质量
总体评价
评定成绩（百分制）
评定小组成员签名
2012年月日
制定人：
王钲旋，单缅审定人：
陈守孔
升级会员 