円周率1万桁の計算（ガウス=ルジャンドル法）

coding UTF-8;

import Time;
import Text;


// 円周率を求める桁数
const int DIGIT = 10000;

// 結果を記載するファイル
const string FILE_PATH = "pi.txt";

// 計算進行状況を表示するかどうか
const bool PRINT_PROGRESS = true;


/*
 * main関数:
 * プログラムの実行開始時に、まずこの関数が呼び出される
 */
void main(){

	// varfloat型（多倍長浮動小数点型）の演算桁数を設定する
	digit( DIGIT * 1.1 ); // 1割くらいマージンを設ける
	if( DIGIT <= 100 ){
		digit( DIGIT + 10 ); // 100桁以下の時は念のため+10桁
	}

	if( PRINT_PROGRESS ){
		println( "計算中…" );
		println( "------------------------------" );
	}

	// 計算開始時刻を記録
	int start = time();


	// !!!!!
	// ここで円周率を計算
	// !!!!!
	varfloat piValue = pi( PRINT_PROGRESS );


	// 結果を見やすい形で出力
	string text = getOutputText( piValue, DIGIT );
	println( "------------------------------" );
	println( text );

	// 計算所要時間（ベンチマークスコア）を表示
	if( PRINT_PROGRESS ){

		// 計算所要時間を求めて表示
		float totaltime = millisecond( time() - start ) / 1000.0;
		println( "------------------------------" );
		println( "計算完了" + EOL );
		println( "SCORE	" + totaltime + " SEC" );
	}

	// 結果をファイルに出力
	bool save = confirm( "結果をファイルに保存しますか ?" );
	if( save ){
		save( FILE_PATH ); // CONSOLEの内容をファイルに出力
		alert( FILE_PATH + " に保存しました。" );
	}
}


/*
 * 平方根 √a を計算する関数（ニュートン法）
 * 収束を速くするため、開始値として近い値をstartに与える
 *
 * ※varfloat型は、多倍長精度の浮動小数点を扱う型
 */
varfloat sqrt( varfloat a, varfloat start ){
	varfloat value = start;
	varfloat valueStock = -1.0vf;
	while( value != valueStock ){
		valueStock = value;
		value = ( value + a / value ) * 0.5vf;
	}
	return value;
}


/*
 * 円周率を計算する関数（ガウス=ルジャンドル法）
 * 引数 printProgress は、trueなら計算進行状況を表示する
 *
 * ※varfloat型は、多倍長精度の浮動小数点を扱う型
 */
varfloat pi( bool printProgress ){
	int start = time(); // 処理開始時間

	// アルゴリズムで使う変数
	varfloat value = 0.0vf;
	varfloat a = 1.0vf;
	varfloat b;
	varfloat aStock = 0.0vf;
	varfloat bStock = 0.0vf;
	varfloat c = 0.25vf;
	varfloat d = 1.0vf;

	// bの初期値を1/√2 = √0.5に設定
	// 2個目の引数の即値は、収束初期値で、結果に近い値
	b = sqrt( 0.5vf, 0.7071067811865475vf );

	// ガウス=ルジャンドル法による
	// 円周率の収束アルゴリズム
	int loopCount = 0;
	while( a != aStock && b != bStock ){

		aStock = a;
		bStock = b;
		a = ( aStock + bStock )*0.5vf;

		// 平方根の計算、収束高速化のため前ループの結果からスタート
		b = sqrt( aStock*bStock, b );

		c = c - d*( a - aStock )*( a - aStock );
		d = 2.0vf*d;

		// 計算進行状況を表示
		if( printProgress ){
			// 現在までの計算所要時間
			float sec = millisecond( time() - start ) / 1000.0;
			loopCount++;
			println( loopCount + " ループ完了	/	" + sec + " SEC" );
		}
	}

	value = (a+b)*(a+b) / (4.0vf*c);
	return value;
}




/*
 * 長い桁数(digit)の円周率を、読みやすい文字列に整形する関数
 */
