圆周率PI

算法设计

⑴选择计算公式
计算圆周率∏式很多，选取收敛速度快且容易操作的计算公式是设计的重要一环。

∏/2=1+1/3+(1•2)/(3•5)+(1•2•3)/(3•5•7)+⋯+(1•2•…n)/(3•5•…•(2n+1) )
=1+1/3(1+2/5(1+...(n-1)/(2n-1)(1+n/(2n+1))...)

⑵确定计算项数
可证明，式中分工第n项之后的所有余项之和Rn<an。因此，只要选取n，满足an<1/pow(10,x+1)即可。即只要使
lg 3 + lg 5/2 + ... + lg (2n+1)/n > x + 1
于是可设置对数累加实现计算到x位所需的项数n。为确保准确，算法可设置计算位数超过x位（例如x+5位），只打印输出x位。


⑶模拟乘除综合运算
设置a数组，下标根据计算位预设5000，必要时可增加。计算的整数值存放在a(0)中，小数点后第i位存放在a(i)(i=1,2,...)中。
依据公式，应用模拟乘除运算进行计算：
数组除以2n+1，乘以n，加上1；再除以2n-1，乘以n－1，加上1；……。这些数组操作设置在j(j=n,n-1,...,1）循环中实施。
按公式实施除法操作：被除数为c，除数d分别取2n+1,2n-1,...,3。商仍存取放在各数组元素(a(i)=c/d)中。余数(c%d)乘10加在后数组元素a(i+1)上，作为后一位的被除数。
按公式实施乘法操作：乘数j分别取n,n-1,...,i。乘积要注意进位，设进位数为b，则对计算的积a(i)=a(i)*j+b，取其十位以上数作为进位数b=a(i)/10，取其个数仍存放在原数组元素a(i)=a(i)%10中。
循环实施除乘操作完成后，按数组元素从高位到你们顺序输出。因计算位数较多，为方便查对，每一行控制打印50位，每10位空一格。

 

// 圆周率PI的程序实现 // 利用公式：∏/2=1+1/3+(1•2)/(3•5)+(1•2•3)/(3•5•7)+⋯+(1•2•…n)/(3•5•…•(2n+1) ) // =1+1/3(1+2/5(1+...(n-1)/(2n-1)(1+n/(2n+1))...) #include <stdio.h> #include <math.h> #define MAXSIZE 5000 int main(void) { double s; int b, x, n, c, i, j, d, a[MAXSIZE]; printf("请输入精确位数："); scanf("%d", &x); // 累加确定计算的项数n for (s = 0, n = 1; n <= MAXSIZE; n++) { s += log10((double)((2 * n + 1) / n)); if (s > x + 1) break; } for (i = 0; i <= x + 5; i++) a[i] = 0; // 按公式分步计算 for (c = 1, j = n; j >= 1; j--) { d = 2 * j + 1; for (i = 0; i <= x + 4; i++) { // 各位实施除2j+1 a[i] = c / d; c = (c % d) * 10 + a[i + 1]; } a[x + 5] = c / d; for (b = 0, i = x + 5; i >= 0; i--) { // 各位实施乘j a[i] = a[i] * j + b; b = a[i] / 10; a[i] %= 10; } a[0]++; c = a[0]; // 整数位加1 } for (b = 0, i = x + 5; i >= 0; i--) { // 按公式各位乘2 a[i] = a[i] * 2 + b; b = a[i] / 10; a[i] %= 10; } printf("/nPI = %d.", a[0]); // 逐位输出计算结果 for (i = 1; i <= x; i++) printf("%d", a[i]); putchar('/n'); return 0; }