All Posts

全ての記事をそのまま表示している。

とても重くなってしまっているかもしれない。

org-mode の概要

    about

org-mode とは emacs 付属の軽量マークアップ言語で、 markdown や ReStructuredText の親戚にあたる。

このページには org-mode の基本的な機能をまとめた。

リンク

これは [[https://www.google.co.jp/][リンク]] です

と書くと

これは リンク です。

となる。

リンクを挿入するコマンドとして C-c C-l (org-insert-link) がある。

ウェブのリンクだけではなく画像のリンクなども同じように貼れる

表

| x |  y |
|--------|
| 1 |  1 |
| 2 |  4 |
| 3 |  9 |

と書くと

x y
1 1
2 4
3 9

となる。

emacs の org-mode エディタには表を編集する機能があり、 “|” を書いてから Tab でいい感じに表をアラインできる。

文字の装飾

通常のテキスト \\
+打ち消し線+ \\
_下線_ \\
*太字* \\
/斜体/ \\
~source code~ \\

と書くと

通常のテキスト
打ち消し線
下線
太字
斜体
source code

となる。

LaTeX

数式 \(e = mc^{2}\) を含んだ文章

とかくと

数式 \(e = mc^{2}\) を含んだ文章

となる。

中央寄せするには

\[ y = f(x) \]

とすれば

\[ y = f(x) \]

となる。

この他に align, array, matrix などの LaTeX の環境も使える。 この場合は \[ ... \] で環境を囲む必要はなく、 テキスト中に begin{env} ... end{env} といきなり書いても LaTeX のコードだと認識される。

文章中の LaTeX コードは C-c C-x C-l (org-toggle-latex-fragment) で画像に変換できる。

Figure 1: org-toggle-latex-fragment を実行した後のエディタ

Figure 1: org-toggle-latex-fragment を実行した後のエディタ

箇条書き

- 項目1

  項目1の説明

- 項目2

  項目3の説明

- 項目3

  項目3の説明

と書くと

  • 項目1

    項目1の説明

  • 項目2

    項目3の説明

  • 項目3

    項目3の説明

となる。

ソースコード

#+BEGIN_SRC python
def fact(n):
    if n == 0:
        return 1
    else:
        return n * fact(n - 1)

# 6! を出力する
print(fact(6))
#+END_SRC

と書くと

def fact(n):
    if n == 0:
        return 1
    else:
        return n * fact(n - 1)

# 6! を出力する
print(fact(6))

となる。

org-src-fontify-natively を t にすればエディタ上でもシンタックスハイライトが行われる。 また C-c ' (org-edit-special) で言語にあったメジャーモードを使ってソースコードを編集できる。



このサイトの全体像

    about

このサイトは

  • org mode で書いた文章を
  • ox-hugo で markdown に変換し
  • markdown から Hugo で静的サイトを作成し
  • そのリポジトリを github pages で公開する

という手順で作られている。

全体像を画像にまとめた。自分が直接編集するファイルを緑で塗った。



BigIntの実装 in C++

    cpp  / コード置き場

多倍長整数を実装できたのでまとめた。C++のいい練習になると思って始めたらかなり時間がかかった。

全体のコードはかなり長くなったので最後に載せた。

使用例

作ったクラスで何かできるかを最初に示す。 BigInt fact(BigInt) と BigInt fib(BigInt) を int と同じ方法で定義し、 operator<< で出力している。

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

BigInt fact(BigInt n) {
  if (n == 1)
    return 1;
  else
    return n * fact(n - 1);
}

BigInt fib(BigInt n) {
    if (n < 2)
        return 1;
    BigInt a = 1;
    BigInt b = 1;
    for (; n > 1; --n) {
        BigInt temp = b;
        b = a + b;
        a = temp;
    }
    return b;
}

int main() {
  std::cout << "fact(100) = " << fact(100) << "\n";
  std::cout << "fib(1000) = " << fib(1000) << "\n";
  return 0;
}

fact(100) = 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000
fib(1000) = 70330367711422815821835254877183549770181269836358732742604905087154537118196933579742249494562611733487750449241765991088186363265450223647106012053374121273867339111198139373125598767690091902245245323403501

実装したこと

使用例を動かすための最低限の関数を実装した。最終的なヘッダは

class BigInt {
 private:
  ...(略)...

 public:
  static std::mt19937_64 random_engine;
  BigInt();
  BigInt(int32_t i);
  BigInt(int64_t i);
  BigInt(uint32_t i, bool sign = false);
  BigInt(uint64_t i, bool sign = false);
  BigInt(const std::string &str);

  static BigInt of_big_uint(BigUInt abs, bool sign = false);
  static BigInt of_vector(std::vector<uint32_t> abs, bool sign = false);
  static BigInt random(unsigned n = 5);
  uint64_t width() const;
  bool sign() const;
  bool is_zero() const;
  bool is_positive() const;
  bool is_negative() const;
  BigInt cshift_left(uint64_t shamt) const;
  BigInt &shift_left(uint64_t shamt);
  BigInt cshift_right(uint64_t shamt) const;
  BigInt &shift_right(uint64_t shamt);
  BigInt pow(BigInt exp) const;
  const uint32_t &operator[](uint64_t i) const;
  uint32_t &operator[](uint64_t i);
  BigInt &operator++();
  BigInt operator++(int);
  BigInt &operator--();
  BigInt operator--(int);
  BigInt operator-() const;
  friend bool operator==(const BigInt &lhs, const BigInt &rhs);
  friend bool operator!=(const BigInt &lhs, const BigInt &rhs);
  friend bool operator<(const BigInt &lhs, const BigInt &rhs);
  friend bool operator>(const BigInt &lhs, const BigInt &rhs);
  friend bool operator<=(const BigInt &lhs, const BigInt &rhs);
  friend bool operator>=(const BigInt &lhs, const BigInt &rhs);
  friend BigInt &operator+=(BigInt &lhs, const BigInt &rhs);
  friend BigInt &operator-=(BigInt &lhs, const BigInt &rhs);
  friend BigInt &operator*=(BigInt &lhs, const BigInt &rhs);
  friend BigInt &operator/=(BigInt &lhs, const BigInt &rhs);
  friend BigInt &operator%=(BigInt &lhs, const BigInt &rhs);
  friend BigInt operator+(const BigInt &lhs, const BigInt &rhs);
  friend BigInt operator-(const BigInt &lhs, const BigInt &rhs);
  friend BigInt operator*(const BigInt &lhs, const BigInt &rhs);
  friend BigInt operator/(const BigInt &lhs, const BigInt &rhs);
  friend BigInt operator%(const BigInt &lhs, const BigInt &rhs);
  friend std::ostream &operator<<(std::ostream &lhs, const BigInt &rhs);
};

となった。

実装する上でつまずいたところを書く。

符号と絶対値を分ける

符号のない多倍長整数を表すクラス BigUInt を作り、符号ありの多倍長整数を絶対値(BigUInt)と符号(bool)で表すことにした。

#include "big_uint.hpp"

class BigInt {
 private:
  bool m_sign;
  BigUInt m_abs;

 public:
  ...(略)...
}

符号と絶対値を分ければ多くの演算は符号による場合分けによって符号なしの演算に帰着できる。例えば加算は

BigInt &operator+=(BigInt &lhs, const BigInt &rhs) {
  // 以下の <, >, +=, -=, + が全てBigUIntのものであることに注意
  if (lhs.m_sign == rhs.m_sign) {
    lhs.m_abs += rhs.m_abs;
  } else if (lhs.m_abs > rhs.m_abs) {
    lhs.m_abs -= rhs.m_abs;
  } else if (lhs.m_abs < rhs.m_abs) {
    lhs.m_sign = !lhs.m_sign;
    lhs.m_abs = rhs.m_abs - lhs.m_abs;
  } else {
    lhs = BigInt(0ul);
  }
  return lhs;
}

となる。

固定長整数の場合には符号付き整数を2の補数で表現することで加算と減算の処理を同じにできる。 しかし、多倍長整数の場合には2項演算子の2つのオペランドのビット長が通常異なるので、減算を加算と同じように行うことができない。そのため2の補数表現は使わなかった。

BigUInt の基数を \(2^{32}\) とする

C++では \(n\) ビット符号なし整数の演算は \(\mod 2^{n}\) 行われるので、 \(2^{n}\)を基数にすると桁同士の演算で都合がいい。

加えて多倍長整数の乗算では桁同士の乗算を桁あふれなく扱う必要があるので、 基数を \(2^{32}\) (型は uint32_t) とした。 uint64_t を使えば uint32_t 同士の乗算を桁あふれなく扱うことができる。

class BigUInt {
 private:
  std::vector<uint32_t> m_vec;
  ...(略)...
};

オーバーフローの検出

多倍長整数の加算と減算ではキャリーを計算するためにオーバーフローを検出する必要がある。

uint32_t の変数 a, b を足したときにオーバーフローが起こるのは

\begin{align*} a + b &\geq 2^{32} \\
\iff \quad \quad a &> (2^{32} - 1) - b \end{align*}

のときである。最後の式の右辺の減算ではオーバーフローは発生しないので、 uint32_t の演算だけを使ってオーバーフローの有無を調べることができる。

実際のコードでは「この2つの数を足したらオーバーフローするか?」ではなく「この2つの数を足した結果オーバーフローが起こったか?」が必要になったので、次の事実を使った。

uint32_t 型の変数 a, b について、 a += b を実行したとする。この加算でオーバーフローが発生したかどうかは

\begin{align*} a < b \end{align*}

を調べれば分かる。

a += b を実行する前の a, b の値を \(A, B\) とする。

オーバーフローが発生した場合には

\begin{align*} a &= A + B \mod 2^{32} \\
&= A + B - 2^{32} \end{align*}

であり、

\begin{align*} b - a = 2^{32} - a > 0 \end{align*}

が成り立つ。 オーバーフローが発生しなかった場合には

\begin{align*} a &= A + B \mod 2^{32} \\
&= A + B \end{align*}

より

\begin{align*} a - b = A \geq 0 \end{align*}

が成り立つ。

コードでは以下のようになった。

BigUInt &operator+=(BigUInt &lhs, const BigUInt &rhs) {
  uint32_t carry = 0;
  uint32_t i;

  lhs.m_vec.resize(std::max(lhs.width(), rhs.width()));
  for (i = 0; i < rhs.width(); ++i) {
    // lhs[i] += carry + lhs[i] を実行し、carryを更新する
    lhs[i] += carry;
    carry = lhs[i] < carry;
    lhs[i] += rhs[i];
    carry |= lhs[i] < rhs[i];
  }

  for (; carry && i < lhs.width(); ++i) {
    // lhs[i] += carry を実行し、carryを更新する
    ++lhs[i];
    carry = lhs[i] < carry;
  }

  if (carry) {
    lhs.m_vec.push_back(1);
  }

  return lhs;
}

減算も似たような関係を使って書いた。

除算アルゴリズム

除算には TAOCP 2巻 の4.3.1項にあるアルゴリズムDを用いた。 このアルゴリズム自体は何も間違っていないが、本に載っている前提条件が普通でないので少しハマった。

具体的には、TAOCPでは \(m + n - 1\) 桁の整数を \(n - 1\) 桁の整数で割った商を \(m\) 桁としているが、 一般には \(m + 1\) 桁になりうる。そして商の最初の桁についてはインデックスの特別扱いが必要になる。 前提条件を読み間違えたせいでメモリリークが発生して大変だった。

効率

効率をあまり意識して書いていないのと、乗算アルゴリズムが \(O(mn)\) なことがあって効率は全然良くない。 冒頭で挙げた使用例は実行に1秒程度かかるが、 pythonインタープリタで同じコードを書いて実行すると一瞬で終わる(悲しい)。

さらなる高速化の手段としては

  • 乗算アルゴリズムの変更 (Schönhage–Strassen algorithm, Karatsuba algorithm)
  • 配列の要素として uint32_t ではなく uint64_t を使う
  • アセンブリを書く(コンディションレジスタを使う / mulhi命令を(あれば)使う)
  • operator<< でバッファを使うようにする、あるいは operator<< を廃止する。

が考えられる。必要になったらやる。

BigUIntクラス

BigUIntクラスのヘッダと実装を載せる。

Code Snippet 1: big_uint.hpp 

#ifndef BIG_UINT_H
#define BIG_UINT_H

#include <cstdint>
#include <ostream>
#include <random>
#include <string>
#include <vector>

class BigUInt {
 private:
  static std::mt19937_64 random_engine;
  std::vector<uint32_t> m_vec;
  BigUInt(std::vector<uint32_t> vec);
  BigUInt &delete_trailing_zeros();

 public:
  BigUInt();
  BigUInt(int32_t i);
  BigUInt(int64_t i);
  BigUInt(uint32_t i);
  BigUInt(uint64_t i);
  BigUInt(const std::string &str);
  static BigUInt of_vector(std::vector<uint32_t> vec);
  static BigUInt random(unsigned n = 5);
  uint64_t width() const;
  bool is_zero() const;
  BigUInt cshift_left(uint64_t shamt) const;
  BigUInt &shift_left(uint64_t shamt);
  BigUInt cshift_right(uint64_t shamt) const;
  BigUInt &shift_right(uint64_t shamt);
  BigUInt pow(BigUInt exp) const;
  const uint32_t &operator[](uint64_t i) const;
  uint32_t &operator[](uint64_t i);
  BigUInt &operator++();
  BigUInt operator++(int);
  BigUInt &operator--();
  BigUInt operator--(int);
  BigUInt operator-();
  friend bool operator==(const BigUInt &lhs, const BigUInt &rhs);
  friend bool operator!=(const BigUInt &lhs, const BigUInt &rhs);
  friend bool operator<(const BigUInt &lhs, const BigUInt &rhs);
  friend bool operator>(const BigUInt &lhs, const BigUInt &rhs);
  friend bool operator<=(const BigUInt &lhs, const BigUInt &rhs);
  friend bool operator>=(const BigUInt &lhs, const BigUInt &rhs);
  friend BigUInt &operator+=(BigUInt &lhs, const BigUInt &rhs);
  friend BigUInt &operator-=(BigUInt &lhs, const BigUInt &rhs);
  friend BigUInt &operator*=(BigUInt &lhs, const BigUInt &rhs);
  friend BigUInt &operator/=(BigUInt &lhs, const BigUInt &rhs);
  friend BigUInt &operator%=(BigUInt &lhs, const BigUInt &rhs);
  friend BigUInt operator+(const BigUInt &lhs, const BigUInt &rhs);
  friend BigUInt operator-(const BigUInt &lhs, const BigUInt &rhs);
  friend BigUInt operator*(const BigUInt &lhs, const BigUInt &rhs);
  friend BigUInt operator/(const BigUInt &lhs, const BigUInt &rhs);
  friend BigUInt operator%(const BigUInt &lhs, const BigUInt &rhs);
  friend std::ostream &operator<<(std::ostream &lhs, const BigUInt &rhs);
};

#endif /* BIG_UINT_H */

Code Snippet 2: big_uint.cpp 

#include "big_uint.hpp"
#include <stdexcept>

template <typename T>
constexpr bool bit_at(T a, unsigned n) {
  static_assert(std::is_arithmetic_v<T>);
  return a & (static_cast<T>(1) << n);
}

// position of the highest non-zero bit.
// returns width of the argument type if a == 0.
template <typename T>
constexpr unsigned hnz(T a) {
  static_assert(std::is_unsigned_v<T>);
  unsigned i = std::numeric_limits<T>::digits - 1;
  for (; i != std::numeric_limits<T>::max(); --i)
    if (bit_at(a, i)) {
      return i;
    }
  return i;
}

std::mt19937_64 BigUInt::random_engine((std::random_device())());

BigUInt::BigUInt(std::vector<uint32_t> vec) : m_vec(vec) {
  delete_trailing_zeros();
}

BigUInt &BigUInt::delete_trailing_zeros() {
  while (width() > 1 && m_vec[width() - 1] == 0) {
    m_vec.pop_back();
  }
  return *this;
}

BigUInt::BigUInt() {}

BigUInt::BigUInt(int32_t i) : m_vec({static_cast<uint32_t>(i)}) {}

BigUInt::BigUInt(int64_t i) : BigUInt(static_cast<uint64_t>(i)) {}

BigUInt::BigUInt(uint32_t i) : m_vec({i}) {}

BigUInt::BigUInt(uint64_t i)
    : m_vec({static_cast<uint32_t>(i), static_cast<uint32_t>(i >> 32)}) {
  delete_trailing_zeros();
}

BigUInt::BigUInt(const std::string &str) : m_vec({0}) {
  BigUInt base(1);
  BigUInt ten(10);

  for (uint64_t i = str.length() - 1; i != UINT64_MAX; --i) {
    *this += std::stoul(str.substr(i, 1)) * base;
    base *= ten;
  }

  delete_trailing_zeros();
}

uint64_t BigUInt::width() const { return m_vec.size(); }

const uint32_t &BigUInt::operator[](uint64_t i) const { return m_vec[i]; }

uint32_t &BigUInt::operator[](uint64_t i) { return m_vec[i]; }

BigUInt BigUInt::of_vector(std::vector<uint32_t> vec) { return BigUInt(vec); }

BigUInt BigUInt::random(unsigned n) {
  std::vector<uint32_t> vec(n);
  for (auto &i : vec) {
    i = static_cast<uint32_t>(random_engine());
  }
  return BigUInt(vec);
}

bool BigUInt::is_zero() const {
  for (uint32_t i = 0; i < width(); ++i) {
    if ((*this)[i] != 0) {
      return false;
    }
  }
  return true;
}

BigUInt BigUInt::cshift_left(uint64_t shamt) const {
  BigUInt temp(*this);
  return temp.shift_left(shamt);
}

BigUInt &BigUInt::shift_left(uint64_t shamt) {
  if (shamt == 0) {
    return *this;
  }

  const uint64_t shamt_quot = shamt / 32;
  const uint64_t shamt_rem = shamt % 32;
  m_vec.resize(width() + shamt_quot + 1);

  uint64_t i = width() - 1;
  if (shamt_rem == 0) {
    // avoid undefined behavior by shifting 0 or 32.
    for (; i != shamt_quot - 1; --i) {
      (*this)[i] = (*this)[i - shamt_quot];
    }
  } else {
    const uint64_t shamt_rem_inv = 32 - shamt_rem;
    for (; i != shamt_quot; --i) {
      auto lower = (*this)[i - shamt_quot - 1] >> shamt_rem_inv;
      auto upper = (*this)[i - shamt_quot] << shamt_rem;
      (*this)[i] = lower + upper;
    }
    (*this)[i--] = (*this)[0] << shamt_rem;
  }
  for (; i != UINT64_MAX; --i) {
    (*this)[i] = 0;
  }
  return delete_trailing_zeros();
}

BigUInt BigUInt::cshift_right(uint64_t shamt) const {
  BigUInt temp(*this);
  return temp.shift_right(shamt);
}

BigUInt &BigUInt::shift_right(uint64_t shamt) {
  if (shamt == 0) {
    return *this;
  }

  const uint64_t shamt_quot = shamt / 32;
  const uint64_t shamt_rem = shamt % 32;

  if (shamt_rem == 0) {
    // avoid undefined behavior by shifting 0 or 32.
    for (uint64_t i = 0; i < width() - shamt_quot; ++i) {
      (*this)[i] = (*this)[i + shamt_quot];
    }
  } else {
    const uint64_t shamt_rem_inv = 32 - shamt_rem;
    for (uint64_t i = 0; i < width() - shamt_quot - 1; ++i) {
      auto lower = (*this)[i + shamt_quot] >> shamt_rem;
      auto upper = (*this)[i + shamt_quot + 1] << shamt_rem_inv;
      (*this)[i] = lower + upper;
    }
    (*this)[width() - shamt_quot - 1] = (*this)[width() - 1] >> shamt_rem;
  }
  m_vec.resize(width() - shamt_quot);
  return delete_trailing_zeros();
}

BigUInt BigUInt::operator-() {
  throw std::invalid_argument("unary - of BigUInt");
}

BigUInt &BigUInt::operator++() {
  *this += 1;
  return *this;
}

BigUInt BigUInt::operator++(int) {
  auto temp(*this);
  *this += 1;
  return temp;
}

BigUInt &BigUInt::operator--() {
  *this -= 1;
  return *this;
}

BigUInt BigUInt::operator--(int) {
  auto temp(*this);
  *this -= 1;
  return temp;
}

bool operator==(const BigUInt &lhs, const BigUInt &rhs) {
  return lhs.m_vec == rhs.m_vec;
}

bool operator!=(const BigUInt &lhs, const BigUInt &rhs) {
  return lhs.m_vec != rhs.m_vec;
}

bool operator<(const BigUInt &lhs, const BigUInt &rhs) {
  if (lhs.width() < rhs.width()) {
    return true;
  }

  if (lhs.width() > rhs.width()) {
    return false;
  }

  for (uint64_t i = lhs.width() - 1; i != UINT64_MAX; --i) {
    if (lhs[i] == rhs[i]) {
      continue;
    }
    return lhs[i] < rhs[i];
  }
  return false;
}

bool operator<=(const BigUInt &lhs, const BigUInt &rhs) {
  return lhs == rhs || lhs < rhs;
}

bool operator>(const BigUInt &lhs, const BigUInt &rhs) { return rhs < lhs; }

bool operator>=(const BigUInt &lhs, const BigUInt &rhs) { return rhs <= lhs; }

BigUInt &operator+=(BigUInt &lhs, const BigUInt &rhs) {
  uint32_t carry = 0;
  uint32_t i;

  lhs.m_vec.resize(std::max(lhs.width(), rhs.width()));
  for (i = 0; i < rhs.width(); ++i) {
    lhs[i] += carry;
    carry = lhs[i] < carry;
    lhs[i] += rhs[i];
    carry |= lhs[i] < rhs[i];
  }

  for (; carry && i < lhs.width(); ++i) {
    ++lhs[i];
    carry = lhs[i] < carry;
  }

  if (carry) {
    lhs.m_vec.push_back(1);
  }

  return lhs;
}

BigUInt operator+(const BigUInt &lhs, const BigUInt &rhs) {
  BigUInt ans(lhs);
  return ans += rhs;
}

// if lhs < rhs, result is undefined.
BigUInt &operator-=(BigUInt &lhs, const BigUInt &rhs) {
  uint32_t i;
  uint32_t carry = 0;
  for (i = 0; i < rhs.width(); i++) {
    lhs[i] -= carry;
    carry = lhs[i] == UINT32_MAX && carry == 1;
    carry |= lhs[i] < rhs[i];
    lhs[i] -= rhs[i];
  }

  for (; carry && i < lhs.width(); ++i) {
    --lhs[i];
    carry = lhs[i] == UINT32_MAX && carry == 1;
  }

  return lhs.delete_trailing_zeros();
}

BigUInt operator-(const BigUInt &lhs, const BigUInt &rhs) {
  BigUInt ans(lhs);
  return ans -= rhs;
}

BigUInt operator*(const BigUInt &lhs, const BigUInt &rhs) {
  BigUInt ans;
  std::vector<uint32_t> lower(lhs.width());
  std::vector<uint32_t> upper(lhs.width() + 1);

  for (uint32_t i = 0; i < rhs.width(); ++i) {
    for (uint32_t j = 0; j < lhs.width(); ++j) {
      upper[0] = 0;
      lower[j] = lhs[j] * rhs[i];
      upper[j + 1] = (static_cast<uint64_t>(lhs[j]) * rhs[i]) >> 32;
    }
    ans += (BigUInt(lower) + BigUInt(upper)).shift_left(32 * i);
  }
  return ans;
}

BigUInt &operator*=(BigUInt &lhs, const BigUInt &rhs) {
  lhs = lhs * rhs;
  return lhs;
}

// knuth's division algorithm D.
// returns quotient, modifing dividend to be remainder.
static BigUInt divmod(BigUInt &dividend, const BigUInt &divisor) {

  if (divisor.is_zero()) {
    throw std::invalid_argument("BigUInt : division by zero.");
  }

  if (dividend < divisor) {
    return 0;
  }

  // normalize dividend and divisor.
  const uint64_t shamt = 31 - hnz(divisor[divisor.width() - 1]);
  dividend.shift_left(shamt);
  const_cast<BigUInt &>(divisor).shift_left(shamt);

  uint64_t div_msb = divisor[divisor.width() - 1];
  uint64_t div_ssb = divisor.width() == 1 ? 0 : divisor[divisor.width() - 2];

  std::vector<uint32_t> ans;
  ans.resize(dividend.width() - divisor.width() + 1);

  for (uint64_t i = ans.size() - 1; i != UINT64_MAX; --i) {
    uint64_t upper = (i + divisor.width() < dividend.width())
                         ? dividend[i + divisor.width()]
                         : 0;
    uint64_t lower = dividend[i + divisor.width() - 1];
    uint64_t num = (upper << 32) + lower;
    uint64_t q = num / div_msb;
    uint64_t r = num % div_msb;
    uint64_t idx = i + divisor.width() - 2;
    while (r < UINT32_MAX &&
           (q >= UINT32_MAX ||
            q * div_ssb >
                (r << 32) + (idx < dividend.width() ? dividend[idx] : 0))) {
      q--;
      r += div_msb;
    }

    BigUInt sub = (q * divisor).cshift_left(i * 32);

    if (dividend < sub) {
      q--;
      sub -= divisor.cshift_left(i * 32);
    }

    ans[i] = static_cast<uint32_t>(q);
    dividend -= sub;
  }

  const_cast<BigUInt &>(divisor).shift_right(shamt);
  dividend.shift_right(shamt);
  return BigUInt::of_vector(ans);
}

BigUInt operator/(const BigUInt &lhs, const BigUInt &rhs) {
  auto temp(lhs);
  return divmod(temp, rhs);
}

BigUInt &operator/=(BigUInt &lhs, const BigUInt &rhs) {
  return lhs = lhs / rhs;
}

BigUInt operator%(const BigUInt &lhs, const BigUInt &rhs) {
  BigUInt temp(lhs);
  return temp %= rhs;
}

BigUInt &operator%=(BigUInt &lhs, const BigUInt &rhs) {
  divmod(lhs, rhs);
  return lhs;
}

BigUInt BigUInt::pow(BigUInt exp) const {
  BigUInt x = 1;
  for (; exp != 0; --exp) {
    x *= *this;
  }
  return x;
}

std::ostream &operator<<(std::ostream &lhs, const BigUInt &rhs) {
  BigUInt base(10);
  const BigUInt ten(10);

  while (rhs >= base) {
    base *= ten;
  }

  auto temp(rhs);
  while (base != 1) {
    temp %= base;
    lhs << (temp / (base /= ten))[0];
  }
  return lhs;
}

BigIntクラス

BigIntクラスのヘッダと実装を載せる。

Code Snippet 3: big_int.hpp 

#ifndef BITINT_H
#define BITINT_H

#include "big_uint.hpp"

class BigInt {
 private:
  bool m_sign;
  BigUInt m_abs;

 public:
  static std::mt19937_64 random_engine;
  BigInt();
  BigInt(int32_t i);
  BigInt(int64_t i);
  BigInt(uint32_t i, bool sign = false);
  BigInt(uint64_t i, bool sign = false);
  BigInt(const std::string &str);

  static BigInt of_big_uint(BigUInt abs, bool sign = false);
  static BigInt of_vector(std::vector<uint32_t> abs, bool sign = false);
  static BigInt random(unsigned n = 5);
  uint64_t width() const;
  bool sign() const;
  bool is_zero() const;
  bool is_positive() const;
  bool is_negative() const;
  BigInt cshift_left(uint64_t shamt) const;
  BigInt &shift_left(uint64_t shamt);
  BigInt cshift_right(uint64_t shamt) const;
  BigInt &shift_right(uint64_t shamt);
  BigInt pow(BigInt exp) const;
  const uint32_t &operator[](uint64_t i) const;
  uint32_t &operator[](uint64_t i);
  BigInt &operator++();
  BigInt operator++(int);
  BigInt &operator--();
  BigInt operator--(int);
  BigInt operator-() const;
  friend bool operator==(const BigInt &lhs, const BigInt &rhs);
  friend bool operator!=(const BigInt &lhs, const BigInt &rhs);
  friend bool operator<(const BigInt &lhs, const BigInt &rhs);
  friend bool operator>(const BigInt &lhs, const BigInt &rhs);
  friend bool operator<=(const BigInt &lhs, const BigInt &rhs);
  friend bool operator>=(const BigInt &lhs, const BigInt &rhs);
  friend BigInt &operator+=(BigInt &lhs, const BigInt &rhs);
  friend BigInt &operator-=(BigInt &lhs, const BigInt &rhs);
  friend BigInt &operator*=(BigInt &lhs, const BigInt &rhs);
  friend BigInt &operator/=(BigInt &lhs, const BigInt &rhs);
  friend BigInt &operator%=(BigInt &lhs, const BigInt &rhs);
  friend BigInt operator+(const BigInt &lhs, const BigInt &rhs);
  friend BigInt operator-(const BigInt &lhs, const BigInt &rhs);
  friend BigInt operator*(const BigInt &lhs, const BigInt &rhs);
  friend BigInt operator/(const BigInt &lhs, const BigInt &rhs);
  friend BigInt operator%(const BigInt &lhs, const BigInt &rhs);
  friend std::ostream &operator<<(std::ostream &lhs, const BigInt &rhs);
};
#endif // BITINT_H

Code Snippet 4: big_int.cpp 

#include "big_int.hpp"

const uint32_t &BigInt::operator[](uint64_t i) const { return m_abs[i]; }

uint32_t &BigInt::operator[](uint64_t i) { return m_abs[i]; }

std::mt19937_64 BigInt::random_engine((std::random_device())());

BigInt::BigInt() {}

BigInt::BigInt(int32_t i)
    : m_sign(i < 0),
      m_abs(m_sign ? static_cast<uint32_t>(-i) : static_cast<uint32_t>(i)) {}

BigInt::BigInt(int64_t i)
    : m_sign(i < 0),
      m_abs(m_sign ? static_cast<uint64_t>(-i) : static_cast<uint64_t>(i)) {}

BigInt::BigInt(uint32_t i, bool sign) : m_sign(sign), m_abs(i) {}

BigInt::BigInt(uint64_t i, bool sign) : m_sign(sign), m_abs(i) {}

BigInt::BigInt(const std::string &str)
    : m_sign(str[0] == '-'), m_abs(str.substr(m_sign, str.length() - m_sign)) {}

BigInt BigInt::of_big_uint(BigUInt abs, bool sign) {
  BigInt temp;
  temp.m_abs = abs;
  temp.m_sign = sign;
  return temp;
}

BigInt BigInt::of_vector(std::vector<uint32_t> abs, bool sign) {
  BigInt temp;
  temp.m_abs = BigUInt::of_vector(abs);
  temp.m_sign = sign;
  return temp;
}

BigInt BigInt::random(unsigned n) {
  return BigInt::of_big_uint(BigUInt::random(n), random_engine() % 2);
}

uint64_t BigInt::width() const { return m_abs.width(); }

bool BigInt::sign() const { return m_sign; }
bool BigInt::is_zero() const { return m_abs.is_zero(); }
bool BigInt::is_positive() const { return !m_sign; }
bool BigInt::is_negative() const { return m_sign; }

BigInt BigInt::cshift_left(uint64_t shamt) const {
  return BigInt::of_big_uint(m_abs.cshift_left(shamt), m_sign);
}

BigInt &BigInt::shift_left(uint64_t shamt) {
  m_abs.shift_left(shamt);
  return *this;
}

BigInt BigInt::cshift_right(uint64_t shamt) const {
  return BigInt::of_big_uint(m_abs.cshift_right(shamt), m_sign);
}

BigInt &BigInt::shift_right(uint64_t shamt) {
  m_abs.shift_right(shamt);
  return *this;
}

bool operator==(const BigInt &lhs, const BigInt &rhs) {
  return lhs.m_sign == rhs.m_sign && lhs.m_abs == rhs.m_abs;
}

bool operator!=(const BigInt &lhs, const BigInt &rhs) {
  return lhs.m_sign != rhs.m_sign || lhs.m_abs != rhs.m_abs;
}

BigInt BigInt::operator-() const { return BigInt::of_big_uint(m_abs, !m_sign); }

bool operator<(const BigInt &lhs, const BigInt &rhs) {
  if (lhs.is_positive() && rhs.is_negative()) {
    return false;
  }
  if (lhs.is_negative() && rhs.is_positive()) {
    return true;
  }
  if (lhs.is_positive() && rhs.is_positive()) {
    return lhs.m_abs < rhs.m_abs;
  }
  return lhs.m_abs > rhs.m_abs;
}

bool operator<=(const BigInt &lhs, const BigInt &rhs) {
  return lhs == rhs || lhs < rhs;
}

bool operator>(const BigInt &lhs, const BigInt &rhs) { return rhs < lhs; }

bool operator>=(const BigInt &lhs, const BigInt &rhs) { return rhs <= lhs; }

BigInt &operator+=(BigInt &lhs, const BigInt &rhs) {
  if (lhs.m_sign == rhs.m_sign) {
    lhs.m_abs += rhs.m_abs;
  } else if (lhs.m_abs > rhs.m_abs) {
    lhs.m_abs -= rhs.m_abs;
  } else if (lhs.m_abs < rhs.m_abs) {
    lhs.m_sign = !lhs.m_sign;
    lhs.m_abs = rhs.m_abs - lhs.m_abs;
  } else {
    lhs = BigInt(0ul);
  }
  return lhs;
}

BigInt operator+(const BigInt &lhs, const BigInt &rhs) {
  BigInt ans = lhs;
  return ans += rhs;
}

BigInt &operator-=(BigInt &lhs, const BigInt &rhs) {
  if (lhs.m_sign != rhs.m_sign) {
    lhs.m_abs += rhs.m_abs;
  } else if (lhs.m_abs > rhs.m_abs) {
    lhs.m_abs -= rhs.m_abs;
  } else if (lhs.m_abs < rhs.m_abs) {
    lhs.m_sign = !lhs.m_sign;
    lhs.m_abs = rhs.m_abs - lhs.m_abs;
  } else {
    lhs = BigInt(0ul);
  }
  return lhs;
}

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

BigInt &operator*=(BigInt &lhs, const BigInt &rhs) {
  lhs.m_abs *= rhs.m_abs;
  lhs.m_sign = lhs.sign() != rhs.sign();
  return lhs;
}

BigInt operator*(const BigInt &lhs, const BigInt &rhs) {
  auto temp(lhs);
  return temp *= rhs;
}

BigInt &operator/=(BigInt &lhs, const BigInt &rhs) {
  lhs.m_abs /= rhs.m_abs;
  lhs.m_sign = lhs.sign() != rhs.sign();
  return lhs;
}

BigInt operator/(const BigInt &lhs, const BigInt &rhs) {
  auto temp(lhs);
  return temp /= rhs;
}

BigInt &operator%=(BigInt &lhs, const BigInt &rhs) {
  lhs.m_abs %= rhs.m_abs;
  lhs.m_sign = lhs.sign() != rhs.sign();
  return lhs;
}

BigInt operator%(const BigInt &lhs, const BigInt &rhs) {
  auto temp(lhs);
  return temp %= rhs;
}

BigInt &BigInt::operator++() {
  *this += 1;
  return *this;
}

BigInt BigInt::operator++(int) {
  auto temp(*this);
  *this += 1;
  return temp;
}

BigInt &BigInt::operator--() {
  *this -= 1;
  return *this;
}

BigInt BigInt::operator--(int) {
  auto temp(*this);
  *this -= 1;
  return temp;
}

BigInt BigInt::pow(BigInt exp) const {
  BigInt x = 1;
  for (; exp > 0; --exp) {
    x *= *this;
  }
  return x;
}

std::ostream &operator<<(std::ostream &lhs, const BigInt &rhs) {
  if (rhs.is_negative()) {
    lhs << "-";
  }
  lhs << rhs.m_abs;
  return lhs;
}


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

    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) の計算には５秒程度の時間がかかった。



有限体の実装 in C++

    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