程序计算精确圆周率π的方法

原创已于 2025-09-03 21:05:06 修改 · 1.6w 阅读

15 ·

CC 4.0 BY-SA版权

文章标签：

#算法

于 2016-10-27 16:24:21 首次发布

c-c++ 同时被 2 个专栏收录

8 篇文章

订阅专栏

Linux开发

4 篇文章

订阅专栏

本文介绍了一种使用简单级数计算圆周率π的有效方法，并通过C语言实现了算法，能够计算到小数点后1000位以上的精度。此外，还展示了如何利用相似原理计算自然底数e。

圆周率π=3.14159265358……[事实上π是超越数], 自古以来，人们就对π的精确计算充满兴趣。从最早的几何方法到现代的计算机算法，计算π的精度逐步提高。本文将介绍快速计算圆周率的算法以及它们在计算机上的应用。

首先给一种极简级的获取方法：python+mpmath库【pip install mpmath】

from mpmath import mp
mp.dps = 1000 
print( mp.pi)   ## 输出1000位圆周率

Python 运行效果
这里提前为大家轻松列举1000位圆周率，其实mpmath内部的算法也是用c语言实现的。

下面我们选一种关于π比较简单, 并且收敛非常快的级数。
π /2= 1 + 1/3 + (1/3)(2/5) + (1/3)(2/5)(3/7) +(1/3)(2/5)(3/7)(4/9) + …
调整下：
π = 2 + 2/3 + (2/3)(2/5)+(2/3)(2/5)(3/7)+(2/3)(2/5)(3/7)(4/9) + …

void main()
{
	double x = 2, z = 2;
	int a = 1, b = 3;
	while (z > 0)
	{   //when b->max, z->0
		z = z * a / b;
		x += z;
		a++;
		b += 2;
	}
	printf("π=%.14f\n", x);
}

执行效果： π=3.1415926535898 ；这个程序如此简单它可以计算到小数点后14位，因受限于double的精度。

感兴趣的可以使用GMP大数库，MPFR高精度库再做如下尝试】

#include <stdio.h>
#include <gmp.h>

int main()
{
    // 初始化高精度变量
    mpf_set_default_prec(256);  // 设置全局默认精度，单位是比特，可以根据需要调整
    mpf_t x, z, temp_a, temp_b; // 定义高精度变量
    mpf_t small_zero;
    
    mpf_init(x);
    mpf_init(z);
    mpf_init(temp_a);
    mpf_init(temp_b);
    mpf_init(small_zero);

    // 初始化变量值
    mpf_set_ui(x, 2);      // x = 2
    mpf_set_ui(z, 2);      // z = 2
    mpf_set_ui(temp_a, 1); // a = 1
    mpf_set_ui(temp_b, 3); // b = 3
    mpf_set_str(small_zero,"0.00000000000000000001",20 ); // small_zero = 0

    int cnt = 0;
    // 循环计算
    while (mpf_cmp(z, small_zero) > 0)
    { // 当 z > 0 时继续循环
        // z = z * a / b
        mpf_mul(z, z, temp_a); // z = z * a
        mpf_div(z, z, temp_b); // z = z / b
        // x += z
        mpf_add(x, x, z);
        // a++
        mpf_add_ui(temp_a, temp_a, 1);
        // b += 2
        mpf_add_ui(temp_b, temp_b, 2);
        gmp_printf("[%d]x=%.50Ff\n", ++cnt, x);

    }

    // 输出结果
    // 不是所有小数点都准确,通过对比上下两位数据就可以判断第几位开始不准确
    // 通过对比目前有效位是26位
    gmp_printf("[--]π=%.26Ff\n", x); 

    // 释放内存
    mpf_clear(x);
    mpf_clear(z);
    mpf_clear(temp_a);
    mpf_clear(temp_b);
    mpf_clear(small_zero);

    return 0;
    //LINUX: gcc main.c -lgmp;./a.out
}

把最开始的那个简单C程序改造下, 将double类型的x，z变成数组保存每位小数，让它精确到小数点后面 1000 位再测试试试！！


void main()
{
	const int ARRSIZE = 1010;
	const int DISPCNT = 1000;	 // 定义数组大小，显示位数
	char x[ARRSIZE], z[ARRSIZE]; // x[0] x[1] . x[2] x[3] x[4] .... x[ARRSIZE-1]
	int a = 1, b = 3, c, d, i;
	int Run = 1, Cnt = 0;
	memset(x, 0, ARRSIZE);
	memset(z, 0, ARRSIZE);
	x[1] = 2;
	z[1] = 2;

	while (Run && (++Cnt < 100000))
	{
		d = 0;
		i = ARRSIZE - 1;
		do
		{
			c = z[i] * a + d;
			z[i] = c % 10;
			d = c / 10;
		} while (i--);

		d = 0;
		i = 0;
		do
		{
			c = z[i] + d * 10;
			z[i] = c / b;
			d = c % b;
		} while (i++ < ARRSIZE);

		Run = 0;
		for (i = ARRSIZE - 1; i > 0; i--)
		{
			c = x[i] + z[i];
			x[i] = c % 10;
			x[i - 1] += c / 10;
			Run |= z[i];
		}
		a++;
		b += 2;
	}
	printf("calc %d times\r\n", Cnt);
	printf("Pai=%d.\r\n", x[1]);
	for (i = 0; i < DISPCNT; i++)
	{
		if (i && ((i % 100) == 0))
			printf("\r\n");
		printf("%d", (int)x[i + 2]);
	}
}

