CPUの浮動小数点演算性能やクロック周波数を測定する

CPUのカタログスペックはメーカーが公表しているので、それを参照することもできますが、本当にそうなのか確かめてみるという話です。 最近のIntelのCPUしか持っていないので、それでしか確かめていないませんが、四種類のCPUで試してみました。

CPUのクロック周波数

CPUのクロック周波数を求めることも、意外と簡単です。依存関係が直列につながっている一連の命令列を何クロックで処理できるかを計測し、それを実時間で割ればクロック周波数が分かります。 それ以外の命令も実行する必要があるので計算が狂ってきそうに思えますが、実際には並列かつアウトオブオーダーに実行されるため、分岐予測が当たる前提のもと、それらを考慮する必要はありません。

#include <iostream>
#include <thread>
#include <chrono>
#include <atomic>

std::atomic<bool> begin;
std::atomic<bool> end;
std::atomic<double> x;

void job() {
  double a = 0.0;
  while( !begin.load() );
  while( !end.load() ) a += 1.0;
  x.store(a);
}

int main() {
  std::thread th( &job );
  const auto startT = std::chrono::system_clock::now();
  begin.store( true );
  system( "sleep 5" );
  end.store( true );
  th.join();
  const auto endT = std::chrono::system_clock::now();
  const auto time = std::chrono::duration_cast<std::chrono::nanoseconds>( endT - startT ).count();
  std::cout << (long long)x << std::endl;
  std::cout << time << std::endl;
  std::cout << x*3 / time << " GHz" << std::endl;
}

コンパイルは、-std=c++11 -pthread -O3 -march=nativeで行います。

直列に実行すべき命令は、a += 1.0;の部分です。この部分は(SandyBridge以降のAVX命令に対応している場合)vaddsd命令となります。 この命令のレイテンシ(結果が出て次の命令で使えるようになるまでの時間)は3クロックなので、加算が行われた回数の三倍がかかったクロック数になります。 これをかかった時間で割ることにより、クロック周波数を求めることができます。

実測してみた結果はこれです。十回行ってもっともよかった結果を採用します。

  • Intel(R) Core(TM) i5 2540M(ターボブースト時3.30GHz)→3.13GHz
  • Intel(R) Core(TM) i5 7200U(ターボブースト時3.10GHz)→3.06GHz*1
  • Intel(R) Xeon(R) E5-2643 v3(ターボブースト時3.70GHz)→3.66GHz
  • Intel(R) Xeon(R) E5-2603 v3(1.60GHz)→1.58GHz

おおむねカタログスペックと同じ値を実験でも得ることができました。

CPUの浮動小数点演算性能

無駄な計算をするだけであれば、CPUの浮動小数点演算性能の理論値の90%程度を出すことは簡単です。 無駄な計算は、コンパイラの最適化が見抜けないように注意して書く必要があります。

#include <iostream>
#include <chrono>
#include <thread>
#include <atomic>
#include <cmath>

constexpr int ND = 32;
constexpr int NThreads = 4;
std::atomic<bool> ready;

void task( int N_loops ) {
    double x[ND];

    for( int i = 0; i < ND; ++i ) {
        x[i] = 1e-9 * i;
    }
    while( !ready.load() );
    for( int i = 0; i < N_loops; ++i ) {
        for( int j = 0; j < ND; ++j ) {
            x[j] = std::fma( x[j], x[j], x[j] );
        }
    }
    for( int i = 0; i < ND; ++i ) {
        std::cout << x[i] << " ";
    }
    std::cout << std::endl;
}

int main( int argc, char* argv[] ) {
    const int N = std::stoi( argv[1] );
    std::thread th[NThreads];
    for( auto& t : th ) {
        t = std::move( std::thread( &task, N ) );
    }
    const auto start = std::chrono::system_clock::now();
    ready.store( true );
    for( auto& t : th ) {
        t.join();
    }
    const auto end = std::chrono::system_clock::now();
    std::cout << N * ND*2. * NThreads / std::chrono::duration_cast<std::chrono::nanoseconds>( end - start ).count() << "GFLOPS" << std::endl;
}

コンパイルは、先ほどと同じオプションで行います。NThreadsは、そのCPUの動かせるスレッド数に応じて適宜変更してください。

std::fmaを計算するには、足し算と掛け算の二演算が必要であるため、2FLOPをこなしたことになります。

理論性能は、AVXの場合、ymmレジスタdoubleが4個入り、二つの命令が同時に実行可能であるので8FLOP/clock・Coreとなります。 また、AVX2の場合は、2FLOPに相当するfma命令が同時に二命令実行可能であるため、16FLOP/clock・Coreとなります。

実験してみた結果はこれです。同様に、十回行ってみて最も良いものを採用します。

  • Intel(R) Core(TM) i5 2540M(8FLOP/clock×3.30GHz×2Core=52.8GFLOPS)→2.09GFLOPS
  • Intel(R) Core(TM) i5 7200U(16FLOP/clock×3.10GHz×2Core=99.2GFLOPS)→94.99GFLOPS

SandyBridge世代のi5 2540MはAVXしか動かせず、AVX2で追加されたFMA命令に対応していないため、std::fmaがライブラリ関数コールとなり、非常に低速になります。 そこで、以下のように計算部分を変更します。

    for( int i = 0; i < ND; ++i ) {
        x[i] = 1 + 1e-9 * i;
    }
    while( !ready.load() );
    for( int i = 0; i < N_loops; ++i ) {
        for( int j = 0; j < ND; ++j ) {
            x[j] = x[j] * x[j];
        }
    }

