seq/p_recursive.hpp

View this file on GitHub
Last update: 2024-12-17 23:15:20+09:00
Include: #include "seq/p_recursive.hpp"

Depends on

Code

#include "linalg/solve_linear.hpp"
#include "mod/all_inverse.hpp"

/*
return: polynomial sequence F[0](x),...,F[m](x)
F[0](n)A[n+m]+...+F[m](n)A[n]=0
何かが失敗したら return empty
*/
template <typename mint>
vvc<mint> find_p_recursive(vc<mint> A, int max_m, int max_deg) {
  int N = (1 + max_m) * (1 + max_deg);
  auto idx = [&](int i, int d) -> int { return (1 + max_deg) * i + d; };

  assert(len(A) >= N + max_m);
  vv(mint, mat, N, N);
  FOR(row, N) {
    FOR(i, max_m + 1) {
      mint pw = 1;
      FOR(j, max_deg + 1) { mat[row][idx(i, j)] = A[row + max_m - i] * pw, pw *= row; }
    }
  }
  auto Xs = solve_linear<mint>(mat, vc<mint>(N, mint(0)));
  Xs.erase(Xs.begin());
  if (Xs.empty()) return {};
  vvc<mint> F(max_m + 1);
  FOR(i, max_m + 1) {
    vc<mint> f(max_deg + 1);
    FOR(j, max_deg + 1) f[j] = Xs[0][idx(i, j)];
    while (len(f) && f.back() == mint(0)) POP(f);
    F[i] = f;
  }
  while (F[0].empty()) F.erase(F.begin());
  int m = len(F) - 1;
  auto poly_eval = [&](vc<mint>& f, mint x) -> mint {
    mint pw = 1, ans = 0;
    FOR(i, len(f)) ans += f[i] * pw, pw *= x;
    return ans;
  };

  FOR(n, len(A) - m) {
    mint sm = 0;
    FOR(i, m + 1) sm += poly_eval(F[i], n) * A[n + m - i];
    if (sm != 0) return {};
  }
  return F;
}

// A[0], ..., A[N]
template <typename mint>
vc<mint> extend_p_recursive(vc<mint> A, vvc<mint> F, int N) {
  auto poly_eval = [&](vc<mint>& f, mint x) -> mint {
    mint pw = 1, ans = 0;
    FOR(i, len(f)) ans += f[i] * pw, pw *= x;
    return ans;
  };
  int M = len(F) - 1;
  if (len(A) >= N + 1) return {A.begin(), A.begin() + N + 1};
  A.resize(N + 1);
  vc<mint> S(N - M + 1);
  FOR(n, N - M + 1) S[n] = poly_eval(F[0], n);
  S = all_inverse<mint>(S);
  FOR(n, N - M + 1) {
    mint ans = 0;
    FOR(i, 1, M + 1) ans += A[n + M - i] * poly_eval(F[i], n);
    A[n + M] = (-ans) * S[n];
  }
  return A;
}

#line 1 "seq/p_recursive.hpp"

