近況とか今年度の反省とか

今週は無理しすぎて風邪をひいてしまいました。つらい。

来年度からはもっと規則正しい生活をするようにしたいです。

それから、「いつまでにやるか明確に決まっていないけれど、いずれやらないといけないこと」は3月に片づけようと思っていたのに、いろいろ自分に言い訳して結局片付きませんでした。 こういうのを先延ばしにする癖を直したいです。

そういえば、今年度はなんだか、Twitterでは知っているけれどあったことのない人間にたくさん会った一年間でした。 そういう機会を増やしていきたいですね。

最大公約数をもっと高速に求める（その4）

先週の記事の続きです。

最大公約数をもっと高速に求める（その3） - よーる

前回示したコードは、以下のようなものでした。

uint64_t gcd_impl( uint64_t n, uint64_t m ) {
        constexpr uint64_t K = 5;
        for( int i = 0; i < 80; ++i ) {
                uint64_t t = n - m;
                uint64_t s = n - m * K;
                bool q = t < m;
                bool p = t < m * K;
                n = q ? m : t;
                m = q ? t : m;
                if( m == 0 ) { return n; }
                n = p ? n : s;
        }
        return gcd_impl( m, n % m );
}

uint64_t gcd_pre( uint64_t n, uint64_t m ) {
    for( int i = 0; i < 4; ++i ) {
        uint64_t t = n - m;
        bool q = t < m;
        n = q ? m : t;
        m = q ? t : m;
        if( m == 0 ) { return n; }
    }
    return gcd_impl( n, m );
}

uint64_t gcd( uint64_t n, uint64_t m ) {
        return n > m ? gcd_pre( n, m ) : gcd_pre( m, n );
}

分岐を外していい場合

（その1）で示したコードでは、互減法と互除法のどちらか適切な方を選ぶ分岐をするより、10回の固定回ループの方が速いことを使って高速化しました。 これは、分岐予測ミスのペナルティ時間で10回くらいのループは実行できるということによっています。 ちなみに、10回というのをちゃんと最適化すると、8回か9回がベストなようです。

今回は、80回の固定回ループとなっているため、それと同じ論理では固定回数ループの方が速いとは言えません。 つまり、互減法と互除法のどちらか適切な方を選ぶ分岐をすることで高速化する可能性が残されています。

nがmの何倍大きいときに互除法が適切だと判断するべきかという問題は、実験してみると32倍あたりに最適解があるようです。 これはループ6回分に相当します。

（その1）で示したコードでは、ループ一周当たり3cycleで8周がベストという結果でした。つまり、分岐ミスペナルティと引き換えに無駄な仕事をする価値があるのは24cycleくらいだということになります。 今回のコードでは、ループ一周当たり4cycleで6周がベストという結果だったので、やはり24cycleくらいだということが支持されます*1

そこで、以下のようなコードとしてみます。

uint64_t gcd_impl( uint64_t n, uint64_t m ) {
    constexpr uint64_t K = 5
    if( m == 0 ) { return n; }
    for( ; n / 32 < m; ) {
        uint64_t t = n - m;
        uint64_t s = n - m * K;
        bool q = t < m;
        bool p = t < m * K;
        n = q ? m : t;
        m = q ? t : m;
        if( m == 0 ) { return n; }
        n = p ? n : s;
    }
    return gcd_impl( m, n % m );
}

このようにすることで、さらに一割ほど高速化しました。

今週はちょっと忙しくて記事を書く時間が取れなかったのでちょっと少ないですが、ここまで！

ここまでで、std::gcd関数の2.4倍ほど速いgcd関数を作れています。

*1:もちろん、分岐ミスペナルティの確率が違う＆引き換えにできる仕事の内容が違うため、単純な比較はできないのですが、20-30cycleというのは一つの目安になります。

最大公約数をもっと高速に求める（その3）

先週の記事の続きです。

最大公約数をもっと高速に求める（その2）【cmova命令は遅い】 - よーる

前回示したコードは、以下のようなものでした。

uint64_t gcd_impl( uint64_t n, uint64_t m ) {
        for( int i = 0; i < 10; ++i ) {
                uint64_t t = n - m;
                bool q = t < m;
                n = q ? m : t;
                m = q ? t : m;
                if( m == 0 ) { return n; }
        }
        return gcd_impl( m, n % m );
}

uint64_t gcd( uint64_t n, uint64_t m ) {
        return n > m ? gcd_impl( n, m ) : gcd_impl( m, n );
}

このコードの機械語命令列

gcd_impl関数のループ一回分の機械語命令列は、以下のようになっています（実際にはループアンローリングされており、レジスタ割り付けだけが異なる命令列が3回周期で現れます）。

sub    rdi, rsi    # t = n - m;
cmp    rdi, rsi    # q = t < m;
mov    rax, rdi
cmovb  rax, rsi    # n = q ? m : t;
cmovnb rdi, rsi    # m = q ? t : m;
test   rdi, rdi    # m == 0
jz     Label       # if( m == 0 )

ループ一回分は7命令になっているわけですが、だからといって7 cycleかかるわけではありません。近年のCPUでは、これくらいの命令列は並列に実行されます。

今回の命令列で、並列実行のボトルネックになるのは依存の連鎖です。 具体的には、t = n - mを計算し、その結果を使ってq = t < mを計算し、その結果によって次ループのnとmが算出される部分です。

qはtが算出されなければ計算できませんし、nとmはqが算出されなければ計算できません。また、次のループのtはnとmが算出されなければ計算できません。

つまり、どんなに並列に実行できるCPUを持ってきても、この部分は各々1cycleかかってしまいます。

よってこのコードでは、ループ一回分に3cycleかかることになります。

なお、mov命令が気になるかもしれませんが、最近のCPUでは実質的に0cycleで実行されると考えてよいようです。

互減法をさらに高速化する

互減法の一ステップでは、以下のように計算が進みます。

GCD(n, m) = if n-m >= m then GCD(n-m, m) else GCD(m, n-m)

これを拡張し、以下のようにしてみます。

GCD(n, m) = if n-Km >= m then GCD(n-Km, m) else if n-m >= m then GCD(n-m, m) else GCD(m, n-m)

もしn-Km >= mであれば、互減法のKステップを一気に進めることができるため、互減法の効率が上がることが予想されます。 また、それにより、アルゴリズム全体に占める除算の割合が減少し、高速化すると予想されます。 ただし、以下の点に注意する必要があります。

  1. 途中結果n-Kmが負になる可能性がある。
  2. Kmがオーバーフローする可能性がある。