この時、計算部分は以下のようにコンパイルされており、最適になっています。

 .LBB0_13:                               # =>This Inner Loop Header: Depth=1
 vmulpd %ymm7, %ymm7, %ymm7
 vmulpd %ymm6, %ymm6, %ymm6
 vmulpd %ymm3, %ymm3, %ymm3
 vmulpd %ymm0, %ymm0, %ymm0
 vmulpd %ymm1, %ymm1, %ymm1
 vmulpd %ymm2, %ymm2, %ymm2
 vmulpd %ymm4, %ymm4, %ymm4
 vmulpd %ymm5, %ymm5, %ymm5
 vmulpd %ymm7, %ymm7, %ymm7
 vmulpd %ymm6, %ymm6, %ymm6
 vmulpd %ymm3, %ymm3, %ymm3
 vmulpd %ymm0, %ymm0, %ymm0
 vmulpd %ymm1, %ymm1, %ymm1
 vmulpd %ymm2, %ymm2, %ymm2
 vmulpd %ymm4, %ymm4, %ymm4
 vmulpd %ymm5, %ymm5, %ymm5
 addl $-2, %ecx
 jne .LBB0_13

結果は、以下の通りです。

  • Intel(R) Core(TM) i5 2540M(8FLOP/clock×3.30GHz×2Core=52.8GFLOPS)→48.11GFLOPS

元のstd::fmaを使ったソースコードで出力されるアセンブリコードは、ここまで品質が良くなく、途中でレジスタ移動が発生しています。 イントリンシックを使って最適な機械語コードになるようにしてみます。

#include <iostream>
#include <chrono>
#include <thread>
#include <atomic>
#include <immintrin.h>

constexpr int NR = 16;
constexpr int ND = 64;
constexpr int NThreads = 12;
std::atomic<bool> ready;

void task( long long N_loops ) {
  double x[ND];
  __m256d ymm[NR];

  for( int i = 0; i < ND; ++i ) {
    x[i] = 1e-9 * i;
  }
  for( int i = 0; i < NR; ++i ) {
    ymm[i] = _mm256_load_pd( &x[i*4] );
  }
  while( !ready.load() );

  for( long long i = 0; i < N_loops; ++i ) {
    for( int j = 0; j < NR; ++j ) {
      ymm[j] = _mm256_fmadd_pd( ymm[j], ymm[j], ymm[j] );
    }
  }
  for( int i = 0; i < NR; ++i ) {
    _mm256_store_pd( &x[i*4], ymm[i] );
  }
  for( int i = 0; i < ND; ++i ) {
    std::cout << x[i] << " ";
  }
  std::cout << std::endl;
}

int main( int argc, char* argv[] ) {
  const long long N = std::stoll( argv[1] );
  std::thread ts[NThreads];
  for( int i = 0; i < NThreads; ++i ) {
    ts[i] = std::move( std::thread( &task, N ) );
  }
  const auto start = std::chrono::system_clock::now();
  ready.store( true );
  for( int i = 0; i < NThreads; ++i ) {
    ts[i].join();
  }
  const auto end = std::chrono::system_clock::now();
  std::cout << N * ND * 2. * NThreads / std::chrono::duration_cast<std::chrono::nanoseconds>( end - start ).count() << " GFLOPS" << std::endl;
}

これを使った時の結果は以下の通りです。

  • Intel(R) Core(TM) i5 7200U(16FLOP/clock×3.10GHz×2Core=99.2GFLOPS)→98.02GFLOPS
  • Intel(R) Xeon(R) E5-2643 v3(16FLOPS/clock×3.70GHz×6Core×2Socket=710.4GFLOPS)→669.8GFLOPS
  • Intel(R) Xeon(R) E5-2603 v3(16FLOPS/clock×1.60GHz×6Core=153.6GFLOPS)→152.2GFLOPS

最高の場合、理論性能の99%程度を出すことに成功しています。E5-2643v3ではそこまで達成できていませんが、全てのコアをターボクロック周波数にできない制約によるもの(つまり括弧内の計算による理論性能は正しくない)だと思われます。

*1:vaddsd命令のレイテンシは4クロックなので該当部分を変更したソースコードで実験

わりざんするアルゴリズム(その2.5)

三週間前の記事(わりざんするアルゴリズム(その2) - よーる)の割り算する回路・アルゴリズムを改良する話です。

わりざんするアルゴリズム(その1) - よーるわりざんするアルゴリズム(その2) - よーるで紹介した方法は、商と剰余を求めるのに常に整数型のビット幅+αサイクルかかっていました。つまり、9 / 7を計算するというように、実際に入っている値(の絶対値)が小さい場合でもそれだけのサイクル数がかかってしまうということです。

回復法の場合

回復法を使った除算法において、実際に立つ商が小さい場合、商の上位桁はずっと0が立ち続けることになります。筆算を行う場合、上位桁の0は普通書きません。つまり、上位桁に0を立て続ける部分は、実質的な計算が行われているわけではなく、無駄であるということになります。

実際にどの位に初めて商として1が立つのかは計算してみないとわからないですが、除数・被除数のleading zeros(最上位桁から連続する0の個数)を調べれば、初めて商として1が立つ桁か、その一つ上の桁を求めることはできます。これにより無駄な計算を極力省くことができます。

非回復法の場合

非回復法の場合、回復法のような単純な方法を使うことはできません。これは、商の表現が、各桁に \overline{1}(上線付きの1であり、-1を表します。以下形の似ているTで代用します)と1しか使えない(つまり、0が使えない)ためです。

しかし、全く解決策がないわけではありません。 商としてT111111111111TT1T == -27が立つ状況を考えてみます。立つ商が負でその絶対値が小さい場合、上位桁はT1111....という形になることが分かります(商が正でその絶対値が小さい時は、上位桁は1TTTT...という形になります)。この部分の計算はやや無駄であり、なんとか計算サイクル数を減らすことを考えます。

回復法での改良のように、途中から計算を始めることを考えます。もし上位桁に0が使えるのであれば、0000000000T1TT1T == -27のように表すことができます。したがって、三十二の位から商を立て始めればよいことになります。

その後、BSD表現から通常の二の補数表現に変換するところに二つの注意が必要です。

二の補数表現における符号拡張

回復法の場合は符号なしで計算していたため、計算していない上位桁は0で埋めればよかったのですが、非回復法の場合は符号付きで計算しているため、計算していない上位桁を何で埋めるのかが問題になります。

BSD表現では、最上位桁にTが立っていれば負、1が立っていれば正なので、それを利用して符号拡張用のマスクを作成します。

ビット演算テクニックによるBSD表現から二の補数表現への変換

もともとの非回復法では、quotient = (quotient << 1) | 1;というビット演算テクニックでBSD表現から二の補数表現へ変換していました。 この時、最上位ビットが左シフトによって暗黙に追い出されること込みでこの変換は成り立っています。

途中から計算を始めた場合、商の長さが本来のビット幅より短くなっているため、未計算の上位桁に追い出されることになります。 この追い出されたビットは不要であるため、適切にマスクを掛けて消すことが必要です。

実際、最上位桁がTの場合、それは商が負であることを示しているため、未計算部分は1で埋める必要があります。このとき、最上位桁のTは0で表されていますが、左シフトで追い出されてきたこの0も含めて1で埋める必要があります。 また、最上位桁が1の場合、それは商が正であることを示しているため、未計算部分は0で埋める必要があります。このとき、左シフトで追い出されてきた最上位桁の1も含めて0で埋める必要があります。

この部分の手順を間違えると答えが完全におかしくなるので注意してください。

以上の二つを踏まえて改良したアルゴリズムC++で記すと以下のようになります

改良されたアルゴリズムC++コード

以下のコードは-std=c++17オプションをつけてコンパイルします。 また、ライブラリとしてGitHub - lpha-z/lbl: header only fixed width integer libraryが必要です。

#include "lbl/xintN_t.hpp"

template<std::size_t N>
constexpr int countSignificantBits( lbl::intN_t<N> x ) {
    bool top_bit = x & 1ull<<(N-1);
    for( int i = N-2; i >= 0; --i ) {
        bool ith_bit = x & 1ull<<i;
        if( top_bit != ith_bit ) { return i+2; }
    }
    return 1;
}

template<std::size_t N>
constexpr auto nonrestoring_divider( lbl::intN_t<N> dividend, lbl::intN_t<N> divisor ) {
    lbl::intN_t<N> remainder = dividend < 0 ? -1 : 0;
    lbl:intN_t<N> quotient = 0; // any value is allowed

    int loop_count = countSignificantBits( dividend );
    for( int i = loop_count-1; i >= 0; --i ) {
        bool new_bit = dividend & 1ull<<i;
        if( (remainder < 0) != (divisor < 0) ) {
            quotient = (quotient<<1) | 0;
            remainder = ((remainder<<1) | new_bit) + divisor;
        } else {
            quotient = (quotient<<1) | 1;
            remainder = ((remainder<<1) | new_bit) - divisor;
        }
    }

    bool positive = quotient & 1ull<<(loop_count-1);
    quotient = (quotient << 1) | 1;
    if( positive ) {
       quotient &= (1ull << loop_count) - 1;
    } else {
        quotient |= ~((1ull << loop_count) - 1);
    }

    if( remainder == 0 ) {
        /* do nothing */
    } else if( remainder == divisor ) {
        assert( divisor < 0 );
        quotient += 1;
        remainder = 0;
    } else if( remainder == -divisor ) {
        assert( divisor > 0 );
        quotient -= 1;
        remainder = 0;
    } else if( (remainder < 0) != (dividend < 0) ) {
        if( (remainder < 0) != (divisor < 0) ) {
            quotient -= 1;
            remainder += divisor;
        } else {
            quotient += 1;
            remainder -= divisor;
        }
    }

    return std::pair {
        quotient,
        remainder
    };
} 

template<std::size_t N>
constexpr auto nonrestoring_divider( lbl::uintN_t<N> dividend, lbl::uintN_t<N> divisor ) {
    auto [quotient, remainder] = nonrestoring_divider<N+1>( static_cast<lbl::intN_t<N+1>>(dividend), static_cast<lbl::intN_t<N+1>>(divisor) );
    return std::pair {
        static_cast<lbl::uintN_t<N>>(quotient),
        static_cast<lbl::uintN_t<N>>(remainder)
    };
}

さらなる改良

先の方法では、商の最上位ビットは左シフトで追い出された後マスクを掛けられて消されます。しかし、符号拡張のための情報として利用されるため、商の最上位ビットが役に立っていないわけではありません。

