競プロ日記

問題
ACコード

#include <bits/stdc++.h>
#define GET_REP(_1, _2, _3, NAME, ...) NAME
#define rep(...) GET_REP(__VA_ARGS__, irep, _rep)(__VA_ARGS__)
#define rep1(...) GET_REP(__VA_ARGS__, irep1, _rep1)(__VA_ARGS__)
#define _rep(i, n) irep (i, 0, n)
#define _rep1(i, n) irep1(i, 1, n)
#define irep(i, a, n) for (int i = a; i < (int)(n); ++i)
#define irep1(i, a, n) for (int i = a; i <= (int)(n); ++i)
#define rrep(i, n) for (int i = (int)(n) - 1; i >= 0; --i)
#define rrep1(i, n) for (int i = (int)(n); i >= 1; --i)
#define allrep(X, x) for (auto &&X : x)
#define all(x) (x).begin(), (x).end()
#ifdef LOCAL
  #include "../../Lib/cout_container.hpp"
  #define debug(x) cerr << #x " => " << (x) << endl
#else
  #define debug(x) 0
#endif
using lint = long long;
constexpr int    INF  = 1 << 29;
constexpr lint   INFL = 1LL << 61;
constexpr int    MOD  = (int)1e9 + 7;
constexpr double EPS  = 1e-9;
using namespace std;
namespace { struct INIT { INIT() { cin.tie(0); ios::sync_with_stdio(false); cout << fixed << setprecision(15); } } INIT; }

class VariableModInt {
  long long x;
  static long long &mod() {
    static long long mod_num = 1;
    return mod_num;
  }

  public:
  VariableModInt(const long long num = 0) : x(num % get_mod()) { if (x < 0) x += get_mod(); }
  bool operator==(const VariableModInt &r) const { return this->x == r.x; }
  bool operator!=(const VariableModInt &r) const { return this->x != r.x; }
  bool operator<(const VariableModInt &r) const { return this->x < r.x; }
  bool operator>(const VariableModInt &r) const { return this->x > r.x; }
  bool operator<=(const VariableModInt &r) const { return this->x <= r.x; }
  bool operator>=(const VariableModInt &r) const { return this->x >= r.x; }
  VariableModInt operator+(const VariableModInt &r) const { return VariableModInt(*this) += r; }
  VariableModInt operator-(const VariableModInt &r) const { return VariableModInt(*this) -= r; }
  VariableModInt operator*(const VariableModInt &r) const { return VariableModInt(*this) *= r; }
  VariableModInt operator/(const VariableModInt &r) const { return VariableModInt(*this) /= r; }
  friend VariableModInt operator+(const long long &l, const VariableModInt &r) { return VariableModInt(l) += r; }
  friend VariableModInt operator-(const long long &l, const VariableModInt &r) { return VariableModInt(l) -= r; }
  friend VariableModInt operator*(const long long &l, const VariableModInt &r) { return VariableModInt(l) *= r; }
  friend VariableModInt operator/(const long long &l, const VariableModInt &r) { return VariableModInt(l) /= r; }
  VariableModInt &operator++() { return *this += 1; }
  VariableModInt &operator--() { return *this -= 1; }
  VariableModInt &operator+=(const VariableModInt &r) {
    x += r.x;
    if (x >= get_mod()) x -= get_mod();
    return *this;
  }
  VariableModInt &operator-=(const VariableModInt &r) {
    if (x < r.x) x += get_mod();
    x -= r.x;
    return *this;
  }
  VariableModInt &operator*=(const VariableModInt &r) {
    x = x * r.x % get_mod();
    return *this;
  }
  VariableModInt &operator/=(VariableModInt r) {
    *this *= r.inverse();
    return *this;
  }
  VariableModInt inverse(void) const {
    long long a = x, b = get_mod(), u = 1, v = 0;
    while (b) {
      long long t = a / b;
      swap(a -= t * b, b);
      swap(u -= t * v, v);
    }
    return VariableModInt(u);
  }
  VariableModInt pow(long long n) const {
    VariableModInt ret(1), mul(x);
    while (n) {
      if (n & 1) ret *= mul;
      mul *= mul;
      n  >>= 1;
    }
    return ret;
  }
  friend istream &operator>>(istream &is, VariableModInt &r) {
    long long n;
    is >> n;
    r = VariableModInt(n);
    return is;
  }
  friend ostream &operator<<(ostream &os, const VariableModInt &r) { return os << r.x; }
  static void set_mod(const long long mod_num) { mod() = mod_num; }
  static long long get_mod() { return mod(); }
  explicit operator int() const { return x; }
  explicit operator long long() const { return x; }
  explicit operator double() const { return x; }
};

template <typename T>
vector<T> lagrange_interpolation(const vector<T> &yvec) {
  const int n = yvec.size() - 1;
  vector<T> ret(n + 1), fac(n + 1), c(n + 2);
  fac[0] = 1;
  for (int i = 1; i <= n; ++i) fac[i] = fac[i - 1] * i;
  c[0] = 0, c[1] = 1;
  for (int i = 1; i <= n; ++i) {
    for (int j = i + 1; j > 0; --j) {
      c[j] = c[j - 1] - c[j] * i;
    }
  }
  for (int i = 0; i <= n; ++i) {
    T c2 = 0, qi = yvec[i] / (fac[n - i] * fac[i]);
    if ((n - i) & 1) qi *= -1;
    for (int j = n; j >= 0; --j) {
      c2 = c[j + 1] + c2 * i;
      ret[j] += c2 * qi;
    }
  }
  return ret;
}