1.については、移項することで回避可能です。 2.については、前処理でmを小さくしておけば回避可能です。

改良されたメインループ

メインループの内部を以下のようにしてみます。メインループの前に前処理を置くことでうまくmは小さくしているとして、m*Kはオーバーフローしないと仮定します。

uint64_t t = n - m;
uint64_t s = n - m * K;
bool q = t < m;
bool p = t < m * K;
n = q ? m : t;
m = q ? t : m;
if( m == 0 ) { return n; }
n = p ? n : s;

n = p ? n : sが追加されたため、さらに依存の連鎖が伸び、このループ一回には最低でも4cycleかかるようになりました。 もし、m*Kの計算に4cycle以上かかる場合、そのせいで依存の連鎖が伸び、ループ一回にかかるサイクル数はもっと長くなります。 あるいは、m*Kの計算が他の命令の実行の裏で並列に実行できない場合も、ループ一回にかかるサイクル数はもっと長くなります。

最近のIntel社製CPUであれば、MUL命令は3cycleで完了するため、m*Kの計算のせいで依存の連鎖が伸びるということはありません（Kの値によってはコンパイラ最適化によりLEA命令を駆使した命令列となり、もっと速く計算されます）。 また、m*Kの計算は、他の命令の実行の裏で行うことができるため、そのせいでループ一回にかかるサイクル数が長くなるということはありません。

元のコードのループ一回には3cycleかかっていたところ、4cycleかけることで、Kステップ一気に実行できる可能性が追加されたことになります。 nとmの大きさに大きな差がなかった場合でも、少なくとも1ステップは進むため、それほど損したことにはなりません。

Kとして何をとるか

Kとして2や3などの小さな数をとれば、n-Km >= mである可能性がそれなりに高いため、改良により互減法のステップを効率よく進めることができるようになります。 Kとして大きな数をとる場合、n-Km >= mである可能性が低く、意味のある仕事をできる確率（pがfalseになり、n = sが実行される確率）が下がります。しかし、意味のある仕事をできた場合、互減法のステップを一気に進めることができます。

このトレードオフを考慮に入れ、いくつかの実験を行った結果、以下にようにすることでstd::gcdより二倍以上高速なgcd関数を作れることがわかりました。

uint64_t gcd_impl( uint64_t n, uint64_t m ) {
        constexpr uint64_t K = 5;
        for( int i = 0; i < 80; ++i ) {
                uint64_t t = n - m;
                uint64_t s = n - m * K;
                bool q = t < m;
                bool p = t < m * K;
                n = q ? m : t;
                m = q ? t : m;
                if( m == 0 ) { return n; }
                n = p ? n : s;
        }
        return gcd_impl( m, n % m );
}

uint64_t gcd_pre( uint64_t n, uint64_t m ) {
    for( int i = 0; i < 4; ++i ) {
        uint64_t t = n - m;
        bool q = t < m;
        n = q ? m : t;
        m = q ? t : m;
        if( m == 0 ) { return n; }
    }
    return gcd_impl( n, m );
}

uint64_t gcd( uint64_t n, uint64_t m ) {
        return n > m ? gcd_pre( n, m ) : gcd_pre( m, n );
}

gcd_pre関数は、m*Kがオーバーフローしないようにするための前処理です。互減法を4回繰り返すことで、mは大きくとも0xffff'ffff'ffff'ffff / 5であると保証できます。

Kとして5をとったことにより、32のような大きな商が立つ場合でも、6回程度ループを回れば互除法の一ステップを行ったことになります。これ以上大きな商が立つ場合は、除算を行った方がよさそうです。 そのようなことの発生頻度が低いため、単純な互減法であれば10回程度で除算に切り替えたほうがよかったところ、80回まで増やすことが可能になりました。 なお、gcd_implのループ回数80回は目安であり、必ずしも最適な値とは限りません。

最大公約数をもっと高速に求める（その2）【cmova命令は遅い】

先週の記事の続きです。

最大公約数をもっと高速に求める（その1） - よーる

前回示したコードは、以下のようなものでした。

uint64_t gcd_impl( uint64_t n, uint64_t m ) {
        for( int i = 0; i < 10; ++i ) {
                uint64_t t = n - m;
                bool q = m > t;
                n = q ? m : t;
                m = q ? t : m;
                if( m == 0 ) { return n; }
        }
        return gcd_impl( m, n % m );
}

uint64_t gcd( uint64_t n, uint64_t m ) {
        return n > m ? gcd_impl( n, m ) : gcd_impl( m, n );
}

コンパイラの気まぐれにつきあう

先ほど示したコードの

                bool q = m > t;

を

                bool q = t < m;

に変えてみましょう。私の環境では、25%ほど高速化しました。わけがわかりません。こんな高速化は誰も思いつけません。

アセンブリを確認しても、レジスタの割り付けは一切変わっておらず、単にcmp命令のオペランドが逆転し、後続のcmov命令の条件が変わるだけです。

他の環境でも試してみましたが、25%などという劇的な高速化ではないものの、10%程度は性能が向上するようです。

Intel社製CPUとつきあう

Intel Architecture Code Analyzerを使って調べてみたところ、Broadwell以降の場合、

  • cmovb命令やcmovnb命令は1つのμOPに変換される。実行は一サイクルで完了する。
  • cmovbe命令やcmovnbe命令は2つのμOPに変換される。実行には合計で二サイクルかかる。
    • cmova命令はcmovnbe命令の別名称。
    • cmovna命令はcmovbe命令の別名称。

というようになっているようです*1。

Code Analyzerなんて持ち出さなくても、インテルの最適化リファレンスマニュアルに書いてあったようです（最後の方、C-17ページ。pdfとしては759ページ目796ページ目（2019/09/15修正））。

https://software.intel.com/sites/default/files/managed/9e/bc/64-ia-32-architectures-optimization-manual.pdf

bool q = m > tの場合、cmovbe命令とcmovnbe命令が使われ、レイテンシが二サイクルとなります。一方、bool q = t < mの場合、cmovb命令とcmovnb命令が使われ、レイテンシは一サイクルとなります。 これにより高速化したということのようでした。

なお、Haswell以前の場合、前者は2μOP、後者は3μOPに分解され、それらが直列に実行されるようなので、その差は小さくなるとはいえ、やはり性能に差が出ます。