そこで、商を立てるときにその処理をしてしまえば、後処理の部分で余計なことをする必要がなくなります。

具体的には、商の最上位ビットにTが立つ場合、それは商が負であることを意味しているため、符号拡張を見越してquotient = -1とします。 商の最上位ビットに1が立つ場合、それは商が正であることを意味しているため、符号拡張を見越してquotient = 0とします。

それを実装したのが以下のC++コードになります。

#include "lbl/xintN_t.hpp"

template<std::size_t N>
constexpr int countSignificantBits( lbl::intN_t<N> x ) {
    bool top_bit = x & 1ull<<(N-1);
    for( int i = N-2; i >= 0; --i ) {
        bool ith_bit = x & 1ull<<i;
        if( top_bit != ith_bit ) { return i+2; }
    }
    return 1;
}

template<std::size_t N>
constexpr auto nonrestoring_divider( lbl::intN_t<N> dividend, lbl::intN_t<N> divisor ) {
    lbl::intN_t<N> remainder = dividend < 0 ? -1 : 0;
    lbl:intN_t<N> quotient = 0; // any value (or non-initiaziled) is allowed

    int loop_count = countSignificantBits( dividend );
    for( int i = loop_count-1; i >= 0; --i ) {
        bool new_bit = dividend & 1ull<<i;
        if( (remainder < 0) != (divisor < 0) ) {
            quotient = i == loop_count - 1 ? lbl::intN_t<N>(-1) : (quotient<<1) | 0;
            remainder = ((remainder<<1) | new_bit) + divisor;
        } else {
            quotient = i == loop_count - 1 ? lbl::intN_t<N>(0) : (quotient<<1) | 1;
            remainder = ((remainder<<1) | new_bit) - divisor;
        }
    }

    quotient = (quotient << 1) | 1;

    if( remainder == 0 ) {
        /* do nothing */
    } else if( remainder == divisor ) {
        assert( divisor < 0 );
        quotient += 1;
        remainder = 0;
    } else if( remainder == -divisor ) {
        assert( divisor > 0 );
        quotient -= 1;
        remainder = 0;
    } else if( (remainder < 0) != (dividend < 0) ) {
        if( (remainder < 0) != (divisor < 0) ) {
            quotient -= 1;
            remainder += divisor;
        } else {
            quotient += 1;
            remainder -= divisor;
        }
    }

    return std::pair {
        quotient,
        remainder
    };
} 

template<std::size_t N>
constexpr auto nonrestoring_divider( lbl::uintN_t<N> dividend, lbl::uintN_t<N> divisor ) {
    auto [quotient, remainder] = nonrestoring_divider<N+1>( static_cast<lbl::intN_t<N+1>>(dividend), static_cast<lbl::intN_t<N+1>>(divisor) );
    return std::pair {
        static_cast<lbl::uintN_t<N>>(quotient),
        static_cast<lbl::uintN_t<N>>(remainder)
    };
}

ビット行列の行列積高速化に関するメモ

行列の要素が0 or 1であり、加法がbit or、乗法がbit andであるような束(加法がxorではないので、ブール環にはなっていません。また、束であることは特に使いません)上での行列積演算を考えます。

これが疎行列、つまり行列の成分がほとんど零であり、非零な要素が \Theta(N)個くらいしかないような行列でありれば、適当な圧縮表現上で取り扱うことで計算量は \Theta(N^{2})となります。

密行列であっても、N ≦256であれば、AVX2のvptest命令(二つのymmレジスタに乗っかっている256bit整数のbitwise andを取り、その結果が0であるかが分かる命令)を使えば、 \Theta(N^{2})になります(Nが有限である時に漸近記法を用いるのは厳密にはおかしいですが、記号の濫用ということでお許しください)。以下、N ≦256の場合のことを考えます。

N = 256の時、一つの行列を保持するのに必要なメモリは8KiBであり、三つの行列を保持しても24KiBであるため、L1キャッシュのサイズが32KiBであれば全ての情報が収まります。 その場合、行列積を速くする種々のテクニックを使う必要はあまりありません。 普通の二重ループ(最内ループはvptest命令一つになるので三重ループにはなりません)を書けば十分です。

行列のべき乗を計算する場合、もともとの行列が疎行列であってもべき乗を繰り返すうちに密行列になって行ってしまうので、どこかで手法を切り替える必要があります。

疎行列の手法と密行列の手法を切り替えるべき境界について調べたところ、行列の非零成分の数が、 N\sqrt{N}個くらいのところだということが分かりました。 つまり、行列の非零成分の数が N\sqrt{N}個よりも少なければ疎行列の手法を使った方が高速であり、それよりも多ければ密行列の手法を使った方が高速であるということです。

疎行列の手法の場合、演算結果が零になるようなビットにはアクセスしないのが高速化につながります。この効果は行列が疎である(非零な成分が少ない)ほど高まります。 一方、非零な要素の数が O(N)を超える場合、計算量は O(N^{2})に収まらないため、非零な要素の数が増えるにしたがって、非零な要素の数によらない密行列の手法が相対的に高速化します。

わりざんするアルゴリズム(その2)

先週の記事(わりざんするアルゴリズム(その1) - よーる)に引き続き、割り算する回路・アルゴリズムの話です。

非回復法(non-restoring division)

回復法の場合、「計算してみて、その結果に応じて……」という、ある意味試行錯誤に近いアルゴリズムでした。 非回復法では、被除数(部分剰余)と除数の符号の関係のみによって商として何が立つかが決まり、そのような試行錯誤的な手順はありません。

