[BZOJ 3601]一个人的数论

danihao123 posted @ 2018年5月29日 14:35 in 题解 with tags BZOJ 伯努利数 莫比乌斯反演 狄利克雷卷积 , 96 阅读
转载请注明出处:http://danihao123.is-programmer.com/

本来想做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;
}

登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter