クヌースの筆算（TAOCP Vol II, Algorithm 4.3.1D）における数十年前のバグ

筆算の物語 2026年8月13日 クヌースの『The Art of Computer Programming』にある定理を、アルゴリズム D の数十年前のバグ（+ llvm の「バグ」）を発見した後に入手しました。 私はクヌースの『The Art of Computer Programming』で有名な筆算アルゴリズム、Algorithm D を実装していて、どうしても見過ごせない問題に遭遇しました。アルゴリズムの正当性は定理 B に依存していましたが、その証明が私を悩ませました。不自然に感じられ、単純な命題を証明するために非常に複雑な経路をたどり、問題とは無関係に見える、コーナーケースではない特別なケースを孤立させていました。何か奇妙な点があったので、自分で定理を証明しようとしましたが、失敗しました。しかし、その失敗は、数十年間正しいとされてきた Algorithm D の反例を私に与え、それとともに私の名前を冠したアルゴリズムの正当性に関する定理をもたらしました。この記事では、背景を説明し、必要であれば筆算を一から解説し、バグがどのように発生し、なぜそれほど長く隠されていたのかについての私の考えを共有し、最後に筆算を実装するより現代的な方法のプレビューで締めくくります。このブログを書く過程で、llvm のこのアルゴリズムの実装にも「バグ」を見つけ、それについても詳しく説明します。バグだけに興味がある場合は、「The bug」に直接ジャンプできます。 目次 どうしてこうなったのか 筆算を一から解説 ハードウェアにおける多倍長整数 筆算を中規模除算に削減 中規模除算を小規模除算に削減 正規化 バグ なぜ数十年間隠されていたのか 悪用できるか llvm の「バグ」 AI は見つけられなかった チェック 少しだけトリット（trit）を増やす より強い境界 商の肢（limb）を倍増させる 定数による除算 どうしてこうなったのか インタビューの準備として、小さなプロジェクトを行うことにしました。素体上の算術演算のための小さなライブラリを構築することです。これには、固定サイズの多倍長整数、算術演算、いくつかの体演算、定数時間、定メモリアクセスが含まれ、現代の暗号プロトコルへの出発点となります。実装を進めるうちに、プロセスは 1 つのルールを持つゲームになりました。それは、あらゆるコストをかけて除算を避けることです。ほとんどのことは除算なしで達成できますし、私の知る限り、暗号ライブラリは実行時に除算命令を実行しません。除算が必要な場合は、通常、乗算といくつかのシャッフリングに置き換えられます。 なぜ私たちは除算を避けるためにこれほど手間をかけるのでしょうか？乗算は非常に単純で、自然数の公理と考えることさえできます。除算ははるかに複雑です。まず、どこでも定義されているわけではありません。ゼロで割ることはできません。しかし、5 を 2 で割ることもできません！実際には「剰余付き除算」があり、これは 2 つの答えを返す、より複雑な操作です。商と剰余、つまり被除数と除数の倍数との間の最小の非負の差です。「最小」が具体的に何を意味するのか、なぜこの特定の定義を選択するのか、そしてなぜサイズの概念がそもそも関係してくるのかという問題が隠されています。別の選択をすることもできます。例えば、ゼロ中心の剰余です。しかし、さらに進んで、異なるサイズ関数を選択することもできます。これにより、まったく異なる除算アルゴリズムが得られます。 これらすべてを念頭に置くと、理論的な複雑さが実践に移行するのも不思議ではありません。現代のマシンでは、乗算命令は 1〜2 サイクルのコストで完全にパイプライン処理できますが、除算は最大 20 サイクルのコストがかかり、通常はパイプライン処理されません。最終的に、多倍長実装で残された唯一のギャップは、多倍長除算アルゴリズムでした。そこで、私は当然のこととして、筆算を実装することにしました。参考資料として、Donald Knuth の『The Art of Computer Programming』第 II 巻、第 3 版、アルゴリズム 4.3.1D を使用しました。筆算が実際にどのように機能するかを見てみましょう。 筆算を一から解説 ハードウェアにおける多倍長整数 暗号用の整数は数百または数千ビットに達し、単一のレジスタをはるかに超えるため、基数 b で格納します。各肢（limb）は 1 つのマシンワードです。x=(xn−1,…,x0)b=∑i=0n−1xibi,0≤xi<b。 多倍長算術アルゴリズムの主な構成要素は、単一/倍精度肢（limb）で動作する 4 つのプリミティブ命令です。 addc: x, y -> s, carry # s = (x+y) mod b, carry = 1 iff overflow subc: x, y -> d, borrow # d = (x−y) mod b, borrow = 1 iff underflow mul: x, y -> (hi, lo) # x·y = hi·b + lo div: (hi, lo), y -> (q1, q0), r # hi·b + lo = q·y + r, 0 ≤ r < y, q = q1·b + q0 最初の 3 つは特筆すべきものではありませんが、除算は異質なものとして残ります。2 つの単一肢乗数（multiplicand）を乗算すると常に 2 つの肢の積が返されますが、2 つの肢の被除数（dividend）を 1 つの肢の除数（divisor）で割った商は、常に単一の肢に収まるとは限らないため、2 つの肢の商を使用します。さらに、剰余が必要な副産物として現れます。ゼロ除算では未定義の動作を想定します。つまり、q と r は任意の値を取ることができますが、一部のアーキテクチャではこのケースを異なって扱います。 筆算を中規模除算に削減 私たちのタスクは、多倍長整数 u を v で割ること、つまり u=q·v+r かつ 0≤r<v となる整数 q と r を見つけることです。一般性を失うことなく、除数 v が非ゼロの最上位肢 vn−1 を持つ n 肢の整数であり、被除数 u を v より厳密に多くの肢を持つように先頭にゼロをパディングすると仮定します。u の肢の数を n+m+1 とします。また、後続の処理で自然になる別の要件があります。u の上位 n 肢を n 肢の整数として読み取ったものが v より厳密に小さくなければなりません。(un+m,…,um+1)b<(vn−1,…,v0)b、同等に⌊ubm+1⌋<v。もしそれが既に満たされていない場合は、さらにゼロを u に追加することで保証されます。この仮定は、商のサイズを m+1 肢に固定します。⌊ubm+1⌋<v⟺u<vbm+1⟺q=⌊uv⌋<bm+1。パディングが施されると、すべてのオペランドは固定された形状を持ちます。q と r を u を v で割った商と剰余とすると、次のようになります。 u = (u_{n+m}, u_{n+m-1}, ..., u_0)_b v = (v_{n-1}, v_{n-2}, ..., v_0)_b 0 < v_{n-1} q = (q_{m}, q_{m-1}, ..., q_0)_b u = q·v + r r = (r_{n-1}, r_{n-2}, ..., r_0)_b 0 ≤ r < v この説明では、3 つの異なる除算アルゴリズムを 3 つのレベルで使用します。私たちが望むのは、n+m+1 を n で割る、「長い」n+m+1/n 除算です。私たちが持っているのは、「短い」2/1 除算命令です。それらを橋渡しするために、「中規模」n+1/n 除算を使用します。q の肢を計算する自然な方法は、上から下へ進むことです。 最上位の肢は q_m = ⌊⌊u/v⌋b^m⌋ = ⌊⌊u/b^m⌋v⌋ = ⌊(u_{n+m},…,u_m)_b / (v_{n-1},…,v_0)_b⌋。 分子 ⌊u/b^m⌋ は、u の上位 n+1 肢にすぎません。したがって、商 q_m の最上位の肢は、n+1/n 除算の商自体です。そして、u の上位 n 肢が v より小さいという仮定により、それは単一の肢に収まります。(u_{n+m},…,u_{m+1})_b < v ⟺ (u_{n+m},…,u_m)_b < bv、したがって q_m = ⌊(u_{n+m},…,u_m) / v⌋ < b。 q_m を計算した後、v の倍数 b^m q_m v を u から減算し、更新された u でアルゴリズムを続行します。この u を更新する操作は、u の上位 n+1 肢を R_m、つまり q_m に対応する剰余（n 肢に収まることがわかっています）に置き換えることに相当します。したがって、更新された u は 1 肢少なくなり、アルゴリズムを続行して残りの m 肢の q を計算します。この折りたたみ技術が正しい答えを計算することが容易ではないかもしれませんが、以下のアルゴリズムを観察し、ループの各エントリおよびエグジットで 2 つの不変条件が満たされていることに注意することで、容易に推測できます。 アルゴリズム 1: 筆算 入力: u = (u_{n+m}, ..., u_0)_b、上位 n 肢が v より小さくなるようにパディング済み v = (v_{n-1}, ..., v_0)_b, v_{n-1} > 0 開始: r = u, q = 0 // r は u から始まり、最終的な剰余に縮小します ループ: for k = m down to 0: u' = (r_{k+n}, ..., r_k)_b // 位置 k における r の上位 n+1 肢 (q_k, R_k) = ⌊u'/v⌋, u' mod v // n+1/n 除算 q += q_k·b^k // 商の k 番目の肢 r -= q_k·v·b^k // (r_{k+n},...,r_k) を (0, R_k) に置き換えます 戻り値: (q, r) 不変条件: 1. u = q·v + r 2. r の上位 n 肢 < v // 肢 k+1 から n+k まで k=0 まで不変条件 2 を実行すると r<v となり、ペア (q,r) が得られます。