Float デモ — 多倍長浮動小数点の威力

Float は仮数部が任意精度の Int で表現された浮動小数点数型である。 double の約 15 桁の壁を超え、数十万桁以上の精度で計算できる。

このページでは 5 つのデモを通じて、多倍長浮動小数点ならではの威力を紹介する。以下の出力はすべて calx で実際に計算したものである。ソースコードはページ末尾に掲載しており、 calx をビルドすれば手元で同じ結果を確認できる。

Demo 1: 数学定数 100 桁

Float は主要な数学定数の計算関数を備えている。 7 つの定数を 100 桁精度で計算し、所要時間を計測する。

sqrt(2) = 1.4142135623730950488016887242096980785696718753769480731766 797379907324784621070388503875343276415727 (81 us) ln 2 = 0.6931471805599453094172321214581765680755001343602552541206 800094933936219696947156058633269964186875 (71 us) ln 10 = 2.3025850929940456840179914546843642076011014886287729760333 279009675726096773524802359972050895982983 (42 us) e = 2.7182818284590452353602874713526624977572470936999595749669 676277240766303535475945713821785251664274 (15 us) gamma = 0.5772156649015328606065120900824024310421593359399235988057 672348848677271199468391476858899497227991 (206 us) G = 0.9159655941772190150546035149323841107741493742816721342664 981196217630197762547694793565129261151062 (187 us) pi = 3.1415926535897932384626433832795028841971693993751058209749 445923078164062862089986280348253421170680 (601 us)

上記の計算時間はすべて初回計算の値である。各定数は thread_local キャッシュされるため、同一精度での 2 回目以降の呼び出しは即座に返る。 G はカタラン定数 $G = \sum_{k=0}^{\infty} \frac{(-1)^k}{(2k+1)^2}$ である。

Demo 2: π の大量桁計算 (Chudnovsky)

Chudnovsky の公式は 1 項あたり約 14 桁を生成し、 binary splitting と NTT 乗算を組み合わせることで大量桁の π を高速に計算できる。

pi (1,000 桁, 先頭 200 桁): 3.1415926535897932384626433832795028841971693993751058209749 445923078164062862089986280348253421170679821480865132823066 470938446095505822317253594081284811174502841027019385211055 5964462294895493038196 ... 1000 桁: 146 us 10000 桁: 1 ms 100000 桁: 21 ms 1000000 桁: 419 ms

100 万桁がわずか 0.4 秒。Chudnovsky + binary splitting + NTT 乗算の組み合わせにより、桁数が 10 倍になっても計算時間は約 20 倍程度に抑えられる。

Demo 3: Ramanujan の「ほぼ整数」

$e^{\pi\sqrt{163}}$ は整数に極めて近い — その差はわずか $7.5 \times 10^{-13}$。この驚くべき現象は Heegner 数と呼ばれる特殊な整数で起こり、 j-不変量と虚数乗法 (CM) 理論で説明される。

e^(pi * sqrt(163)) の高精度計算: pi = 3.141592653589793238462643383280... sqrt(163)= 12.767145334803704661710952009781... e^(pi * sqrt(163)) = 262537412640768743.9999999999992500725971981856888793538563 他の Heegner 数による類似の現象: e^(pi*sqrt(19 )) = 885479.77768015431949753789 e^(pi*sqrt(43 )) = 884736743.99977746603490666194 e^(pi*sqrt(67 )) = 147197952743.99999866245422450683 e^(pi*sqrt(163)) = 262537412640768743.99999999999925007260 d が大きいほど整数に近づくことがわかる。

$$e^{\pi\sqrt{163}} \approx 262537412640768744 - 7.5 \times 10^{-13}$$

Heegner 数は類数 1 の虚二次体 $\mathbb{Q}(\sqrt{-d})$ の判別式に対応する。 $d = 1, 2, 3, 7, 11, 19, 43, 67, 163$ の 9 個しか存在しない (Stark-Heegner の定理)。

