[BZOJ 4025]二分图

只会做水题力……qwq(搞错大小写,自裁请

每条边存在的时间事一个区间,因此考虑线段树分治……首先当然要写一个可撤销并查集。

然后每个点维护他到父亲的边的边权膜2,合并啥的和食物链就很类似力……然后判定已联通两点间边数的奇偶性就直接把两个点到根的值加起来就行了,因为LCA到根的部分会算两遍,对答案无影响。

代码:

/**************************************************************
    Problem: 4025
    User: danihao123
    Language: C++
    Result: Accepted
    Time:9928 ms
    Memory:34840 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <utility>
#include <vector>
#include <stack>
const int maxn = 100005;
const int maxno = maxn << 2;
const int maxm = 200005;
int n, m, T;
 
int p[maxn], siz[maxn], d[maxn];
void init_dsu() {
  for(int i = 1; i <= n; i ++) {
    p[i] = i; siz[i] = 1; d[i] = 0;
  }
}
int get_fa(int x) {
  if(p[x] == x) {
    return x;
  } else {
    return get_fa(p[x]);
  }
}
int get_d(int x) {
  if(p[x] == x) {
    return 0;
  } else {
    return (d[x] + get_d(p[x])) % 2;
  }
}
typedef std::pair<int, int> upd;
upd link_set(int x, int y, int v) {
  if(siz[x] > siz[y]) std::swap(x, y);
  p[x] = y; siz[y] += siz[x]; d[x] = v;
  return std::make_pair(x, y);
}
upd merge_set(int x, int y) {
  int v = (get_d(x) + get_d(y) + 1) % 2;
  return link_set(get_fa(x), get_fa(y), v);
}
void unlink_set(const upd &pr) {
  int x = pr.first, y = pr.second;
  p[x] = x; siz[y] -= siz[x]; d[x] = 0;
}
bool is_same(int x, int y) {
  return (get_fa(x) == get_fa(y));
}
 
std::vector<upd> event[maxno];
std::stack<std::pair<int, upd> > S;
void update(int o, int L, int R, int ql, int qr, upd v) {
  if(ql <= L && R <= qr) {
    event[o].push_back(v);
  } else {
    int M = (L + R) / 2;
    if(ql <= M) update(o << 1, L, M, ql, qr, v);
    if(qr > M) update(o << 1 | 1, M + 1, R, ql, qr, v);
  }
}
int val[maxno];
void solve(int o, int L, int R) {
  for(int i = 0; i < event[o].size(); i ++) {
    const upd &pr = event[o][i];
    int x = pr.first, y = pr.second;
    if(is_same(x, y)) {
      if((get_d(x) + get_d(y)) % 2 == 0) {
        val[o] = 0;
        break;
      }
    } else {
      S.push(std::make_pair(o, merge_set(x, y)));
    }
  }
  if(L == R) {
    if(val[o] != 0) val[o] = 1;
  } else {
    if(val[o] != 0) {
      int M = (L + R) / 2;
      solve(o << 1, L, M);
      solve(o << 1 | 1, M + 1, R);
    }
  }
  while((!S.empty()) && S.top().first >= o) {
    upd v = S.top().second; S.pop();
    unlink_set(v);
  }
}
void print(int o, int L, int R) {
  if(L == R) {
    if(val[o]) {
      puts("Yes");
    } else {
      puts("No");
    }
  } else {
    if(val[o] != -1) {
      val[o << 1] = val[o];
      val[o << 1 | 1] = val[o];
    }
    int M = (L + R) / 2;
    print(o << 1, L, M);
    print(o << 1 | 1, M + 1, R);
  }
}
 
int main() {
  memset(val, -1, sizeof(val));
  scanf("%d%d%d", &n, &m, &T);
  for(int i = 1; i <= m; i ++) {
    int u, v, s, t; scanf("%d%d%d%d", &u, &v, &s, &t);
    s ++; if(s <= t) update(1, 1, T, s, t, std::make_pair(u, v));
  }
  init_dsu(); solve(1, 1, T); print(1, 1, T);
  return 0;
}

[BZOJ 2654]tree

APIO的时候听了一下wqs二分,做了这题,然后现在才写题解……

wqs二分的思路大抵就是你要求必须选\(k\)个,那就二分每个操作的一个“额外代价”,然后进行没有限制的选择。然后最后选出来的个数事和你二分的代价正相关/反相关的。

这道题的话,就二分选择黑边的代价,然后跑一般的最小生成树(有相同边权时要选择黑边!)。当然我们会遇到一个问题,就是二分到\(x\)的时候选的比\(k\)多,到\(x + 1\)的时候又比\(k\)少了。这道题的处理方法,事考虑代价为\(x\)时,其实存在选\(k\)个的最优解(如果说大于\(x\)那就没了),因此我们钦点代价为\(x\)跑一遍限制只能选\(k\)条黑边的最短路即可。

代码:

/**************************************************************
    Problem: 2654
    User: danihao123
    Language: C++
    Result: Accepted
    Time:5756 ms
    Memory:5856 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
#include <vector>
const int maxn = 50005;
const int maxm = 100005;
struct Edge {
  int u, v, d;
  int typ;
  bool operator <(const Edge &res) const {
    if(d == res.d) {
      return typ < res.typ;
    } else {
      return d < res.d;
    }
  }
};
 
Edge E[maxm];
int n, m, need;
 
int par[maxn], siz[maxn];
void init_dsu() {
  for(int i = 1; i <= n; i ++) {
    par[i] = i;
    siz[i] = 1;
  }
}
int get_fa(int x) {
  if(par[x] == x) {
    return x;
  } else {
    return (par[x] = get_fa(par[x]));
  }
}
void link_set(int x, int y) {
  if(siz[x] > siz[y]) std::swap(x, y);
  siz[y] += siz[x];
  par[x] = y;
}
void merge_set(int x, int y) {
  return link_set(get_fa(x), get_fa(y));
}
bool is_same(int x, int y) {
  return (get_fa(x) == get_fa(y));
}
 
typedef long long ll;
int ans = 0x7fffffff;
bool kruskal(int delta) {
  std::vector<Edge> vec;
  for(int i = 1; i <= m; i ++) {
    Edge e = E[i];
    if(e.typ == 0) {
      e.d += delta;
    }
    vec.push_back(e);
  }
  std::sort(vec.begin(), vec.end());
  int used = 0; ll tot = -(ll(delta)) * (ll(need));
  init_dsu();
  for(int i = 0; i < m; i ++) {
    const Edge &e = vec[i];
    int u = e.u, v = e.v;
    if(!is_same(u, v)) {
      if(used == need && e.typ == 0) continue;
      merge_set(u, v); tot += e.d;
      if(e.typ == 0) used ++;
    }
  }
  if(used == need) {
    ans = std::min(ans, (int(tot)));
    return true;
  } else {
    return false;
  }
}
 
int main() {
  scanf("%d%d%d", &n, &m, &need);
  for(int i = 1; i <= m; i ++) {
    scanf("%d%d%d%d", &E[i].u, &E[i].v, &E[i].d, &E[i].typ);
    E[i].u ++; E[i].v ++;
  }
  int L = -10000001, R = 10000001;
  while(true) {
#ifdef LOCAL
    printf("Range (%d, %d)\n", L, R);
    fflush(stdout);
#endif
    if(R - L <= 3) {
      for(int i = R; i >= L; i --) {
        if(kruskal(i)) {
          break;
        }
      }
      break;
    }
    int M = L + (R - L) / 2;
    if(kruskal(M)) {
      L = M;
    } else {
      R = M;
    }
  }
  printf("%d\n", ans);
  return 0;
}

[BZOJ 1004][HNOI2008]Cards

肉肉肉肉,,,

碰置换群啥的不是第一次了……这个题就是给你一堆置换(不要忘记还有一个幺元啊),然后限制颜色数量,求本质不同解个数。那么考虑使用Burnside引理,接下来考虑怎么计算每个置换的不动点数量,这个要求每个循环的颜色一致(不就事Polya定理了吗),所以说可以用背包DP搞一搞。

代码:

/**************************************************************
    Problem: 1004
    User: danihao123
    Language: C++
    Result: Accepted
    Time:156 ms
    Memory:3172 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <utility>
#include <vector>
int sr, sb, sg, m, p;
int pow_mod(int a, int b) {
  int ans = 1, res = a;
  while(b) {
    if(1 & b) ans = (ans * res) % p;
    res = (res * res) % p; b >>= 1;
  }
  return ans;
}
int inv(int x) {
  return pow_mod(x % p, p - 2);
}
 
int d[65][21][21][21];
std::vector<int> len;
int dp() {
  int n = len.size();
  d[0][0][0][0] = 1;
  for(int i = 1; i <= n; i ++) {
    int l = len[i - 1];
    for(int j = 0; j <= sr; j ++) {
      for(int k = 0; k <= sb; k ++) {
        for(int t = 0; t <= sg; t ++) {
          d[i][j][k][t] = 0;
          if(j >= l) d[i][j][k][t] += d[i - 1][j - l][k][t];
          if(k >= l) d[i][j][k][t] += d[i - 1][j][k - l][t];
          if(t >= l) d[i][j][k][t] += d[i - 1][j][k][t - l];
          d[i][j][k][t] %= p;
        }
      }
    }
  }
  return d[n][sr][sb][sg];
}
 
int next[65]; bool vis[65];
int main() {
  scanf("%d%d%d%d%d", &sr, &sb, &sg, &m, &p);
  int n = sr + sb + sg;
  int ans = 0;
  for(int i = 1; i <= n; i ++) {
    len.push_back(1);
  }
  ans = dp();
  for(int i = 1; i <= m; i ++) {
    memset(vis, 0, sizeof(vis));
    for(int i = 1; i <= n; i ++) {
      scanf("%d", &next[i]);
    }
    len.clear();
    for(int i = 1; i <= n; i ++) {
      if(!vis[i]) {
        int p = i, cnt = 0;
        do {
          vis[p] = true; cnt ++;
          p = next[p];
        } while(p != i);
        len.push_back(cnt);
      }
    }
    ans = (ans + dp()) % p;
  }
  ans = (ans * inv(m + 1)) % p;
  printf("%d\n", ans);
  return 0;
}

[BZOJ 3640]JC的小苹果

逆矩阵文明,,,

很显然我们可以定义一个状态\(f[i][j]\)表示当前血量为\(i\)走到\(j\)的概率,然后肥肠爆歉这个东西没法DP(可能会有的点的伤害为0,这样可以凿出来环)。考虑高消,这个东西有个好处是不可能从血量低的到血量高的状态,所以可以从大到小枚举血量,这样各层事独立的,复杂度比直接高消降低了很少。可惜复杂度为\(O(sn^3)\)(设血量为\(s\)),会T掉。

考虑转移的过程,转移时等价于解这样一个方程:

\[
\begin{aligned}
a_{11}x_{1} + a_{12}x_{2} + \ldots + a_{1n}x_{n} &= c_1\\
a_{21}x_{1} + a_{22}x_{2} + \ldots + a_{2n}x_{n} &= c_2\\
&\ldots\\
a_{n1}x_{1} + a_{n2}x_{2} + \ldots + a_{nn}x_{n} &= c_n
\end{aligned}
\]

其中的未知数\(x\)事我们要求的东西,\(c\)表示从高血量状态转移过来的概率(这个可以视作常数)。根据友矩阵那一套理论,这一系列方程等价于以下等式:

\[
\begin{bmatrix}
a_{11} & a_{12} & \ldots & a_{1n}\\
a_{21} & a_{22} & \ldots & a_{2n}\\
\vdots & \vdots & \ddots & \vdots\\
a_{n1} & a_{n2} & \ldots & a_{nn}
\end{bmatrix}
\begin{bmatrix}
x_{1}\\ x_{2} \\ \vdots \\ x_{n}
\end{bmatrix}
=
\begin{bmatrix}
c_1\\ c_2\\ \vdots\\ c_n
\end{bmatrix}
\]

不妨将之简记为\(AB = C\),其中我们只有\(B\)不知道,要求的也是\(B\)。我们除一下就可以得到\(B = A^{-1}C\),然后我们还发现每一层的\(A\)都是一样的,所以每一层的\(A^{-1}\)也都是一样的,预处理即可。这样转移部分的复杂度就变成了\(O(sn^2)\)(矩阵乘法在这里事方阵乘列向量)。

至于逆矩阵的求法,我们知道对矩阵做初等变化也就等价于乘上另一个矩阵。因此,我们将一个矩阵\(A\)用类似于高消的手段消为单位阵\(I\),所做的初等变换也就等价于乘上\(A^{-1}\)。我们对一个单位阵\(I\)作用上一样的操作,也就等于给这个单位阵乘上了\(A^{-1}\),这样我们就得到了\(A^{-1}\)。

代码:

/**************************************************************
    Problem: 3640
    User: danihao123
    Language: C++
    Result: Accepted
    Time:8708 ms
    Memory:13416 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cassert>
#include <cmath>
#include <algorithm>
#include <functional>
#include <utility>
#include <vector>
#include <queue>
#include <set>
#include <map>
int n, m, hp;
const int maxn = 151, maxm = 5005;
const int maxh = 10005;
typedef double R;
typedef R Mat[maxn][maxn];
Mat D, F;
void calc_inv() {
  for(int i = 1; i <= n; i ++) {
    F[i][i] = 1;
  }
  for(int i = 1; i <= n; i ++) {
    int r = i;
    for(int j = i + 1; j <= n; j ++) {
      if(fabs(D[j][i]) > fabs(D[r][i])) {
        r = j;
      }
    }
    // assert(r > -1);
    if(r != i) {
      for(int j = 1; j <= n; j ++) {
        std::swap(D[r][j], D[i][j]);
        std::swap(F[r][j], F[i][j]);
      }
    }
    R bs = D[i][i];
    for(int j = 1; j <= n; j ++) {
      D[i][j] /= bs; F[i][j] /= bs;
    }
    for(int k = 1; k <= n; k ++) {
      if(k != i) {
        R rate = D[k][i];
        for(int j = 1; j <= n; j ++) {
          D[k][j] -= rate * D[i][j];
          F[k][j] -= rate * F[i][j];
        }
      }
    }
  }
}
void matrix_mul(const Mat &A, const Mat &B, int a, int b, int c, Mat &res) {
  static Mat C;
  for(int i = 1; i <= a; i ++) {
    for(int j = 1; j <= c; j ++) {
      C[i][j] = 0;
    }
  }
  for(int i = 1; i <= a; i ++) {
    for(int j = 1; j <= c; j ++) {
      for(int k = 1; k <= b; k ++) {
        C[i][j] += A[i][k] * B[k][j];
      }
    }
  }
  for(int i = 1; i <= a; i ++) {
    for(int j = 1; j <= c; j ++) {
      res[i][j] = C[i][j];
    }
  }
}
 
int first[maxn], deg[maxn];
int next[maxm << 1], to[maxm << 1];
int gcnt = 0;
void add_edge(int u, int v) {
  gcnt ++;
  next[gcnt] = first[u]; first[u] = gcnt;
  to[gcnt] = v;
}
void ins_edge(int u, int v) {
  deg[u] ++;
  add_edge(u, v);
  if(u != v) {
    deg[v] ++;
    add_edge(v, u);
  }
}
 
int atk[maxn];
R f[maxh][maxn];
R solve() {
  for(int i = 1; i <= n; i ++) {
    D[i][i] = 1.00;
  }
  for(int i = 1; i < n; i ++) {
    for(int j = first[i]; j; j = next[j]) {
      int v = to[j];
      if(!atk[v]) {
        D[v][i] -= 1.00 / (R(deg[i]));
      }
    }
  }
  calc_inv();
  R ans = 0; static Mat C;
  f[hp][1] = 1.00;
  for(int i = hp; i >= 1; i --) {
    for(int j = 1; j <= n; j ++) {
      C[j][1] = f[i][j];
    }
    matrix_mul(F, C, n, n, 1, C);
    for(int j = 1; j < n; j ++) {
      for(int e = first[j]; e; e = next[e]) {
        int v = to[e];
        if(atk[v] > 0 && i - atk[v] > 0) {
          f[i - atk[v]][v] += C[j][1] / (R(deg[j]));
        }
      }
    }
    ans += C[n][1];
  }
  return ans;
}
 
int main() {
  scanf("%d%d%d", &n, &m, &hp);
  for(int i = 1; i <= n; i ++) {
    scanf("%d", &atk[i]);
  }
  for(int i = 1; i <= m; i ++) {
    int u, v; scanf("%d%d", &u, &v);
    ins_edge(u, v);
  }
  printf("%.8lf\n", solve());
  return 0;
}

[BZOJ 5358][Lydsy1805月赛]口算训练

后几页有我会做的题……很好,我很安详,,,

考虑将所有数质因数分解。那么询问\([l, r]\)中所有数的积是否是\(d\)的倍数的本质就是对于\(d\)的每一类质因子,询问区间中该类质因子的指数之和是否不小于\(d\)中的。

考虑到数的范围都很小(不超过100000),我们可以先线性筛预处理,这样一次分解质因数的复杂度就降为了\(O(\log n)\)。至于维护区间每类质因子的指数和这种事……就用主席树处理吧。

代码:

/**************************************************************
    Problem: 5358
    User: danihao123
    Language: C++
    Result: Accepted
    Time:2804 ms
    Memory:68408 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <utility>
const int maxn = 100005;
int minp[maxn], mint[maxn], minh[maxn];
int prm[maxn], pcnt;
bool vis[maxn];
void sieve() {
  const int N = 100000;
  vis[1] = true; pcnt = 0;
  for(int i = 2; i <= N; i ++) {
    if(!vis[i]) {
      prm[++ pcnt] = i;
      minp[i] = pcnt; mint[i] = 1;
      minh[i] = i;
    }
    for(int j = 1; j <= pcnt && i * prm[j] <= N; j ++) {
      int v = i * prm[j];
      vis[v] = true; minp[v] = j;
      if(i % prm[j] == 0) {
        mint[v] = mint[i] + 1; minh[v] = minh[i] * prm[j];
        break;
      } else {
        mint[v] = 1; minh[v] = prm[j];
      }
    }
  }
}
 
const int bufsiz = 64 * 1024 * 1024;
char buf[bufsiz]; char *cur;
void init_buf() {
  cur = buf;
}
void *alloc(size_t size) {
  if(cur + size - buf > bufsiz) {
    return malloc(size);
  } else {
    char *ret = cur; cur += size;
    return ret;
  }
}
 
struct Node {
  int v; Node *lc, *rc;
};
Node *build_tree(int L, int R) {
  Node *ret = (Node*)alloc(sizeof(Node));
  ret -> v = 0; 
  if(L == R) {
    ret -> lc = ret -> rc = NULL;
  } else {
    int M = (L + R) / 2;
    ret -> lc = build_tree(L, M);
    ret -> rc = build_tree(M + 1, R);
  }
  return ret;
}
Node *add_p(Node *o, int L, int R, int p, int v) {
  Node *ret = (Node*)alloc(sizeof(Node));
  ret -> v = o -> v;
  ret -> lc = o -> lc; ret -> rc = o -> rc;
  if(L == R) {
    ret -> v += v;
  } else {
    int M = (L + R) / 2;
    if(p <= M) {
      ret -> lc = add_p(ret -> lc, L, M, p, v);
    } else {
      ret -> rc = add_p(ret -> rc, M + 1, R, p, v);
    }
  }
  return ret;
}
int query(Node *o, int L, int R, int p) {
  if(L == R) {
    return o -> v;
  } else {
    int M = (L + R) / 2;
    if(p <= M) {
      return query(o -> lc, L, M, p);
    } else {
      return query(o -> rc, M + 1, R, p);
    }
  }
}
 
Node *rt[maxn];
void solve() {
  init_buf();
  int n, m; scanf("%d%d", &n, &m);
  rt[0] = build_tree(1, pcnt);
  for(int i = 1; i <= n; i ++) {
    rt[i] = rt[i - 1];
    int x; scanf("%d", &x);
    while(x > 1) {
      int p = minp[x], t = mint[x];
      rt[i] = add_p(rt[i], 1, pcnt, p, t);
      x /= minh[x];
    }
  }
  while(m --) {
    int l, r, d; scanf("%d%d%d", &l, &r, &d);
    bool ok = true;
    while(d > 1) {
      int p = minp[d], t = mint[d];
      int st = query(rt[r], 1, pcnt, p) - query(rt[l - 1], 1, pcnt, p);
      if(st < t) {
        ok = false; break;
      }
      d /= minh[d];
    }
    if(ok) {
      puts("Yes");
    } else {
      puts("No");
    }
  }
}
 
int main() {
  sieve();
  int T; scanf("%d", &T);
  while(T --) {
    solve();
  }
  return 0;
}

[BZOJ 4403]序列统计

老早做的题……

一看单调不降就想去用差分……但差分不好推(比下面的颓法要多一步……)。其实我们发现,只要给\([L, R]\)里每种整数分配出现次数,原序列就可以唯一确定了。

因此我们把\([L, R]\)中每个整数的出现次数当做一个变量,他们加起来应该等于一个\([1, n]\)中的整数。用隔板法很容易退出来式子是(令\(l = R - L + 1\)):

\[\sum_{i = 1}^n \binom{i + l - 1}{l - 1}\]

看起来玩不动了……但是我们给式子加上一个\(\binom{l}{l}\)(其实就是1),然后我们会发现式子可以用杨辉三角的递推式合并成一个组合数,即\(\binom{n + l}{l}\)。然后求这个组合数还需要Lucas定理……

代码:

/**************************************************************
    Problem: 4403
    User: danihao123
    Language: C++
    Result: Accepted
    Time:908 ms
    Memory:16448 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
typedef long long ll;
const ll ha = 1000003LL;
const int maxn = 1000003;
ll pow_mod(ll a, ll b) {
  ll ans = 1LL, res = a % ha;
  while(b) {
    if(1LL & b) ans = (ans * res) % ha;
    res = (res * res) % ha;
    b >>= 1;
  }
  return ans;
}
ll f[maxn], g[maxn];
void process() {
  f[0] = 1LL;
  for(int i = 1; i < maxn; i ++) {
    f[i] = (f[i - 1] * (ll(i))) % ha;
  }
  g[maxn - 1] = pow_mod(f[maxn - 1], ha - 2LL);
  for(int i = maxn - 2; i >= 0; i --) {
    g[i] = (g[i + 1] * (ll(i + 1))) % ha;
  }
}
 
ll C(int a, int b) {
  if(a < b) return 0LL;
  return (((f[a] * g[b]) % ha) * g[a - b]) % ha;
}
ll calc(int a, int b) {
  if(a < b) return 0LL;
  if(b == 0) return 1LL;
  if(a < ha && b < ha) {
    return C(a, b);
  } else {
    int as = a % ha, bs = b % ha;
    return (C(as, bs) * calc(a / ha, b / ha)) % ha;
  }
}
 
int main() {
  process();
  int T; scanf("%d", &T);
  while(T --) {
    int n, l, r; scanf("%d%d%d", &n, &l, &r);
    int len = r - l + 1;
    printf("%lld\n", (calc(n + len, len) - 1LL + ha) % ha);
  }
  return 0;
}

[BZOJ 3601]一个人的数论

本来想做JZPKIL的……但卡常太恶心了

上来先颓式子:

\[
\DeclareMathOperator{\id}{id}
\begin{aligned}
f_d(n) &= \sum_{i = 1}^n [(i, n) = 1]i^d\\
&= \sum_{i = 1}^n i^d\sum_{k | n, k | i}\mu(k)\\
&= \sum_{k | n}\mu(k)\sum_{i = 1}^{\frac{n}{k}} (ik)^d\\
&= \sum_{k | n}\mu(k)k^d h_d(\frac{n}{d})\\
&= ((\mu\cdot\id^k)\ast h_d)
\end{aligned}
\]

其中\(h_d\)表示指数为\(d\)时的等幂求和函数。然后我们可能会想,如果这个函数是积性函数,那么我们可以对质因数分解的每一种质因数都求一遍答案,而由于\(\mu\)的存在(对于质数的幂\(p^k\)的因数,只有\(1\)和\(p\)的莫比乌斯函数值不为0),所以需要处理的情况只有两种。

不过很可惜,\(h_d\)显然不是积性函数,所以最后的整个卷积也不是积性函数。但是我们注意到\(h_d\)事一个\(d + 1\)次多项式,所以我们可以把多项式每一项都当成单独的函数,然后每个单独函数卷出来都是一个积性函数。至于求多项式每一项的系数,用伯努利数即可。

代码:

/**************************************************************
    Problem: 3601
    User: danihao123
    Language: C++
    Result: Accepted
    Time:220 ms
    Memory:948 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <utility>
typedef long long ll;
ll ha = 1000000007LL;
ll pow_mod(ll a, ll b) {
  ll ans = 1LL, res = a;
  while(b) {
    if(1LL & b) ans = (ans * res) % ha;
    res = (res * res) % ha;
    b >>= 1;
  }
  return ans;
}
ll inv(ll x) {
  return pow_mod(x, ha - 2LL);
}
 
const int maxn = 105;
ll C[maxn][maxn];
void process_C() {
  C[0][0] = 1LL;
  for(int i = 1; i < maxn; i ++) {
    C[i][0] = C[i][i] = 1LL;
    for(int j = 1; j < i; j ++) {
      C[i][j] = (C[i - 1][j] + C[i - 1][j - 1]) % ha;
    }
  }
}
ll B[maxn];
void process_B() {
  B[0] = 1LL;
  for(int i = 1; i <= 101; i ++) {
    B[i] = 1LL;
    for(int j = 0; j < i; j ++) {
      ll d = (((C[i][j] * B[j]) % ha) * inv(i - j + 1)) % ha;
      B[i] = (B[i] - d + ha) % ha;
    }
  }
}
 
const int maxw = 1005;
typedef std::pair<ll, ll> pii;
pii P[maxw]; ll dv[maxw][3];
int d, w;
void process_dv() {
  for(int i = 1; i <= w; i ++) {
    ll p = P[i].first, t = P[i].second;
    p %= ha;
    dv[i][0] = pow_mod(p, t);
    dv[i][1] = pow_mod(p, t - 1LL);
    dv[i][2] = pow_mod(p, d);
  }
}
ll calc(ll k, int z) {
  ll ans = k;
  for(int i = 1; i <= w; i ++) {
    ll p = P[i].first, t = P[i].second;
    p %= ha;
    ll f1 = pow_mod(dv[i][0], z), f2 = pow_mod(dv[i][1], z);
    ll v1 = (f1 - (dv[i][2] * f2) % ha + ha) % ha;
    ans = (ans * v1) % ha;
  }
#ifdef LOCAL
  printf("Case (%lld, %d) : %lld\n", k, z, ans);
#endif
  return ans;
}
 
int main() {
  process_C(); process_B();
  scanf("%d%d", &d, &w);
  for(int i = 1; i <= w; i ++) {
    scanf("%lld%lld", &P[i].first, &P[i].second);
  }
  process_dv();
  ll ans = 0LL;
  for(int i = 0; i <= d; i ++) {
    ll g = calc((((C[d + 1][i] * B[i]) % ha) * inv(d + 1)) % ha, d + 1 - i);
    ans = (ans + g) % ha;
  }
  printf("%lld\n", ans);
  return 0;
}

[BZOJ 2321][BeiJing2011集训]星器

物理题可还行(其实我是想学势能方法,然后误入了……

给每个星星(假设他在的位置是\((i, j)\))定义势能\(E_p = i^2 + j^2\),定义势函数\(\Phi(S)\)来表示状态\(S\)时所有星星势能的和。

然后我们发现,当两个星星互相吸引时(假设他们同行,坐标分别为\((i, j)\)和\((i, k)\),且$j<k$),他们的坐标会变为\((i, j + 1)\)和\((i, k - 1)\),势能的总减少量(也就是势函数的减小量)为\(2(k - j - 1)\)。

因此,整个过程中势函数的总减小量,就是答案的两倍。因此算出操作前后的势能做差即可……

代码:

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <utility>
typedef long long ll;
int main() {
  int n, m; scanf("%d%d", &n, &m);
  ll E = 0LL;
  for(int i = 1; i <= n; i ++) {
    for(int j = 1; j <= m; j ++) {
      ll delta = i * i + j * j; ll x;
      scanf("%lld", &x); delta *= x;
      E += delta;
    }
  }
  for(int i = 1; i <= n; i ++) {
    for(int j = 1; j <= m; j ++) {
      ll delta = i * i + j * j; ll x;
      scanf("%lld", &x); delta *= x;
      E -= delta;
    }
  }
  E /= 2LL;
  printf("%lld\n", E);
  return 0;
}

[BZOJ 3456]城市规划

aji多项式求逆毁我青春,,,

设\(f_i\)表示\(i\)个点的有标号无向联通图,考虑所有可能的图(记\(F_i\)为\(i\)个点的有标号无向图,显然\(F_i = 2^{\binom{i}{2}}\))和\(f\)的关系(使用图计数的经典套路:枚举1所在的联通块大小):

\[F_n = \sum_{i = 1}^n \binom{n-1}{i-1} f_i F_{n - i}\]

看起来事卷积?但是这个卷积没有办法直接用FFT/NTT求(当然分离一下项啊,移下项就可以分治NTT力)。

考虑进一步化简柿子。完全展开后会发现右边有个非常碍眼的\((n-1)!\),所以两边除一下:

\[\frac{F_n}{(n-1)!} =\sum_{i = 1}^n \frac{f_i}{(i-1)!}\cdot \frac{F_{n - i}}{(n - i)!}\]

然后这个问题就很毒瘤了:我们要求的答案多项式(除上那个阶乘)和一个已知多项式做卷积,可以得到另一个已知多项式……这样就需要多项式除法了,于是乎多项式逆元派上了用场。

代码:

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
typedef long long ll;
const int maxn = 1020010;
const ll ha = 1004535809LL;
const ll bs = 3LL;
ll pow_mod(ll a, ll b) {
  ll ans = 1LL, res = a % ha;
  while(b) {
    if(1LL & b) ans = (ans * res) % ha;
    res = (res * res) % ha;
    b >>= 1;
  }
  return ans;
}
ll inv(ll x) {
  return pow_mod(x,  ha - 2LL);
}
 
int flip(int bi, int x) {
  int ans = 0;
  for(int i = 0; i < bi; i ++) {
    if((1 << i) & x) {
      ans += (1 << (bi - i - 1));
    }
  }
  return ans;
}
void NTT(ll *A, int bi, bool flag = false) {
  int n = 1 << bi;
  for(int i = 0; i < n; i ++) {
    int v = flip(bi, i);
    if(v < i) std::swap(A[v], A[i]);
  }
  for(int L = 1; L < n; L <<= 1) {
    ll xi = pow_mod(3LL, (ha - 1LL) / (ll(L << 1)));
    if(flag) xi = inv(xi);
    for(int i = 0; i < n; i += (L << 1)) {
      ll w = 1LL;
      for(int j = i; j < i + L; j ++) {
        ll v1 = A[j], v2 = A[j + L];
        A[j] = (v1 + (w * v2) % ha) % ha;
        A[j + L] = (v1 - (w * v2) % ha + ha) % ha;
        w = (w * xi) % ha;
      }
    }
  }
}
void poly_mul(ll *A, ll *B, int bi, ll *C) {
  static ll T1[maxn], T2[maxn];
  int n = (1 << bi);
  std::copy(A, A + n, T1);
  std::copy(B, B + n, T2);
  NTT(T1, bi); NTT(T2, bi);
#ifdef LOCAL
  puts("poly_mul :");
  for(int i = 0; i < (n); i ++) {
    printf("%lld ", A[i]);
  }
  puts("");
  for(int i = 0; i < (n); i ++) {
    printf("%lld ", B[i]);
  }
  puts("");
  for(int i = 0; i < (n); i ++) {
    printf("%lld ", T1[i]);
  }
  puts("");
  for(int i = 0; i < (n); i ++) {
    printf("%lld ", T2[i]);
  }
  puts("");
#endif 
  for(int i = 0; i < n; i ++) {
    T1[i] = (T1[i] * T2[i]) % ha;
  }
#ifdef LOCAL
  for(int i = 0; i < (n); i ++) {
    printf("%lld ", T1[i]);
  }
  puts("");
#endif 
  NTT(T1, bi, true);
  ll inv_n = inv(n);
  for(int i = 0; i < n; i ++) {
    T1[i] = (T1[i] * inv_n) % ha;
  }
  std::copy(T1, T1 + n, C);
#ifdef LOCAL
  for(int i = 0; i < (n); i ++) {
    printf("%lld ", C[i]);
  }
  puts("");
#endif 
}
 
void poly_inv(int mod, ll *B, ll *BB) {
  if(mod == 1) {
    BB[0] = inv(B[0]);
  } else {
    poly_inv((mod + 1) >> 1, B, BB);
    int bi = 0, sz = 1;
    while(sz <= ((mod * 2) + 1)) {
      bi ++; sz <<= 1;
    }
    ll inv_sz = inv(sz);
    static ll tmp[maxn];
    std::copy(B, B + mod, tmp);
    std::fill(tmp + mod, tmp + sz, 0LL);
    NTT(tmp, bi); NTT(BB, bi);
    for(int i = 0; i < sz; i ++) {
      tmp[i] = (tmp[i] * BB[i]) % ha;
      tmp[i] = (tmp[i] * (ha - 1LL)) % ha;
      tmp[i] = (tmp[i] + 2LL) % ha;
      tmp[i] = (tmp[i] * BB[i]) % ha;
    }
    NTT(tmp, bi, true);
    for(int i = 0; i < sz; i ++) {
      tmp[i] = (tmp[i] * inv_sz) % ha;
    }
    std::copy(tmp, tmp + mod, BB);
    std::fill(BB + mod, BB + sz, 0LL);
  }
}
 
int main() {
  static ll fac[maxn], A[maxn], B[maxn], BB[maxn];
  int n; scanf("%d", &n);
  int bi = 0, sz = 1;
  while(sz <= n + 1) {
    bi ++; sz <<= 1;
  }
  fac[0] = 1LL;
  for(int i = 1; i <= n; i ++) {
    fac[i] = (fac[i - 1] * (ll(i))) % ha;
  }
  B[0] = 1LL;
  for(int i = 1; i <= n; i ++) {
    B[i] = pow_mod(2LL, (ll(i)) * (ll(i - 1)) / 2LL);
    B[i] = (B[i] * inv(fac[i])) % ha;
  }
  for(int i = 1; i <= n; i ++) {
    A[i] = pow_mod(2LL, (ll(i)) * (ll(i - 1)) / 2LL);
    A[i] = (A[i] * inv(fac[i - 1])) % ha;
  }
  poly_inv(n + 1, B, BB);
  poly_mul(A, BB, bi, A);
  printf("%lld\n", (A[n] * fac[n - 1]) % ha);
  return 0;
}

[BZOJ 3328]PYXFIB

又学了个新套路呢qwq

题面要求的其实是:

\[\sum_{i = 0}^n [k|i]\binom{n}{i} F_i\]

不考虑那个\([k|i]\),式子是非常像一个二项式展开的……

考虑构造矩阵的二项式展开,可以发现(其中\(A\)是斐波那契数列的二项式展开):

\[(I + A)^n = \sum_{i = 0}^n \binom{n}{i} A^i\]

然后\(k=1\)的情况我们就会做辣!接下来的问题是,如何取出一个生成函数中次数为\(k\)倍数的项求和?

然后我们发现单位根有个非常优良的性质:

\[\frac{1}{k} \sum_{i=1}^{k} \xi_k^{ni} = [k|n]\]

这里的\(n\)是一个常数,但是其实可以对应到原式里的某一项的次数。至于证明……满足\(k|n\)的情况左边显然是一堆1取平均值,其他情况可以通过等比数列求和证明左边为0。

于是可以构造多项式\((I+xA)^n\),把所有\(k\)次单位根做为\(x\)(求这个的话,可以利用原根)带入,对结果取平均值即可。

代码:

/**************************************************************
    Problem: 3328
    User: danihao123
    Language: C++
    Result: Accepted
    Time:9080 ms
    Memory:832 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <utility>
#include <cmath>
#include <vector>
typedef long long ll;
typedef ll Mat[2][2];
ll p; int k;
ll pow_mod(ll a, ll b) {
  ll ans = 1LL, res = a % p;
  while(b) {
    if(1LL & b) ans = (ans * res) % p;
    res = (res * res) % p;
    b /= 2LL;
  }
  return ans;
}
ll inv(ll x) {
  return pow_mod(x, p - 2LL);
}
 
void factor(int x, std::vector<int> &V) {
  int bd = sqrt(x + 0.5);
  for(int i = 2; i <= bd; i ++) {
    if(x % i == 0) {
      V.push_back(i);
      while(x % i == 0) x /= i;
    }
  }
  if(x > 1) V.push_back(x);
}
int get_phi() {
  if(p == 2LL) return 1LL;
  std::vector<int> V;
  factor(p - 1LL, V);
  for(int i = 2; i <= (p - 1); i ++) {
    bool flag = true;
    for(int j = 0; j < V.size(); j ++) {
      int up = (p - 1) / V[j];
      if(pow_mod(i, up) == 1LL) {
        flag = false; break;
      }
    }
#ifdef LOCAL
    if(flag) printf("xi : %d\n", i);
#endif
    if(flag) return i;
  }
}
 
void mat_mul(const Mat &A, const Mat &B, int n, int m, int t, Mat &C) {
  Mat D; memset(D, 0, sizeof(D));
  for(int i = 0; i < n; i ++) {
    for(int j = 0; j < t; j ++) {
      for(int k = 0; k < m; k ++) {
        D[i][j] = (D[i][j] + (A[i][k] * B[k][j]) % p) % p;
      }
    }
  }
  memcpy(C, D, sizeof(D));
}
void mat_pow(const Mat &A, int n, ll b, Mat &ret) {
  Mat C, res;
  memset(C, 0, sizeof(C));
  for(int i = 0; i < n; i ++) C[i][i] = 1LL;
  memset(res, 0, sizeof(res));
  memcpy(res, A, sizeof(res));
  while(b) {
    if(1LL & b) mat_mul(C, res, n, n, n, C);
    mat_mul(res, res, n, n, n, res);
    b /= 2LL;
  }
  memcpy(ret, C, sizeof(C));
}
 
int main() {
  int T; scanf("%d", &T);
  while(T --) {
    ll n;
    scanf("%lld%d%lld", &n, &k, &p);
    ll xi = pow_mod(get_phi(), (p - 1) / k);
    ll w = 1LL;
    ll ans = 0;
    for(int i = 0; i < k; i ++) {
      Mat X;
      memset(X, 0, sizeof(X));
      X[0][0] = X[1][0] = X[0][1] = w;
      X[0][0] = (X[0][0] + 1LL) % p;
      X[1][1] = (X[1][1] + 1LL) % p;
      mat_pow(X, 2, n, X);
#ifdef LOCAL
      puts("X : ");
      for(int i = 0; i < 2; i ++) {
        for(int j = 0; j < 2; j ++) {
          printf("%lld ", X[i][j]);
        }
        puts("");
      }
#endif
      ans = (ans + X[0][0]) % p;
      w = (w * xi) % p;
    }
    ans = (ans * inv(k)) % p;
    printf("%lld\n", ans);
  }
  return 0;
}