bzoj3601 - 一个人的数论（莫比乌斯函数，伯努利数）

题目

给定n（以唯一分解的形式给出）和d，求[1, n]中与n互质的数的d次方和。
结果模1e9+7。

题解

由于n非常大，是以质数幂次方给出。因此应该是要整出一个积性函数，然后拆开求解。

推式子......

\[\sum\limits_{i=1}^n{i^d[\gcd(i, n)=1]} \]

\[\sum\limits_{i=1}^n{i^d\sum\limits_{d_1|\gcd(i, j)}{\mu(d_1)}} \]

\[\sum\limits_{d_1|n}{\mu(d_1)\sum\limits_{i=1}^n{i^d[d_1|i]}} \]

\[\sum\limits_{d_1|n}{\mu(d_1)d_1^d\sum\limits_{i=1}^{\frac{n}{d_1}}{i^d}} \]

化到这里发现到头了，但\(\sum\limits_{i=1}^{\frac{n}{d_1}}{i^d}\)不是积性函数。故还要进行进一步的转化。

\(\sum\limits_{i=1}^{\frac{n}{d_1}}{i^d}\)是d次幂和，可以化为一个关于\(\frac{n}{d_1}\)的多项式。即

\[\sum\limits_{d_1|n}{\mu(d_1)d_1^d\sum\limits_{i=1}^{k}{a_i(\frac{n}{d_1})^i}} \]

\[\sum\limits_{i=1}^{k}{a_i\sum\limits_{d_1|n}{\mu(d_1)d_1^d(\frac{n}{d_1})^i}} \]

现在，\(\sum\limits_{d_1|n}{\mu(d_1)d_1^d(\frac{n}{d_1})^i}\)就是积性函数了，可以愉快地拆开求解。

至于多项式系数\(a_i\)怎么求，见伯努利数，当然也可以用高斯消元或者拉格朗日插值法。

#include <bits/stdc++.h>

#define endl '\n'
#define IOS std::ios::sync_with_stdio(0); cin.tie(0); cout.tie(0)
#define mp make_pair
#define seteps(N) fixed << setprecision(N) 
typedef long long ll;

using namespace std;
/*-----------------------------------------------------------------*/

ll gcd(ll a, ll b) {return b ? gcd(b, a % b) : a;}
#define INF 0x3f3f3f3f

const int N = 1000 + 10;
const int M = 1e9 + 7;
const double eps = 1e-5;

ll B[N];
ll C[N][N];
ll inv[N];

void init() {
    // 预处理组合数
    for (int i = 0; i < N; i++) {
      C[i][0] = C[i][i] = 1;
      for (int k = 1; k < i; k++) {
        C[i][k] = (C[i - 1][k] % M + C[i - 1][k - 1] % M) % M;
      }
    }
    // 预处理逆元
    inv[1] = 1;
    for (int i = 2; i < N; i++) {
      inv[i] = (M - M / i) * inv[M % i] % M;
    }
    // 预处理伯努利数
    B[0] = 1;
    for (int i = 1; i < N; i++) {
        ll ans = 0;
        if (i == N - 1) break;
        for (int k = 0; k < i; k++) {
          ans += C[i + 1][k] * B[k];
          ans %= M;
        }
        ans = (ans * (-inv[i + 1]) % M + M) % M;
        B[i] = ans;
    }
}


ll aval(int n, int p) {
    int k = n + 1 - p;
    return inv[n + 1] * C[n + 1][k] % M * B[k] % M;
}

int p[N], a[N];
int rp[N];

inline ll qpow(ll a, ll b, ll m) {
    ll res = 1;
    while(b) {
        if(b & 1) res = (res * a) % m;
        a = (a * a) % m;
        b = b >> 1;
    }
    return res;
}

int main() {
    IOS;
    init();
    int d, w;
    cin >> d >> w;
    ll n = 1;
    for(int i = 1; i <= w; i++) {
        cin >> p[i] >> a[i];
        n = n * qpow(p[i], a[i], M) % M;
        rp[i] = qpow(p[i], M - 2, M);
    }
    ll ans = 0;
    ll pn = n;
    for(int i = 1; i <= d + 1; i++) {
        ll f = 1;
        for(int j = 1; j <= w; j++) {                      
            f = f * (1 - qpow(p[j], d, M) % M * qpow(rp[j], i, M) % M) % M;
        }
        ans += pn * aval(d, i) % M * f % M;
        ans %= M;
        pn = pn * n % M;
    }
    cout << (ans + M) % M << endl;
    
}