Demo 4: 浮動小数点精度の限界を突破する

double は約 15 桁の有効数字しか持たない。 桁落ち、打ち切り誤差、恒等式の検証の 3 つの例で Float の精度を実証する。

桁落ち (catastrophic cancellation)

--- 桁落ち (catastrophic cancellation) --- x = 1e-15 (1 + x) - 1: double: 1.1102230246251565404e-15 Float: 0.000000000000001000000000000000 厳密値: 0.000000000000001

Basel 問題

$$\sum_{k=1}^{\infty} \frac{1}{k^2} = \frac{\pi^2}{6}$$ 1735 年にオイラーが解決した有名な問題である。有限和からの収束は遅い。

--- Basel 問題: Σ(1/k²) → pi²/6 --- pi^2 / 6 (厳密) = 1.644934066848226436472415166646 Σ(1/k², k=1..100000) [double] = 1.64492406689824 Σ(1/k², k=1..1000) [Float] = 1.643934566681559803139058023822 (有限和なので厳密値より小さい — 収束の遅さがわかる)

恒等式の検証: $\exp(\log(x)) = x$

--- 恒等式の検証: exp(log(x)) = x --- x = 3.14159265358979323846264338327950288419716939937510 exp(log(x)) = 3.14159265358979323846264338327950288419716939937510 |x - exp(log(x))| = 0.00000000000000000000000000000000000000000000000000 (ここでは 50 桁精度で計算したので、内部誤差は 10^-50 以下)

double では $(1 + 10^{-15}) - 1$ の結果に 11% の誤差が生じる。 Float は任意の精度で正確な値を返す。

Demo 5: Machin 型公式で π を独立に検算

2 つの異なる arctan 公式で π を計算し、Chudnovsky の結果と 100 桁で完全一致することを確認する。独立したアルゴリズムが同じ値を返すことは、実装の正しさの強力な証拠である。

Machin (1706): $\dfrac{\pi}{4} = 4\arctan\dfrac{1}{5} - \arctan\dfrac{1}{239}$

Machin の公式: pi/4 = 4*arctan(1/5) - arctan(1/239) Machin: 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170680 Chudnovsky: 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170680 差: 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 (Machin: 41 us, Chudnovsky: 28 us)

Størmer: $\dfrac{\pi}{4} = 6\arctan\dfrac{1}{8} + 2\arctan\dfrac{1}{57} + \arctan\dfrac{1}{239}$

Stormer の公式: pi/4 = 6*atan(1/8) + 2*atan(1/57) + atan(1/239) Stormer: 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170680 Chudnovsky: 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170680 差: 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 (Stormer: 51 us)

完全に独立した 3 つのアルゴリズム (Chudnovsky, Machin, Størmer) が 100 桁で完全一致する。calx の exp, log, atan, pi がすべて正しく動作している証拠である。

ソースコードと実行方法

このページの全出力は以下の C++ プログラムから生成された。 calx をビルドすれば、手元で同じ結果を確認できる。

example_float_demo.cpp (クリックで展開)

// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// example_float_demo.cpp
// 多倍長浮動小数点デモ: Float の威力を5つのシナリオで実演

#include <math/core/mp/Float.hpp>
#include <math/core/mp/Float/FloatMath.hpp>
#include <iostream>
#include <iomanip>
#include <string>
#include <chrono>
#include <cmath>

using namespace calx;