さらに、以下の特徴があります。

  • 二の補数で表された符号付き整数を取り扱うことができる。
  • BSD表現という、特殊な表現で商を求める(最後に通常の二の補数表現に変換する)。
  • この手法だけでは商と余りを正しく求められないことがあるため、最後に補正をするステップが必要。

Binary signed digit (BSD)表現

コンピュータ上では、多くの場合二進法、つまり二を底とする位取り記数法を使って数値を表現しています。これは、各位が0か1で表され、その重みが2nとなっている記数法です。

BSD表現は、重みが2nである点は同じですが、各位の数が異なり、各位は-1か1で表されます。-1は \overline{1}のように、上線付きの1で表すこともあります。以降、見た目の似ている"T"で代用します。

まず、BSD表現は普通の意味での記数法ではない点に注意します。そもそも、奇数しか表せず、偶数が表せません。 しかし、表せるならその表現方法は一意です(冗長な記数法ではありません)。

一方、陽に符号のようなものを考えなくても負の数を表せる記数法である点は優れています(二の補数表現もそのような記数法の一つです)。 たとえばT111 (=-8+4+2+1=-1)など、最上位桁がTである数は負の数を表しています。

また、各位の表現として0が使えないため、「記されていない上位桁や小数点以下の部分は0」という暗黙の取り決めが適用できません。 二の補数表現の場合*1、「記されていない上位桁は、最上位桁と同じ数字が入る」という"符号拡張"を行うことで、任意の桁数に拡張することができます。 BSD表現の場合、以下のように最上位のTT11...1に、11TT...Tに、それぞれすることで任意の桁数に拡張することができます。

      T11T1 (=-16+8+4-2+1=-5)
T11111111T1 (=-1024+512+256+128+64+32+16+8+4-2+1=-5)

     1T1TTT (=32-16+8-4-2-1=17)
1TTTTTT1TTT (=1024-512-256-128-64-32-16+8-4-2-1=17)

C++コード

以下のコードは-std=c++17オプションをつけてコンパイルします。

#include "lbl/xintN_t.hpp"

template<std::size_t N>
constexpr auto nonrestoring_divider( lbl::intN_t<N> dividend, lbl::intN_t<N> divisor ) {
    lbl::intN_t<N> remainder = dividend < 0 ? -1 : 0;
    lbl:intN_t<N> quotient = 0; // any value is allowed
    for( int i = N-1; i >= 0; --i ) {
        bool new_bit = dividend & 1ull<<i;
        if( (remainder < 0) != (divisor < 0) ) {
            quotient = (quotient<<1) | 0;
            remainder = ((remainder<<1) | new_bit) + divisor;
        } else {
            quotient = (quotient<<1) | 1;
            remainder = ((remainder<<1) | new_bit) - divisor;
        }
    }

    quotient = (quotient << 1) | 1;
    if( remainder == 0 ) {
        /* do nothing */
    } else if( remainder == divisor ) {
        assert( divisor < 0 );
        quotient += 1;
        remainder = 0;
    } else if( remainder == -divisor ) {
        assert( divisor > 0 );
        quotient -= 1;
        remainder = 0;
    } else if( (remainder < 0) != (dividend < 0) ) {
        if( (remainder < 0) != (divisor < 0) ) {
            quotient -= 1;
            remainder += divisor;
        } else {
            quotient += 1;
            remainder -= divisor;
        }
    }

    return std::pair {
        quotient,
        remainder
    };
} 

template<std::size_t N>
constexpr auto nonrestoring_divider( lbl::uintN_t<N> dividend, lbl::uintN_t<N> divisor ) {
    auto [quotient, remainder] = nonrestoring_divider<N+1>( static_cast<lbl::intN_t<N+1>>(dividend), static_cast<lbl::intN_t<N+1>>(divisor) );
    return std::pair {
        static_cast<lbl::uintN_t<N>>(quotient),
        static_cast<lbl::uintN_t<N>>(remainder)
    };
}

解説

非回復法の場合、回復法とは逆に符号付き整数を取り扱うのが基本であるため、符号なし整数の場合はN+1ビットの符号付き整数用割り算器を呼び出します。Nビットの符号付き整数は、N+1ビットの符号付き整数で表せるため、それを利用します。

非回復法の割り算器は大きく分けて二つのステージに分かれます。for文で書かれた一ビットずつ商を求める部分と、if-else文で書かれた商と余りを補正する部分です。

remainderの初期化は、以下でremainderの符号を確認する部分がありますので、それに備えて符号拡張された表現になるよう、0-1を入れておきます。 quotientの上位ビットは全く使われず、かつ以下の手順ですべてのビットが左シフトで追い出されるため何が入っていても構わないのですが、C++的には演算前に初期化しておかないと未定義動作ですので適当な値で初期化しました。ハードウェアに実装する場合は、変なゴミが入っていても問題なく動作するということです。

一ビットずつ商を求める部分では、余りと除数の符号を比較し、一致するか否かでそのサイクルでの動作が決まります。つまり、余りと除数の符号が同じであれば商として1を立てて余りから除数を引く、余りと除数の符号が異なれば商としてTを立てて余りに除数を足す、という動作です。Tはquotient変数の中では0として表現しておきます。

