有限体の実装 in C++

2018-02-10 cpp / コード置き場

※ 位数が素数である場合に限る

概要

C++の練習として、有限体を実装した。

数学

そんなに高度なことはやっていないので、事実だけ書く。

\(p\) を素数とすると、 \(\text{mod} \ p\) の演算について \(\{0,1,2,\ldots p-1\}\) は体をなす。
体 \(\{0,1,2,\ldots p-1\}\) の加算、減算、乗算の定義は明らか。
除算は逆元との乗算として定義する。
0以外の任意の元の逆元は拡張されたユークリッドの互除法を使って求めることができる。

実装

実装では次のことに注意した。

逆元を除算のたびに計算するのは無駄なので、逆元の表をクラスの static 変数として保持する
逆元の表を毎回の実行の冒頭で求めるのも無駄なので、 constexpr を使ってコンパイル時に求める
体の位数によってサイズの変わる逆元の表を constexpr にするために体の位数をテンプレート化する。
体の位数は素数でなことを static_assert で確かめる
適切なコンストラクタとキャスト演算子を整備することで体の元と数リテラル間の自然な演算ができるようにする

Code Snippet 1: finite_field.hpp

#include <array>
#include <cassert>
#include <ostream>

// x * x > p である x を返す
// あまり賢くない
constexpr int isqrt(int p) {
  int i = 2;
  for (; i * i < p; ++i) {
  }
  return i;
}

// p は素数か?
constexpr bool is_prime(int p) {
  for (int i = 2; i < isqrt(i); ++i) {
    if (p % i == 0) {
      return false;
    }
  }
  return true;
}

// i mod a を返す
constexpr int mod_pos(int i, int a) {
  int rem = i % a; // i < 0 だと rem < 0
  if (rem >= 0) {
    return rem;
  } else {
    return rem + a;
  }
}

// 拡張ユークリッド法でax + by = gcd(p, a) を満たす(x, y)を求める
constexpr std::pair<int, int> ext_euclid(int a, int b) {
  int r1 = a;
  int r2 = b;
  int x1 = 1;
  int x2 = 0;
  int y1 = 0;
  int y2 = 1;
  while (r2 != 0) {
    int r3 = r1 % r2;
    int q = r1 / r2;
    int x3 = x1 - q * x2;
    int y3 = y1 - q * y2;
    r1 = r2;
    r2 = r3;
    x1 = x2;
    x2 = x3;
    y1 = y2;
    y2 = y3;
  }
  return {x1, y1};
}

// arr[i] * i = 1 (mod P) となる arr を返す
template <int P>
constexpr std::array<int, P> construct_inverses() {
  std::array<int, P> arr = {};
  for (int i = 0; i < P; ++i) {
    arr[i] = ext_euclid(i, P).first;
  }
  return arr;
}

// 位数Pの有限体を表す (Pは素数)
template <int P>
class FiniteField {
  // Pが素数なことを確認
  static_assert(P > 1, "P in not prime");
  static_assert(is_prime(P), "P is not prime.");

 private:
  int m_val;
  // コンパイル時に計算される逆元の表
  constexpr static std::array<int, P> inverses = construct_inverses<P>();

 public:
  FiniteField() : m_val(0) {}
  FiniteField(const int &val) : m_val(mod_pos(val, P)) {}

  // このexplicitを抜くと恐ろしい量のエラーが出る
  explicit operator int() const { return m_val; }

  FiniteField<P> inv() const {
    if (m_val == 0) {
      throw std::runtime_error("division by zero");
    }
    return inverses[m_val];
  }

  // 以下は演算子のオーバーロード
  FiniteField<P> operator-() const { return FiniteField<P>(P - m_val); }

  FiniteField<P> &operator+=(const FiniteField<P> &other) {
    m_val = mod_pos(m_val + other.m_val, P);
    return *this;
  }

  FiniteField<P> &operator-=(const FiniteField<P> &other) {
    m_val = mod_pos(m_val - other.m_val, P);
    return *this;
  }

  FiniteField<P> &operator*=(const FiniteField<P> &other) {
    m_val = mod_pos(m_val * other.m_val, P);
    return *this;
  }

  FiniteField<P> &operator/=(const FiniteField<P> &other) {
    return *this *= other.inv();
  }

  // friendにして型の昇格を暗黙に行う
  friend FiniteField<P> operator+(const FiniteField<P> &lhs,
                                  const FiniteField<P> &rhs) {
    return FiniteField<P>(lhs.m_val + rhs.m_val);
  }

  friend FiniteField<P> operator-(const FiniteField<P> &lhs,
                                  const FiniteField<P> &rhs) {
    return FiniteField<P>(lhs.m_val - rhs.m_val);
  }

  friend FiniteField<P> operator*(const FiniteField<P> &lhs,
                                  const FiniteField<P> &rhs) {
    return FiniteField<P>(lhs.m_val * rhs.m_val);
  }

  friend FiniteField<P> operator/(const FiniteField<P> &lhs,
                                  const FiniteField<P> &rhs) {
    return lhs * rhs.inv();
  }

  friend bool operator==(const FiniteField<P> &lhs, const FiniteField<P> &rhs) {
    return rhs.m_val == lhs.m_val;
  }

  friend bool operator!=(const FiniteField<P> &lhs, const FiniteField<P> &rhs) {
    return rhs.m_val != lhs.m_val;
  }

  friend std::ostream &operator<<(std::ostream &rhs,
                                  const FiniteField<P> &lhs) {
    rhs << lhs.m_val;
    return rhs;
  }
};

以下のように使える。なお constexpr で複雑なことを行うので、コンパイルには c++17 の機能が必要になる。

#include "finite_field.hpp"
#include <iostream>

using F17 = FiniteField<17>;
using F8191 = FiniteField<8191>;

int main() {
  F17 x = 4, y = 20;

  // 四則演算
  assert(x + y == 7); // 右辺を F17(7) とする必要はない
  assert(x - y == 1);
  assert(x * y == 12);
  assert(x / y == 7);

  // リテラルとの四則演算
  assert(x + 20 == 7); // 型変換(昇格)は暗黙に行われる
  assert(4 - y == 1);
  assert(4 * 20 == F17(12));

  // 逆元
  for (int i = 1; i < 8191; ++i) {
    F8191 x = i;
    assert(x * x.inv() == 1); // この演算は定数時間で終わる
  }

  std::cout << "OK" << std::endl;
  return 0;
}

OK