// ============================================================================
// Demo 1: 数学定数 100 桁
//   — π, e, γ, log2, √2, Catalan を高精度で計算
// ============================================================================
static void demo_constants() {
    std::cout << "============================================================\n";
    std::cout << " Demo 1: Mathematical Constants — 100 Digits\n";
    std::cout << "============================================================\n\n";

    int prec = 110;  // 余裕を持って計算
    int show = 100;

    struct {
        const char* name;
        const char* symbol;
        Float (*func)(int);
    } constants[] = {
        {"pi",      "pi",       Float::pi},
        {"e",       "e",        Float::e},
        {"gamma",   "gamma",    Float::euler},
        {"log 2",   "ln 2",     Float::log2},
        {"log 10",  "ln 10",    Float::log10},
        {"sqrt 2",  "sqrt(2)",  Float::sqrt2},
        {"catalan", "G",        Float::catalan},
    };

    for (auto& [name, symbol, func] : constants) {
        auto start = std::chrono::high_resolution_clock::now();
        Float val = func(prec);
        auto end = std::chrono::high_resolution_clock::now();
        auto us = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();

        std::string s = val.toDecimalString(show);
        // 80桁ごとに改行
        std::cout << "  " << std::setw(8) << std::left << symbol << " = ";
        for (size_t i = 0; i < s.size(); i += 60) {
            if (i > 0) std::cout << std::string(13, ' ');
            size_t len = std::min<size_t>(60, s.size() - i);
            std::cout << s.substr(i, len) << "\n";
        }
        std::cout << std::string(13, ' ') << "(" << us << " us)\n\n";
    }
}

// ============================================================================
// Demo 2: π 10,000 桁 (Chudnovsky)
//   — 計算速度の実力を示す
// ============================================================================
static void demo_pi_digits() {
    std::cout << "============================================================\n";
    std::cout << " Demo 2: Pi — 1,000 and 10,000 Digits (Chudnovsky)\n";
    std::cout << "============================================================\n\n";

    // 桁数を変えて計算速度を比較
    int benchmarks[] = {1000, 10000, 100000, 1000000};
    for (int digits : benchmarks) {
        auto start = std::chrono::high_resolution_clock::now();
        Float pi = Float::pi(digits);
        auto end = std::chrono::high_resolution_clock::now();
        auto ms = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();

        if (digits == 1000) {
            // 1000 桁は先頭200桁を表示
            std::string s = pi.toDecimalString(200);
            std::cout << "  pi (1,000 桁, 先頭 200 桁):\n\n";
            for (size_t i = 0; i < s.size(); i += 60) {
                size_t len = std::min<size_t>(60, s.size() - i);
                std::cout << "    " << s.substr(i, len) << "\n";
            }
            std::cout << "    ...\n";
        }

        if (ms >= 1000) {
            std::cout << "  " << std::setw(10) << digits << " 桁: "
                      << (ms / 1000) << " ms\n";
        } else {
            std::cout << "  " << std::setw(10) << digits << " 桁: "
                      << ms << " us\n";
        }
    }

    std::cout << std::endl;
}

// ============================================================================
// Demo 3: Ramanujan の「ほぼ整数」
//   — e^(π√163) は整数に極めて近い
// ============================================================================
static void demo_ramanujan() {
    std::cout << "============================================================\n";
    std::cout << " Demo 3: Ramanujan's \"Almost Integer\"\n";
    std::cout << "============================================================\n\n";

    int prec = 60;

    // e^(π√163) ≈ 262537412640768744 (整数に極めて近い)
    Float pi_val = Float::pi(prec);
    Float sqrt163 = sqrt(Float(163), prec);
    Float result = exp(pi_val * sqrt163, prec);

    std::string s = result.toDecimalString(40);

    std::cout << "  e^(pi * sqrt(163)) の高精度計算:\n\n";
    std::cout << "    pi       = " << pi_val.toDecimalString(30) << "...\n";
    std::cout << "    sqrt(163)= " << sqrt163.toDecimalString(30) << "...\n";
    std::cout << "\n    e^(pi * sqrt(163)) =\n";
    std::cout << "      " << s << "\n";

    std::cout << "\n  小数部分は 0.9999999999...75 で、\n";
    std::cout << "  整数 262537412640768744 との差は約 7.5 * 10^-13。\n";
    std::cout << "  この驚くべき現象は j-不変量と虚数乗法の理論で説明される。\n";

    // 他の「ほぼ整数」も計算
    std::cout << "\n  他の Heegner 数による類似の現象:\n";
    for (int d : {19, 43, 67, 163}) {
        Float sq = sqrt(Float(d), prec);
        Float val = exp(pi_val * sq, prec);
        std::string vs = val.toDecimalString(20);
        std::cout << "    e^(pi*sqrt(" << std::setw(3) << d << ")) = " << vs << "\n";
    }

    std::cout << "\n  d が大きいほど整数に近づくことがわかる。\n";
    std::cout << std::endl;
}