この動作は、その桁までで最近接丸めを行った商を求める動作になっています(回復法の場合、その桁までで切り捨てた商を求める動作になっているのと比較してください)。最近接丸めの場合、ぴったり真ん中だったときの動作が問題になります。それがどういう時に発生するかというと、その桁まででの余りが0になった場合に発生します。その場合でも何かの商を立てないといけないわけですが、0は二の補数表現における"符号"が"正"の側に属していますので、除数が正であれば1、除数が負であればTが立つことになります。

このように変な商を立てるのはまずいことになりそうですが、実際は安全です。ループ内では、 \rm -|divisor| \le remainder \lt |divisor|が不変条件として常に成立しています。

quotient = (quotient << 1) | 1;の部分は、BSD表現を通常の二の補数表現に変換する部分です。技巧的ですが、以下のように解釈することができます。BSD表現における、1になっている桁の寄与はquotient分です。また、Tになっている桁の寄与は、-(~quotient)分になります。つまり、BSD表現でのquotientは、二の補数表現ではquotient + (-(~quotient))に対応する値であるということが分かります。ここで、ビット演算テクニックを考えると、-(~q) == q+1であることが分かります。よって、quotient + (-(~quotient)) == quotient + quotient +1となります。これを単純化すると(quotient << 1) | 1になります。

最後に商と余りを補正する部分を見ていきます。この時点で、商は最近接奇数に丸めた値、余りはそれに対応する(商が奇数という制約の上での)絶対値最小剰余になっています。これをC99言語流の商と余りの制約に補正していきます。

  • 余りが0、つまり最終桁でぴったり割り切れた場合には補正は不要です。
  • 余りがぴったり除数と一致するというのは、本来商が偶数で割り切れるのが正しい時、かつ除数が負の場合に発生します。割り切れた次の位での動作では、商として0が立てられず、最近接丸めの動作から、除数が負の時にはTが立ちます。それ以下の桁では、常に余りと除数が一致し、商としてずっとTが立ちます(ループ内不変条件の等号が成立する状況です)。つまり、...T11111のような商が立っていることになります。この場合、1を足すことで正しい商と余りに戻します。
  • 余りがぴったり除数に負号をつけたものに一致するというのは、本来商が偶数で割り切れるのが正しい時、かつ除数が正の場合に発生します。その原理は上述のものと同一です。この場合1を引くことで正しい商と余りに戻します。
  • それ以外の場合、被除数と余りの符号が一致しなければ一致する側に戻します。あるいは、小数点以下の計算を行うことで、最近接偶数を求めているという解釈もできます(割り切れず、かつ商として最近接奇数が不適切ならば、最近接偶数は一意でそれが適切な商となります)。BSD表現では偶数は表せませんが、小数点以下の計算を無限に行えば、循環小数として表現可能です。偶数であれば小数点以下はすべて同じ数字が表れます*2。そのため、小数点以下第一位の桁に立つ商を求めればよいことになります。そのため、(remainder < 0) != (divisor < 0)という、商を求めるforループで使った条件式と同じものが再出現します。

特性

回復法と同様、一ビットずつ商を求めている特性から、ビット長Nに対して少なくともNクロックサイクルはかかる除算器になります。ただし、補正の部分を商を求めているサイクルと同一サイクルで行うのは無理があるので、普通はN+1クロックサイクルかかるはずです。また、符号なし整数の場合、N+1ビットの符号付整数としてこの除算器に入力するため、N+2クロックサイクルかかることになります。

FPGAに実装する場合、ループ部分の回路は典型的なFPGAの論理セルと非常に相性が良いです。典型的なFPGAの論理セルは、LUT→FA→FFという構成になっています。ここで、LUTはLook-up-tableの略で、組み合わせ回路を実現する素子、FAはFull adderの略で加算器、FFはFlip-flopの略で記憶素子です。if( (remainder) < 0) != (divisor < 0) )の部分がLUTに対応、その後の加算がFFに対応、その後に実行される結果を代入する部分がFFに対応、となっているためです。

*1:二の補数表現では各位の表現として0が使えますが、最上位桁の重みだけ -2^{N-1}と一様でないため、任意の桁数に拡張する際にはやはり暗黙のルールが適用できず、特別な操作が必要です。

*2:実数を小数点以下が無限に続く位取り記数法で表現すると1.0000...=0.9999...のように二通りの表現がありえますが、BSD表現でもそのようになります(二の補数表現でももちろん同様です)。補正前の商が正しい商より小さい奇数だったか大きい奇数だったかにより、その二つの表現のどちらになるかが決まります。

わりざんするアルゴリズム(その1)

なんか最近割り算する回路を書いていろいろ知ったのでまとめていきます。

以下、割り算とは以下のものをさすことにします。

  1. 商と余りを両方求める
  2. 符号付きの場合、商はゼロへ丸め、余りは被除数と同じ符号とする

2は、数学的な定義とは異なります。C99で定義された除算の定義に一致します。

回路の説明というよりは、先週作った任意幅整数型lbl::uintN_t<N>lbl::intN_t<N>GitHub - lpha-z/lbl: header only fixed width integer library)を使ってC++アルゴリズムを書き起こしたものの紹介です。

回復法(restoring division)

小学校で習うような、筆算で求める方法です。

「回復」の名は、いったん除数を引いてみて負になってしまうようであれば除数を足すことで元の状態に"戻す"、というところからきているようです。 ソフトウェア的に実装するならともかく、ハードウェアで実現するのであれば、引く前の値を持っておけばよいだけなので、この名前はやや違和感があります。

#include "lbl/xintN_t.hpp"

