JZOJ 7405 诡秘的爪钩中心

首先考虑一个显然重要的子问题：
假定已有 ki 条定向的长为 i 的链 (1≤i≤n)，将其任意连接为若干环，求环大小的乘积的和。

事实上我在场上使用多元拉格朗日反演得出了结论，不过在此不表。
考虑一个经典的组合意义：转化为从每个环中选出一个元素的方案数。
枚举 ci 表示长度为 i 的链被选择多少次，则其同时被钦定不在同一环内（这个方案数可以直接通过简单的组合意义得到）： ∑0≤ci≤ki(n∏j=1(kjcj)jcj)(∑jkj−1)!(∑jcj−1)!

也就是 (n∑j=1kj−1)!∑i≥11(i−1)![xi]n∏j=1(1+jx)kj

我们把其放到容斥中，那么看起来我们务必要再用一元计量 ∑jkj。
也就是从环上每断出长 l 的链会贡献 (−1)l−1t(1+lx)。为了快速计算，我们再用一元 u 计量环长。当然这里是当做链来处理了，之后每一项要乘环长再除以链数（也就是 t 的次数）。

注意到 ∑l≥0(−1)l−1t(1+lx)ul=tu1+u+xtu(1+u)2=tu(1+x+u)(1+u)2

拼成链 11−tu(1+x+u)(1+u)2

其任意一项系数是容易提取的（可以按照 t,x,u 的顺序提取），于是分治并做二维卷积计算系数即可。

为了偷懒我直接写了映射到一维来实现二维的卷积，不过这样会有个问题就是卷积长度达到了 222，板子里的某个优化就用不了了（不过 DIT DIF 当然还是能写的）。
所以大概常数大了点。

代码：

#include <queue>
#include <cstdio>
#include <vector>
#include <cstring>
#include <algorithm>

using ll = long long;

using namespace std;

const int mod = 998244353;
inline int norm(int x) {
  return x >= mod ? x - mod : x;
}
inline int reduce(int x) {
  return x < 0 ? x + mod : x;
}
inline int neg(int x) {
  return x ? mod - x : 0;
}
inline void add(int &x, int y) {
  if ((x += y - mod) < 0)
    x += mod;
}
inline void sub(int &x, int y) {
  if ((x -= y) < 0)
    x += mod;
}
inline void fam(int &x, int y, int z) {
  x = (x + (ll)y * z) % mod;
}
inline int qpow(int a, int b) {
  int ret = 1;
  for (; b; b >>= 1)
    (b & 1) && (ret = (ll)ret * a % mod),
    a = (ll)a * a % mod;
  return ret;
}

namespace Poly {
  const int LG = 21;
  const int N = 1 << LG + 1;
  const int G = 3;
  
  int lg2[N + 5];
  int fac[N + 5], ifac[N + 5], inv[N + 5];
  int rt[N + 5];
  
  inline void init() {
    for (int i = 2; i <= N; ++i)
      lg2[i] = lg2[i >> 1] + 1;
    int w = qpow(G, (mod - 1) >> LG + 1);
    rt[N >> 1] = 1;
    for (int i = (N >> 1) + 1; i <= N; ++i)
      rt[i] = (ll)rt[i - 1] * w % mod;
    for (int i = (N >> 1) - 1; i; --i)
      rt[i] = rt[i << 1];
    fac[0] = 1;
    for (int i = 1; i <= N; ++i)
      fac[i] = (ll)fac[i - 1] * i % mod;
    ifac[N] = qpow(fac[N], mod - 2);
    for (int i = N; i; --i)
      ifac[i - 1] = (ll)ifac[i] * i % mod;
    for (int i = 1; i <= N; ++i)
      inv[i] = (ll)ifac[i] * fac[i - 1] % mod;
  }
  
  struct poly {
    
    vector<int> a;
    inline poly(int x = 0) {
      if (x)
        a.push_back(x);
    }
    inline poly(const vector<int> &o) {
      a = o;
    }
    inline poly(const poly &o) {
      a = o.a;
    }
    
    inline int size() const {
      return a.size();
    }
    inline bool empty() const {
      return a.empty();
    }
    inline void resize(int x) {
      a.resize(x);
    }
    inline int operator[](int x) const {
      if (x < 0 || x >= size())
        return 0;
      return a[x];
    }
    inline void clear() {
      vector<int>().swap(a);
    }
    inline poly modxn(int n) const {
      if (a.empty())
        return poly();
      n = min(n, size());
      return poly(vector<int>(a.begin(), a.begin() + n));
    }
    inline poly rever() const {
      return poly(vector<int>(a.rbegin(), a.rend()));
    }
    