// ============================================================================
// Demo 4: 浮動小数点精度の限界を突破する
//   — double が破綻する計算を多倍長で正確に
// ============================================================================
static void demo_precision_limit() {
    std::cout << "============================================================\n";
    std::cout << " Demo 4: Breaking the Precision Barrier\n";
    std::cout << "============================================================\n\n";

    // 桁落ち: (1 + x) - 1  (x が極めて小さいとき)
    std::cout << "  --- 桁落ち (catastrophic cancellation) ---\n\n";
    {
        // x = 1e-15
        double xd = 1e-15;
        double rd = (1.0 + xd) - 1.0;

        int prec = 60;
        Float one(1);
        one.setPrecision(prec);
        Float xf = pow(Float(10), -15, prec);
        Float rf = (one + xf) - one;

        std::cout << "    x = 1e-15\n";
        std::cout << "    (1 + x) - 1:\n";
        std::cout << std::setprecision(20);
        std::cout << "      double: " << rd << "\n";
        std::cout << "      Float:  " << rf.toDecimalString(30) << "\n";
        std::cout << "      厳密値: 0.000000000000001\n\n";
    }

    // Basel 問題: Σ(1/k², k=1..N) → π²/6
    std::cout << "  --- Basel 問題: Σ(1/k²) → pi²/6 ---\n\n";
    {
        int prec = 50;
        Float pi_val = Float::pi(prec);
        Float exact = sqr(pi_val, prec) / Float(6);

        // double: N=100000 で加算
        double sum_d = 0.0;
        for (int k = 1; k <= 100000; ++k)
            sum_d += 1.0 / (static_cast<double>(k) * k);

        // Float: N=100000 で加算
        Float sum_f = Float::zero();
        for (int k = 1; k <= 1000; ++k) {
            Float kf(k);
            sum_f = sum_f + Float(1) / (kf * kf);
        }

        std::cout << "    pi^2 / 6 (厳密) = " << exact.toDecimalString(30) << "\n";
        std::cout << std::setprecision(15);
        std::cout << "    Σ(1/k², k=1..100000) [double] = " << sum_d << "\n";
        std::cout << "    Σ(1/k², k=1..1000)   [Float]  = " << sum_f.toDecimalString(30) << "\n";
        std::cout << "    (有限和なので厳密値より小さい — 収束の遅さがわかる)\n";
    }

    // e^(ln x) = x の高精度検証
    std::cout << "\n  --- 恒等式の検証: exp(log(x)) = x ---\n\n";
    {
        int prec = 100;
        Float x("3.14159265358979323846264338327950288419716939937510");
        x.setPrecision(prec);

        Float lx = log(x, prec);
        Float elx = exp(lx, prec);

        std::cout << "    x              = " << x.toDecimalString(50) << "\n";
        std::cout << "    exp(log(x))    = " << elx.toDecimalString(50) << "\n";

        Float diff = abs(x - elx);
        std::cout << "    |x - exp(log(x))| = " << diff.toDecimalString(50) << "\n";
        std::cout << "    (ここでは 50 桁精度で計算したので、内部誤差は 10^-50 以下)\n";
    }

    std::cout << std::endl;
}

