BigInteger.GreatestCommonDivisor() の高速化の試み【失敗 orz】

1. 概要

かの有名な D. E. Knuth 氏の著書「The Art of Computer Programming」（邦題: 準数値算法 サイエンス社）で最大公約数 (GCD) の求め方についていくつか解説されていて、その中に非常に興味深いものがありました。

その方法は基本的にビットシフトと減算を使用しており、商や剰余の計算は使用しない、というものでした。
ユークリッドの互除法では剰余の計算を繰り返しているので、「これを使ったら .NET の BigInteger.GreatestCommonDivisor() より高速に計算できるんじゃないか」とか思ってしまったわけです。

本稿では、そのネタになったアルゴリズムの紹介と、BigInteger.GreatestCommonDivisor() との性能比較について述べます。

なお、結論を先に言ってしまうと、

「すべてのケースでBigInteger.GreatestCommonDivisor()の代わりに使いたい」というほどではない。

ということになりました。無念…

詳細については、https://github.com/rougemeilland/Experiment.GreatestCommonDivisor/blob/main/docs/README_ja.md を参照してください。

2. アルゴリズムの紹介

実装を試みたアルゴリズムは、いわゆる 「互減法」 に近いものです。(互除法ではありません)
アルゴリズムの名前やら考案者が誰かとかは、本を紛失してしまったので不明です。申し訳ありません…orz

説明がやや面倒なので、まずは uint 型整数の最大公約数を計算するサンプルコードをご覧ください。

なお、「準数値算法」で紹介されていたアルゴリズムでは「偶数である限り1ビット右シフトし続ける」という部分がありましたが、本稿ではその代わりに u >>= BitOperations.TrailingZeroCount(u) とかを使用しております。

public static uint GreatestCommonDivisor(uint u, uint v)
{
    if (u == 0)
    {
        if (v == 0)
        {
            // If both u and v are 0, the greatest common divisor is undefined.
            throw new ArgumentException($"Both {nameof(u)} and {nameof(v)} are 0. The greatest common divisor of 0 and 0 is undefined.");
        }

        // If u is 0, then the greatest common divisor is equal to v
        return v;
    }

    if (v == 0)
    {
        // If v is 0, then the greatest common divisor is equal to u
        return u;
    }

    // Both u and v are non-zero.
    System.Diagnostics.Debug.Assert(u > 0 && v > 0);

    // Find the largest integer k such that both u and v are divisible by 2^k.
    // Then divide u and v by 2^k.
    var k = 0;
    {
        var zeroBitCount = Math.Min(BitOperations.TrailingZeroCount(u), BitOperations.TrailingZeroCount(v));
        u >>= zeroBitCount;
        v >>= zeroBitCount;
        k += zeroBitCount;
    }

    // Either u or v is odd.
    System.Diagnostics.Debug.Assert(u > 0 && v > 0 && ((u & 1) != 0 || (v & 1) != 0));

    u >>= BitOperations.TrailingZeroCount(u); // If it is even, shift it right until it becomes odd.
    v >>= BitOperations.TrailingZeroCount(v); // If it is even, shift it right until it becomes odd.

    while (true)
    {
        // Both u and v are odd.
        System.Diagnostics.Debug.Assert(u > 0 && v > 0 && (u & 1) != 0 && (v & 1) != 0);

        // If u == v , return u * 2^k as the greatest common divisor.
        if (u == v)
            return u << k;

        // If u < v, swap the values ​​of u and v so that u > v.
        if (u < v)
            (v, u) = (u, v);

        // Both u and v are odd and u > v.
        System.Diagnostics.Debug.Assert(u > 0 && v > 0 && u > v && (u & 1) != 0 && (v & 1) != 0);

        // Subtract v from u, so that u is always a positive even number.
        u -= v;

        // u is an even number other than 0
        System.Diagnostics.Debug.Assert(u > 0 && (u & 1) == 0);

        u >>= BitOperations.TrailingZeroCount(u); // If it is even, shift it right until it becomes odd.
    }
}

このアルゴリズムでは、GCD の定義から導かれる以下の性質を利用しています。