string getOutputText( varfloat piValue, int digit ){
	int lineColumnN = 10; // 行あたり列数
	int columnWidth = 10; // 1列あたり文字数
	int currentIndex = 0; // 現在の文字位置を控える
	int blockLineN = 10; //1ブロックごとの行数
	int lineN = digit / (lineColumnN*columnWidth); //総行数（自動計算）

	string outputText = "3." + EOL; // 結果の出力用文字列

	string piText = (string)piValue; // 円周率の文字列

	// 先頭の "3." を除外
	piText = substring( piText, 2, (lengthOf(piText)-1) );

	for( int i=0; i<lineN; i++ ){

		// 一行の内容を書き出す
		for( int j=0; j<lineColumnN; j++ ){

			// 切り取り始点インデックスを求める
			int subindex = i*lineColumnN*columnWidth + j*columnWidth;

			// 始点インデックスから、1列あたり文字数の分だけpiTextを切り取り、outputTextに追加
			outputText += substring( piText, subindex, subindex+columnWidth );
			outputText += " ";

		}

		// 行の折り返し
		outputText += EOL;

		// ブロックの区切り行ごとに、現在の行数を表示する
		if( (i+1)%blockLineN==0 ){
			outputText += EOL + " - " + (i+1)*lineColumnN*columnWidth + " DIGIT -" + EOL + EOL;
		}
	}

	return outputText;
}

coding UTF-8;

import Time;
import Text;

// 円周率を求める桁数
const int DIGIT = 10000;

// 結果を記載するファイル
const string FILE_PATH = "pi.txt";

// 計算進行状況を表示するかどうか
const bool PRINT_PROGRESS = true;


/*
 * main関数:
 * プログラムの実行開始時に、まずこの関数が呼び出される
 */
void main(){

	// varfloat型（多倍長浮動小数点型）の演算桁数を設定する
	digit( DIGIT * 1.1 ); // 1割くらいマージンを設ける
	if( DIGIT <= 100 ){
		digit( DIGIT + 10 ); // 100桁以下の時は念のため+10桁
	}

	if( PRINT_PROGRESS ){
		println( "計算中…" );
		println( "------------------------------" );
	}

	// 計算開始時刻を記録
	int start = time();


	// !!!!!
	// ここで円周率を計算
	// !!!!!
	varfloat piValue = pi( PRINT_PROGRESS );


	// 結果を見やすい形で出力
	string text = getOutputText( piValue, DIGIT );
	println( "------------------------------" );
	println( text );

	// 計算所要時間（ベンチマークスコア）を表示
	if( PRINT_PROGRESS ){

		// 計算所要時間を求めて表示
		float totaltime = millisecond( time() - start ) / 1000.0;
		println( "------------------------------" );
		println( "計算完了" + EOL );
		println( "SCORE	" + totaltime + " SEC" );
	}

	// 結果をファイルに出力
	bool save = confirm( "結果をファイルに保存しますか ?" );
	if( save ){
		save( FILE_PATH ); // CONSOLEの内容をファイルに出力
		alert( FILE_PATH + " に保存しました。" );
	}
}

/*
 * 平方根 √a を計算する関数（ニュートン法）
 * 収束を速くするため、開始値として近い値をstartに与える
 *
 * ※varfloat型は、多倍長精度の浮動小数点を扱う型
 */
varfloat sqrt( varfloat a, varfloat start ){
	varfloat value = start;
	varfloat valueStock = -1.0vf;
	while( value != valueStock ){
		valueStock = value;
		value = ( value + a / value ) * 0.5vf;
	}
	return value;
}

/*
 * 円周率を計算する関数（ガウス=ルジャンドル法）
 * 引数 printProgress は、trueなら計算進行状況を表示する
 *
 * ※varfloat型は、多倍長精度の浮動小数点を扱う型
 */
varfloat pi( bool printProgress ){
	int start = time(); // 処理開始時間

	// アルゴリズムで使う変数
	varfloat value = 0.0vf;
	varfloat a = 1.0vf;
	varfloat b;
	varfloat aStock = 0.0vf;
	varfloat bStock = 0.0vf;
	varfloat c = 0.25vf;
	varfloat d = 1.0vf;

	// bの初期値を1/√2 = √0.5に設定
	// 2個目の引数の即値は、収束初期値で、結果に近い値
	b = sqrt( 0.5vf, 0.7071067811865475vf );

	// ガウス=ルジャンドル法による
	// 円周率の収束アルゴリズム
	int loopCount = 0;
	while( a != aStock && b != bStock ){

		aStock = a;
		bStock = b;
		a = ( aStock + bStock )*0.5vf;

		// 平方根の計算、収束高速化のため前ループの結果からスタート
		b = sqrt( aStock*bStock, b );

		c = c - d*( a - aStock )*( a - aStock );
		d = 2.0vf*d;

		// 計算進行状況を表示
		if( printProgress ){
			// 現在までの計算所要時間
			float sec = millisecond( time() - start ) / 1000.0;
			loopCount++;
			println( loopCount + " ループ完了	/	" + sec + " SEC" );
		}
	}

	value = (a+b)*(a+b) / (4.0vf*c);
	return value;
}