using mint = VariableModInt;
int main(void) {
  int p;
  cin >> p;
  mint::set_mod(p);
  vector<mint> a(p);
  rep (i, p) cin >> a[i];
  auto ans = lagrange_interpolation(a);
  rep (i, ans.size()) cout << ans[i] << (i != ans.size() - 1 ? " " : "\n");
  return 0;
}

解説が少ない！！！
PDF解法のフェルマーの小定理を使う解法はかなり頭の良い解法らしいし、そもそもフェルマーの小定理なんて知らないので、ラグランジュ補間で頑張ろうと思ってライブラリの作成から始めた。
ラグランジュ補間自体は理解できるのだが、係数を復元するというのでかなり悩まされた。
おそらく丁寧に書かれているであろう解説記事はあったものの、解説を読んでもコードを読んでもとにかく理解できない。
数学力の無さが原因なのかと少し寂しい気持ちになる。
総じて大卒程度の数学力はあると思うが、実際は大卒ではないし、受けた数学の講義も数学Iや応用数学等で、幾何学や代数学等とジャンル分けして詳しく学べる講義もなかった。
数学とは違うところでFFTなんかは出てきたのだが、数学的な演習はしていないのでどういうところで使えるか等は分からない。
競プロやっていると少し大学行っておけば良かったなという気持ちになる。

アルゴリズムは理解してACできたものの、解説記事に関しては未だに理解していない。
分からない部分は、結局のところ(x-xi)の総乗の係数の計算をやっているだけだと思うので、解説記事の方法は置いておいて自分流になんとか完成させた。
オーダーもO(N^2)に収まっているので問題ないだろう。
こういう他の解説は完全に理解できなかったが、自分なりに実装できた問題のためにブログを書いていると言っても過言ではない。
後から理解できなかったところが理解できて、過去の自分の解法と比較できるようになりたい。

コードで計算した順番に内容を書いておく。
求める多項式をN次多項式として、入力の点をN+1点とする。
・0!~N!まで階乗を計算する。
DP的に計算するだけでO(N)。
階乗の計算は高速化のために事前計算として行っているのだが、オーダー的には改善されない。
係数の復元ではなく、任意のaに対してf(a)を求める場合にオーダーレベルで高速化される。

・N+1次式{(x-xi)の総乗}の係数(C0~CN+1)を計算する。
解説記事はおそらくここを解説しているであろう部分が理解できなかった。
自分なりに愚直に計算した。
各桁ごとに計算し、(x-xi)をかけるごとに係数の数が1つ増えていくので、全体としてはO(N^2)。

・ラグランジュ補間多項式の各iについて直接計算できる部分Qiを計算する。
xの係数に関係なく、xiに代入することで直接計算できる部分をQiとおく。
具体的には、ラグランジュ補間多項式
\[L(x) = \sum_{i=0}^{N+1} y_i \frac{\prod_{j=0, j \neq i}^{N+1} (x - x_j)}{\prod_{j=0, j \neq i}^{N+1} (x_i - x_j)}\]のxがかからない部分。
\[Q(x) = \frac{y_i} {\prod_{j=0, j \neq i}^{N+1} (x_i - x_j)}\]
xiが0~Nまで連続しているので、(N-i)!*i!*-1^(N-i)で計算できる。
事前に階乗を計算しているのでO(1)。

・i=jを含まないN次式{(xi-xj)の総乗}の係数(C20~C2N)を各iについて計算する。
C2=C/(xi-xi)なので、c2はcを使って上位の桁から各桁ごとにO(1)で求められ、全桁でO(N)。
各iについて計算するのでO(N^2)。

・C2*Qiを各iについて全て足し合わせる。
ラグランジュ補間の公式に従って全て足し合わせる。
今回は係数について考えているので、係数ごとにQiをかけて足し合わせる。
QiがO(N)、C2がO(N^2)なので、全体としてO(N^2)

ラグランジュ補間の公式の見た目があまりにも複雑で、通常の書き方ではとても理解できそうになかったので、急遽数式の表示スクリプトを導入した。
しかしLATEXの書き方微塵も覚えていなくて笑った。
5,6年前？に学校で少し習ったのだが、それから自分含めクラスメイトが本当に誰一人として使っていないんじゃないかというくらい使わなかったLATEX形式。
論文は知る限りみんなWord。
当時は難しく感じたが、かなり簡単に書けてこんなところでも成長を感じてしまう。
最も当時は段落や章等も含んだ1つの文書として書けるようになる講義だったので、さすがに今回ほど簡単とはいかないだろうが、プログラムを触ったこともなかった時代と比べたらさすがに成長しているだろう。

ついでに今回、シンタックスハイライトで実体参照していないがためにtemplateあたりが崩れているのに気付いた。
変換しなくても正常に見えるじゃーんと思っていたが実はダメだった・・・。
わりと問題ない部分ではあるが、一応今回から実体参照にして修正することにした。
どんどんブログ書くのが面倒臭くなっていく。