大規模なxに対するexp(x)/2の効率的な計算方法

問題点

IEEE 754 binary64（倍精度浮動小数点数）において、sinh(x)関数を実装する際に問題が発生します。xが約709.783から710.476の範囲にある場合、直接exp(x)/2を計算しようとするとオーバーフローが発生します。

倍精度浮動小数点数の最大値DBL_MAXは約1.798×10³⁰⁸であり、以下の制約があります：

• 有限のexp(x)が得られる最大のx：約709.783
• 有限のsinh(x)およびcosh(x)が得られる最大のx：約710.476

このため、exp(710.0)はオーバーフローし、後続の/2演算では有限の結果を得ることができません。

// 問題のある実装例
double problematic_sinh(double x) {
    return exp(x) / 2;  // xが大きいとオーバーフロー
}

解決策

基本アプローチ：数学的恒等式の活用

オーバーフローを回避するために、数学的恒等式を利用します：

eˣ/2 = e^(x - v) × eᵛ/2

適切に選択された定数vを使用することで、減算x - vを正確に行い、スケーリング係数eᵛ/2を事前に計算しておくことができます。

最適な解決策

方法1：1024 log(2)を利用した高精度計算（Mark Dickinson氏の解法）

double metd_sinh(double x) {
    // c1 + c2は1024 * log(2)の高精度近似
    static const double c1 = 0x1.62e42fefa39efp+9;
    static const double c2 = 0x1.acp-46;
    return exp(x - c1 - c2) * 0x1p+1023;
}

特徴：

• ULP誤差：約0.612
• 高速な計算
• 減算と乗算は正確に行われる
• エラーの主要因はexp()関数の実装精度のみ

方法2：改良版定数利用アプローチ

double optimized_sinh(double x) {
    static const double v = 0x1.62e4307c58800p-1;
    static const double scale = 0x1.000000465a709000p+1 / 2;
    return exp(x - v) * scale;
}

特徴：

• ULP誤差：約1.008
• vの選択により減算が正確
• scale定数はeᵛ/2を高精度で表現

方法3：シンプルな近似解法

double simple_sinh(double x) {
    double s = exp(x / 2);
    return 0.5 * s * s;
}

特徴：

• ULP誤差：約1.918
• 特別な定数が不要
• 実装が簡潔

注意点

0.5 / exp(-x)というアプローチもありますが、非正規化数がサポートされていない環境では精度低下や完全な値の損失が発生する可能性があります。

パフォーマンス比較

各種手法の性能比較（テスト環境：gcc 12.4.0, Cygwin）：

手法	最大ULP誤差	標準偏差	平均実行時間(μs)
metd_sinh (1024 log(2)利用)	0.612	0.289	0.0147
改良版定数利用	1.007	0.408	0.0199
シンプルな近似	1.918	0.751	0.0169
標準ライブラリsinh()	1.918	0.751	0.0234

技術的詳細

正確な減算の重要性

x - vの減算を正確に行うことが重要です。適切に選択されたvは：

1. xとの減算が正確に行える（Sterbenzの補題）
2. 数学的なeᵛとその倍精度表現が非常に近い（0.001 ULP以内）

// 減算が正確に行われるvの例
static const double v = 0x1.62e4307c58800p-1;  // 10下位ビットがゼロ

スケーリング係数の最適化

scale = eᵛ/2の定数を高精度で計算することで、最終結果の精度を大幅に向上できます。理想的なvを選択すると、eᵛが0x1.000000xxxxxxxp+1の形式になり、倍精度表現での誤差を最小化できます。

実装のベストプラクティス

高精度な定数の計算

定数の計算には多倍長演算や高精度数学ライブラリを使用します：

# Pythonでの高精度計算例
from decimal import Decimal, getcontext

getcontext().prec = 100  # 100桁の精度
v = Decimal('0.69314719694034466')
scale = (v.exp() / 2).to_eng_string()

環境適応型実装

環境に応じて最適な手法を選択できます：

double adaptive_sinh(double x) {
    if (x < 710.12928648366403) {
        // 1024 * log(2)を利用
        static const double c1 = 0x1.62e42fefa39efp+9;
        static const double c2 = 0x1.abc9e3b39803fp-46;
        return exp(x - c1 - c2) * 0x1p+1023;
    } else {
        // 1025 * log(2)を利用
        static const double c1 = 0x1.633ce8fb9f87ep+9;
        static const double c2 = -0x1.3ae594e9bd8bp-45;
        return exp(x - c1 - c2) * 0x1p+1023 * 2;
    }
}

結論

大規模なxに対するexp(x)/2の計算には、以下のいずれかの方法を推奨します：

1. 最高精度が求められる場合：1024 log(2)を利用したMark Dickinson氏の解法（ULP誤差約0.612）
2. バランスの良い選択：改良版定数利用アプローチ（ULP誤差約1.008）
3. 実装簡便性が優先の場合：シンプルな近似解法（ULP誤差約1.918）

いずれの方法も、単純なexp(x)/2の計算よりもはるかに優れた結果を提供し、オーバーフローを効果的に回避できます。

パフォーマンスチップ

ほとんどの現代的なプロセッサでは、exp()関数の計算コストが支配的です。定数回数の追加演算はパフォーマンスにほとんど影響しません。