ガウス=ルジャンドルのアルゴリズム

未分類
ガウス=ルジャンドルのアルゴリズム; Gauss Legendre Algorithm

概要

ガウス=ルジャンドルのアルゴリズムは円周率を計算する際に用いられる反復計算である。円周率に収束するスピードは他の収束計算に比べかなり早く、1回計算するごとに有効桁が約2倍になる。

アルゴリズム

・初期値
\(\quad \displaystyle a_0 = 1, \ b_0 = \frac{1}{2} \sqrt{2}, \ t_0 = \frac{1}{4}, \ p_0 = 1 \)

・反復式
\(\quad \displaystyle a_{n+1} = \frac{a_n + b_n}{2}\)
\(\quad \displaystyle b_{n+1} = \sqrt{a_n b_n} \)
\(\quad \displaystyle t_{n+1}= t_{n} - p_{n} (a_{n}-a_{n+1})^2 \)
\(\quad \displaystyle p_{n+1} = 2 p_{n} \)

・\(\pi \)の算出
\(\quad \displaystyle \pi_{n} = \frac{(a_{n}+b_{n})^2}{4 t_{n}}\)
\(\quad \displaystyle \lim_{n \to \infty} \pi_{n} = \pi\)

計算結果

仮数部が1024bitある浮動小数点型でも7回の反復計算で有効桁数に達してしまう。
\( \begin{eqnarray} 1 : 3.&&140... \\ 2 : 3.&&14159264... \\ 3 : 3.&&1415926535897932382... \\ 4 : 3.&&14159265358979323846264338327950288419711... \\ 5 : 3.&&14159265358979323846264338327950288419716939937510 \\ &&5820974944592307816406286208998625... \\ 6 : 3.&&14159265358979323846264338327950288419716939937510 \\ &&58209749445923078164062862089986280348253421170679 \\ &&82148086513282306647093844609550582231725359408128 \\ &&481117450284102701936... \\ 7 : 3.&&14159265358979323846264338327950288419716939937510 \\ &&58209749445923078164062862089986280348253421170679 \\ &&82148086513282306647093844609550582231725359408128 \\ &&48111745028410270193852110555964462294895493038196 \\ &&44288109756659334461284756482337867831652712019091 \\ &&45648566923460348610454326648213393607260249141273 \\ &&72458702... \\ \end{eqnarray} \)

ソースコード


namespace MultiPrecision {
    class Program {
        static void Main(string[] args) {

            string pi_true =
                "3.14159265358979323846264338327950288419716939937510" +
                  "58209749445923078164062862089986280348253421170679" +
                  "82148086513282306647093844609550582231725359408128" +
                  "48111745028410270193852110555964462294895493038196" +
                  "44288109756659334461284756482337867831652712019091" +
                  "45648566923460348610454326648213393607260249141273" +
                  "72458700660631558817488152092096282925409171536436" +
                  "78925903600113305305488204665213841469519415116094";

            using(StreamWriter sw = new StreamWriter("gauss_legendre.txt")) {
                for(int i = 1; i < 8; i++) {
                    Precision1024 pi = GaussLegendreAlgorithm(i);
                    string pi_approx = $"{pi:F400}";

                    int match_index;
                    for(match_index = 0; match_index < 400; match_index++) {
                        if(pi_true[match_index] != pi_approx[match_index]) {
                            break;
                        }
                    }

                    sw.WriteLine($"{i} : {pi_approx.Substring(0, match_index + 1)}...");
                }
            }

            Console.WriteLine("END");
            Console.Read();
        }

        static Precision1024 GaussLegendreAlgorithm(int precision_level) {
            Precision1024 a = 1, b = Precision1024.Sqrt(2) / 2, t = 0.25, p = 1;

            while(precision_level > 0) {
                Precision1024 a_next = (a + b) / 2;
                Precision1024 b_next = Precision1024.Sqrt(a * b);
                Precision1024 t_next = t - p * (a - a_next) * (a - a_next);
                Precision1024 p_next = 2 * p;

                a = a_next;
                b = b_next;
                t = t_next;
                p = p_next;

                precision_level--;
            }

            return ((a + b) * (a + b)) / (4 * t);
        }
    }
}

関連項目

高精度浮動小数点