template<std::size_t N>
constexpr lbl::uintN_t<N> safe_abs( lbl::intN_t<N> x ) {
    return x < 0 ? -x : x;
}

template<std::size_t N>
constexpr auto divider( const lbl::uintN_t<N> dividend, const lbl::uintN_t<N> divisor ) {
    lbl::uintN_t<N> remainder = 0;
    lbl::uintN_t<N> quotient = 0;
    for( int i = N-1; i >= 0; --i ) {
        quotient <<= 1;
        remainder <<= 1;
        remainder |= dividend & 1ull<<i ? 1 : 0;

        const lbl::uintN_t<N+1> test_sub = static_cast<lbl::uintN_t<N+1>>(remainder) - static_cast<lbl::uintN_t<N+1>>(divisor);

        const bool success = static_cast<lbl::intN_t<N+1>>(test_sub) >= 0;

        quotient |= success ? 1 : 0;
        remainder = success ? static_cast<lbl::uintN_t<N>>(test_sub) : remainder;
    }
    return std::pair {
        quotient,
        remainder
    };
}

template<std::size_t N>
constexpr auto divider( const lbl::intN_t<N> dividend, const lbl::intN_t<N> divisor ) {
    const bool quotient_sign = dividend < 0 ^ divisor < 0;
    const bool remainder_sign = dividend < 0;
    const auto [quotient, remainder] = divider_impl( safe_abs(dividend), safe_abs(divisor) );
    return std::pair {
        quotient_sign ? -static_cast<lbl::intN_t<N>>(quotient) : static_cast<lbl::intN_t<N>>(quotient),
        remainder_sign ? -static_cast<lbl::intN_t<N>>(remainder) : static_cast<lbl::intN_t<N>>(remainder)
    };
}


reminder |= dividend & 1ull<<i ? 1 : 0;の部分は、筆算で「次の桁をおろしてくる」部分に相当します。dividendは他の部分で使われないので破壊しても問題なく、dividend自体をシフトすることによって、ハードウェアは簡素化しますが、見やすさを優先しました。そのほか、最初の方はremainderの上位桁が使われていないことを利用すると、実はdividendと共存可能で、記憶素子をケチることができます。

要するに、引いてみても正だったら商として1を立て、そうでなければ0を立てる、というだけのアルゴリズムです。

引いてみて正かどうかを判定する部分は、オーバーフローしないように、結果を一桁多く求める必要があります(ボローが発生するかの確認が必要です)。reminderdivisor自体はN桁のまま計算してよいのですが、C++ではうまく書けなかったのでこのような形になっています。

for文で書いてある部分は、ハードウェアでは普通、複数サイクルかけて計算を行うことになるでしょう。 for文を回る回数はN回であり、各ループ内で一回加算器*1を使っています。つまり、1クロックに加算器が一回しか使えないような状況*2では、割り算を終わらせるのにNクロックかかるということになります。


符号付きの除算を行う場合、符号がついていると面倒なことになりますが、被除数・除数ともに絶対値を取って計算をした後、最後の段階で符号を正しく付け直すと簡単になります。

絶対値を取ったり負号をつけたりするのにも加算器が必要なことを考えると、素直に作るとN+2クロックかかるはずです。少し工夫すると、最終サイクルでの符号調整はループの中に組み込めるので、ループの最後の回では処理を変えることにより、余分な2クロックのうち1クロックは除去可能です。

また、ループの最初の回で立つ商が1であるパターンは限られているので、そのパターンを書き下せば、加算器を使わずに商を立てることができます。それと並行して絶対値を取る処理を行えば、余分なもう1クロックを削ることができます。そこまでする必要があるかは疑問ですが……。


回復法で除算を行う場合、加算器の結果がマルチプレクサの入力になりますが、典型的なFPGA (Field programmable gate array)に実装する場合、このことが大きな弱点になります。 典型的なFGPAの論理セルでは、組み合わせ回路を実現する LUT (Look-up table) の直後にフルアダーが配置してあります。マルチプレクサはLUTで実現、加算器はフルアダーで実現することになりますが、順序が逆なため二つの論理セルを使わないと実現できないことになります。

*1:加算器と言っていますが、専ら減算に使っていますね……。

*2:長大な桁数の加算器は(繰り上がりが連鎖するような、最悪な場合に)遅延が大きく、直列に複数回使うと1クロック内での最大遅延が大きくなりすぎ、クロック周波数が下がる原因となりえます。

半端なbit長の、固定幅の整数型

C99/C++11で追加された、固定幅の整数型(uint8_tint32_tのようなものたち)は特定のビット幅で計算したいときに大変便利です。何ビットの整数が要求されているかを明示することで移植性が高まります。

また、uint16_tは必ずオーバーフローせず、結果の上位ビットが無視されるだけの結果になるなど、ハードウェアに近いことをやりたいときにも便利です。しかし、この用途ではint13_tuint42_tが欲しくなることはないでしょうか。

大きめのbit幅の整数型を用い、適宜ビットマスクを掛けることで所望の動作をシミュレートすることはできますが、ユーザーコードで行うのは不便です。 当然そのようなライブラリが求められますが、少し探した感じではそのようなライブラリは見当たりませんでした。

そこで自作してみたというのが今日のお話です。

GitHubで公開してありますので、興味があれば使っていただけると嬉しいです。間違いなどの指摘も歓迎です。

github.com

