Day 160-Q5 — 2026-09-25 — 浮動小数の厳密な総和（Shewchuk の partials と math.fsum の half-even 補正）

問題

$N$ 個の倍精度浮動小数 $x_1,\dots,x_N$ が与えられる（入力の各トークンを最も近い IEEE 754 倍精度に丸めた値。Python の float(token) と同じ）。実数としての厳密な和 $\sum x_i$ を、最近接偶数丸め（round-half-even）で倍精度に 1 回だけ丸めた値 を求め、Python の repr の形式（その値に戻る最短の 10 進表記）で出力せよ。

入力形式

N
x_1
...
x_N

制約

入出力例

3
1e100
1.0
-1e100

1.0

3
1.0
1.1102230246251565e-16
1.232595164407831e-32

1.0000000000000002

（例1: 左から足すと $10^{100}+1$ で 1 が消え 0.0 になる。例2: $x_2=2^{-53}$ は 1.0 の ulp のちょうど半分、$x_3=2^{-106}$。厳密な和は「ちょうど半分」をわずかに超えるので 1 つ上の $1+2^{-52}$ に丸めるのが正しい。左から足すと $1+2^{-53}$ が同点の偶数丸めで 1.0 になり、その後 $2^{-106}$ は消える。**Python 3.12 以降の組み込み sum() は Neumaier 補正で例1 は 1.0 を返すが、例2 は 1.0（誤り）を返す**——実測）

本問は Shewchuk「Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates」（Discrete & Computational Geometry 18, 1997）の展開（expansion）を使った総和で、CPython の math.fsum の実装（Modules/mathmodule.c）と同じ手順である。Python 3.12 の変更点は What's New（https://docs.python.org/3/whatsnew/3.12.html 、2026-09-25 参照）に「sum() now uses Neumaier summation」とある。

概念図：partials と half-even の同点

ヒント（クリックで段階的に開示）

for x in xs:
    i = 0
    for y in partials:
        if abs(x) < abs(y): x, y = y, x
        hi = x + y; lo = y - (hi - x)
        if lo: partials[i] = lo; i += 1
        x = hi
    partials[i:] = [x]
# 最後: 上から hi に足していき lo != 0 で止める
# lo がちょうど半 ulp で、次の partial が lo と同符号なら hi を 1 ulp 動かす

模範解答 (Python)

📋 コピーimport sys


def exact_sum(xs):
    partials = []                       # 重なりのない double の列（絶対値の昇順）。総和が厳密に保たれる
    for x in xs:
        i = 0
        for y in partials:
            if abs(x) < abs(y):
                x, y = y, x
            hi = x + y                  # 丸めた和
            lo = y - (hi - x)           # 丸めで失った分（Fast2Sum: |x|>=|y| なら厳密）
            if lo:
                partials[i] = lo
                i += 1
            x = hi
        partials[i:] = [x]
    # 上から足し、初めて不正確になったところで止める
    n = len(partials)
    hi = 0.0
    if n:
        n -= 1
        hi = partials[n]
        lo = 0.0
        while n:
            x = hi
            n -= 1
            y = partials[n]
            hi = x + y
            lo = y - (hi - x)
            if lo:
                break
        # 丁度中間（half-even の同点）に見えて、残りの partial が同じ向きなら切り上げる
        if n and ((lo < 0 and partials[n - 1] < 0) or (lo > 0 and partials[n - 1] > 0)):
            y = lo * 2.0
            x = hi + y
            if y == x - hi:
                hi = x
    return hi + 0.0                     # 厳密な和が 0 のとき -0.0 を 0.0 にそろえる（実数の 0）


def main():
    data = sys.stdin.buffer.read().split()
    n = int(data[0])
    xs = [float(t) for t in data[1:1 + n]]
    print(repr(exact_sum(xs)))


main()

Step-by-Step 解説

迷った他案と採らなかった理由

前提が変わったら: 集計（金額など）で「正しく丸めた 1 つの値」が要るなら math.fsum を使う（自作しない）。GPU／並列の総和で再現性（実行順に依らない結果）が要るなら、厳密和か固定小数点のアキュムレータが必要——補正付き総和は順序で結果が変わる。

よくあるミス

次のステップ

• 発展1: 同じ partials を使って、厳密な和の 符号 だけを高速に判定する（Shewchuk の幾何述語 orient2d の考え方）
• 発展2: $\sum x_iy_i$（内積）の厳密計算。積の誤差は FMA（math.fma, Python 3.13+）で $x y - \mathrm{fl}(xy)$ として取り出せる（TwoProduct）
• 扱わなかった範囲: 途中のオーバーフロー（$|x_i|$ が $10^{308}$ 近くで $hi$ が inf になる。math.fsum は特別扱いしている）、inf／nan の入力、並列化

自己評価

自分の回答

気づき・メモ