luogu P6271 [湖北省队互测2014]一个人的数论

https://www.luogu.com.cn/problem/P6271

设\(G(n)=\sum\limits_{i=1}^n [i \perp n] i ^k\)

按照套路显然有

\(G(n)=\sum\limits_{i=1}^n \sum\limits_{d|n,d|i}\mu(d) i^k\)
\(= \sum\limits_{d|n}\mu(d) \sum\limits_{i=1}^{n/d}(id)^k\)

\(= \sum\limits_{d|n}\mu(d) d^k \sum\limits_{i=1}^{n/d}i^k\)

后面那个就是个自然数幂和，设\(F(n)=\sum\limits_{i=1}^ni^k\)
\(= \sum\limits_{d|n}\mu(d) d^k F(n/d)\)

根据常识，我们知道，\(F(n)\)为一个\(k+1\)次的多项式
\(= \sum\limits_{d|n}\mu(d) d^k \sum\limits_{i=0}^{k+1}f_i(\frac{n}{d})^i\)

交换求和顺序
\(= \sum\limits_{i=0}^{k+1}f_i n^i \sum\limits_{d|n}\mu(d) d^{k-i}\)

我们来看看后面那个东西怎么计算

可以发现显然是个积性函数，所以有
\(\sum\limits_{d|n}\mu(d) d^{k-i}=\prod\limits_{p|n}\sum\limits_{j=0}^{a[p]}\mu(p^j) p^{j(k-i)}\)

然而我们发现当\(j>1\)的时候\(\mu(p_j)=0\)
所以这条式子实际上是
\(\prod\limits_{p|n}(1-p^{k-i})\)

带回上式可以得到

\(= \sum\limits_{i=0}^{k+1}f_i n^i\prod\limits_{p|n}(1-p^{k-i})\)

\(f_i\)直接用拉格朗日插值大力吧多项式搞出来即可

后面的可以直接计算
时间复杂度\(O(k^3+nk)\)

注意\(k+1\)次多项式要带\(k+2\)个值进去
code:

#include<bits/stdc++.h>
#define mod 1000000007
#define ll long long
#define N 3050
using namespace std;
ll qpow(ll x, ll y) { 
    if(y < 0) y += mod - 1;
    ll ret = 1;
    for(; y; y >>= 1, x = x * x % mod) if(y & 1) ret = ret * x % mod;
    return ret;
}

int len, m, k;
ll g[N], f[N], x[N], y[N], p[N];
void MUL(int a, int b) {
    len ++;
    for(int i = len; i >= 1; i --) {
        g[i] = (g[i - 1] * a % mod + g[i] * b % mod) % mod;
    }
    g[0] = g[0] * b % mod;
}
void MULv(int a) {
    for(int i = 0; i <= len; i ++) g[i] = g[i] * a % mod;
}
void ADD() {
    for(int i = 0; i <= len; i ++) f[i] = (f[i] + g[i]) % mod, g[i] = 0;
}

void cha(int n) {
    for(int i = 1; i <= n; i ++) x[i] = i;
    for(int i = 1; i <= n; i ++) y[i] = (y[i - 1] + qpow(x[i], k)) % mod;

    // for(int i = 1; i <= n; i ++) printf("%lld ", x[i]); printf("\n");
    // for(int i = 1; i <= n; i ++) printf("%lld ", y[i]); printf("\n");

    for(int i = 1; i <= n; i ++) {
        len = 0; g[0] = 1;
        ll M = y[i];
        for(int j = 1; j <= n; j ++) {
            if(i == j) continue;
            MUL(1, mod - x[j]);
            M = M * qpow((x[i] - x[j] + mod) % mod, mod - 2) % mod;
        } MULv(M);
        ADD();
    }    
}
int main() {
    scanf("%d%d", &k, &m); 
    ll n = 1;
    for(int i = 1; i <= m; i ++) {
        int mi;
        scanf("%lld%d", &p[i], &mi);
        n = n * qpow(p[i], mi) % mod;
    }

    cha(k + 2);

    //for(int i = 0; i <= k + 2; i ++) printf("%lld ", f[i]); printf("\n");

//    int ss = 0;
//     for(int i = 0; i <= k + 2; i ++) ss = (ss + f[i] * qpow(8, i) % mod) % mod;
//     printf("%d\n", ss);
//     printf("%lld \n", n);
    ll ans = 0;
    for(int i = 0; i <= k + 2; i ++) {
        int ret = f[i] * qpow(n, i) % mod;
        for(int j = 1; j <= m; j ++) {
            ret = ret *  (1 - qpow(p[j], k - i) + mod) % mod;
        }
        ans = (ans + ret) % mod;
    //    printf("%lld\n", ans);
    }
    printf("%lld", ans);
    return 0;
}