[51Nod 1355]斐波那契的最小公倍数

MIN-MAX容斥第一题……

众所周知斐波那契数有个非常好的性质:

\[\gcd(\{f_i | i\in T\}) = f_{\gcd\{T\}}\]

但很不幸的事情是,lcm没有类似的性质。那我们就考虑怎么把lcm改成gcd。

根据MIN-MAX容斥,可以得到这个式子(可以理解为对指数的容斥。这里设\(S\)为全集):

\[\prod_{T\subseteq S, T\neq\emptyset} f_{\gcd\{T\}}^{(-1)^{|T| + 1}}\]

gcd套在上面非常烦人,考虑用狄利克雷卷积去化解它。考虑有个函数\(g\)满足\(f_n = \prod_{d|n} g(d)\),然后xjb化简一下之后柿子变成:

\[\prod_{d} g(d)^{\sum_{T\subseteq S, T\neq\emptyset, \forall x\in T, d | x}(-1)^{|T| + 1}}\]

再去考虑那个指数。当且仅当\(d\in S\)时,\(g(d)\)才会对答案做贡献。并且无论\(S\)中\(d\)的约数有多少个(只要大于),最后指数肯定都是1。因为其实大于零的情况那个指数也是一个所有最小值都是1的MIN-MAX容斥……这样求出来最大值肯定是1了,那个指数自然就是1。

于是乎柿子变成了:

\[\prod_{\exists x\in S, d | x} g(d)\]

然后搞出来\(S\)中所有数的约数很容易。当务之急就是求\(g\)了(这里其实就是莫比乌斯反演了)。不过球莫比乌斯函数反而让问题复杂化了……移下项可以发现:

\[g(n) = f_n\cdot (\prod_{d | n, d\neq n} g(d))^{-1}\]

然后就很愉快的搞出来了(逃

代码:

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <cctype>
#include <algorithm>
#include <utility>
#include <queue>
const int maxn = 1000005;
using ll = long long;
const 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 /= 2LL;
  }
  return ans;
}
ll inv(ll x) {
  return pow_mod(x, ha - 2LL);
}

ll f[maxn], g[maxn];
void process_f() {
  int N = 1000000;
  f[0] = 0LL; f[1] = 1LL;
  for(int i = 2; i <= N; i ++) {
    f[i] = (f[i - 1] + f[i - 2]) % ha;
  }
  std::fill(g + 1, g + 1 + N, 1LL);
  for(int i = 1; i <= N; i ++) {
    g[i] = (f[i] * inv(g[i])) % ha;
    for(int j = i * 2; j <= N; j += i) {
      g[j] = (g[j] * g[i]) % ha;
    }
  }
#ifdef LOCAL
  puts("g has been processed!");
#endif
}

bool vis[maxn];
void process_vis() {
  int N = 1000000;
  for(int i = 1; i <= N; i ++) {
    for(int j = i * 2; j <= N; j += i) {
      vis[i] = (vis[i] || vis[j]);
    }
  }
}

int main() {
  process_f();
  int n; scanf("%d", &n);
  for(int i = 1; i <= n; i ++) {
    int a; scanf("%d", &a);
    vis[a] = true;
  }
  process_vis();
  ll ans = 1LL;
  for(int i = 1; i <= 1000000; i ++) {
    if(vis[i]) {
      ans = (ans * g[i]) % ha;
    }
  }
  printf("%lld\n", ans);
  return 0;
}

[BZOJ 4259]残缺的字符串

终于肝掉了这道提……

考虑用编辑距离函数处理这个问题。考虑每个可能的结尾\(i\),定义\(f_i\)(这里把\(*\)看做\(0\)把):

\[f_i = \sum_{j = 0}^{m - 1} A_j B_{i - m + 1 + j}\,(A_j - B_{i - m + 1 + j})^2\]

考虑翻转\(A\),式子变成:

\[f_i = \sum_{j = 0}^{m - 1} A_{m - 1 - j}\,B_{i - m + 1 + j}\,(A_j - B_{i - m + 1 + j})^2\]

式子给人一种卷积的即视感……再化简一下,得到:

\[f_i = \sum_{j = 0}^{m - 1} A_{m - 1 - j}^3\,B_{i - m + 1 + j}\,- 2A_{m - 1 - j}^2\,B_{i - m + 1 + j}^2\,+ A_{m - 1 - j}\,B_{i - m + 1 + j}^3\]

这是三个卷积……跑三遍FFT即可得到\(f_i\),当且仅当\(f_i = 0\)时,以\(i\)结尾的长度为\(m\)的子串才能匹配。

顺便提个几个坑……这个题用long double + std::complex会因为常数巨大而T掉。并且注意eps不要调那么细啊……这个东西精度误差极大,我直接调成1了。

代码:

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <cmath>
#include <algorithm>
#include <utility>
#include <queue>
#include <complex>
constexpr int maxn = 300005;
using R = double;
using C = std::complex<R>;
constexpr R pi = acos(-1);
inline int flip(int bi, int x) {
  int ans = 0;
  for(int i = 0; i < bi; i ++) {
    if(x & (1 << i)) {
      ans += (1 << (bi - i - 1));
    }
  }
  return ans;
}
inline void FFT(C *A, int bi, int n, R flag = 1.00) {
  for(int i = 0; i < n; i ++) {
    int fl = flip(bi, i);
    if(i < fl) std::swap(A[fl], A[i]);
  }
  for(int L = 1; L < n; L <<= 1) {
    R rd = pi / (R(L));
    C xi(cos(rd), sin(flag * rd));
    for(int i = 0; i < n; i += (L << 1)) {
      C w(1.00, 0.00);
      for(int j = i; j < i + L; j ++) {
        C v1 = A[j], v2 = A[j + L];
        A[j] = v1 + v2 * w;
        A[j + L] = v1 - v2 * w;
        w = w * xi;
      }
    }
  }
}