#line 1 "linalg/solve_linear.hpp"
/*
0 行目に解のひとつ。
1～行目に解空間の基底が行ベクトルとして入る。
解なし = empty
*/
template <typename T>
vc<vc<T>> solve_linear(vc<vc<T>> a, vc<T> b, int n = -1, int m = -1) {
  if (n == -1) {
    n = len(a);
    assert(n > 0);
    m = len(a[0]);
  }
  assert(n == len(a) && n == len(b));
  int rk = 0;
  FOR(j, m) {
    if (rk == n) break;
    FOR(i, rk, n) if (a[i][j] != 0) {
      swap(a[rk], a[i]);
      swap(b[rk], b[i]);
      break;
    }
    if (a[rk][j] == 0) continue;
    T c = T(1) / a[rk][j];
    for (auto&& x: a[rk]) x *= c;
    b[rk] *= c;
    FOR(i, n) if (i != rk) {
      T c = a[i][j];
      if (c == T(0)) continue;
      b[i] -= b[rk] * c;
      FOR(k, j, m) { a[i][k] -= a[rk][k] * c; }
    }
    ++rk;
  }
  FOR(i, rk, n) if (b[i] != 0) return {};
  vc<vc<T>> res(1, vc<T>(m));
  vc<int> pivot(m, -1);
  int p = 0;
  FOR(i, rk) {
    while (a[i][p] == 0) ++p;
    res[0][p] = b[i];
    pivot[p] = i;
  }
  FOR(j, m) if (pivot[j] == -1) {
    vc<T> x(m);
    x[j] = -1;
    FOR(k, j) if (pivot[k] != -1) x[k] = a[pivot[k]][j];
    res.eb(x);
  }
  return res;
}
#line 2 "mod/all_inverse.hpp"
template <typename mint>
vc<mint> all_inverse(vc<mint>& X) {
  for (auto&& x: X) assert(x != mint(0));
  int N = len(X);
  vc<mint> res(N + 1);
  res[0] = mint(1);
  FOR(i, N) res[i + 1] = res[i] * X[i];
  mint t = res.back().inverse();
  res.pop_back();
  FOR_R(i, N) {
    res[i] *= t;
    t *= X[i];
  }
  return res;
}
#line 4 "seq/p_recursive.hpp"

/*
return: polynomial sequence F[0](x),...,F[m](x)
F[0](n)A[n+m]+...+F[m](n)A[n]=0
何かが失敗したら return empty
*/
template <typename mint>
vvc<mint> find_p_recursive(vc<mint> A, int max_m, int max_deg) {
  int N = (1 + max_m) * (1 + max_deg);
  auto idx = [&](int i, int d) -> int { return (1 + max_deg) * i + d; };

  assert(len(A) >= N + max_m);
  vv(mint, mat, N, N);
  FOR(row, N) {
    FOR(i, max_m + 1) {
      mint pw = 1;
      FOR(j, max_deg + 1) { mat[row][idx(i, j)] = A[row + max_m - i] * pw, pw *= row; }
    }
  }
  auto Xs = solve_linear<mint>(mat, vc<mint>(N, mint(0)));
  Xs.erase(Xs.begin());
  if (Xs.empty()) return {};
  vvc<mint> F(max_m + 1);
  FOR(i, max_m + 1) {
    vc<mint> f(max_deg + 1);
    FOR(j, max_deg + 1) f[j] = Xs[0][idx(i, j)];
    while (len(f) && f.back() == mint(0)) POP(f);
    F[i] = f;
  }
  while (F[0].empty()) F.erase(F.begin());
  int m = len(F) - 1;
  auto poly_eval = [&](vc<mint>& f, mint x) -> mint {
    mint pw = 1, ans = 0;
    FOR(i, len(f)) ans += f[i] * pw, pw *= x;
    return ans;
  };

  FOR(n, len(A) - m) {
    mint sm = 0;
    FOR(i, m + 1) sm += poly_eval(F[i], n) * A[n + m - i];
    if (sm != 0) return {};
  }
  return F;
}

// A[0], ..., A[N]
template <typename mint>
vc<mint> extend_p_recursive(vc<mint> A, vvc<mint> F, int N) {
  auto poly_eval = [&](vc<mint>& f, mint x) -> mint {
    mint pw = 1, ans = 0;
    FOR(i, len(f)) ans += f[i] * pw, pw *= x;
    return ans;
  };
  int M = len(F) - 1;
  if (len(A) >= N + 1) return {A.begin(), A.begin() + N + 1};
  A.resize(N + 1);
  vc<mint> S(N - M + 1);
  FOR(n, N - M + 1) S[n] = poly_eval(F[0], n);
  S = all_inverse<mint>(S);
  FOR(n, N - M + 1) {
    mint ans = 0;
    FOR(i, 1, M + 1) ans += A[n + M - i] * poly_eval(F[i], n);
    A[n + M] = (-ans) * S[n];
  }
  return A;
}