// +++++++++++++++++++++++++++++++++++++++++++++++++++++
// +++ 自然対数の底ｅ と 円周率π を求めるプログラム +++
// +++++++++++++++++++++++++++++++++++++++++++++++++++++
//
//	original maked by Y.S.
//	transfered(fortran to C) by G-HAL
//		1995 / 09 / 16
//	transfered(C to C++) by G-HAL
//		1995 / 09 / 17
//	transfered(DECIMAL to BINALY) by G-HAL
//		1995 / 09 / 17
//		1995 / 11 / 26	DEBUG
//		1995 / 11 / 27

#include "EPPP_C.h"
#include <stdlib.h>
#include <stdio.h>
#include <ctype.h>
#include <math.h>

//	（デフォルトの）求める精度（桁数）を KETAM とする。
signed long KETAM = 1000;

//
//	ｅを求める
//
void calE(void)
{
	unsigned long I;
	LONGNUM ANS, WORK0, WORK1;

	if( ANS.INIT(KETAM) || WORK0.INIT(KETAM) || WORK1.INIT(KETAM) )
	{
		fputs("ERROR : out of memory.", stderr);
		exit(-1);
	}

	ANS   = 1;
	WORK0 = 1;

	for( I=1; ; I++)
	{
		WORK0 /= I;
		if( WORK0 == 0 || 65534l<I)	break;		// 終了判定
		ANS += WORK0;
	}

	printf("現在の最大除数は %lu。\n", I);
	puts("現在の有効桁数は、");
	ANS.KETA(WORK0);
	puts("自然対数の底の値は、");
	ANS.HYOUZI(WORK0, WORK1);
	return;
}

//
//	円周率を求める
//
void calPI(void)
{
	signed long I, COUNT;
	LONGNUM ANS, WORK0, WORK1, WORK2;

	if( ANS.INIT(KETAM) || WORK0.INIT(KETAM) 
		|| WORK1.INIT(KETAM) || WORK2.INIT(KETAM) )
	{
		fputs("ERROR : out of memory.", stderr);
		exit(-1);
	}

//	   変数をリセット
	ANS	= 0;
	WORK0 = 0;
	WORK1 = 16;
	WORK2 = 4;
	COUNT = 1;
	I     = 3;
//												   //１項目を計算する。
//	  １項目の計算									  16/5 を計算して足す。
//													  4/239を計算して引く。
	WORK1 /= 5;
	ANS += WORK1;
	WORK2 /= 239;
	ANS -= WORK2;

//	  Ｌoop
	while(1){
//										  // WORK11 = WORK11 / 5 / 5
//											 WORK00 = WORK11 （WORK11の計算結果
//											 はまたループして使うのでWORK00にコ
//	  16 arctan (1/5) の count項目 の 計算	 ピーしてつぎの計算をする）
//											 WORK00 = WORK00 / I
//											 ANS = ANS - WORK00
		WORK1 /= 25;
		WORK0 = WORK1;
		WORK0 /= I;
		ANS -= WORK0;

//											  //　上と同様に計算
//												 WORK21 = WORK21 / 239 / 239
//	  -4 arctan (1/239) の count項目 の 計算	 WORK00 = WORK21
//												 WORK00 = WORK00 / I
//												 ANS = ANS + WORK00

		WORK2 /= 239;
		WORK2 /= 239;
		WORK0 = WORK2;
		WORK0 /= I;
		ANS += WORK0;

		I += 2;
		COUNT++;
//													 //上の計算をもう１回して
//													   からループする。
//	 上と同様に２つのarctanを計算					  （＋−が１項ずつ逆転する
//														 ため）
		WORK1 /= 25;
		WORK0 = WORK1;
		WORK0 /= I;
		ANS += WORK0;
//													//足す数が0なら終了。
//	　　  終了判定									 遅く0になるarctan(1/5)で
//													 判断している。
		if( WORK0 == 0 || 65528l < I)	break;
		WORK2 /= 239;
		WORK2 /= 239;
		WORK0 = WORK2;
		WORK0 /= I;
		ANS -= WORK0;
		I += 2;
		COUNT++;
	}

//		終了
	printf("現在の最大除数は %lu。 \n", I);
	puts("現在の有効桁数は、");
	ANS.KETA(WORK0);
	puts("円周率の値は、");
	ANS.HYOUZI(WORK0, WORK1);
	return;
}


int main( const unsigned int argc, unsigned char *argv[])
{
	if(argc <= 1 || 3 < argc)
	{
		fputs("want to get ｅ>EPPP E [length]\n", stderr);
		fputs("want to get π>EPPP P [length]\n", stderr);
		fputs("want to get max calclately ｅ>EPPP GE [最大除数]\n", stderr);
		fputs("want to get max calclately π>EPPP GP [最大除数]\n", stderr);
		exit(-1);
	}
	if(argc == 3)
	{
		KETAM  = atol(argv[2]);
		if(!KETAM)
		{
			fputs("option error.\n", stderr);
			exit(-1);
		}
	}else
		puts("default 1000.");

	argv[1][0] = toupper(argv[1][0]);
	if(argv[1][0] == 'E')
		calE();
	else if(argv[1][0] == 'P')
		calPI();
	else if(argv[1][0] == 'G'){
		double X=0;
		unsigned long I;
		argv[1][1] = toupper(argv[1][1]);
		if(argv[1][1] == 'E'){
			for(I=2; I<KETAM; I++)
				X += log10(I);
		}else if(argv[1][1] == 'P'){
			X = KETAM*log10(25) + log10(KETAM);
		}
		printf("%lf",X);
	}

	exit(0);
}