円周率πを計算で求めるその1

2019.06.21

はじめに

円周率 πを計算で求める方法を紹介します。

ここで紹介するアルゴリズムは、単位円を小さく分割した台形で近似し、それら台形の面積の合計で円周率 πを求めるものです。

（半径×半径×円周率π)

で求められるので、単位円(半径1の円)の面積は、

（1×1×円周率π)＝円周率π

となります。

よって、単位円の面積の計算をすることで円周率 πを求めることができます。

$単位円の面積＝円周率$

ここからは、単位円の面積を計算する方法を説明していきます。

単位円の右上１／４を台形（一番上は三角形）に分割します。このとき台形の高さを全て同じにします。

下図の例は、単位円の右上１／４を２つの台形に分割した例です。

オレンジ色で塗りつぶされている部分が台形の範囲です。２つの台形に分割しているので、それぞれの台形の高さは0.5です。

$単位円の分割(y方向に２分割)$

以下の例は、単位円の右上１／４を４つの台形に分割した例です。

４つの台形に分割しているので、それぞれの台形の高さは1を4で割った0.25になります。

$単位円の分割(y方向に４分割)$

このオレンジ色で塗られた部分の面積を計算し、それを４倍した値を円周率 πとします。

分割する数を増やすとオレンジ色で塗られた部分の形が、単位円に近づいていきます。

そして、それは台形の面積の合計値が、単位円の１／４の面積に近づくことを意味しています。

Javaソースコード

以下は、円周率 πを計算するJava ソースコードです。

CalcPI1.java

ダウンロード

public class CalcPI1 {
	// πを計算するメソッド
	public static double getPI_1( int yn )
	{
		double upperbase;	// 台形の上底
		double lowerbase;	// 台形の下底
		double height;		// 台形の高さ
		double lowery;		// 台形の下底の座標
		double area;		// 円の面積
		double aread;		// 台形の面積

		// 台形の高さ(円の半径をynで割る)
		height = 1.0 /  (double)yn;

		// 上底の初期値
		upperbase = 0.0;

		// 下底のy座標の初期値
		lowery = 1.0 - height;

		// 計算開始
		area = 0.0;
		for ( int i = 0; i < yn; ++ i ) {
			// 下底の計算
			lowerbase = Math.sqrt( 1.0 - lowery * lowery );

			// 台形の面積計算
			aread = ( upperbase + lowerbase ) * height / 2.0;

			// 円の面積に台形の面積を足す
			area += aread;

			// 下底を上底に代入
			upperbase = lowerbase;

			// 下底のy座標からheightを引く
			lowery -= height;
		}

		// 1/4円の面積を4倍
		area = area * 4.0;

		// 結果を戻す
		return area;
	}


	public static void main( String[] args ) {
		// yの分割数
		int yn;

		// 1000分割
		yn = 1000;
		System.out.println( "分割数:" + yn+ "     π=" + getPI_1( yn ) );

		// 10000分割
		yn = 10000;
		System.out.println( "分割数:" + yn+ "    π=" + getPI_1( yn ) );

		// 100000分割
		yn = 100000;
		System.out.println( "分割数:" + yn+ "   π=" + getPI_1( yn ) );

		// 1000000分割
		yn = 1000000;
		System.out.println( "分割数:" + yn+ "  π=" + getPI_1( yn ) );

		// 10000000分割
		yn = 10000000;
		System.out.println( "分割数:" + yn+ " π=" + getPI_1( yn ) );
	}
}

コンパイルソースコードが「ANSI」の場合

C:\talavax\javasample>javac -encoding sjis CalcPI1.java

コンパイルソースコードが「UTF-8」の場合

C:\talavax\javasample>javac CalcPI1.java

実行

C:\talavax\javasample>java CalcPI1

出力結果

分割数:1000     π=3.1415554669110293
分割数:10000    π=3.1415914776112435
分割数:100000   π=3.141592616398606
分割数:1000000  π=3.1415926524339723
分割数:10000000 π=3.1415926531488605

分割数を大きくするとπの計算値が、実際のπの値に近づいています。

実際の円周率 πは、

３．１４１５９２６５３５８９７９３２３８４６２・・・

です。分割数を1000000にすると小数点以下第９桁まで一致しました。

Javaソースコードの解説

public class CalcPI1 {

クラス名を、CalcPI1としています。

002
003

	// πを計算するメソッド
	public static double getPI_1( int yn )

円周率 πを計算するメソッドです。int型の変数ynは分割数で、計算する台形の個数を表します。

		double upperbase;	// 台形の上底
		double lowerbase;	// 台形の下底
		double height;		// 台形の高さ
		double lowery;		// 台形の下底の座標
		double area;		// 円の面積
		double aread;		// 台形の面積

計算に使用する変数を宣言しています。

012
013

		// 台形の高さ(円の半径をynで割る)
		height = 1.0 /  (double)yn;

半径1.0を変数ynで割って台形の高さを計算し、double型の変数heightに代入しています。

015
016

		// 上底の初期値
		upperbase = 0.0;

最初に計算する台形の上底に0.0を代入しています。このプログラムでは、上から下に向けて台形の面積を計算していきます。

018
019

		// 下底のy座標の初期値
		lowery = 1.0 - height;

初期の下底のy座標を変数loweryに代入しています。上底のy座標は1.0なので、下底のy座標は1.0から高さheightを引いた値になります。

021
022
023

		// 計算開始
		area = 0.0;
		for ( int i = 0; i < yn; ++ i ) {

ここから台形の面積を、指定した分割数yn回計算していきます。そして、変数areaに台形の面積を足して行きます。変数areaの初期値は0.0です。

024
025

			// 下底の計算
			lowerbase = Math.sqrt( 1.0 - lowery * lowery );

下底のy座標loweryから、下底を計算しています。

下底はy座標loweryを持つ単位円上のプラスのx座標になるので、1²=lowery²+lowerbase²を満たすlowerbaseが下底となります。

027
028

			// 台形の面積計算
			aread = ( upperbase + lowerbase ) * height / 2.0;

上底と下底と高さで台形の面積を計算して変数areadに代入しています。

030
031

			// 円の面積に台形の面積を足す
			area += aread;

変数areaに台形の面積areadを足しています。

033
034

			// 下底を上底に代入
			upperbase = lowerbase;

下底を上底に代入しています。計算した台形の下底が、次に計算する台形の上底になるためです。

036
037

			// 下底のy座標からheightを引く
			lowery -= height;

下底はy座標loweryから高さheightを引いています。

040
041

		// 1/4円の面積を4倍
		area = area * 4.0;

for文を抜けると、変数areaには単位円の四分の一の面積が格納されています。変数areaに４を掛けて単位円全体の面積を求めます。

043
044

		// 結果を戻す
		return area;

return文で、単位円の面積areaを戻します。

	public static void main( String[] args ) {
		// yの分割数
		int yn;

		// 1000分割
		yn = 1000;
		System.out.println( "分割数:" + yn+ "     π=" + getPI_1( yn ) );

		// 10000分割
		yn = 10000;
		System.out.println( "分割数:" + yn+ "    π=" + getPI_1( yn ) );

		// 100000分割
		yn = 100000;
		System.out.println( "分割数:" + yn+ "   π=" + getPI_1( yn ) );

		// 1000000分割
		yn = 1000000;
		System.out.println( "分割数:" + yn+ "  π=" + getPI_1( yn ) );

		// 10000000分割
		yn = 10000000;
		System.out.println( "分割数:" + yn+ " π=" + getPI_1( yn ) );
	}