C++で標数2の体

概要

　本記事では標数2の体をC++で定義する際の要点について概説する．異なる拡大による体を型で静的に区別できるような仕組みを作り，C++における自然で美しい体の実装を目指す．

動機

　夏休み中，学友がC++を学ぶために勉強会を開いた．勉強会は様々な数学的対象をC++で表現する内容だった．本記事は勉強会の一回で扱った内容の焼き直しである．

準備

群

　ある集合と演算が次の1~3の条件を満たすとき，集合と演算の組を群という．また，1~4の条件を満たすとき，可換群またはアーベル群という．考えている演算が明らかなとき単にを群と呼ぶ．

1. 演算の結合性 :
2. 単位元の存在 :
3. 逆元の存在 :
4. 交換法則 :

　演算が加法的である場合，の逆元をと書き，演算が乗法的である場合，をと書く．

環

　ある集合と加法の組が可換群であり，かつ以下の1~3の条件を満たすとき，と乗法の組を環という．また，1~4の条件を満たすとき，可換環という．

1. 乗法の結合性 :
2. 乗法単位元の存在 :
3. 分配律 : かつ
4. 交換法則 :

　環の乗法単位元について，となるが存在するとき，その最小の値をの標数とする．また，そのようなが存在しないとき，の標数をとする．

イデアル

　環の部分集合が以下の1~3の条件を満たすときをイデアルという．1~4の条件を満たすときを素イデアルという．

1. 
2. 
3. 
4. 

　例えば任意の整数の倍数全体の集合は，整数全体と通常の加法，乗法からなる環のイデアルとなっている．また，任意の素数の倍数全体の集合は素イデアルとなっている．可換環にイデアルを求めることができれば，をで割った余りからなる新たな環を定義できる．について，ならばとする．この同値関係による同値類をと書くと加法 : と乗法 : によって剰余環が定義できる．

多項式環

　可換環や体を係数とする多項式環を定義する．を変数とし，可換環を係数とする一変数多項式環をと書く．多項式の各項のの次数のうち最高の数をとし，の次数がである，またはは次多項式であるという．の次の項の係数がである多項式をモニック多項式であるという．定数でない多項式が，定数でない多項式の二つの積で表せないときを既約多項式という．

体

　ある集合が可換環であり，かつが乗法に関して可換群をなすとき，を体という．つまり体は四則演算を自由に行える代数系である．環の素イデアルで割った剰余環は体となる．

体の代数拡大

　体の拡大とは，体との組のことをいう．がの根であるとき，を上代数的であるという．また，の拡大体の任意の元が代数的であるとき，を代数拡大という．代数拡大の例にがある．はの根をに追加した体である．また，代数拡大でない拡大の例にがある．などは上代数的でない．

　の既約多項式の根をに加えるとを拡大できる．そのように拡大した体の任意の元はとを用いてと表され，はのベクトル空間となる．のベクトル空間としてのの次元数を拡大次数という．

標数2の体のビット表現

　整数環の部分集合はの素イデアルである．剰余環は標数の体になる．の拡大体の元は，上既約な次多項式の根とによって以下のように表される．

　ここで上式のをビットの整数の各ビットに対応させるとの元を2進数として表現できる．また，を変数に置き換えればの多項式も2進数として表現できる．

実装

型の定義

　型は値に意味を与える．あるビット表現が単なる整数なのか，あるいは標数2の体の元を表現しているのか，あるいは標数2の体の元を表現しているとしてどの体の元なのかを適切に表現し計算するために，異なる体には異なる型を与える必要がある．C++ではテンプレート引数に整数引数を取れるため，体を拡大する原始多項式のビット表現をテンプレート引数で与えることができる．本記事では簡単のために拡大次数がののみ扱うこととするが，任意次数の拡大であっても理論は同じである．

　原始多項式のビット表現Polynomialをテンプレート引数で受け取れば，異なる原始多項式による異なる拡大体を異なる型で表現することができる．また，原始多項式のビット表現では最高次数の項がビット数によって常に一定なので省略されている．例えばT = unsigned charなら，Polynomialで直接表現されているのは以下の項のみで，の項は常に係数がである．

template<class T, std::enable_if_t<std::is_unsigned_v<T>, T> Polynomial>
class gf {
    T value_;
};

