数値アルゴリズム

要約

この単元のこれまでのレッスンでは、ループの書き方と、不変条件と変量を使ってループについて推論する方法を学んだ。このレッスンでは、その推論を、コンピュータよりも古く、いまも使われている四つの古典的な数値アルゴリズムに当てはめる。試し割りによる素数判定。これは証明できる理由によって $\sqrt{n}$ で止まり、不注意に書くと停止条件がオーバーフローする。最大公約数を求めるユークリッドの互除法。その不変条件、停止性、負の数にまつわる落とし穴、そして std::gcd との比較を扱う。逐次除算による基数変換とホーナー法。ホーナー法は数をその桁から組み立て直し、同じ図式で多項式の値を計算する。そして二分法。区間を二つに分けながら根を近似する方法で、中間値の定理に保証され、必要な周回数は前もって計算できる。どのアルゴリズムも数学的な対象として示される。すなわち、言うとおりのことをするという証明、停止するという証明、そしてステップ数の最初の見積もりを伴って示される。

学習目標

$\sqrt{n}$ まで割る素数判定を実装し、 $\sqrt{n}$ まで調べれば十分である理由を説明する。
ユークリッドの互除法を実装し、その不変条件（最大公約数は変わらない）を述べ、std::gcd と比較する。
逐次除算によって数を別の基数に変換し、位取り表現や多項式の値をホーナー法で計算する。
二分法で根を許容誤差まで近似し、必要な周回数を計算する。

約数と素数

試し割り

整数 $n \geq 2$ は、正の約数が 1 と $n$ だけのとき素数といい、ほかの約数をもつとき合成数という。0 と 1 はどちらでもない。この定義からすぐにアルゴリズムが思い浮かぶ。各候補 $d$ を 2 から $n - 1$ まで試し、どれかが割り切ったらただちにその数を合成数と判定するのである。これはレッスンループの最小の約数を求めるループであり、正しい。しかし、十桁の素数に対しては百億回近い除算を行い、そのほとんどすべてが無駄である。理由は次のとおりである。

命題

\sqrt{n}

まで試せば十分

$n \geq 2$ が合成数ならば、約数 $d$ で $2 \leq d \leq \sqrt{n}$ を満たすものをもつ。

証明

$n = a \cdot b, \quad 2 \leq a \leq b$
$n$ が合成数ならば、1 にも $n$ にも等しくない二つの因数の積として書ける。二つのうち小さいほうを $a$ と名づける。
$a^2 \leq a \cdot b = n$
$a \leq b$ の両辺に正の数 $a$ を掛けると $a^2 \leq ab = n$ が得られ、そこから $a \leq \sqrt{n}$ となる。つまり、どの因数分解でも小さいほうの因数は $\sqrt{n}$ の左側にある。
$d \leq \sqrt{n} \iff \frac{n}{d} \geq \sqrt{n}$
したがって、約数は $(d, n/d)$ という組で、 $\sqrt{n}$ を中心に対称に現れる。 $n = 36$ のとき、組は $(2, 18)$ 、 $(3, 12)$ 、 $(4, 9)$ 、 $(6, 6)$ である。 $\sqrt{n}$ までのどの候補も $n$ を割り切らなければ、その先のどの候補も割り切らない。その相手がもっと前に現れているはずだからである。

ここから、除算の回数はおよそ $n$ 回からおよそ $\sqrt{n}$ 回に減ることがわかる。 $n$ が二十億に近いとき、二十億回から四万六千回に減る。プログラムはレッスン正しいループの設計のフラグのパターンを使う。isPrime は真から始まり、ループが約数を見つければそれを下ろす。

C++ のプログラム

#include <iostream>#include <print>int main(){    int n{};    std::cin >> n;    bool isPrime{n >= 2};    for (int d{2}; isPrime && d <= n / d; ++d) {        if (n % d == 0) {            isPrime = false;        }    }    if (isPrime) {        std::println("{} は素数である", n);    } else {        std::println("{} は素数ではない", n);    }}

プログラムへの入力

2147483647

出力

2147483647 は素数である

三つの判断について述べておく。フラグは n >= 2 で初期化される。そのため 0、1、負の数はループに入らずに素数でないと判定される。この場合を忘れるのが典型的な誤りである。 $n = 1$ のときループは一度も回らず、プログラムは 1 を素数と判定してしまうからである。for の条件は二つの部分からなる。最初の約数でループを止めるフラグと、上限である。そして上限は d <= n / d と書かれており、これは正の整数については $d^2 \leq n$ と同値である（整数除算は切り捨てを行うが、 $d$ が切り捨てた商を超えないのは、ちょうど $n/d$ を超えないときである。 $d$ は整数だからである）。入力の数 $2^{31} - 1$ は int の最大値であり、素数である。プログラムはこれを $46\,339$ 回の剰余、すなわち 2 から $46\,340$ までの候補についての剰余で判定する。

ループの不変条件は、2 から $d - 1$ までのどの整数も $n$ を割り切らないということであり、これはフラグが真である間について述べたものである。入るときには、2 から 1 までの範囲が空なので成り立つ。 $d$ が $n$ を割り切らなければ範囲に割り切らない数が一つ加わるだけなので、維持される。そしてフラグが真のまま出るときには、不変条件と偽になった条件 $d^2 > n$ とが合わさって、 $\sqrt{n}$ までのどの候補も $n$ を割り切らないことを述べる。これがまさに命題の要求である。変量は、まだ試していない候補の個数、つまり $d$ と $\sqrt{n}$ の間の整数の個数である。これは一周ごとに一つ減り、0 を下回らない。

オーバーフローする条件

上限を書く自然な方法は d * d <= n であり、これは誤りである。同じ入力では、ループは $d = 46\,341$ に達し、そのとき積 $46\,341^2 = 2\,147\,488\,281$ は int の最大値を超える。

C++ のプログラム

#include <iostream>#include <print>int main(){    int n{};    std::cin >> n;    bool isPrime{n >= 2};    for (int d{2}; isPrime && d * d <= n; ++d) {        if (n % d == 0) {            isPrime = false;        }    }    std::println("素数？{}", isPrime);}

プログラムへの入力

2147483647

出力

main.cpp:9:33: runtime error: signed integer overflow: 46341 * 46341 cannot be represented in type 'int'

プログラムはコード 1 で終了しました。

未定義動作サニタイザは積のところでプログラムを止める。サニタイザがなければ、符号付きオーバーフローは未定義動作であり、規格は何が起こるかを定めていない。Linux 上の GCC で、最適化の有無にかかわらず観察されることは教訓的である。積は桁あふれして一周し、負になるか $n$ より小さくなる。そしてどの int も最大値を超えないので、条件 d * d <= n は二度と偽にならない。ループは二十億周ほど候補を試し続け、 $d = n$ （これは $n$ を割り切る）に達して、プログラムは素数を合成数と判定してしまう。オーバーフローの誤りはオーバーフローとしてではなく、誤った答えとして現れる。d <= n / d という形はオーバーフローしようがない。掛けるのではなく割るからである。レッスンビットと数学ライブラリは別の出口、すなわち std::sqrt で上限を一度だけ計算する方法を示したが、その場合は上限の正確さが浮動小数点の丸めに左右される。整数除算にはその問題がない。

コース「代数と幾何 I」のレッスン「素数と算術の基本定理」は、1 より大きいすべての整数がただ一通りに素数の積として書けることを証明している。試し割りは、見つけた約数で $n$ を割り、商について続けることによって、その積を求める最も直接的な方法である。