    inline void dif() {
      int n = size();
      for (int len = n >> 1; len; len >>= 1)
        for (int j = 0; j < n; j += len << 1)
          for (int k = j, *w = rt + len; k < j + len; ++k, ++w) {
            int R = norm(a[k] + a[k + len]);
            a[k + len] = (ll)*w * (a[k] - a[k + len] + mod) % mod,
            a[k] = R;
          }
    }
    inline void dit() {
      int n = size();
      for (int len = 1; len < n; len <<= 1)
        for (int j = 0; j < n; j += len << 1)
          for (int k = j, *w = rt + len; k < j + len; ++k, ++w) {
            int R = (ll)*w * a[k + len] % mod;
            a[k + len] = reduce(a[k] - R),
            add(a[k], R);
          }
      reverse(a.begin() + 1, a.end());
      for (int i = 0; i < n; ++i)
        a[i] = (ll)a[i] * inv[n] % mod;
    }
    inline void ntt(int type = 1) {
      type == 1 ? dif() : dit();
    }
  };

  struct poly2D {
    int n, m;
    vector<int> a;
    inline poly2D(int r = 0, int s = 0): n(r), m(s) {
      a.resize(n * m);
    }
    inline poly2D(int r, int s, vector<int> vec): n(r), m(s), a(vec) {}
      
    inline int size() const {
      return a.size();
    }
    inline bool empty() const {
      return a.empty();
    }
    inline void resize(int r, int s) {
      n = r, m = s, a.resize(n * m);
    }
    inline int operator[](int x) const {
      if (x < 0 || x >= size())
        return 0;
      return a[x];
    }
    
    friend inline poly2D operator*(const poly2D &a, const poly2D &b) {
      int n = a.n + b.n - 1, m = a.m + b.m - 1;
      poly aBuf, bBuf, resBuf;
      poly2D ret(n, m);
      int tot = n * m, lim = 1;
      for (; lim < tot; lim <<= 1);
      aBuf.resize(lim), bBuf.resize(lim), resBuf.resize(lim);
      for (int i = 0; i < a.n * a.m; ++i) {
        int x = i / a.m, y = i % a.m;
        aBuf.a[x * m + y] = a[i];
      }
      for (int i = 0; i < b.n * b.m; ++i) {
        int x = i / b.m, y = i % b.m;
        bBuf.a[x * m + y] = b[i];
      }
      aBuf.ntt(), bBuf.ntt();
      for (int i = 0; i < lim; ++i)
        fam(resBuf.a[i], aBuf[i], bBuf[i]);
      resBuf.ntt(-1);
      for (int i = 0; i < tot; ++i)
        ret.a[i] = resBuf[i];
      return ret;
    }
  };
}
using Poly::fac;
using Poly::ifac;
using Poly::inv;
using Poly::init;
using Poly::poly;
using Poly::poly2D;

inline int binom(int n, int m) {
  return n < m || m < 0 ? 0 : (ll)fac[n] * ifac[m] % mod * ifac[n - m] % mod;
}

const int N = 2e3;

int n;
int p[N + 5];

struct UnionFind {
  int fa[N + 5], size[N + 5];
  inline UnionFind() {
    for (int i = 1; i <= N; ++i)
      size[i] = 1;
  }
  inline bool isRoot(int x) {
    return !fa[x];
  }
  inline int find(int x) {
    return isRoot(x) ? x : fa[x] = find(fa[x]);
  }
  inline void merge(int x, int y) {
    int fx = find(x), fy = find(y);
    if (fx != fy)
      fa[fx] = fy, size[fy] += size[fx];
  }
} uf;

struct comparer {
  inline bool operator()(const poly2D &x, const poly2D &y) {
    return x.size() > y.size();
  }
};
priority_queue<poly2D, vector<poly2D>, comparer> q;

int ans = 1;

int main() {
  freopen("C.in", "r", stdin), freopen("C.out", "w", stdout);
  init();
  
  scanf("%d", &n);
  for (int i = 1; i <= n; ++i)
    scanf("%d", p + i), uf.merge(i, p[i]);
  for (int i = 1; i <= n; ++i)
    if (uf.isRoot(i)) {
      int size = uf.size[i];
      ans = (ll)ans * size % mod;
      poly2D buf(size + 1, size + 1);
      for (int k = 1; k <= size; ++k)
        for (int t = 0; t <= k; ++t) {
          buf.a[k * (size + 1) + t] = (ll)binom(size + t - 1, t + k - 1) * binom(k, t) % mod * size % mod * inv[k] % mod;
          if ((size - k) & 1)
            buf.a[k * (size + 1) + t] = neg(buf[k * (size + 1) + t]);
        }
      buf.a[0] = size & 1 ? neg(size) : size;
      q.push(buf);
    }
  if (n & 1)
    ans = neg(ans);
  while (q.size() > 1) {
    poly2D x = q.top();
    q.pop(), x = x * q.top(), q.pop(), q.push(x);
  }
  poly2D f = q.top();
  for (int k = 1; k <= n; ++k)
    for (int t = 1; t <= k; ++t)
      ans = (ans + (ll)fac[k - 1] * ifac[t - 1] % mod * f[k * (n + 1) + t]) % mod;
  printf("%d\n", ans);
}