{
	pi.kn (for Kuin 1.00):
		Last Modified: 2013/09/23 10:38:07.
		(以前のコードをKuin 1.00で動作するように修正しました。)
		Created by @tatt61880

	下記のサイトを参考に作成しました。
		http://mail2.nara-edu.ac.jp/~asait/c_program/sample0/pi.htm
}
const M: int :: 2502
const N: int :: 2078

func Main()
	var x: []byte32 :: #[@N]byte32
	var y: []byte32 :: #[@N]byte32

	do @Atan(x, 5)
	do @Mul(x, 4$ byte32)

	do @Atan(y, 239)

	do @Sub(x, y)
	do @Mul(x, 4$ byte32)

	do @Display(x)
end func

func Atan(x: []byte32, m: int)
	var e: int
	var m2: byte32
	var y: []byte32 :: #[@N]byte32
	var z: []byte32 :: #[@N]byte32
	for i(0, @N - 1)
		do x[i] :: 0$ byte32
		do y[i] :: 0$ byte32
	end for
	do y[0] :: 1$ byte32

	do @Div(y, m$ byte32)
	do @Dup(y, x)
	do m2 :: (m * m)$ byte32
	var i: int :: 1
	while (e = -1, skip)
		do @Div(y, m2)
		do @Dup(y, z)
		do @Div(z, (2 * i + 1)$ byte32)
		if ((i$ byte32).And(1$ byte32) = 0$ byte32 )
			do @Add(x, z)
		else
			do @Sub(x, z)
		end if
		do i :+ 1
		do e :: @if_zero(z)
	end while
end func

