BZOJ3601 一个人的数论

Description

定义

[f_k(n)=sum_{substack{1leq ileq n\gcd(i,n)=1}}i^k ]

给出(n=prod_{i=1}^w p_i^{a_i})，求(f_k(n))。(1leq wleq 1000, 1leq q_i,a_ileq 10^9)。保证(p_i)都为质数且互不相同。

Solution

令(g_k(n)=sum_{i=1}^ni^k)，则

[egin{aligned} g_k(n)&=sum_{i=1}^ni^k\ &=sum_{d|n}sum_{substack{1leq ileq n\gcd(i,n)=d}}i^k\ &=sum_{d|n}d^ksum_{substack{1leq ileq frac nd\gcdleft(i,frac nd ight)=d}}i^k\ &=sum_{d|n}d^kf_kleft(frac nd ight) end{aligned} ]

也即(g_k(n) = (n^k) * f_k(n))（(*)表示狄利克雷卷积）。

由于((n^k)*(mu(n)n^k)=[n=1])，所以(f_k(n)=(mu(n)n^k)*g_k(n))。

显然(g_k(n))可以写成(sum_{i=1}^{k+1}a_in^i)，那么

[egin{aligned} f_k(n)&=sum_{d|n}mu(d)d^ksum_{i=1}^{k+1}a_ileft(frac nd ight)^i\ &=sum_{i=1}^{k+1}a_isum_{d|n}mu(d)d^kleft(frac nd ight)^i\ &=sum_{i=1}^{k+1}a_in^isum_{d|n}mu(d)d^{k-i}\ end{aligned} ]

(sum_{d|n}mu(d)d^{k-i})显然是积性函数，所以对每个质因子(O(1))求出之后乘起来即可。

所有(a_i)提前高消出来就行了。

Code

#include <algorithm>
#include <cstdio>
#include <cstring>
const int D = 105;
const int mod = 1000000007;
typedef long long LL;
inline LL pow_mod(LL a, LL b) {
  LL ans = 1;
  for (a %= mod, (b += mod - 1) %= mod - 1; b; b >>= 1, a = a * a % mod)
    if (b & 1) ans = ans * a % mod;
  return ans;
}
inline LL inv(LL a) { return pow_mod(a, -1); }
int d;
LL a[D];
void solve() {
  static LL A[D][D];
  LL t = 0;
  for (int i = 0; i <= d; ++i) {
    LL j = 1;
    for (int k = 0; k <= d; ++k) A[i][k] = j = j * (i + 1) % mod;
    A[i][d + 1] = ((t += d ? A[i][d - 1] : 1) %= mod);
  }
  for (int i = 0; i <= d; ++i) {
    int j = i;
    while (!A[j][i]) ++j;
    for (int k = i; k <= d + 1; ++k) std::swap(A[i][k], A[j][k]);
    LL inv1 = inv(A[i][i]);
    for (int j = i; j <= d + 1; ++j)
      A[i][j] = A[i][j] * inv1 % mod;
    for (int j = i + 1; j <= d; ++j)
      for (int k = d + 1; k >= i; --k)
        A[j][k] = (A[j][k] - A[j][i] * A[i][k] % mod) % mod;
  }
  for (int i = d; ~i; --i) {
    a[i + 1] = A[i][d + 1];
    for (int j = i - 1; ~j; --j)
      A[j][d + 1] = (A[j][d + 1] - A[j][i] * A[i][d + 1] % mod) % mod;
  }
}
const int N = 1050;
int p[N], q[N];
int main() {
  int w;
  scanf("%d%d", &d, &w);
  solve();
  LL ans = 0, n = 1;
  for (int i = 0; i < w; ++i) {
    scanf("%d%d", &p[i], &q[i]);
    n = n * pow_mod(p[i], q[i]) % mod;
  }
  LL y = 1;
  for (int i = 1; i <= d + 1; ++i) {
    y = y * n % mod;
    LL t = a[i] * y % mod;
    for (int j = 0; j < w; ++j)
      t = t * (1 - pow_mod(p[j], d - i)) % mod;
    ans = (ans + t) % mod;
  }
  printf("%d
", (int)((ans + mod) % mod));
  return 0;
}