// ============================================================================
// Demo 5: Machin 型公式で π を独立に検算
//   — π/4 = 4·arctan(1/5) - arctan(1/239)
// ============================================================================
static void demo_machin() {
    std::cout << "============================================================\n";
    std::cout << " Demo 5: Machin's Formula — Independent Verification of Pi\n";
    std::cout << "============================================================\n\n";

    int prec = 150;
    int show = 100;

    // 除算前に精度を設定（デフォルト精度で割ると桁が失われる）
    auto makeRecip = [&](int d) {
        Float num(1); num.setPrecision(prec);
        Float den(d); den.setPrecision(prec);
        return num / den;
    };

    // Machin の公式: π/4 = 4·arctan(1/5) - arctan(1/239)
    std::cout << "  Machin の公式: pi/4 = 4*arctan(1/5) - arctan(1/239)\n\n";
    {
        Float one_5th = makeRecip(5);
        Float one_239th = makeRecip(239);

        auto start = std::chrono::high_resolution_clock::now();
        Float atan_1_5 = atan(one_5th, prec);
        Float atan_1_239 = atan(one_239th, prec);
        Float pi_machin = (atan_1_5 * Float(4) - atan_1_239) * Float(4);
        auto end = std::chrono::high_resolution_clock::now();
        auto us = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();

        auto start2 = std::chrono::high_resolution_clock::now();
        Float pi_chud = Float::pi(prec + 10);  // キャッシュ回避
        auto end2 = std::chrono::high_resolution_clock::now();
        auto us2 = std::chrono::duration_cast<std::chrono::microseconds>(end2 - start2).count();

        std::cout << "    Machin:     " << pi_machin.toDecimalString(show) << "\n";
        std::cout << "    Chudnovsky: " << pi_chud.toDecimalString(show) << "\n";

        Float diff = abs(pi_machin - pi_chud);
        std::cout << "\n    差: " << diff.toDecimalString(show) << "\n";
        std::cout << "    (Machin: " << us << " us, Chudnovsky: " << us2 << " us)\n";
    }

    // Størmer の公式: π/4 = 6·arctan(1/8) + 2·arctan(1/57) + arctan(1/239)
    std::cout << "\n  Stormer の公式: pi/4 = 6*atan(1/8) + 2*atan(1/57) + atan(1/239)\n\n";
    {
        Float a8 = atan(makeRecip(8), prec);
        Float a57 = atan(makeRecip(57), prec);
        Float a239 = atan(makeRecip(239), prec);

        Float pi_stormer = (a8 * Float(6) + a57 * Float(2) + a239) * Float(4);
        Float pi_chud = Float::pi(prec);

        std::cout << "    Stormer:    " << pi_stormer.toDecimalString(show) << "\n";
        std::cout << "    Chudnovsky: " << pi_chud.toDecimalString(show) << "\n";

        Float diff = abs(pi_stormer - pi_chud);
        std::cout << "\n    差: " << diff.toDecimalString(show) << "\n";
    }

    std::cout << "\n  異なるアルゴリズムで同一の結果が得られることは、\n";
    std::cout << "  calx の多倍長演算 (exp, log, atan) の正しさを示す強力な証拠である。\n";
    std::cout << std::endl;
}

// ============================================================================
// main
// ============================================================================
int main() {
    std::cout << "###########################################################\n";
    std::cout << "#  calx Float Demo — 多倍長浮動小数点の威力              #\n";
    std::cout << "###########################################################\n\n";

    demo_constants();
    demo_pi_digits();
    demo_ramanujan();
    demo_precision_limit();
    demo_machin();

    std::cout << "###########################################################\n";
    std::cout << "#  全デモ完了                                             #\n";
    std::cout << "###########################################################\n";

    return 0;
}

API の詳細は Float API リファレンスを参照のこと。

ビルドと実行

cd calx
mkdir build && cd build
cmake .. -G "Visual Studio 17 2022" -A x64
cmake --build . --config Release --target example-float-demo
examples\Release\example-float-demo.exe

計測環境

このページの計算時間は以下の環境で計測した参考値である。

CPU: AMD Ryzen Threadripper PRO 5995WX (Zen 3)
OS: Windows 11 Pro
コンパイラ: MSVC 17.x (Visual Studio 2022), Release x64