func Add(a: []byte32, b: []byte32)
	var j: int
	var x: byte32
	for i(0, @N - 1)
		do x :: a[i] + b[i]
		if(x <= 16#ffff$ byte32)
			do a[i] :: x
		else
			do a[i] :: x.And(16#ffff$ byte32)
			do j :: i - 1
			while(a[j] = 16#ffff$ byte32)
				do a[j] :: 0$ byte32
				do j :- 1
			end while
			do a[j] :+ 1$ byte32
		end if
	end for
end func

func Sub(a: []byte32, b: []byte32)
	var j: int
	for i(0, @N - 1)
		if (a[i] >= b[i])
			do a[i] :: a[i] - b[i]
		else
			do a[i] :: 16#10000$ byte32 + a[i] - b[i]
			do j :: i - 1
			while(a[j] = 0$ byte32)
				do a[j] :: 16#ffff$ byte32
				do j :- 1
			end while
			do a[j] :- 1$ byte32
		end if
	end for
end func

func Div(a: []byte32, d: byte32)
	var x: byte32
	var q: byte32
	var res: byte32
	do res :: 0$ byte32
	for i(0, @N - 1)
		do res :: res.Shl(16)
		do x :: a[i] + res
		do q :: x / d
		do a[i] :: q
		do res :: x - q * d
	end for
end func

func Mul(a: []byte32, d: byte32)
	var i: int
	var x: byte32
	var q: byte32
	do q :: 0$ byte32
	for i(@N - 1, 0, -1)
		do x :: a[i] * d + q
		do a[i] :: x.And(16#ffff$ byte32)
		do q :: x.Shr(16)
	end for
end func

func Dup(a: []byte32, b: []byte32)
	for i(0, @N - 1)
		do b[i] :: a[i]
	end for
end func

func if_zero(a: []byte32): int
	var j: int
	do j :: 0
	for i(0, @N - 1)
		if (a[i] <> 0$ byte32)
			do j :: -1
			break i
		end if
	end for
	return j
end func

func Display(a: []byte32)
	var j: int
	var w: []byte32 :: #[@M]byte32
	do @to_decimal(a, w)
	var str: []char :: ""
	do str :~ (w[0]$ int).ToStrF("4d") ~ "."
	for i(1, @M - 1)
		do str :~ (w[i]$ int).ToStrF("04d") ~ " "
		do j :: i % 12
		if(j=0)
			do str :~ "\n     "
		end if
	end for
	do str :~ "\n"
	do Dbg@Log(str)
end func

func to_decimal(a: []byte32, w: []byte32)
	var b: []byte32 :: #[@N]byte32
	for i(0, @N - 1)
		do b[i] :: a[i]
	end for
	do w[0] :: b[0]
	do b[0] :: 0$ byte32
	for i(1, @M - 1)
		do @Mul(b, 10000$ byte32)
		do w[i] :: b[0]
		do b[0] :: 0$ byte32
	end for
end func
	

/*
	pi.c
		下記のサイトのコードを少し改変しました。
		http://mail2.nara-edu.ac.jp/~asait/c_program/sample0/pi.htm
*/
#include <stdio.h>
#define M 2502
#define N 2078

void Add(unsigned a[], unsigned b[]);
void Sub(unsigned a[], unsigned b[]);
void Div(unsigned a[], unsigned d);
void Mul(unsigned a[], unsigned d);
void Dup(unsigned a[], unsigned b[]);
void Atan(unsigned x[], int m);
int if_zero(unsigned a[]);
void Display(unsigned a[]);
void to_decimal(unsigned a[], unsigned w[]);

void main(){
	unsigned x[N], y[N];

	Atan(x, 5);
	Mul(x, 4);

	Atan(y, 239);

	Sub(x, y);
	Mul(x, 4);

	Display(x);
}

void Atan(unsigned x[], int m){
	int i, e;
	unsigned m2;
	unsigned y[N], z[N];
	for(i = 0; i < N; i++){
		x[i] = 0;
		y[i] = 0;
	}
	y[0] = 1;

	Div(y, m);
	Dup(y, x);
	m2 = m*m;
	i = 1;
	do {
		Div(y, m2);
		Dup(y, z);
		Div(z, 2*i+1);
		if( (i & 1) == 0 ){
			Add(x, z);
		} else {
			Sub(x, z);
		}
		i++;
		e = if_zero(z);
	} while (e == -1);
}

void Add(unsigned a[], unsigned b[]){
	int i,j;
	unsigned x;
	for(i = 0; i < N; i++){
		x = a[i] + b[i];
		if(x <= 0xffff){
			a[i] = x;
		} else {
			a[i] = x & 0xffff;
			j = i-1;
			while(a[j] == 0xffff){
				a[j] = 0;
				j--;
			}
			a[j]++;
		}
	}
}

void Sub(unsigned a[], unsigned b[]){
	int i,j;
	for(i = 0; i < N; i++){
		if (a[i] >= b[i]){
			a[i] = a[i] - b[i];
		} else {
			a[i] = 0x10000 + a[i] - b[i];
			j = i-1;
			while(a[j] == 0){
				a[j] = 0xffff;
				j--;
			}
			a[j]--;
		}
	}
}

void Div(unsigned a[], unsigned d){
	int i;
	unsigned x, q, res;
	res = 0;
	for(i = 0; i < N; i++){
		res = res << 16;
		x = a[i] + res;
		q = x/d;
		a[i] = q;
		res = x - q*d;
	}
}

void Mul(unsigned a[], unsigned d){
	int i;
	unsigned x, q;
	q = 0;
	for(i = N-1; i >= 0; i--){
		x = a[i]*d + q;
		a[i] = x & 0xffff;
		q = x >> 16;
	}
}

void Dup(unsigned a[], unsigned b[]){
	int i;
	for(i = 0; i < N; i++){
		b[i] = a[i];
	}
}

int if_zero(unsigned a[]){
	int i, j;
	j = 0;
	for(i = 0; i < N; i++){
		if (a[i] != 0) {
			j = -1;
			break;
		}
	}
	return j;
}

void Display(unsigned a[]){
	int i,j;
	unsigned w[M];
	to_decimal(a, w);
	printf("%4.1u.",w[0]);
	for(i = 1; i < M; i++){
		printf("%4.4u ", w[i]);
		j = i%12;
		if(j==0) printf("\n     ");
	}
	printf("\n");
}

void to_decimal(unsigned a[], unsigned w[]){
	int i;
	unsigned b[N];
	for(i = 0; i < N; i++){
		b[i] = a[i];
	}
	w[0] = b[0];
	b[0] = 0;
	for(i = 1; i < M; i++){
		Mul(b, 10000);
		w[i] = b[0];
		b[0] = 0;
	}
}