　ここに元のビット表現を格納するメンバ変数を定義し，基本的な演算を行う演算子を追加していく．

加法

　の加法はと同様に排他的論理和演算で表すことができる．以下の表のようにの加算をビット毎 (=元の式の項ごと) に行えばよい．加法の逆演算についてはであるため，加法と全く同じ演算で逆演算を行える．

乗法

　の乗法はのを原始元に置き換えた乗法と変わらないが，元の多項式基底表現の次数がを超える元は原始多項式で割って余りをとらなければコンピューターで表現できる桁数を超えてしまう．乗法の最も素朴な実装は，筆算をするように次数の低い項から掛けていき，ビット表現が型でサポートされる最大値を上回るとき，すなわち計算途中でを加算するとき，その代わりに，原始多項式の最高次数の項を移行した残りの辺のビット表現を加算するものである．例えばの乗法を考える．ただし，原始多項式をとする．だが，このままビット表現にすると101になり，2ビットで表せる数を超える．このままでは決まったサイズで表すことができないので，を等しい値で置き換える．なので，よりが求まる．

template<class T, std::enable_if_t<std::is_unsigned_v<T>, T> Polynomial>
class gf {
public:
    gf() = default;
    gf(const gf&) = default;
    explicit gf(T value)
        : value_{value}
    {}
    friend gf operator+(gf a, gf b) noexcept {
        return gf{a.value_ ^ b.value_};
    }
    friend gf operator-(gf a, gf b) noexcept { return a + b; }
    friend gf operator*(gf a, gf b) noexcept {
        gf res{};
        while(b.value_) {
            res.value_ ^= (b.value_ & 1 ? a.value_ : 0);
            a.value_ = !(a.value_ & (1 << (8*sizeof(T) - 1))) ? a.value_ << 1 : ((a.value_ << 1) ^ Polynomial);
            b.value_ >>= 1;
        }
        return res;
    }

private:
    T value_;
};

カラツバ法やSchönhage–Strassen法といった乗法の高速化手法もある．

除法

　除法は割る数の逆元を掛けることで実現する．有限群では任意の元についてとなる．体は零元を除いて乗法に関して群をなすので，においても任意の元の乗はその元の逆元となる．を愚直に計算してもよいが，大きな体では逆元計算に時間がかかるようになる．べき乗は次数の2進展開をもちいて高速化できる．例えばならば，25を2進展開し1 1001とし，となる．8乗はと計算すれば掛け算は3回で済む．C++で書けば以下のようになる．

template<class T>
constexpr T pow(T a, unsigned int e) {
    T r{1};
    while (e) {
        if (e & 1) {
            r *= a;
            --e;
        }
        else {
            a *= a;
            e >>= 1;
        }
    }
    return r;
}

　最終的なソースコードは以下のようになる．

#include <type_traits>
#include <limits>

template<class T>
constexpr T pow(T a, unsigned int e) {
{
    T r{1};
    while (e) {
        if (e & 1) {
            r *= a;
            --e;
        }
        else {
            a *= a;
            e >>= 1;
        }
    }
    return r;
}

template<class T, std::enable_if_t<std::is_unsigned_v<T>, T> Polynomial>
class gf {
public:
    gf() = default;
    gf(const gf&) = default;
    explicit gf(T value)
        : value_{value}
    {}
    constexpr gf inv() const noexcept {
        return pow(*this, std::numeric_limits<T>::max() - 1);
    }
    constexpr explicit operator T() const noexcept { return value_; }
    friend gf operator+(gf a, gf b) noexcept {
        return gf{a.value_ ^ b.value_};
    }
    friend gf operator-(gf a, gf b) noexcept { return a + b; }
    friend gf operator*(gf a, gf b) noexcept {
        gf res{};
        while(b.value_) {
            res.value_ ^= (b.value_ & 1 ? a.value_ : 0);
            a.value_ = !(a.value_ & (1 << (8*sizeof(T) - 1))) ? a.value_ << 1 : ((a.value_ << 1) ^ Polynomial);
            b.value_ >>= 1;
        }
        return res;
    }
    friend gf operator/(gf a, gf b) noexcept {
        return a * b.inv();
    }
private:
    T value_;
};