constexpr R eps = 1.00;
inline int sign(R x) {
  if(fabs(x) < eps) {
    return 0;
  } else {
    if(x < 0.00) {
      return -1;
    } else {
      return 1;
    }
  }
}
int main() {
  static char A[maxn], B[maxn];
  static int V1[maxn << 2], V2[maxn << 2];
  static C P1[maxn << 2], P2[maxn << 2], P3[maxn << 2];
  static R F[maxn << 2];
  int m, n; scanf("%d%d", &m, &n);
  if(m > n) {
    puts("0"); return 0;
  }
  scanf("%s%s", A, B);
  std::reverse(A, A + m);
  for(int i = 0; i < m; i ++) {
    V1[i] = (A[i] == '*') ? 0 : (A[i] - 'a' + 1);
  }
  for(int i = 0; i < n; i ++) {
    V2[i] = (B[i] == '*') ? 0 : (B[i] - 'a' + 1);
  }
  int sz = 1, bi = 0;
  while(sz <= (m + n)) {
    sz *= 2; bi ++;
  }
  std::fill(P1, P1 + sz, C(0, 0));
  std::fill(P2, P2 + sz, C(0, 0));
  for(int i = 0; i < m; i ++) {
    P1[i] = V1[i] * V1[i] * V1[i];
  }
  for(int i = 0; i < n; i ++) {
    P2[i] = V2[i];
  }
  FFT(P1, bi, sz); FFT(P2, bi, sz);
  for(int i = 0; i < sz; i ++) {
    P3[i] = P1[i] * P2[i];
  }
  FFT(P3, bi, sz, -1.00);
  for(int i = 0; i < sz; i ++) {
    F[i] += P3[i].real() / (R(sz));
  }
  std::fill(P1, P1 + sz, C(0, 0));
  std::fill(P2, P2 + sz, C(0, 0));
  for(int i = 0; i < m; i ++) {
    P1[i] = V1[i] * V1[i];
  }
  for(int i = 0; i < n; i ++) {
    P2[i] = V2[i] * V2[i];
  }
  FFT(P1, bi, sz); FFT(P2, bi, sz);
  for(int i = 0; i < sz; i ++) {
    P3[i] = P1[i] * P2[i];
  }
  FFT(P3, bi, sz, -1.00);
  for(int i = 0; i < sz; i ++) {
    F[i] += (R(-2)) * P3[i].real() / (R(sz));
  }
  std::fill(P1, P1 + sz, C(0, 0));
  std::fill(P2, P2 + sz, C(0, 0));
  for(int i = 0; i < m; i ++) {
    P1[i] = V1[i];
  }
  for(int i = 0; i < n; i ++) {
    P2[i] = V2[i] * V2[i] * V2[i];
  }
  FFT(P1, bi, sz); FFT(P2, bi, sz);
  for(int i = 0; i < sz; i ++) {
    P3[i] = P1[i] * P2[i];
  }
  FFT(P3, bi, sz, -1.00);
  for(int i = 0; i < sz; i ++) {
    F[i] += P3[i].real() / (R(sz));
  }
  static int ans[maxn]; int cnt = 0;
  for(int i = 0; i < n; i ++) {
#ifdef LOCAL
    printf("F[%d] : %.10Lf\n", i, F[i]);
#endif
    if(sign(F[i]) == 0 && i >= m - 1) {
      ans[cnt ++] = i - m + 1 + 1;
    }
  }
  printf("%d\n", cnt);
  for(int i = 0; i < cnt; i ++) {
    if(i) putchar(' ');
    printf("%d", ans[i]);
  }
  putchar('\n');
  return 0;
}

[BZOJ 3622]已经没有什么好害怕的了

已经……没有什么好害怕的了

OI给了我足够的欢乐了,OI让我看到了新的世界,让我有了人生中最好的几位朋友,我没有什么可以遗憾的了……就算可能要面临最残酷的结局吧……

有点跑题了呢……还是说这个题吧

显然可以看出来这个题要求你的从大到小的匹配(姑且称之为好的匹配)恰好要有\(\frac{n + k}{2}\)个,剩下的全都不是好的匹配。

首先把糖果(记为\(A\))和药片(记为\(B\))分别排序,对于每一个\(A\)中元素就很容易得到\(B\)中有多少比它小的。考虑设计状态\(f[i][j]\)表示对(排序后的)\(A\)的前\(i\)项,已经确定了的好的匹配有至少\(j\)个。这个DP转移显然。

然后发现\(f[n][j]\times (n - j)!\)就是对于所有元素的完美匹配中,有至少\(j\)个是好的匹配的方案数。这正是广义容斥原理的用武之地。考虑记\(t_j = f[n][j]\times (n - j)!\),那么答案就是(令\(\frac{n + k}{2} = l\)):

\[\sum_{i = l}^n (-1)^{n - i}\binom{i}{l} t_i\]

代码:

