C++经典算法题-PI算法

cwl_java

发布于 2022-11-30 08:40:44

7620

发布于 2022-11-30 08:40:44

文章被收录于专栏：cwl_Javacwl_Java

17.Algorithm Gossip: 长 PI

说明

圆周率后的小数位数是无止境的，如何使用电脑来计算这无止境的小数是一些数学家与程式设计师所感兴趣的，在这边介绍一个公式配合大数运算，可以计算指定位数的圆周率。

解法

首先介绍J.Marchin的圆周率公式：

PI = [16/5 - 16 / (3*53) + 16 / (5*55) - 16 / (7*57) + ......] -

[4/239 - 4/(3*2393) + 4/(5*2395) - 4/(7*2397) + ......]

可以将这个公式整理为：

PI = [16/5 - 4/239] - [16/(53) - 4/(2393)]/3+ [16/(55) - 4/(2395)]/5 + ......

也就是说第n项，若为奇数则为正数，为偶数则为负数，而项数表示方式为： [16/52*n-1 - 4/2392*n-1] / (2*n-1)

如果我们要计算圆周率至10的负L次方，由于[16/52*n-1 - 4/2392*n-1]中16/52*n-1比4/2392*n-1来的大，具有决定性，所以表示至少必须计算至第n项：

[16/52*n-1 ] / (2*n-1) = 10-L

将上面的等式取log并经过化简，我们可以求得： n = L / (2log5) = L / 1.39794

所以若要求精确度至小数后L位数，则只要求至公式的第n项，其中n等于： n = [L/1.39794] + 1

在上式中[]为高斯符号，也就是取至整数（不大于L/1.39794的整数）；为了计简方便，可以在程式中使用下面这个公式来计简第n项：

[W -1/52- V -1 / (2392)] / (2*n-1)

这个公式的演算法配合大数运算函式的演算法为： div(w, 25, w); div(v, 239, v);

div(v,  239, v);
sub(w, v, q);
div(q, 2*k-1, q)

至于大数运算的演算法，请参考之前的文章，必须注意的是在输出时，由于是输出阵列中的整数值，如果阵列中整数位数不满四位，则必须补上0，在C语言中只要使用格式指定字%04d，使得不足位数部份自动补上0再输出，至于Java的部份，使用 NumberFormat来作格式化。

代码示例

#include <stdio.h>
#define L 1000 
#define N L/4+1

	// L 为位数，N是array长度

    void add(int*, int*, int*); void sub(int*, int*, int*); void div(int*, int, int*);

    int main(void) {
        int s[N+3] = {0};
        int w[N+3] = {0};
        int v[N+3] = {0};
        int q[N+3] = {0};
        int n = (int)(L/1.39793 + 1); int k;

        w[0] = 16*5; v[0] = 4*239;

        for(k = 1; k <= n; k++) {
			// 套用公式
            div(w, 25, w);
            div(v, 239, v);
            div(v, 239, v);
            sub(w, v, q);
            div(q, 2*k-1, q);

            if(k%2) // 奇数项
                add(s, q, s);
            else	// 偶数项
                sub(s, q, s);
        }

        printf("%d.", s[0]); for(k = 1; k < N; k++)
            printf("%04d", s[k]); printf("\n");
        return 0;

    }

    void add(int *a, int *b, int *c) { int i, carry = 0;

        for(i = N+1; i >= 0; i--) {
            c[i] = a[i] + b[i] + carry; if(c[i] < 10000)
                carry = 0; else { // 进位
                c[i] = c[i] - 10000;
                carry = 1;
            }
        }
    }

    void sub(int *a, int *b, int *c) { int i, borrow = 0;
        for(i = N+1; i >= 0; i--) {
            c[i] = a[i] - b[i] - borrow; if(c[i] >= 0)
                borrow = 0; else { // 借位
                c[i] = c[i] + 10000;
                borrow = 1;
            }
        }
    }

    void div(int *a, int b, int *c) {	// b 为除数
        int i, tmp, remain = 0; for(i = 0; i <= N+1; i++) {
            tmp = a[i] + remain; c[i] = tmp / b;
            remain = (tmp % b) * 10000;
        }
    }

本文参与腾讯云自媒体分享计划，分享自作者个人站点/博客。

原始发表：2020-01-17，如有侵权请联系 cloudcommunity@tencent.com 删除

java

algorithm

本文分享自作者个人站点/博客前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

java

algorithm

登录后参与评论

0 条评论

热度