1. 整数 u と v がともに $2^k$ の倍数である場合、GCD(u >> k, v >> k) << k == GCD(u, v)
2. 整数 u が $2^n$ の倍数であり、整数 v が奇数である場合、GCD(u >> n, v) == GCD(u, v)
3. 任意の整数 u と v に対して、GCD(u - v, v) == GCD(u, v)
4. 任意の整数 u に対して、GCD(u, u) == u

1～3の性質を利用して 与えられた2つの整数 u と v をどんどん小さくしていき、最終的に u == v になったら 4 の性質を利用してゴール、というアルゴリズムになります。

一般的に言われる互減法と異なるのは以下の点です。

• 最初の段階で、与えられた u と v が 共に $2^k$ の倍数であれば、その $k$ を覚えておいて、以降の計算では u >> k と v >> k の最大公約数を求める。
  その結果を $k$ ビット分だけ左シフトして最終結果とする。
• u と v の片方だけが偶数であれば、偶数である方を奇数になるまで右シフトする。
• u - v は明らかに偶数になるので、減算結果は奇数になるまで右シフトする。

3. 性能の比較と分析

.NET のBigInteger.GreatestCommonDivisor()の実装と性能を比較してみた結果が以下になります。

データの ビット数	BigInteger の場合 [msec]	本稿紹介の アルゴリズム [msec]
8	0.000000	0.000000
16	0.000000	0.000000
32	0.000100	0.000000
64	0.000200	0.000100
128	0.002000	0.000900
256	0.004400	0.003100
512	0.009800	0.009300
1024	0.022800	0.029800

8ビットとか16ビットの短いデータも測定してみたのですが、あまりに短時間すぎるのか、測定不能でした。

64ビット～512ビットぐらいだと、本稿のアルゴリズムの方が高速です。
1024ビットになると、逆に BigInteger の方が高速です。

実を言うと、BigInteger.GreatestCommonDivisor() で採用されているアルゴリズムは純粋な「ユークリッドの互除法」ではありません。

ユークリッドの互除法にしろ、互減法にしろ、2つの整数を何とかして小さくしながら最終的に解を求める、という繰り返しの計算であるのは同じです。
その繰り返しの最中で、どうやら計算中の2つの整数の桁数の差が小さいケースが大半を占めるようなのですが、桁数の差が小さいと値がなかなか小さくなってくれません。

その点を改良したのが Lehmer の GCD アルゴリズムというものらしく、BigInteger.GreatestCommonDivisor() ではそれを採用しているようです。
その辺が長いビット長での性能差に効いてきているのでしょうか。

一応、Lehmer の GCD アルゴリズム についても調べてみたのですが、私のおつむではさっぱりでした…orz

4. まとめ

• ほぼビットシフトと減算のみで可能な GCD 計算アルゴリズムを紹介した。
• 本稿紹介のアルゴリズムの性能は、BigInteger.GreatestCommonDivisor() で採用されているアルゴリズムに比べて、比較的短い整数では高速であるものの、512ビット～1024ビット程度以上の整数のGCDの計算では BigInteger.GreatestCommonDivisor() の方が高速である。

まぁ、32ビット整数とか 64ビット整数とかの最大公約数を求めたいけどBigInteger といういかにも重たそうなもののお世話にはなりたくない、という場合には便利かもしれません。

5. 蛇足

2. アルゴリズムの紹介 で紹介したサンプルコードは uint 対応ですが、コード中の uint を ulong に書き換えればそのまま 64 ビットの GCD 計算ができます。
当然、BigInteger.GreatestCommonDivisor() を使うより高速です。

同じように、UInt128 に置き換えれば 128 ビット整数の GCD 計算も一応可能ではありますが、実際に計算してみると、何故か BigInteger.GreatestCommonDivisor() より時間がかかります。
x64 の機械語にはあまり詳しくないのですが、x64 の CPU に 128 ビット整数を扱える命令がそう豊富にあるとは思えないので、単純に UInt128 の四則演算やらビット演算やらの基本的な演算が低速なんだろうなぁ、と思っています。