読者です 読者をやめる 読者になる 読者になる

altebute.hatenablog.com

犬も歩けば規格にあたる

Boost.Multiprecision と BBP formula を用いて整数演算のみで円周率を近似する。

前回の記事で boost を導入したので、 Boost.Multiprecision と BBP formula を用いて整数演算のみで円周率を近似してみた。

整数演算のみなので、如何なる実行環境でも同じ結果が得られる筈。
cpp_rationalcpp_intが内部的に浮動小数点演算を使用していたら同じ結果は得られないが、そんなことは無いだろう。たぶん。

#include <iostream>
#include <boost/multiprecision/cpp_int.hpp>

boost::multiprecision::cpp_int power(int m, unsigned int k)
{
    boost::multiprecision::cpp_int result = 1;

    while (k--)
    {
        result *= m;
    }

    return result;
}

boost::multiprecision::cpp_rational bbp_formula(unsigned int precision)
{
    using cpp_rational = boost::multiprecision::cpp_rational;

    auto f = [](int k) -> cpp_rational
    {
        cpp_rational a(4, 8 * k + 1);
        cpp_rational b(2, 8 * k + 4);
        cpp_rational c(1, 8 * k + 5);
        cpp_rational d(1, 8 * k + 6);

        return (a - b - c - d) / power(16, k);
    };

    cpp_rational result = 0;

    for (unsigned int i = 0; i < precision; ++i)
    {
        result += f(i);
    }

    return result;
}

int main()
{
    using namespace std;

    for (int i = 0; i < 10; ++i)
    {
        auto result = bbp_formula(i + 1);
        cout << fixed << setprecision( 16 ) << result.convert_to< double >() << "(" << result << ")" << endl;
    }
}
3.1333333333333333(47/15)
3.1414224664224664(102913/32760)
3.1415873903465816(615863723/196035840)
3.1415924575674357(357201535487/113700787200)
3.1415926454603360(16071212445820879/5115625817702400)
3.1415926532280878(40413742330349316707/12864093722915635200)
3.1415926535728809(4318127540987083098959311/1374502686106089789849600)
3.1415926535889724(16331158360096799798177512637/5198369158853231585211187200)
3.1415926535897523(2090388270092909093421859664191/665391252333213642907031961600)
3.1415926535897913(64251934196540737654784844866951/20452025861189303550405613977600)