まとめ

cmov命令は、使う条件によっては遅いことがあります。（はじめて知りました！） そのため、不等号の方向を変えるだけという一見等価な変換で高速化したり、遅くなったりすることがありえるようです。

今回の記事はもっとアルゴリズム的なところを掘り下げようと思っていたのですが、Intel社製CPUに最適化するだけの回になってしまいました。（でも25%高速化！）

まだまだ続きます。

*1:似たような命令なのに動作原理が違うあたり、x86はやはりCISCだなぁという感じです

最大公約数をもっと高速に求める（その1）

最大公約数を求める非常に有名なアルゴリズムとして、ユークリッドの互除法が知られています。 このアルゴリズムは、二つの入力を十進法で表したとき、小さいほうの桁数の高々五倍の回数の剰余演算で最大公約数が求まることが知られています（ラメの定理）*1。二つの入力が符号なし64bit整数であれば、高々92回です。これは最悪ケースであり、多くの場合はかなり少ない回数で済みます。実験してみると、平均的には37回程度のようです。

しかし、剰余演算は非常に遅い演算です。Intelの最近のCPUでは、21cycle～83cycleかかるとされています。ここから計算すると、二つの符号なし64bit整数の最大公約数を求めるには、1μ秒弱の時間がかかることになります。

実際、以下のコードを実行してみると、乱数生成を除いて7秒ほどかかります。つまり一回当たり0.42μ秒ということで、だいたい計算と合います。

#include <iostream>
#include <random>
#include <chrono>