以下是运行的效果，可以看到经过3343轮计算后前面的999位数字都是对的！
在这里插入图片描述
这下心里就有底了, 是不是改变数组大小就可以计算更多位数呢？答案是肯定的。
如果把定义数组大小和显示位数改为:
const int ARRSIZE = 10100;
const int DISPCNT = 10000;
执行的结果表明计算33538轮后圆周率精度就可以达到惊人的10000位！
初步结论，这个算法可以通过计算大概3.3轮就能输出一位正确的数值。

趣味扩展
基于以上算法，其实我们就可以对π的一些属性展开研究，如统计π每个位数数字出现的概率，等等

算法提高精度的原理:
以上程序的原理是利用数组把计算结果保存起来, 其中数组每一项保存10进制数的一位,
小数点定位在数组第1个数和第2个数之间, 即小数点前面2位整数, 其余都是小数位。
利用电脑模拟四则运算的笔算方法来实现高精度的数据计算,没想到最原始的方法竟然是精度最高的。

利用这个原理我们就可以扩展到其他常量的计算如自然底数e=2.71828…任意位数

#include <stdio.h>
#include <string.h>

// e = 1+1/1!+1/2!+1/3!+…+1/n! = 1+1/2*(1+1/3*(1+…*(1+1/n))…)));
int main()
{
#define TIMES 100 /// 输出e有效位数
	int n, e, d, i, f[TIMES + 1];
	i = TIMES;
	while (i--)
		f[i] = 10;
	printf("e=2.");
	i = 0;
	e = 0;
	while (i++ < TIMES)
	{
		d = 0;
		for (n = TIMES; n > 0; n--)
		{
			d = d + f[n] * 10;
			f[n] = d % n; // 余数
			d = d / (n);  // 整数
		}
		printf("%d", (e + d / 10) % 10);
		e = d % 10;
	}
	sleep(1);
	printf("\n\n");
	return 0;
}

用python+mpmath打印下答案对比

import mpmath
mpmath.mp.dps = 1000
# 计算自然对数的底数e的值，1000位
print(mpmath.e)

#####代码赏析#########

Bellard在2000年的国际C语言代码混淆竞赛获奖作品，使用NTT数论变换在较短时间内计算出已知当时最大的梅森素数(2^6972593-1)
它是第 38 个梅森素数，于 1999 年由 Nayan Hajratwala 使用 GIMPS（Great Internet Mersenne Prime Search）项目发现。
它的十进制表示有 2,098,960 位，是一个极大的素数。

#include <stdio.h>

int m = 754974721, N, t[1 << 22], a, *p, i, e = 1 << 22, j, s, b, c, U;
f(d)
{
    for (s = 1 << 23; s; s /= 2, d = d * 1LL * d % m)
        if (s < N)
            for (p = t; p < t + N; p += s)
                for (i = s, c = 1; i; i--)
                    b = *p + p[s], p[s] = (m + *p - p[s]) * 1LL * c % m, *p++ = b % m, c = c * 1LL * d % m;
    for (j = 0; i < N - 1;)
    {
        for (s = N / 2; !((j ^= s) & s); s /= 2)
            ;
        if (++i < j)
            a = t[i], t[i] = t[j], t[j] = a;
    }
}

main()
{
    *t = 2;
    U = N = 1;
    while (e /= 2)
    {
        N *= 2;
        U = U * 1LL * (m + 1) / 2 % m;
        f(362);
        for (p = t; p < t + N;)
            *p++ = (*p * 1LL * *p % m) * U % m;
        f(415027540);
        for (a = 0, p = t; p < t + N;)
            a += (6972593 & e ? 2 : 1) * *p, *p++ = a % 10, a /= 10;
    }
    while (!*--p)
        ;
    t[0]--;
    while (p >= t)
        printf("%d", *p--);
}

作为文章结束，给出一个某年Obfuscated C Contest佳作选录: 800位精度圆周率

#include <stdio.h>
long a=10000, b, c=2800, d, e, f[2801], g; 
main(){
for(;b-c;)f[b++]=a/5; 
for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a) 
for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);
}