一年生の夢の実装 in C++

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

式

\begin{align*} (x + y)^{p} = x^{p} + y^{p} \end{align*}

を一年生の夢 (Freshman’s dream) という。

普通の多項式としてはこの式は間違っているが、有限体上の多項式としてみれば合っている。つまり、 \(p\) を素数とすると

\begin{align*} \forall x, y \in \mathbb{F}_{p}, \ (x + y)^{p} = x^{p} + y^{p} \end{align*}

が成り立つ。

前回の記事で有限体を実装したので、その上に多項式を実装してこの事実を確かめるコードを書いた。

多項式

まず多項式を実装した。特に難しい部分はなかったが、次数が大きくなったときと係数が0か1のときに operator<< が破綻しないような工夫をした。また operator+ と operator+= などでコードが使いまわせる部分についてはなるべく使いまわした。

Code Snippet 1: polynomial.hpp

#include <algorithm>
#include <initializer_list>
#include <ostream>
#include <vector>

template <typename Coeff>
class Polynomial {
 private:
  std::vector<Coeff> m_coeffs;
  // m_vec末尾の0を削除する
  Polynomial &delete_trailing_zeros() {
    while (m_coeffs.size() > 1 && m_coeffs[m_coeffs.size() - 1] == 0) {
      m_coeffs.pop_back();
    }
    return *this;
  }

 public:
  static bool print_zero;

  Polynomial(std::vector<Coeff> coeffs)
      : m_coeffs(coeffs.size() == 0 ? std::vector<Coeff>{0} : coeffs) {
    delete_trailing_zeros();
  }

  Polynomial() : m_coeffs(std::vector<Coeff>{0}) {}

  Polynomial(const Coeff &c) : Polynomial(std::vector<Coeff>{c}) {}

  Polynomial(std::initializer_list<Coeff> coeffs)
      : Polynomial(std::vector<Coeff>(coeffs)) {}

  size_t degree() const { return m_coeffs.size() - 1; }

  Coeff operator[](size_t i) const { return m_coeffs[i]; }

  friend bool operator==(const Polynomial &lhs, const Polynomial &rhs) {
    return lhs.m_coeffs == rhs.m_coeffs;
  }

  friend bool operator!=(const Polynomial &lhs, const Polynomial &rhs) {
    return lhs.m_coeffs != rhs.m_coeffs;
  }

  Polynomial &operator+=(const Polynomial &other) {
    m_coeffs.resize(1 + std::max(degree(), other.degree()));
    for (size_t i = 0; i <= other.degree(); ++i) {
      m_coeffs[i] += other[i];
    }
    return delete_trailing_zeros();
  }

  Polynomial &operator-=(const Polynomial &other) {
    m_coeffs.resize(1 + std::max(degree(), other.degree()));
    for (size_t i = 0; i <= other.degree(); ++i) {
      m_coeffs[i] -= other[i];
    }
    return delete_trailing_zeros();
  }

  Polynomial &operator*=(const Polynomial &other) {
    *this = *this * other;
    return *this;
  }

  friend Polynomial operator+(const Polynomial &lhs, const Polynomial &rhs) {
    auto temp(lhs);
    return temp += rhs;
  }

  friend Polynomial operator-(const Polynomial &lhs, const Polynomial &rhs) {
    auto temp(lhs);
    return temp -= rhs;
  }

  friend Polynomial operator*(const Polynomial &lhs, const Polynomial &rhs) {
    std::vector<Coeff> answer(lhs.degree() + rhs.degree() + 1);
    for (size_t i = 0; i <= rhs.degree(); ++i) {
      for (size_t j = 0; j <= lhs.degree(); ++j) {
        answer[i + j] += rhs[i] * lhs[j];
      }
    }
    return {answer};
  }

  // このイテレータは今回の例では使っていない
  typename std::vector<Coeff>::const_iterator begin() const {
    return m_coeffs.cbegin();
  }

  typename std::vector<Coeff>::const_iterator end() const {
    return m_coeffs.cend();
  }

  friend std::ostream &operator<<(std::ostream &rhs, const Polynomial &lhs) {
    // 次数が 0 or 1, 係数が 0 or 1 の場合でそれぞれ場合分けをしている
    // print_zero が true の場合にのみ係数の 0 が出力される

    if (lhs.degree() == 0) {
      rhs << lhs[0];
      return rhs;
    }

    rhs << lhs[0];
    if (lhs[1] != 0 || print_zero) {
      rhs << " + ";
      if (lhs[1] != 1) {
        rhs << lhs[1];
      }
      rhs << "x";
    }
    if (lhs.degree() == 1) {
      return rhs;
    }

    for (unsigned i = 2; i <= lhs.degree(); ++i) {
      if (lhs[i] != 0 || print_zero) {
        rhs << " + ";
        if (lhs[i] != 1) {
          rhs << lhs[i];
        }
        rhs << "x^" << i;
      }
    }
    return rhs;
  }

  Polynomial &pow(unsigned n) {
    Polynomial p{1};
    for (unsigned i = 0; i < n; ++i) {
      p *= *this;
    }
    *this = p;
    return *this;
  }
};

template <typename T>
bool Polynomial<T>::print_zero = false;

確認

上の多項式クラスを使って以下のようにして１年生の夢が成り立つことを確かめた。

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

int main() {
  // 位数が 2, 3, 11, 101, 8192 の場合について確認する
  using P2 = Polynomial<FiniteField<2>>;
  std::cout << "In F2, (1 + x)^2 = " << P2{1, 1}.pow(2) << "\n";

  using P3 = Polynomial<FiniteField<3>>;
  std::cout << "In F3, (1 + x)^3 = " << P3{1, 1}.pow(3) << "\n";

  using P11 = Polynomial<FiniteField<11>>;
  std::cout << "In F11, (1 + x)^11 = " << P11{1, 1}.pow(11) << "\n";

  using P101 = Polynomial<FiniteField<101>>;
  std::cout << "In F101, (1 + x)^101 = " << P101{1, 1}.pow(101) << "\n";

  using P8291 = Polynomial<FiniteField<8291>>;
  std::cout << "In F8291, (1 + x)^8291 = " << P8291{1, 1}.pow(8291) << "\n";

  return 0;
}

In F2, (1 + x)^2 = 1 + x^2
In F3, (1 + x)^3 = 1 + x^3
In F11, (1 + x)^11 = 1 + x^11
In F101, (1 + x)^101 = 1 + x^101
In F8291, (1 + x)^8291 = 1 + x^8291

P8291{1, 1}.pow(8291) の計算には５秒程度の時間がかかった。