// 16M
std::pair<uint64_t, uint64_t> arr[0x100'0000];

int main() {
        std::mt19937 mt;
        std::uniform_int_distribution<uint64_t> dist( 1, 0xffff'ffff'ffff'ffff );
        for( auto& [e1, e2] : arr ) {
                e1 = dist( mt );
                e2 = dist( mt );
        }

        const auto start = std::chrono::system_clock::now();

        for( auto& [e1, e2] : arr ) {
                std::gcd( e1, e2 );
        }

        const auto finish = std::chrono::system_clock::now();
        std::cout << std::chrono::duration_cast<std::chrono::nanoseconds>(finish - start).count() / 1e9 << " sec" << std::endl;
}

コンパイル：clang++ -std=c++17 -march=native -O3

環境：Windows subsystem for Linux、clang5.0.1、Intel(R) Core(TM) i5-7200U@2.50GHz（ターボブースト時3.1GHz）

最大公約数を求めるだけの問題に出会うことはそれほどありませんが、アルゴリズムの内部で最大公約数を求める必要があり、ユークリッドの互除法を部品的に使うアルゴリズムはそれなりにあります。そういったアルゴリズム群の（定数倍）高速化のために、gcd関数を高速化することを考えます。

互除法の仕組み

まずは、ユークリッドの互除法で最大公約数（greatest common divisor: GCD）が求まる理由を復習します。以下では、GCD(n, m)と書いたらnとmの（数学的事実としての）最大公約数を意味することとし、gcd(n, m)と書いた場合は、（帰納的な定義を持つ）最大公約数を求める関数の定義か、それにnとmを代入したものを意味することとします。

まず、GCD(n, m) = GCD(n-m, m)であることが知られています（これの証明は大学入試問題として有名ですが、ここでは割愛します*2）。そのため、gcd(n, m) = if n < m then gcd(m-n, n) else gcd(n-m, m); gcd(0, m) = mといった形で最大公約数を求める関数を、帰納的に定義することができます。これは互減法というアルゴリズムであり、これを使うとコードサイズを小さくできる可能性があります（下記参照。golfer向けであり、実際にプログラムサイズを小さくするのに役に立つとは限りません）。しかし、最悪の入力、つまりgcd(99999,1)などに対して99999回の引き算が発生するため、低速です。

http://golf.shinh.org/reveal.rb?Modular+Inverse/tails_1484205185&c

2019/03/04 01:50 追記：（注意）上記コードは、互減法そのものではありません。モジュラ逆数を求めよという問題に対する解答であり、いわゆる拡張ユークリッドの互除法を使うのが一般的であるところ、互除法ではなく互減法を用いてコードを縮めた例です。

このコードを提出したtailsさんから、混乱を招く恐れがあるとの指摘をいただいたため、追記しました。

GCD(n, m) = GCD(n-m, m)を繰り返し適用すると、GCD(n, m) = GCD(n-m, m) = GCD(n-2m, m) = ... = GCD(mod(n, m), m)となることがわかります。nとmは対称なので、n≧mと仮定しても一般性を失いません。この条件のもと、gcd関数は、gcd(n, m) = gcd(m, mod(n, m)); gcd(n, 0) = nのように定義することができます。これがユークリッドの互除法です。

互除法の各ステップの解析

いくつか手で実験をしてみればわかりますが、[1, N]の区間から一様ランダムに取り出した二つの整数n, mに対する互除法の場合、商は1が立ち剰余はn-mになる、というパターンが多いことがわかります。

これの理由は単純で、商として2以上が立つには、片方がもう一方の二倍以上である必要があるからです。そのようになる確率は、二数の選び方が一様ランダムであれば1/2です。 互除法を進めていくと確率は片寄っていきますが、商として1が立つパターンが頻発することは変わりません。また、商が1でなくとも、商は小さいことが多く、大きな商の出現は稀です。一つ一つの数に対して確率が小さいのはもちろんのこと、それらすべての確率の総和も小さいです。

N=264-1として実験してみると、1が立つのは41%、2が立つのは17%、3が立つのは9%、4が立つのは6%、5が立つのは4%、6が立つのは3%、までで八割を占めています。そのあとを追ってみると、13以下が90%、26以下が95%、128以下が99%です。一様ランダムな場合、商はその逆数に比例した確率で出現します。互除法を進めていくにつれ、一様ランダム性が崩れる（区間が減少するのに加え、偶然小さい数が出現する確率が大きくなる）ため解析は複雑になります。しかし、小さい数が出現すると互除法が速やかに終了するため、全体的な確率には大きな影響を及ぼさないようにも思われます。結局のところ、商の出現確率はその逆数に比例していると考えても、大きく外れてはいないでしょう。

剰余演算は遅いので避ける

剰余演算を使うことで単純な引き算と比べて高速化できる、というのが互除法でした。しかし、剰余演算は非常に遅い演算であり、できる限り避けたいところです。

互除法の各ステップにおいて、商として1が立つ場合、数十サイクルもかけて除算を行うのは非効率的であり、1cycleで終わる引き算を行えばよいことになります。 先ほどの解析によれば、商として1が立つことは多い割合で発生するため、数十サイクルのロスを頻繁に行っていることがうかがえます。

商として何が立つかは実際に割り算を行わないといけないため、そのような効率化はできないように思えるかもしれません。

しかし、m != 0の時、`n >= m && m > n - mがn / m == 1と等価なため、以下のようにgcd関数を作ることができます。

uint64_t gcd_impl( uint64_t n, uint64_t m ) {
        if( m == 0 ) { return n; }
        uint64_t t = n - m;
        if( m > t ) {
                return gcd_impl( m, t );
        } else {
                return gcd_impl( m, n % m );
        }
}

uint64_t gcd( uint64_t n, uint64_t m ) {
        return n > m ? gcd_impl( n, m ) : gcd_impl( m, n );
}

しかし、これを使っても速度が数%改善される程度にとどまります。これは、商として1が出るパターンの出現頻度は四割であり、かつそれは予測不可能に出現するため、分岐予測が全く当たらないのが原因です。

当てられない分岐は遅いので避ける

高速化のテクニックに、当てられない分岐を取り除くというものがあります。

m == 0も予測不可能に出現しますが、これなしというわけにはいかないこと、出現は最後に一回だけであり影響は軽微であることから、この分岐を取り除くことは考えません。 n > mは、入力の最初に50%の確率で発生し、これは絶対に当てられません。

よって、m > tの分岐を取り除くことを考えます。その前にちょっと準備をします。

予測できない分岐なしに互減法を実装する

予測できない分岐を取り除くため、x86の命令セットには、cmov (conditional move)命令が用意されています。これは分岐なしにcond ? a : b;のような計算をすることができる命令です。昔のIntelはさぼっていたのかこの命令のレイテンシが2cycleでしたが、最近では1cycleとなっているため、さらに気軽に使えるようになりました。cmov命令を使う欠点は、分岐すれば計算する必要がなかったかもしれないのに、aもbも計算しないといけなくなる点です。

互減法は以下のように書けます。

while( m != 0 ) {
        uint64_t t = n - m;
        bool q = m > t;
        n = q ? m : t;
        m = q ? t : m;
}

このコードを書いたとしても、必ずしもcmov命令が使われるとは限りません。実際g++でコンパイルすると条件分岐を使ったコードにコンパイルされてしまいます。一方clang++でコンパイルするとcmov命令が使われました。

互除法と互減法を組み合わせる

先ほどの実験では、商として1が現れることが多く、特に128以下の商が立つ状況は99%の確率で発生することを明らかにしました。そこで、互除法ではなく互減法の方が速いのではないか、と思いたくなりますが、これは誤りです。一様ランダムに選ばれた符号なし64bit整数同士の互除法は実験してみると平均37ステップかかることがわかります。つまり、三回に一回は128を超える商がどこかのステップで出現することとなり、そのような場合に引き算を続けるのは非常に非効率的です。

一方、多くの場合に商が1であり、そのような場合は除算より引き算の方が一桁速いことを考えると、互減法が優れている局面は確かに存在することもわかります。

よって、うまく互除法と互減法を切り替える、または組み合わせるなどの手法が有効そうであることがわかります。

ここで、切り替えるために条件分岐を使ってしまうのは、やはり当てられない分岐となるためによくなさそうです。実際、先ほどの失敗は当てられない分岐を作ってしまったことによるものでした。

分岐せずにうまく切り替える方法はあるのでしょうか？

キーとなるアイデアは、「分岐予測に失敗して時間を浪費するくらいなら、その時間に無駄になるかもしれない仕事をしておけば、その仕事が無駄ではなかった場合に高速化する」というものです。このアイデアはニュートン・ラフソン法を用いている平方根関数を高速化するときにも使いました。これは、「収束判定という予測できない分岐をするかわりに、どんな入力に対しても収束する回数だけ回す」というものです。すでに収束している場合は無駄な加速を行っていることになりますが、それでも分岐予測に失敗するよりは速いのです。

このアイデアを使うと、以下のようなコードになります。

uint64_t gcd_impl( uint64_t n, uint64_t m ) {
        for( int i = 0; i < 10; ++i ) {
                uint64_t t = n - m;
                bool q = m > t;
                n = q ? m : t;
                m = q ? t : m;
                if( m == 0 ) { return n; }
        }
        return gcd_impl( m, n % m );
}

uint64_t gcd( uint64_t n, uint64_t m ) {
        return n > m ? gcd_impl( n, m ) : gcd_impl( m, n );
}

ここで、ループの回数10回は適当であり、最適化の余地があります。このコードは、互減法を10回行った後、互除法的な除算を一回やる、というのを繰り返すものになっています。商が小さければ互減法の部分でnとmが小さくなってくれるため、トータルの除算回数が減ることになり、高速化します。一方、商が大きくなる場合でも、10回しか引き算せずに除算を行うようになっているため、延々と引き算をし続ける、といった非効率的な状況にはまることを避けられます。

なお、10回のループに関してはループアンローリングが行われることを期待しています。もしそうでなくても、少なくとも当てられる分岐にはなっています。よって、この関数の中には、if( m == 0 ) { return n; }という避けようのない分岐を除いて、予測不可能な分岐はありません。

先ほどのアイデアがどう使われているかを別の言葉でも説明してみます。

もし商が大きい場合、10回の引き算しても意味はなく、なるべく早く除算モードに切り替えることが高速化につながります。しかし、商が大きいという予測不可能な分岐をしてモードを切り替えるくらいなら、10回程度の引き算という無駄な仕事をした方が高速です。

逆に（運よく）商が小さい場合、引き算は無駄な仕事とならず、高速化に貢献します。

まとめ

最大公約数を求める有名なアルゴリズム、ユークリッドの互除法には除算が含まれており、一般的なCPUで計算すると遅くなりがちです。互減法を部分的に組み合わせて高速化を図ろうとしましたが、予測不可能な分岐が追加されたため、それほど高速化はしませんでした。そこで、無駄かもしれない仕事と引き換えに予測不可能な分岐を除去する方法を考えました。その手法を使うと、元のナイーブな互除法に比べて1.5倍ほど高速に動作するプログラムができました。

*1:実際の比例係数は、黄金比phiを底とする対数を用いてlog_phi(10)とあらわせて、これは4.8ほどです。

*2:というより忘れて思い出せません……。

完全精度 expf 関数の作り方

C言語では、expf関数などの超越関数の精度は一切保証されないことになっています*1。 そうは言っても、数学的な結果に近い結果が返ってくることが、常識的に信じられています。

数学的な結果に最も近いfloatやdoubleが返ってくることを、完全精度（丸め誤差0.5ulp保証*2）といいます。また、数学的な結果に最も近い、または二番目に近い値が返ってくることを、丸め誤差1ulp保証といいます。

一般に、超越関数の結果は、丸め境界に非常に近い値になり得るため、完全精度を実現するのは困難です（テーブルメーカーのジレンマ）。一方、1ulp保証ならば現実的な計算量で実現可能です。

ちょっと調べてみたところ、手元でビルドしたg++8.2では、double std::exp(double)は完全精度っぽいですが、float std::exp(float)は単にdouble std::exp(double)関数の結果を丸めた結果*3を返しているような挙動を示しており、完全精度になっていません。 そこで、自分で作ってみます。

完全精度expf関数の作り方

ここでは、高速化については考えず、とにかく完全精度を追求するものとします。また、テーブルに頼ることはできるだけ避けます。

正しい方針

実は、expfであれば完全精度を達成することは意外と簡単です。

最も丸め境界に近いのはexpf(-0x1.d2259ap+3) = 0x1.fa6636p-22 であり、真の値は0x1.fa66350000001p-22付近にあります。この場合であっても、double精度で正確に計算し最後に一回だけ丸めればぎりぎり有効数字が足ります。

しかし、そのためには簡易的なdouble-double演算*4が必要になってくるため、やや計算量が高くつきます。

どのくらい雑にやってよいのか

実際には丸め境界に近い値はそれほど多くありません。double精度の下2bitを正確に求める必要があるのは、上記の例に加えて

  • expf(-0x1.e1dbe2p-8) = 0x1.fc3fd2p-1 真の値は0x1.fc3fd10000002p-1付近
  • expf(-0x1.c1c4b8p-10) = 0x1.ff1f4ep-1 真の値は0x1.ff1f4effffffcp-1付近
  • expf(0x1.fdff02p-17) = 0x1.000100p+0 真の値は0x1.0000ff0000003p+0付近
  • expf(0x1.cd3982p-14) = 0x1.000734p+0 真の値は0x1.000734ffffffcp+0付近
  • expf(0x1.8d7cb6p-12) = 0x1.0018dap+0 真の値は0x1.0018d90000004p+0付近

くらいしか存在しないため、運が悪くなければ慎重にdoubleで計算するだけで丸め誤差の問題を発生させずに済むかもしれません。

形式的な証拠とはなりえませんが、定義域のすべての値について正しく丸められていることを確認すれば、完全精度を達成していること自体は確かめられます。 expfの定義域は-104.0fから88.0fほどであり、22.4億通り（231よりちょっと多いくらい）で正しく丸められているかを確認すれば目的を達成できます。 これは、最近の*5パソコンであれば、特に工夫をしなくても1時間程度で終わります。

実際にとった方針

基本的にはテーラー展開を使って求めることになりますが、引数の値が大きい場合、収束に必要な項数が多くなります。これは誤差の増大も意味するため、避けるべきです。

方針としては、 \exp(s+t) = \exp(s) \times \exp(t) をうまく活用していきます。 まず、 \exp(s)が簡単に計算できるように、かつtの絶対値が小さくなるように、うまくsを選びます。 こうすると、 \exp(t)をテーラー展開で求めるときに必要な項数が少なくなります。

\exp(s) の計算法

\exp(s) を簡単に計算する手法として、メジャーなのは以下の二つです。

  • s を \ln(2) の整数倍になるようにとる。すると単に二のべき乗となり、std::ldexp関数で効率よく計算できる。tの絶対値は0.5以下となる。
  • \exp(s) のテーブルを持つ。tの絶対値を小さくすることができるため、上の方法より収束が速い。収束の速さはテーブルサイズとのトレードオフとなる。

しかし、今回は難しいことを考えず、sとして整数をとりました。 \exp(s)は繰り返し二乗法により求めました。これによりいくらか誤差が乗りますが、最終的な結果への影響はないようです。

\exp(t) の計算法

tの絶対値が0.5以下の時、| t15 / 15! | < 2^{-55} であるため、double精度を出すためには14次までのテーラー近似を利用すれば十分です。

なお、テーラー展開を13次までとした場合は最終的な結果に誤差が発生しました。

テーラー近似は展開中心点から離れるほど近似度が悪くなることが知られています。この問題を解決する方法として、minimax近似というものあります。 minimax近似は、特定の区間内での最大誤差が最も小さくなる近似多項式です。しかしその係数は複雑になります。

今回はややこしいことを考えず、単純なテーラー近似を利用しました。

ここで、テーラー近似を求める段階で注意が必要です。

一般に多項式の値を求める場合、乗算回数が少なくなる*6ホーナー法が用いられます。

ホーナー法とは、 a_0 + x \times ( a_1 + x \times ( a_2 + x \times a_3 ) ) のように計算する多項式計算法のことです。

この方法は、係数の誤差が蓄積する性質があるため、限られた精度で正確な値を求めるのに向いていないようです（むしろ乗算回数が少なくなることから、多倍長演算向けだと思われます）。

今回は愚直に、 \sum_{i=0}^{14} c_i t^{i} を求めることにしました。ここで、絶対値の小さいほうから足していくことで計算誤差を最小化します。

実際のコード

constexpr std::array<double, 15> make_exp_coeff() {
        std::array<double, 15> arr {};
        double fact = 1.0;
        for( int i = 1; i < 15; ++i ) {
                fact *= i;
                arr[i] = 1. / fact;
        }
        return arr;
}

double power( double x, unsigned int n ) {
        double ret = 1.0;
        for( ; n; n >>= 1 ) {
                if( n & 1 ) { ret *= x; }
                x *= x;
        }
        return ret;
}

// |t| <= 0.5  --> | t^15/15! | < 2^-55
double exp_taylor14( double t1 ) {
        constexpr auto coeff = make_exp_coeff();
        const double t2 = t1 * t1;
        const double t4 = t2 * t2;
        const double t8 = t4 * t4;
        double ret = 0.0;
        ret += t8 * t4 * t2      * coeff[14];
        ret += t8 * t4      * t1 * coeff[13];
        ret += t8 * t4           * coeff[12];
        ret += t8      * t2 * t1 * coeff[11];
        ret += t8      * t2      * coeff[10];
        ret += t8           * t1 * coeff[9];
        ret += t8                * coeff[8];
        ret +=      t4 * t2 * t1 * coeff[7];
        ret +=      t4 * t2      * coeff[6];
        ret +=      t4      * t1 * coeff[5];
        ret +=      t4           * coeff[4];
        ret +=           t2 * t1 * coeff[3];
        ret +=           t2      * coeff[2];
        ret +=                t1 * coeff[1];
        return ret + 1.0;
}

float myExp( float x ) {
        constexpr double exp_plus_one  = 0x1.5bf0a8b145769p+1;
        constexpr double exp_minus_one = 0x1.78b56362cef38p-2;
        const int s = std::round( x ); // round to nearest(ties to max magnitude)
        const double t = x - s; // |t| <= 0.5

        const double exp_s = s > 0 ? power( exp_plus_one, s ) : power( exp_minus_one, -s );
        const double exp_t = exp_taylor14( t );
        return static_cast<float>( exp_s * exp_t );
}

検証コード

#include <iostream>
#include <cstdint>
#include <cstring>
#include <cassert>

#include <kv/interval.hpp>
#include <kv/dd.hpp>
#include <kv/rdd.hpp>

using idd = kv::interval<kv::dd>;

template<class To, class From>
To bit_cast(const From& from) noexcept {
        To to;
        static_assert( sizeof to == sizeof from );
        std::memcpy(&to, &from, sizeof to);
        return to;
}


float next( float x ) {
        assert( !std::signbit(x) );
        return bit_cast<float>( bit_cast<uint32_t>(x) + 1 );
}

float before( float x ) {
        assert( !std::signbit(x) );
        if( x == 0 ) { return -0.0f; }
        return bit_cast<float>( bit_cast<uint32_t>(x) - 1 );
}

idd intervalFloat( float x ) {
        double y = next(x);
        double z = before(x);
        return idd( (x+z)/2, (x+y)/2 );
}

float nearestFloat( idd ia ) {
        // 二重丸め
        float x = (float)(double)mid(ia);
        float y = next(x);
        float z = before(x);
        idd ix = intervalFloat(x);
        idd iy = intervalFloat(y);
        idd iz = intervalFloat(z);
        if( subset(ia, ix) ) { return x; }
        if( subset(ia, iy) ) { return y; }
        if( subset(ia, iz) ) { return z; }
        std::cerr.precision(50);
        std::cerr << x << std::endl;
        std::cerr << y << std::endl;
        std::cerr << z << std::endl;
        std::cerr << ia.lower() << std::endl;
        std::cerr << ia.upper() << std::endl;
        throw nullptr;
}

#include "my_exact_expf.hpp"

void check( float x ) {
        idd ix(x);
        auto exact_exp_range = exp(ix);
        float nearest = nearestFloat(exact_exp_range);
        float my = myExp(x);

        if( nearest != my ) {
                std::cerr << "x     " << x << std::endl;
                std::cerr << "x     " << std::hexfloat << x << std::defaultfloat << std::endl;
                std::cerr << "my    " << my << std::endl;
                std::cerr << "exact " << exact_exp_range << std::endl;
                std::cerr << abs(mid(exact_exp_range) - my) / (nearest - before(nearest)) << " ulp" << std::endl;
        }
}

int main() {
        std::cerr.precision(50);
        const uint32_t minimum = bit_cast<uint32_t>(-104.0f);
        const uint32_t maximum = bit_cast<uint32_t>(88.0f);
        for( uint32_t i = bit_cast<uint32_t>(-0.0f); i < minimum; ++i ) {
                check( bit_cast<float>(i) );
        }
        for( uint32_t i = bit_cast<uint32_t>(0.0f); i < maximum; ++i ) {
                check( bit_cast<float>(i) );
        }
        return 0;
}

この検証コードのコンパイルには、区間演算ライブラリkvが必要です。

github.com

なお、そのままコンパイルするにはboostが必要ですが、このライブラリはC++03対応を維持するためにC++11で導入された標準関数を使っていないようです。 そのため、以下の強引な置換により、boostへの依存は解消されます。

sed -i s/boost/std/g
sed -i s/enable_if_c/enable_if/g

コンパイルオプションは、高速化のため、g++-8 -std=c++17 -march=native -O3 -DKV_FASTROUND -DKV_USE_TPFMA -I ../../kv/を指定しました。使ったCPUはIntel(R) Core(TM) i5-7200U CPU @ 2.50GHzです。かかった時間はほぼ一時間ぴったりでした。

区間演算によりexp(x)の真値を含む区間が求まります。正しくfloatに丸めた場合の可能性が二つ以上ある場合はthrow nullptrされるはずですが、そのようなことは起こりませんでした。 また、(float)(double)mid(ia)の部分はdoubleに丸めてからfloatに丸めているので正しく丸められていない可能性がありますが、依然1ulp保証ではあるので前後のfloat値を確認しています。 本当は正しく求める方法もあるのですが、手を抜きました。

このプログラムが何も出力せずに終了したことにより、実装したmyExp関数は完全精度であることを確かめました。0.0付近の非正規化数であっても0.5ulp保証されています。

ただし、inf付近でどちらに丸めるべきかについては、よくわからなかったので、その部分について0.5ulp（？？）保証されているかはわかりません。

g++のよくわからない挙動

g++では、なぜかstd::expf関数が使えません。::expf関数はあるらしいです。::expf関数はおしいですが完全精度になっていません（この記事で書いたコードが意味のないものとならなくてよかったです）。また、type detectionイディオムを使った確認により、std::exp(float)オーバーロード自体は存在し、返り値の型はfloatであることを確認しました。

#include <cmath>

template<class T>
struct TD;

int main() {
  TD<decltype(std::exp(1.0f))> td;
}

なぜstd::exp(double)と::expf(float)は完全精度であるのにstd::exp(float)だけ精度が低いのでしょうか？

::expf 関数の精度

以下の98個の入力について0.5ulpを超える誤差があることを確認しました。最大の誤差は0.5000031997 ulp 未満でした。

  • -0x1.c1c4b8p-10 （上で紹介した入力のうちの一つ）
  • -0x1.abb574p-9
  • -0x1.f5723ep-9
  • -0x1.f9c67ep-9
  • -0x1.fd2ddcp-9
  • -0x1.021d74p-8
  • -0x1.2e4ab4p-8
  • -0x1.400432p-8
  • -0x1.4b0f34p-8
  • -0x1.4be042p-8
  • -0x1.4d8260p-8
  • -0x1.4dcdc2p-8
  • -0x1.58ef7ap-8
  • -0x1.db8546p-7
  • -0x1.f691dap-7
  • -0x1.f99ab6p-7
  • -0x1.02f812p-6
  • -0x1.074b54p-6 （最大誤差となる入力。真値は0x1.f7d67affff94ap-1付近ですが、0x1.f7d67cp-1が返ってきます）
  • -0x1.a59e4ep-6
  • -0x1.b17ee0p-6
  • -0x1.b1db4ap-6
  • -0x1.8e02c8p-5
  • -0x1.cd9502p-4
  • -0x1.d0a5b6p-4
  • -0x1.fdcbf6p-4
  • -0x1.2b6a50p-3
  • -0x1.0faa34p-2
  • -0x1.19662ep-2
  • -0x1.d66d74p-2
  • -0x1.074b34p-1
  • -0x1.c2aaeap-1
  • -0x1.df1c00p-1
  • -0x1.efdceap-1
  • -0x1.08abe2p+0
  • -0x1.74a396p+0
  • -0x1.7f4296p+0
  • -0x1.b75242p+0
  • -0x1.ff2602p+0
  • -0x1.7dd9aap+1
  • -0x1.d2445ap+1
  • -0x1.f66118p+1
  • -0x1.2b0b04p+2
  • -0x1.90ced0p+2
  • -0x1.bc933ep+2
  • -0x1.e4c6c0p+2
  • -0x1.e8ebf0p+2
  • -0x1.1089b0p+3
  • -0x1.287d08p+3
  • -0x1.93e788p+3
  • -0x1.d2259ap+3 （上で紹介した、最も丸め境界に近くなる入力）
  • -0x1.de6f7cp+3
  • -0x1.1f8dc4p+4
  • -0x1.433aaep+4
  • -0x1.b43edcp+4
  • -0x1.c95efap+4
  • -0x1.23134ap+5
  • -0x1.444328p+5
  • -0x1.2219dep+6
  • -0x1.2c4726p+6
  • -0x1.5ce26ap+6
  • 0x1.6438b6p-8
  • 0x1.9785f6p-8
  • 0x1.1c2c32p-6
  • 0x1.2dc3f0p-6
  • 0x1.2e3554p-6
  • 0x1.c649eap-6
  • 0x1.3a7784p-5
  • 0x1.980654p-5
  • 0x1.eaddb4p-5
  • 0x1.eff8cep-5
  • 0x1.4feb94p-4
  • 0x1.599d64p-3
  • 0x1.72793ap-3
  • 0x1.b28da4p-3
  • 0x1.f3518cp-3
  • 0x1.54f076p-2
  • 0x1.0d5f92p-1
  • 0x1.39a300p-1
  • 0x1.70f4e4p-1
  • 0x1.b94234p-1
  • 0x1.bf4b50p-1
  • 0x1.13f37ap+0
  • 0x1.c8264ap+0
  • 0x1.4c1e18p+1
  • 0x1.69512ep+1
  • 0x1.172096p+2
  • 0x1.24b172p+3
  • 0x1.3983a8p+3
  • 0x1.56f768p+3
  • 0x1.ddc228p+3
  • 0x1.5d4bd4p+4
  • 0x1.6255f4p+5
  • 0x1.91fceap+5
  • 0x1.d8afbep+5
  • 0x1.097f6ap+6
  • 0x1.15298ap+6
  • 0x1.1ddcb8p+6
  • 0x1.3829eep+6

*1:超越関数だけでなく、sqrt関数以外の無理関数は精度が保証されません。

*2:ulpはunit in the last placeの略で、最後の桁の重みのことです。数学的な結果に最も近い値に丸めた場合、その丸め誤差は0.5ulpになります。

*3:これでは二回丸めが発生し、正しくありません。なぜ二回丸めると正しくないかというと、0.495→0.5→1.0のような間違いが発生するためです。

*4:doubleを二個使う多倍長表現のことをdouble-doubleと呼びます。疑似四倍精度とも呼ばれます。

*5:FMA命令に対応しているHaswell以降のCPUでは、double-double演算が高速化します。FMA命令は計算結果を一回だけ丸めるため、演算誤差の補正が簡単になるためです。

*6:必ずしも乗算回数が最小となるわけではありません。さらに乗算回数が少ない方法として、Knuthのadaptation of coefficientsという方法があるようです。

C++におけるinline指定子の誤解

C++にはinline指定子があります。しかし、このinline指定子は、その名前から誤解されがちです。

(誤解)inline指定子を付けると、その関数はインライン展開される

これは誤解であり、正しくありません。正しくは、以下のような意味になります。

  • inline指定子を付けた場合、複数回の定義が許される。また、複数の翻訳単位*1で定義されたとしても実体は一つになることが保証される*2。
  • (C++14まで)inline指定子を付けた場合、コンパイラはその関数をインライン展開するヒントとすることができる。

「inline指定子を付けた場合、コンパイラはその関数をインライン展開するヒントとすることができる。」という記述はC++17以降で無くなりました。 ただし、コンパイラの最適化は見える結果が変わらない範囲で自由に行えるため、inline指定子をヒントにすることは問題ありません。（2020/09/02 なくなっていないとの指摘を受けましたので削除）

ちなみに、テンプレート関数、constexpr関数、クラス定義内で定義される関数、はいずれも暗黙にinline指定されたものとみなされます。default定義されたメンバ関数も暗黙にinline指定されたものとみなされます。

ところで、inline指定されている場合（これには上述の暗黙指定の場合を含みます）、複数回の定義が許されますが、その場合でも全ての定義は字面上同一かつ意味上同一でないといけません（one definition rule, ODR）。

外部リンケージを持つ（static指定されておらず、無名名前空間に囲まれていない）関数がODRを守れていない場合、コンパイラは特に指摘してくれません（診断不要：no diagnostics required）。このようなことが発生すると、リンクの順番によってプログラムの実行結果が変わるなどの意味不明なバグにつながり、非常にデバッグが困難です。つまり、うっかり同名の関数や変数を書いてしまってもコンパイラが指摘してくれず、バグが埋め込まれるということになります。そのため、特にグローバル名前空間にinlineなものを置くことは忌避されるべきです。

うっかりを防ぐためには、static inlineとすべきです。ただし、実行効率まで考えると常に最善とは限りません（後述します）。

(誤解)inline指定子を付けると、その関数はインライン展開されやすくなる

先ほどの誤解よりは少しましになりましたが、依然誤解です。「コンパイラはその関数をインライン展開するヒントとすることができる」をうがって読むと、「コンパイラはその関数をインライン展開 しない ことの根拠としても使える」ことがわかります。そんな例は実在するのかという感じですが、以下に実例を示します。

inc.h

#include <array>
#include <cstdint>

#define M1(i) s += (i)*(i), arr[s%SIZE] = i;
#define M2(i) M1(i) M1(i+1)
#define M4(i) M2(i) M2(i+2)
#define M8(i) M4(i) M4(i+4)
#define M16(i) M8(i) M8(i+8)
#define M32(i) M16(i) M16(i+16)
#define M64(i) M32(i) M32(i+32)
#define M128(i) M64(i) M64(i+64)


DECL_SPEC std::array<std::size_t, 128> LargeFunction(std::size_t s) {
    std::array<std::size_t, 128> arr {};
    for( int t = 0; t < 128; ++t ) {
        M128(0);
    }
    return arr;
}

source1.cpp

#include "inc.h"

std::array<std::size_t, 128> callLargeFunction(std::size_t s) { return LargeFunction(s); }

source2.cpp

#include "inc.h"

std::array<std::size_t, 128> anotherCallLargeFunction(std::size_t s) { return LargeFunction(s); }

main.cpp

#include <iostream>
#include <array>
#include <cstdint>

std::array<std::size_t, 128> callLargeFunction(std::size_t s);
std::array<std::size_t, 128> anotherCallLargeFunction(std::size_t s);

int main() {
    std::cout << callLargeFunction(1)[3] << std::endl;
    std::cout << anotherCallLargeFunction(2)[4] << std::endl;
}

以下のコンパイルオプションで分割コンパイルし、リンクした場合、バイナリのサイズは13120バイトとなります。コンパイルオプションとして-Osの代わりに-O1、-O2、-O3を使う場合、17176バイトとなります。-fltoオプションを付ける場合、四つの最適化オプション-O1、-O2、-O3、-Osのいずれでも12936バイトとなります。

g++-8 -std=c++17 -Os -DDECL_SPEC=inline source1.cpp -c -o source1_inline.o
g++-8 -std=c++17 -Os -DDECL_SPEC=inline source2.cpp -c -o source2_inline.o
g++-8 -std=c++17 -Os main.cpp -c -o main_inline.o

g++-8  source1_inline.o source2_inline.o main_inline.o -o inline

一方、-DDECL_SPEC=inlineの代わりに-DDECL_SPEC=staticや-DDECL_SPEC='static inline'を指定した場合、最適化オプションによらず17176バイトとなります。また、-fltoオプションを付けた場合、-O1の場合のみ12896バイト、その他の三つの最適化オプション-O2、-O3、-Osのいずれでも12936バイトとなります。

つまり、 staticやstatic inlineの代わりにinlineとすることでインライン展開を抑制することができる 場合がある（具体的には-Os指定かつ-flto指定なし）ということです。

なぜこのようなことが起こったのでしょうか？

-fltoオプションを付けた場合、他の翻訳単位で使われるか否かを見てからインライン展開をするかどうかを判断することができるようにコンパイルすることになります。よって、基本的にはインライン展開せずにコンパイルし、リンク時にインライン展開を行うか判断することになります。そのため約13KBの小さなバイナリとなります。

一方、-fltoオプションを付けない場合、他の翻訳単位で使われるか否かはわからないままインライン展開するべきかの判断を下す必要があります。インライン展開をすべきと判断した場合、二つの関数にインライン展開されるためバイナリサイズが増大し、約17KBのバイナリとなります。

ここで、static指定子がついている場合（-DDECL_SPEC=staticと-DDECL_SPEC='static inline'とした場合）を考えてみると、他の翻訳単位からその関数が呼ばれることはない*3ため、実体を作る必要はありません。そのため、その翻訳単位で一回しか呼ばれないのならインライン展開してしまうのが最善となります。

一方、単独でinlineとしている場合、「複数の翻訳単位で定義されたとしても実体は一つになることが保証される」ことから、インライン展開しないで実体を作っておけば、他の翻訳単位で作った実体とリンク時にマージできる可能性があり、バイナリサイズを小さくする役に立てることができます。そのため、むしろインライン展開を抑制する働きを持つことがあります。

まとめ

inline指定子はもはやコンパイラにインライン展開を指示するものではありません。また、インライン展開を促進する効果があるとも限らず、むしろインライン展開を抑制する働きを持つことさえあります。この効果はstatic inlineとした場合には働かないため、ODR違反に細心の注意を払いつつstatic inlineとしないことが有効な局面もあり得るということを示しました。

補遺

inlineはむしろインライン展開を抑制する働きを持つことさえある、と書きましたが、これはあくまで特殊例です。現にコンパイラはinlineキーワードがついている関数についてインライン展開閾値を下げるということを行っているようです。

また、関数のマージはinlineを付けない限り発生しないわけではなく、コンパイラの最適化によって勝手に行われることがあるようです。しかも、関数ポインタをとっているのに勝手にマージされて同じアドレスになって困るなどの現象も発生することがあったらしいです（VC9で発生した例→RTTI を使わずに Boost.Any を実装する方法 - melpon日記 - HaskellもC++もまともに扱えないへたれのページ）。

*1:難しい用語ですが、分割コンパイルするときの一つのcppファイルだと理解すれば大体あっています。

*2:平たく言うと、ヘッダファイルで定義してソースファイルにインクルードする場合、inlineを付けないと実体が複数になって困りますが、inlineを付ければそういうことが起こらないということです。つまり、ヘッダファイルに定義を書くことができるようになるという利点があります。

*3:"その関数"自体が呼ばれることはあり得ますが、staticがついている以上、実体としては別の関数として扱われます。