（转）圆周率PI的高精度计算(C/C++)_c++计算派高精度-CSDN博客

牛人实在看不懂。。。有人讲解一下吗？

#include <stdlib.h>
#include <stdio.h>

#define BITS 10000
int a = 10000, b, c = BITS * 7 / 2, d, e, f[BITS * 7 / 2 + 1], g;

int 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 );
    //getchar();
    return 0;
}

某次碰到pi，想用编程打印出它的比较多的有效位（至少比背的要多）。

开始考虑到 pi/4 = arctan(1)

arctan(x)展成多项式 arctan(x) = (1/1!)x - (1/3)(x^3) + (1/5)(x^5) - ....

所以有 pi/4 = 1 - 1/3 + 1/5 - 1/7 + .....

但是上式后面的式子收敛太慢了，编程很难求到很多的有效位，

而后查到Machin公式 pi/4 = 4arctan(1/5) - arctan(1/239)，这个公式可以自己证明一下（忘的差不多了，涂了半天才算对）

这样泰勒展开式收敛很快，我就是照着这个公式编程的，代码如下

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

#define KS_MIN(a, b) ((a)<(b)?(a):(b))
#define KS_MAX(a, b) ((a)>(b)?(a):(b))

#define BITS 10000

typedef struct BigInt
{
    int L;
    int* d;
    BigInt()
    {
        L = 1;
        d = ( int* )malloc( ( BITS + 16 ) * sizeof( int ) );
        d[0] = 0;
    }
    ~BigInt()
    {
        free( d );
    }

} BigInt;

// out = in / x;
void division( BigInt& out, BigInt& in, int x )
{
    out.L = in.L;
    memset( out.d, 0, out.L * sizeof( int ) );
    int i, s = 0;
    for ( i = in.L - 1; i >= 0; i-- )
    {
        s = s * 10 + in.d[i];
        out.d[i] = s / x;
        s %= x;
    }
    for ( i = out.L - 1; i >= 0; i-- )
    {
        if ( out.d[i] ) break;
    }
    out.L = i + 1;
    if ( out.L == 0 ) out.L++;
}

// a += b;
void addequal( BigInt& a, BigInt& b )
{
    int len = KS_MAX( a.L, b.L ) + 1;
    memset( a.d + a.L, 0, ( len - a.L )*sizeof( int ) );
    a.L = len;
    int i;
    for ( i = 0; i < b.L; i++ )
        a.d[i] += b.d[i];
    for ( i = 0; i < a.L - 1; i++ )
    {
        if ( a.d[i] > 9 )
        {
            a.d[i + 1]++;
            a.d[i] -= 10;
        }
    }
    for ( i = a.L - 1; i >= 0; i-- )
    {
        if ( a.d[i] ) break;
    }
    a.L = i + 1;
    if ( a.L == 0 ) a.L++;
}

// a -= b;
void subequal( BigInt& a, BigInt& b )
{
    int i;
    for ( i = 0; i < b.L; i++ )
        a.d[i] -= b.d[i];
    for ( i = 0; i < a.L - 1; i++ )
    {
        if ( a.d[i] < 0 )
        {
            a.d[i + 1]--;
            a.d[i] += 10;
        }
    }
    for ( i = a.L - 1; i >= 0; i-- )
    {
        if ( a.d[i] ) break;
    }
    a.L = i + 1;
    if ( a.L == 0 ) a.L++;
}

// a *= k
void mulequal( BigInt& a, int k )
{
    int i;
    for ( i = 0; i < a.L; i++ )
        a.d[i] *= k;
    memset( a.d + a.L, 0, 12 * sizeof( int ) );
    a.L += 12;
    for ( i = 0; i < a.L - 1; i++ )
    {
        if ( a.d[i] > 9 )
        {
            a.d[i + 1] += a.d[i] / 10;
            a.d[i] %= 10;
        }
    }
    for ( i = a.L - 1; i >= 0; i-- )
    {
        if ( a.d[i] ) break;
    }
    a.L = i + 1;
    if ( a.L == 0 ) a.L++;
}

void BigIntCopy( BigInt& dst, BigInt& src )
{
    dst.L = src.L;
    memcpy( dst.d, src.d, dst.L * sizeof( int ) );
}

void argtanx( BigInt& out, int x, int digit )
{
    BigInt a, b;
    a.L = digit + 1;
    memset( a.d, 0, ( a.L - 1 )*sizeof( int ) );
    a.d[a.L - 1] = 1;
    division( b, a, x );
    BigIntCopy( a, b );
    BigIntCopy( out, b );
    int i = 3, xx = x * x, p = -1;
    while ( 1 )
    {
        division( b, a, xx );
        BigIntCopy( a, b );
        division( b, a, i );
        if ( b.L == 1 && b.d[0] == 0 )
            break;
        if ( p < 0 ) subequal( out, b );
        else addequal( out, b );
        p *= -1;
        i += 2;
    }
}

void print( BigInt& a )
{
    printf( "%d.", a.d[a.L - 1] );
    int i;
    int low = a.L > BITS ? a.L - BITS : 0; // 不打算全部打印出来
    for ( i = a.L - 2; i >= low; i-- )
        printf( "%d", a.d[i] );
    printf( "\n" );
}

int main( int argc, char* argv[] )
{
    BigInt a, b;
    argtanx( a, 5, BITS );
    argtanx( b, 239, BITS );
    mulequal( a, 16 );
    mulequal( b, 4 );
    subequal( a, b );
    print( a );
    return 0;
}

///* pi/4 = 4*arctan(1/5) - arctan(1/239) */

https://blog.csdn.net/himulakensin/article/details/8478224

测试，第一段代码最准确，第二段代码最后2、3个数有误！反正我都看不懂。唉。。。