以下、このライブラリを作るにあたって、移植性と効率性を両立するために工夫した点です。

下位bitを使ってシミュレートするのではなく、上位bitを使ってシミュレートする

単純な発想では、uint13_tをシミュレートするためには、uint16_tの下位13bitを用いる方法が考えられます。つまり、演算の後に常に0x1fffとのビットごと論理積を取りその結果を保持するという実装です。

この方法は、標準にある整数型に変換するときのコストが全くかからないというメリットがありますが、毎演算ごとにビットごと論理積が発生するというデメリットがあります(コンパイラがどのくらい見抜いてくれるかは確認していないですが……)。

そうではなく、あえてuint16_tの上位13bitを使うという方法を使っています。この方法を使うと、上位ビットが無視される状況が何の追加演算もなしに得られます。また、表現は単に定数倍(uint13_tuint16_tでシミュレートする場合は、3bit分なので8倍)されているだけであるので、大小比較は特に工夫せずにそのまま行って問題ありません。

代わりに、右シフトの場合は下位ビットが0になるように少しだけ補正が必要です。また、乗除算はビットマスクを掛ける処理の代わりにオフセットをずらす分の処理が少し増えます(ライブラリの実装では、左右のオペランドのどちらもオフセットを取り除き通常の表現にした後計算し、オフセットを戻すという処理を行っています。オフセットを取り除く処理とオフセットを戻す処理は打ち消す合うはずなので冗長ですが、コーナーケースが怖く、安全な実装に逃げました)。

ナイーブな実装と比べ、標準にある整数型との変換コストがかかる点だけは劣りますが、シフト演算一回分であり、頻繁な変換は普通は発生しないと考えられるため、問題ないものと考えました。

safe_abs

safe_absC++で安全に絶対値を求める - よーる で紹介した、未定義動作であるオーバーフローを起こさずに絶対値を取る関数です(詳細は以前の記事を参照ください)。

bit_cast

bit_castは移植性のある方法で、符号なし整数型を符号あり整数型(二の補数表現)に変換する(ビット表現を解釈しなおす)関数です。単にstatic_castを使う場合、変換後の値が負になる場合、処理系定義の動作となります。std::memcpyC++20で入る予定のstd::bit_castを用いても同じ動作が実現でき、そのような問題はなくなりますが、constexpr化できなくなります。

この関数では、uintN_tintN_tの変換において値が変わらない場合は移植性が保たれることを利用し、

  • 値が変わらない範囲(つまり、0x0000...から0x7fff...の範囲)にある場合はそのまま変換
  • 値が変わってしまう範囲(つまり、0x8000...から0xffff...の範囲)にある場合は、0x8000...を引いた値を求め、その値を変換(0x0000...から0x7fff...の範囲にあるので変換は移植性がある)したのち、0x8000...に相当する値(intN_tの最小値にあたる値)を加算する(この加算ではオーバーフローは発生しえない)

という方法で移植性を担保しています。最適化コンパイラを用いた場合、この処理は消失します。

算術シフト

符号あり整数型の右シフトの挙動は処理系定義となっています。x86にはsar (Shift arithmetic right) 命令が存在するため(他の多くのアーキテクチャにも算術シフト命令はあります)、多くの処理系では算術シフトとなるはずですが、やはり移植性のあるものではありません。 このライブラリを用いた場合、必ず算術シフトとなります。

なるべくオーバーヘッドを減らしたいため、少なくともx86のg++やclang++ではsar命令にコンパイルされるようなコードを書きたいところですが、うまくいきませんでした。

移植性よく、しかし既存のx86処理系にsar命令を推論させる方法は、試した限りでは存在しませんでした。コンパイラsar命令を出力するのは、

  1. 符号あり整数型を右シフトした場合
  2. 定数での除算の場合

に限られるようです。特に、以下のようなコードであっても最適化されてsar命令になることはありませんでした。

constexpr uint64_t sar1( uint64_t x ) noexcept {
    const uint64_t sign = x & 0x8000'0000'0000'0000;
    return sign | (x>>1);
}

1の場合を移植性よく使うことはできませんので、2の場合をうまく活用する方法を考えたいところですが、右算術シフトと除算では、被除数が負の時の丸め方向が異なります(右算術シフトは負の無限大への丸め、除算はゼロへの丸め)。そのため、正負で場合分けするコードが出現してしまいます。

この定数除算がsar命令を使ったコードに変換される動作は、コンパイラバックエンドで行われているようです。つまり、コンパイラフロントエンドでの情報(特に、未定義動作が起こる値が来ているかどうか)の情報が失われてしまっている点が厄介です。

定数除算がsar命令に変換されるのはほぼ定型文として変換されているようで、正の数が来ないはず、などの情報を活用できていません。コンパイラフロントエンドでの情報が残っていれば、定型文として変換した後にあり得ない分岐が取り除かれる最適化がかかるはずですが、失われているためにそれができていません。

おわりに

このライブラリはあまりビット効率のことを考えていないので、uintN_t<3>などとしても64bit使ってしまっています。8bitで十分なので、そういった点は改良の余地がありますが、あまり気にしないことにしました。

また、通常の整数と混ぜて使いたいためexplicitがほとんどつけていません。意図せぬ暗黙変換が発生してわけのわからないことになることが発生しそうですが、利便性を優先しました。この辺は好みが分かれるところですが、あまり新設ではないかもしれません。

再びになりますが、興味があればぜひ使ってみてください。間違いなどの指摘も歓迎です。