/**************************************************************
    Problem: 3622
    User: danihao123
    Language: C++
    Result: Accepted
    Time:1844 ms
    Memory:63680 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
typedef long long ll;
const ll ha = 1000000009LL;
const int maxn = 2005;
ll C[maxn][maxn], F[maxn];
inline void process() {
  int N = 2000;
  C[0][0] = 1LL;
  for(int i = 1; i <= N; 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;
    }
  }
  F[0] = 1LL;
  for(int i = 1; i <= N; i ++) {
    F[i] = (F[i - 1] * (ll(i))) % ha;
  }
}
 
int A[maxn], B[maxn];
ll cp[maxn]; int n;
void gen_cp() {
  std::sort(A + 1, A + 1 + n); std::sort(B + 1, B + 1 + n);
  for(int i = 1; i <= n; i ++) {
    cp[i] = n - (B + n + 1 - std::upper_bound(B + 1, B + 1 + n, A[i]));
#ifdef LOCAL
    printf("cp[%d] : %lld\n", i, cp[i]);
#endif
  }
}
 
ll f[maxn][maxn];
ll dp(int dec) {
  process(); gen_cp();
  f[0][0] = 1LL;
  for(int i = 1; i <= n; i ++) {
    for(int j = 0; j <= n; j ++) {
      f[i][j] = f[i - 1][j];
      if(j > 0) {
        f[i][j] += (f[i - 1][j - 1] * std::max(0LL, cp[i] - ll(j - 1))) % ha;
        if(f[i][j] >= ha) f[i][j] -= ha;
      }
#ifdef LOCAL
      printf("f[%d][%d] : %lld\n", i, j, f[i][j]);
#endif
    }
  }
  for(int j = 0; j <= n; j ++) {
    f[n][j] = (f[n][j] * F[n - j]) % ha;
  }
  ll ans = 0;
  for(int j = dec; j <= n; j ++) {
    ll a = ((j - dec) % 2 == 0) ? 1LL : (ha - 1LL);
    ll delta = (C[j][dec] * f[n][j]) % ha;
    ans += (a * delta) % ha;
    if(ans >= ha) ans -= ha;
  }
  return ans;
}
 
int main() {
  int k; scanf("%d%d", &n, &k);
  if((n - k) & 1) {
    puts("0"); return 0;
  }
  for(int i = 1; i <= n; i ++) {
    scanf("%d", &A[i]);
  }
  for(int i = 1; i <= n; i ++) {
    scanf("%d", &B[i]);
  }
  printf("%lld\n", dp((n + k) / 2));
  return 0;
}

[HDU 4625]JZPTREE

第一次用到第二类斯特林数……

首先众所周知第二类斯特林数有一个很好康的式子:

\[x^k = \sum_{i = 0}^k S(k, i) x^{\underline{i}}\]

然后带进去试试?一番化简之后得到:

\[E_u = \sum_{i = 0}^k S(k, i) \sum_{v = 1}^n dist(u, v)^{\underline{i}}\]

那个下降幂很恶心……把它强行搞成组合数后发现:

\[E_u = \sum_{i = 0}^k S(k, i) i! \sum_{v = 1}^n \binom{dist(u, v)}{i}\]

然后考虑每个点预处理后面那个和式……显然可以树形DP苟。用经典的树形DP套路,设两个状态分别表示一个点的子树内部的贡献以及外面的点给他的贡献。至于转移,由于那是个组合数,所以可以考虑用Pascal定理苟……

代码:

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <cmath>
#include <algorithm>
#include <utility>
const int maxn = 50005;
const int maxk = 505;
const int ha = 10007;
int first[maxn];
int next[maxn << 1], to[maxn << 1];
int cnt;
void add_edge(int u, int v) {
  cnt ++;
  next[cnt] = first[u]; first[u] = cnt;
  to[cnt] = v;
}

int n, k;
int down[maxn][maxk];
void dfs_1(int x, int fa = 0) {
  down[x][0] = 1;
  for(int i = first[x]; i; i = next[i]) {
    int v = to[i];
    if(v != fa) {
      dfs_1(v, x);
      down[x][0] = (down[x][0] + down[v][0]) % ha;
      for(int j = 1; j <= k; j ++) {
        down[x][j] = (down[x][j] + down[v][j] + down[v][j - 1]) % ha;
      }
    }
  }
}
int up[maxn][maxk];
void dfs_2(int x, int fa = 0) {
  if(fa) {
    up[x][0] = n - down[x][0];
    for(int j = 1; j <= k; j ++) {
      up[x][j] = (up[fa][j] + up[fa][j - 1] + down[fa][j] + down[fa][j - 1]) % ha;
      up[x][j] = (up[x][j] - (2 * down[x][j - 1]) % ha + ha) % ha;
      up[x][j] = (up[x][j] - down[x][j] + ha) % ha;
      if(j > 1) up[x][j] = (up[x][j] - down[x][j - 2] + ha) % ha;
    }
  }
  for(int i = first[x]; i; i = next[i]) {
    int v = to[i];
    if(v != fa) {
      dfs_2(v, x);
    }
  }
}

int S[maxk][maxk];
int fac[maxn];
void process() {
  S[1][1] = 1;
  for(int i = 2; i <= k; i ++) {
    S[i][0] = 0;
    for(int j = 1; j < i; j ++) {
      S[i][j] = ((j * S[i - 1][j]) % ha + S[i - 1][j - 1]) % ha;
    }
    S[i][i] = 1;
  }
  fac[0] = 1;
  for(int i = 1; i <= n; i ++) {
    fac[i] = (fac[i - 1] * (i % ha)) % ha;
  }
}

int main() {
  int T; scanf("%d", &T);
  n = 50000; k = 500; process();
  while(T --) {
    scanf("%d%d", &n, &k);
    cnt = 0; memset(first, 0, sizeof(first));
    for(int i = 0; i < n - 1; i ++) {
      int u, v; scanf("%d%d", &u, &v);
      add_edge(u, v); add_edge(v, u);
    }
    memset(down, 0, sizeof(down));
    memset(up, 0, sizeof(up));
    dfs_1(1); dfs_2(1);
    for(int i = 1; i <= n; i ++) {
      int ans = 0;
      for(int j = 1; j <= k; j ++) {
        int f_v = (down[i][j] + up[i][j]) % ha;
        ans = (ans + (S[k][j] * ((f_v * fac[j]) % ha)) % ha) % ha;
      }
      printf("%d\n", ans);
    }
  }
  return 0;
}

[BZOJ 5091][Lydsy0711月赛]摘苹果

秒,秒啊.jpg

首先根据期望线性性,每个点的期望可以分开算。不妨设\(f_{i, j}\)表示\(j\)在某个第\(i\)步被走到的概率。那么每个点\(j\)的期望访问次数就是\(\sum_{i = 0}^k f_{i, j}\)。

然后考虑去求那个\(f_{i, j}\)。显然\(f_{0, j} = \frac{d_j}{2m}\)。但是考虑\(f_{1, j}\)的转移方程:

\[f_{1, j} = \sum \frac{f_{0, u}}{d_u}\]

然后展开之后发现每个\(\frac{f_{0, u}}{d_u}\)都是\(\frac{1}{2m}\),于是乎\(f_{1, j} = \frac{d_j}{2m}\)。

如此一来,对于任意\(i\),都有\(d_{i, j} = \frac{d_j}{2m}\),然后就很好做了……

代码:

/**************************************************************
    Problem: 5091
    User: danihao123
    Language: C++
    Result: Accepted
    Time:532 ms
    Memory:2384 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
typedef long long ll;
const ll ha = 1000000007LL;
ll pow_mod(ll a, ll b) {
  if(!b) return 1LL;
  ll ans = pow_mod(a, b >> 1);
  ans = (ans * ans) % ha;
  if(1LL & b) ans = (ans * a) % ha;
  return ans;
}
ll inv(ll v) {
  return pow_mod(v, ha - 2LL);
}
 
const int maxn = 100005;
ll a[maxn], d[maxn];
int main() {
  int n, m, k; scanf("%d%d%d", &n, &m, &k);
  for(int i = 1; i <= n; i ++) {
    scanf("%lld", &a[i]);
  }
  for(int i = 1; i <= m; i ++) {
    int u, v; scanf("%d%d", &u, &v);
    d[u] ++; d[v] ++;
  }
  ll ans = 0;
  ll inv_2m = inv(2 * m);
  for(int i = 1; i <= n; i ++) {
    ans += (((d[i] * inv_2m) % ha) * a[i]) % ha;
    if(ans >= ha) ans -= ha;
  }
  ans = (ans * (ll(k))) % ha;
  printf("%lld\n", ans);
  return 0;
}

[SPOJ]LCS2

和LCS那题很像,但是现在有十个串了……

还是考虑构建一个串的SAM。然后考虑SAM上每个点,剩余若干个串都可能会覆盖一下这个点,我们要取最短的那个。

然后我们就直接把其他几个串都在SAM上跑一遍(方法类似于LCS)就行了,并且注意每个点对自己的祖先也是有贡献的。

代码:

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
#include <queue>
#include <unordered_set>
const int maxn = 100005;
const int INF = 0x3f3f3f3f;
const int maxno = maxn << 1;
struct Node {
  Node *fa; int len; int ans[12];
  Node *ch[26];
  Node *lc, *rb;
};
Node pool[maxno];
Node *rt, *last;
Node *cnt;
void init_pool() {
  rt = last = cnt = pool;
  rt -> fa = NULL; rt -> len = 0;
  rt -> ans[0] = 0;
  // memset(rt -> ans, 0, sizeof(rt -> ans));
}
Node *alloc_node(int len = 0, Node *fa = NULL) {
  Node *ret = ++ cnt;
  ret -> len = len; ret -> fa = fa;
  // memset(ret -> ans, 0, sizeof(ret -> ans));
  ret -> ans[0] = len;
  return ret;
}

int idx(char c) {
  return c - 'a';
}
void insert(char cc) {
  int c = idx(cc);
  Node *np = alloc_node(last -> len + 1);
  Node *p = last;
  while(p != NULL && p -> ch[c] == NULL) {
    p -> ch[c] = np;
    p = p -> fa;
  }
  if(!p) {
    np -> fa = rt;
  } else {
    Node *q = p -> ch[c];
    if(q -> len == p -> len + 1) {
      np -> fa = q;
    } else {
      Node *nq = alloc_node(p -> len + 1, q -> fa);
      memcpy(nq -> ch, q -> ch, sizeof(q -> ch));
      q -> fa = np -> fa = nq;
      while(p != NULL && p -> ch[c] == q) {
        p -> ch[c] = nq;
        p = p -> fa;
      }
    }
  }
  last = np;
}
void dfs(int id, Node *u) {
  for(Node *c = u -> lc; c; c = c -> rb) {
    dfs(id, c);
    u -> ans[id] = std::max(u -> ans[id], c -> ans[id]);
  }
}
void run(int id, char *S) {
  int n = strlen(S);
  Node *u = rt;
  int len = 0;
  for(int i = 0; i < n; i ++) {
    int c = idx(S[i]);
    if(u -> ch[c]) {
      len ++; u = u -> ch[c];
    } else {
      Node *p = u;
      while(p != NULL && p -> ch[c] == NULL) p = p -> fa;
      if(p == NULL) {
        len = 0; u = rt;
      } else {
        len = p -> len + 1; u = p -> ch[c];
      }
    }
    u -> ans[id] = std::max(u -> ans[id], len);
  }
  dfs(id, rt);
}
int tot = 0;
int dfs_2(Node *u) {
  int ans = *std::min_element(u -> ans, u -> ans + tot + 1);
  for(Node *c = u -> lc; c; c = c -> rb) {
    ans = std::max(ans, dfs_2(c));
  }
  return ans;
}

int main() {
  static char buf[maxn];
  scanf("%s", buf);
  init_pool();
  int n = strlen(buf);
  for(int i = 0; i < n; i ++) {
    insert(buf[i]);
  }
  for(int i = 1; i <= (cnt - pool); i ++) {
    Node *u = pool + i;
    Node *f = u -> fa;
    u -> rb = f -> lc;
    f -> lc = u;
  }
  while(scanf("%s", buf) == 1) {
#ifdef LOCAL
    if(buf[0] == '-') break;
#endif
    run(++ tot, buf);
    // memset(buf, 0, sizeof(buf));
  }
  printf("%d\n", dfs_2(rt));
  return 0;
}

[LibreOJ 6024]XLkxc

拉格朗日插值的模板题qwq

考虑\(f(n) = \sum_{i = 1}^n i^k\),这显然是一个\(k + 1\)次多项式(通过差分之类的手段可以证明),然后我们发现\(g(n) = \sum_{i = 1}^n f(n)\)可以用类似手段证明是\(k + 2\)次多项式。类似的,我们会发现,答案是一个\(k + 3\)次多项式。

分别对\(g\)以及答案插值即可。

代码:

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <cctype>
#include <algorithm>
#include <utility>
#include <queue>
#include <climits>
typedef long long ll;
const ll ha = 1234567891LL;
int k;
ll a, n, d;
ll pf[205], pg[205], fac[205];
ll pow_mod(const ll &a, ll b) {
  if(!b) return 1LL;
  ll ans = pow_mod(a, b >> 1);
  ans = (ans * ans) % ha;
  if(b & 1LL) ans = (ans * a) % ha;
  return ans;
}
ll inv(ll v) {
  return pow_mod(v, ha - 2LL);
}

void process() {
  fac[0] = 1LL;
  for(int i = 1; i <= k + 5; i ++) {
    fac[i] = (fac[i - 1] * inv(i)) % ha;
#ifdef LOCAL
    printf("fac[%d] : %lld\n", i, fac[i]);
#endif
  }
  for(int i = 1; i <= k + 3; i ++) {
    pf[i] = (pf[i - 1] + pow_mod(i, k)) % ha;
#ifdef LOCAL
    printf("pf[%d] : %lld\n", i, pf[i]);
#endif
  }
  for(int i = 1; i <= k + 3; i ++) {
    pg[i] = (pg[i - 1] + pf[i]) % ha;
#ifdef LOCAL
    printf("pg[%d] : %lld\n", i, pg[i]);
#endif
  }
#ifdef LOCAL
  puts("Processing finished!");
#endif
}

ll g(ll x) {
  static ll sim_f[205], sim_g[205];
  sim_f[0] = 1LL;
  for(int i = 1; i <= k + 3; i ++) {
    sim_f[i] = (x - (ll(i)) + ha) % ha;
    sim_f[i] = (sim_f[i] * sim_f[i - 1]) % ha;
  }
  sim_g[k + 4] = 1LL;
  for(int i = k + 3; i >= 1; i --) {
    sim_g[i] = (x - (ll(i)) + ha) % ha;
    sim_g[i] = (sim_g[i] * sim_g[i + 1]) % ha;
  }
  ll ret = 0;
  for(int i = 1; i <= k + 3; i ++) {
    ll ans = (((pg[i] * sim_f[i - 1]) % ha) * sim_g[i + 1]) % ha;
    ans = (ans * fac[i - 1]) % ha;
    ans = (ans * fac[k + 3 - i]) % ha;
    if(1 & (k + 3 - i)) {
      ret = (ret - ans + ha) % ha;
    } else {
      ret = (ret + ans) % ha;
    }
  }
#ifdef LOCAL
  printf("g(%lld) : %lld\n", x, ret);
#endif
  return ret;
}
ll h(ll x) {
  static ll ph[205];
  ph[0] = g(a);
  for(int i = 1; i <= k + 4; i ++) {
    ph[i] = (ph[i - 1] + g(a + (ll(i)) * d)) % ha;
  }
  static ll sim_f[205], sim_g[205];
  sim_f[0] = 1LL;
  for(int i = 1; i <= k + 4; i ++) {
    sim_f[i] = (x - (ll(i)) + ha) % ha;
    sim_f[i] = (sim_f[i] * sim_f[i - 1]) % ha;
  }
  sim_g[k + 5] = 1LL;
  for(int i = k + 4; i >= 1; i --) {
    sim_g[i] = (x - (ll(i)) + ha) % ha;
    sim_g[i] = (sim_g[i] * sim_g[i + 1]) % ha;
  }
  ll ret = 0;
  for(int i = 1; i <= k + 4; i ++) {
    ll ans = (((ph[i] * sim_f[i - 1]) % ha) * sim_g[i + 1]) % ha;
    ans = (ans * fac[i - 1]) % ha;
    ans = (ans * fac[k + 4 - i]) % ha;
    if(1 & (k + 4 - i)) {
      ret = (ret - ans + ha) % ha;
    } else {
      ret = (ret + ans) % ha;
    }
  }
  return ret;
}

int main() {
  int T; scanf("%d", &T);
  while(T --) {
    scanf("%d%lld%lld%lld", &k, &a, &n, &d);
    process();
    printf("%lld\n", h(n));
  }
  return 0;
}

[LibreOJ 2249][NOI2014]购票

在紧张刺激的等待之后终于肝掉了这道题……

本题的DP方程长成这样(其中\(a\)指\(v\)的某个满足距离限制的祖先,\(d_v\)指\(v\)到根的路径长):

\[f_v = min(f_a + p_v(d_v - d_a) + q_v)\]

化简之后发现:

\[f_v = q_v + p_v d_v + min(f_a - p_v d_a)\]

利用\(min\)中那一块很容易发现是截距式……但是问题在于,我们的转移来源是树上的连续一段祖先,怎样维护他们的凸包?

答案很狂暴啊……用树链剖分套上向量集那题的线段树套凸包,然后完了……

(注意一点细节:本题因为数据范围过大,故可能存在两个向量叉乘爆long long,所以在求凸包时如果直接用叉积判断是否需要删点会炸掉,建议用斜率判断)

代码:

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
#include <vector>
#include <queue>
#include <deque>
#include <cmath>
#include <set>
#include <climits>
using ll = long long;
using T = ll;
using R = long double;
const R eps = 1e-8;
int sign(R x) {
  if(fabsl(x) < eps) {
    return 0;
  } else {
    if(x > (R(0.00))) {
      return 1;
    } else {
      return -1;
    }
  }
}

struct Point {
  T x, y;
  Point(T qx = 0LL, T qy = 0LL) {
    x = qx; y = qy;
  }
};
using Vector = Point;
Vector operator +(const Vector &a, const Vector &b) {
  return Vector(a.x + b.x, a.y + b.y);
}
Vector operator -(const Point &a, const Point &b) {
  return Vector(a.x - b.x, a.y - b.y);
}
Vector operator *(const Vector &a, T lam) {
  return Vector(a.x * lam, a.y * lam);
}
Vector operator *(T lam, const Vector &a) {
  return Vector(a.x * lam, a.y * lam);
}
inline T dot(const Vector &a, const Vector &b) {
  return (a.x * b.x + a.y * b.y);
}
inline T times(const Vector &a, const Vector &b) {
  return (a.x * b.y - a.y * b.x);
}
inline bool cmp(const Point &a, const Point &b) {
  if(a.x == b.x) {
    return a.y < b.y;
  } else {
    return a.x < b.x;
  }
}
inline R slope(const Vector &a) {
  R dx = a.x, dy = a.y;
  return (dy / dx);
}

inline void andrew(Point *P, int L, std::vector<Point> &bot, std::vector<Point> &top) {
  std::sort(P + 1, P + 1 + L, cmp);
  for(int i = 1; i <= L; i ++) {
    if(i != 1 && (P[i].x == P[i - 1].x)) continue;
    while(bot.size() > 1 && sign(slope(P[i] - bot.back()) - slope(bot.back() - bot[bot.size() - 2])) <= 0) {
      bot.pop_back();
    }
    bot.push_back(P[i]);
  }
  for(int i = L; i >= 1; i --) {
    if(i != L && (P[i].x == P[i + 1].x)) continue;
    while(top.size() > 1 && sign(slope(P[i] - top.back()) - slope(top.back() - top[top.size() - 2])) <= 0) {
      top.pop_back();
    }
    top.push_back(P[i]);
  }
  std::reverse(top.begin(), top.end());
}

const int maxn = 200005;
const int maxno = maxn << 2;
const int N = 200000;
bool zen[maxno];
std::vector<Point> bot[maxno], top[maxno];
Point P[maxn];
inline void maintain(int o, int L, int R) {
  static Point tmp[maxn];
  const int lc = o << 1, rc = o << 1 | 1;
  const bool used = zen[o];
  zen[o] = (zen[lc] && zen[rc]);
  if(zen[o] != used) {
    std::copy(P + L, P + R + 1, tmp + 1);
    int len = R - L + 1;
    andrew(tmp, len, bot[o], top[o]);
  }
}
void modify(int o, int L, int R, const int &p, const Point &v) {
  if(L == R) {
    zen[o] = true;
    P[L] = v;
    bot[o].push_back(v); top[o].push_back(v);
  } else {
    const int M = (L + R) / 2;
    if(p <= M) {
      modify(o << 1, L, M, p, v);
    } else {
      modify(o << 1 | 1, M + 1, R, p, v);
    }
    maintain(o, L, R);
  }
}
inline T calc_ans(T k, const Point &v) {
  return v.y - k * v.x;
}
inline T search(const std::vector<Point> &vec, const T &k) {
  int l = 0, r = vec.size() - 1;
  while(r - l > 2) {
    int lm = (l * 2 + r) / 3, rm = (2 * r + l) / 3;
    if((calc_ans(k, vec[lm]) > calc_ans(k, vec[rm]))) {
      l = lm;
    } else {
      r = rm;
    }
  }
  T ans = LLONG_MAX;
  for(int i = l; i <= r; i ++) {
    ans = std::min(ans, calc_ans(k, vec[i]));
  }
  return ans;
}
T query(int o, int L, int R, const int &ql, const int &qr, const T &k) {
  if(ql <= L && R <= qr) {
    return search(bot[o], k);
  } else {
    int M = (L + R) / 2;
    T ans = LLONG_MAX;
    if(ql <= M) {
      ans = std::min(ans, query(o << 1, L, M, ql, qr, k));
    }
    if(qr > M) {
      ans = std::min(ans, query(o << 1 | 1, M + 1, R, ql, qr, k));
    }
    return ans;
  }
}

int first[maxn];
int next[maxn << 1], to[maxn << 1];
ll dist[maxn << 1];
void add_edge(int u, int v, ll d) {
  static int cnt = 0;
  cnt ++;
  next[cnt] = first[u];
  first[u] = cnt;
  to[cnt] = v;
  dist[cnt] = d;
}

int fa[maxn], dep[maxn], hson[maxn];
ll d[maxn];
int siz[maxn];
int bs[maxn][18];
void dfs_1(int x, int f = -1, int depth = 1) {
  fa[x] = bs[x][0] = f; dep[x] = depth;
  siz[x] = 1;
  int max_siz = 0;
  for(int i = first[x]; i; i = next[i]) {
    int v = to[i];
    if(v != f) {
      d[v] = d[x] + dist[i];
      dfs_1(v, x, depth + 1);
      siz[x] += siz[v];
      if(siz[v] > max_siz) {
        hson[x] = v; max_siz = siz[v];
      }
    }
  }
}
int dfn[maxn], tid[maxn], up[maxn];
void dfs_2(int x, int a = 1, int f = 0) {
  static int cnt = 0;
  dfn[x] = ++ cnt; tid[cnt] = x;
  up[x] = a;
  if(hson[x]) {
    dfs_2(hson[x], a, x);
  } else {
    return;
  }
  for(int i = first[x]; i; i = next[i]) {
    int v = to[i];
    if(v != f && v != hson[x]) {
      dfs_2(v, v, x);
    }
  }
}
int k_anc(int x, ll k) {
  int yx = x;
  for(int j = 17; j >= 0; j --) {
    int a = bs[x][j];
    if(a != -1 && d[yx] - d[a] <= k) {
      x = a;
    }
  }
#ifdef LOCAL
  printf("%d's %lld-th anc : %d\n", yx, k, x);
#endif
  return x;
}
int n;
ll get_up(int x, int anc, ll k) {
  ll ans = LLONG_MAX;
  while(up[x] != up[anc]) {
    ans = std::min(ans, query(1, 1, n, dfn[up[x]], dfn[x], k));
    x = fa[up[x]];
  }
  return std::min(ans, query(1, 1, n, dfn[anc], dfn[x], k));
}

ll p[maxn], q[maxn], l[maxn];
ll f[maxn];
void dp(int x) {
#ifdef LOCAL
  printf("processing %d...\n", x);
  printf("d : %lld\n", d[x]);
#endif
  if(x != 1) {
#ifdef LOCAL
    printf("b : %lld\n", get_up(fa[x], k_anc(x, l[x]), p[x]));
#endif
    f[x] = get_up(fa[x], k_anc(x, l[x]), p[x]) + d[x] * p[x] + q[x];
  } else {
    f[x] = 0;
  }
#ifdef LOCAL
  printf("ans : %lld\n", f[x]);
#endif
  modify(1, 1, n, dfn[x], Point(d[x], f[x]));
  for(int i = first[x]; i; i = next[i]) {
    int v = to[i];
    dp(v);
  }
}

int main() {
  int t; scanf("%d%d", &n, &t);
  for(int i = 2; i <= n; i ++) {
    int father; T s;
    scanf("%d%lld%lld%lld%lld", &father, &s, &p[i], &q[i], &l[i]);
    add_edge(father, i, s);
  }
  memset(bs, -1, sizeof(bs));
  dfs_1(1); dfs_2(1);
  for(int j = 1; (1 << j) < n; j ++) {
    for(int i = 1; i <= n; i ++) {
      int a = bs[i][j - 1];
      if(a != -1) {
        bs[i][j] = bs[a][j - 1];
      }
    }
  }
  dp(1);
  for(int i = 2; i <= n; i ++) {
    printf("%lld\n", f[i]);
  }
  return 0;
}

[LibreOJ 2353][NOI2007]货币兑换

emmm做了一下这道神题……(我可能是少有的用动态凸包苟的?

首先DP方程长这样:

\[f_i = max(f_{i - 1}, f_j\cdot\frac{A_iR_j+B_i}{A_jR_j+B_j})\]

然后这个方程炒鸡复杂……首先\(f_{i - 1}\)不要管了,然后设\(a_i = \frac{f_i}{A_iR_i + B_i}\)。在xjb推了一番之后我们终于得到了截距式……

\[-a_j R_j \frac{A_i}{B_i} + \frac{f_i}{B_i} = a_j\]

但是这玩意太毒瘤了……斜率不可能单调的,这还好,在凸壳上二分/三分一下即可。但问题在于,横坐标也不单调……

这个时候就需要动态维护凸包了(其实是我不会CDQ),我直接把我向量集那题的二进制分组线段树搬了过来……(逃

代码:

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
#include <vector>
#include <cmath>
#include <climits>
#include <deque>
#include <cassert>
using R = double;
const R eps = 1e-8;
int sign(R x) {
  if(fabs(x) < eps) {
    return 0;
  } else {
    if(x > 0.00) {
      return 1;
    } else {
      return -1;
    }
  }
}

struct Point {
  R x, y;
  Point(R qx = 0, R qy = 0) {
    x = qx; y = qy;
  }
};
using Vector = Point;
Vector operator +(const Vector &a, const Vector &b) {
  return Vector(a.x + b.x, a.y + b.y);
}
Vector operator -(const Point &a, const Point &b) {
  return Vector(b.x - a.x, b.y - a.y);
}
Vector operator *(const Vector &a, R lam) {
  return Vector(a.x * lam, a.y * lam);
}
Vector operator *(R lam, const Vector &a) {
  return Vector(a.x * lam, a.y * lam);
}
inline R dot(const Vector &a, const Vector &b) {
  return (a.x * b.x + a.y * b.y);
}
inline R times(const Vector &a, const Vector &b) {
  return (a.x * b.y - a.y * b.x);
}
inline bool cmp(const Point &a, const Point &b) {
  if(sign(a.x - b.x) == 0) {
    return a.y < b.y;
  } else {
    return a.x < b.x;
  }
}
inline void andrew(Point *P, int L, std::vector<Point> &bot, std::vector<Point> &top) {
  std::sort(P + 1, P + 1 + L, cmp);
  for(int i = 1; i <= L; i ++) {
    if(i != 1 && sign(P[i].x - P[i - 1].x) == 0) continue;
    while(bot.size() > 1 && sign(times(P[i] - bot.back(), bot.back() - bot[bot.size() - 2])) >= 0) {
      bot.pop_back();
    }
    bot.push_back(P[i]);
  }
  for(int i = L; i >= 1; i --) {
    if(i != L && sign(P[i].x - P[i + 1].x) == 0) continue;
    while(top.size() > 1 && sign(times(P[i] - top.back(), top.back() - top[top.size() - 2])) >= 0) {
      top.pop_back();
    }
    top.push_back(P[i]);
  }
  std::reverse(top.begin(), top.end());
}

const int maxn = 1000005;
const int N = 1000000;
const int maxno = maxn << 2;
bool zen[maxno];
std::vector<Point> bot[maxno], top[maxno];
Point P[maxn];
inline void maintain(int o, int L, int R) {
  static Point tmp[maxn];
  const int lc = o << 1, rc = o << 1 | 1;
  const bool used = zen[o];
  zen[o] = (zen[lc] && zen[rc]);
  if(zen[o] != used) {
    std::copy(P + L, P + R + 1, tmp + 1);
    int len = R - L + 1;
    andrew(tmp, len, bot[o], top[o]);
  }
}
void modify(int o, int L, int R, const int &p, const Point &v) {
  if(L == R) {
    zen[o] = true;
    P[L] = v;
    bot[o].push_back(v); top[o].push_back(v);
  } else {
    const int M = (L + R) / 2;
    if(p <= M) {
      modify(o << 1, L, M, p, v);
    } else {
      modify(o << 1 | 1, M + 1, R, p, v);
    }
    maintain(o, L, R);
  }
}
inline R calc_ans(R k, const Point &v) {
  return v.y - k * v.x;
}
inline R search(const std::vector<Point> &vec, const R &k) {
  int l = 0, r = vec.size() - 1;
  while(r - l > 2) {
    int lm = (l * 2 + r) / 3, rm = (2 * r + l) / 3;
    if(sign(calc_ans(k, vec[lm]) - calc_ans(k, vec[rm])) == 1) {
      r = rm;
    } else {
      l = lm;
    }
  }
  R ans = INT_MIN;
  for(int i = l; i <= r; i ++) {
    ans = std::max(ans, calc_ans(k, vec[i]));
  }
  return ans;
}
R query(int o, int L, int R, const int &ql, const int &qr, const double &k) {
  if(ql <= L && R <= qr) {
    return search(top[o], k);
  } else {
    int M = (L + R) / 2;
    double ans = INT_MIN;
    if(ql <= M) {
      ans = std::max(ans, query(o << 1, L, M, ql, qr, k));
    }
    if(qr > M) {
      ans = std::max(ans, query(o << 1 | 1, M + 1, R, ql, qr, k));
    }
    return ans;
  }
}

int n, s;
R A[maxn], B[maxn], Rate[maxn];
R f[maxn];
R dp() {
  static double a[maxn];
  f[0] = s; f[1] = s; a[1] = f[1] / (A[1] * Rate[1] + B[1]);
  modify(1, 1, n, 1, Point(a[1] * Rate[1], a[1]));
  for(int i = 2; i <= n; i ++) {
    f[i] = query(1, 1, n, 1, i - 1, -A[i] / B[i]) * B[i];
    f[i] = std::max(f[i], f[i - 1]);
    a[i] = f[i] / (A[i] * Rate[i] + B[i]);
    if(i < n) modify(1, 1, n, i, Point(a[i] * Rate[i], a[i]));
  }
  return f[n];
}

int main() {
  scanf("%d%d", &n, &s);
  for(int i = 1; i <= n; i ++) {
    scanf("%lf%lf%lf", &A[i], &B[i], &Rate[i]);
  }
  printf("%.3lf\n", dp());
  return 0;
}

[LibreOJ 2035][SDOI2016]征途

又做了一道适合我这种沙茶做的题qwq

考虑方差,它满足这个式子:

\[Var(X) = E[X^2] - (E[X])^2\]

对于这个题展开,发现后半部分是一个常量(\((\frac{s_n}{m})^2\))。最小化的就是前半部分,前半部分相当于要求你求一个序列划分,使得各部分平方和的平均值最小。这个值乘上\(m^2\)之后就变成了各部分平方和乘上\(m\)。

然后就很简单了……DP方程化简之后是这样的:

\[f_{t, i} = ms_i^2 + max(f_{t - 1, j} + ms_j^2 - 2ms_i s_j)\]

求截距式形式,得:

\[2ms_i s_j + d_i = f_j + ms_j^2\]

然后用类似于上上题的方法搞就行。还有我想想啥时候补一下斜率优化的tutorial……

代码:

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
#include <deque>
#include <cmath>
#include <climits>
typedef long long ll;
typedef ll T;
struct Point {
  T x, y;
  Point(T qx = 0LL, T qy = 0LL) {
    x = qx; y = qy;
  }
};
typedef Point Vector;
Vector operator +(const Vector &a, const Vector &b) {
  return Vector(a.x + b.x, a.y + b.y);
}
Vector operator -(const Point &a, const Point &b) {
  return Vector(a.x - b.x, a.y - b.y);
}
Vector operator *(const Vector &a, T lam) {
  return Vector(a.x * lam, a.y * lam);
}
Vector operator *(T lam, const Vector &a) {
  return Vector(a.x * lam, a.y * lam);
}
inline T dot(const Vector &a, const Vector &b) {
  return (a.x * b.x + a.y * b.y);
}
inline T times(const Vector &a, const Vector &b) {
  return (a.x * b.y - a.y * b.x);
}
 
const int maxn = 3005;
T f[maxn][maxn];
T S[maxn];
int n; T m;
void dp() {
  for(int t = 1; t <= m; t ++) {
    std::deque<Point> Q;
    Q.push_back(Point(0LL, 0LL));
    for(int i = 1; i <= n; i ++) {
      ll k = 2 * m * S[i];
      Vector st(1LL, k);
      while(Q.size() > 1 && times(Q[1] - Q[0], st) > 0LL) {
        Q.pop_front();
      }
      f[t][i] = m * S[i] * S[i];
      f[t][i] += Q.front().y - Q.front().x * k;
#ifdef LOCAL
      printf("f[%d][%d] : %lld\n", t, i, f[t][i]);
#endif
      if(t == 1) continue;
      Vector ins(S[i], f[t - 1][i] + m * S[i] * S[i]);
      while(Q.size() > 1 && times(ins - Q.back(), Q.back() - Q[Q.size() - 2]) > 0LL) {
        Q.pop_back();
      }
      Q.push_back(ins);
    }
  }
}

int main() {
  scanf("%d%lld", &n, &m);
  for(int i = 1; i <= n; i ++) {
    scanf("%lld", &S[i]);
  }
  for(int i = 1; i <= n; i ++) {
    S[i] += S[i - 1];
  }
  dp();
  printf("%lld\n", f[m][n] - S[n] * S[n]);
  return 0;
}