# 【长 PI】

/*

PI = [16/5 - 16 / (3*5^3 ) + 16 / (5*5^5) - 16 / (7*5^7) + ] ......] -
[4/239 - 4/(3*239^3) + 4/(5*239^5) - 4/(7*239^7) + ......]

PI = [16/5 - 4/239] - [16/(5^3)- 4/(239^3)]/3+  [16/(5^5)- 4/(239^5)]/5 + ......

[16/5^(2*n-1)- 4/239^(2*n-1)] / (2*n-1)

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

n = L / (2log5) = L /  1.39794

n = [L/1.39794] + 1

[W(n) -1/5^2 - V(n) - 1 / (239^2)] / (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)

*/

#include <stdio.h>

#define L 1000                //L为位数，N是array的长度
#define N L/4 + 1

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)
{
}
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)
{
int i, tmp, remain = 0;
for(i = 0; i <= N + 1; i++)
{
tmp = a[i] + remain;
c[i] = tmp / b;
remain = (tmp % b) * 10000;
}
}

posted @ 2017-01-20 10:31  天秤libra  阅读(801)  评论(0编辑  收藏