【bzoj3601】一个人的数论（莫比乌斯反演+拉格朗日插值）

传送门

题意：
求

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

思路：
我们对上面的式子进行变换，有：

\[\begin{aligned} &\sum_{i=1}^{n}i[gcd(i,n)=1]\\ =&\sum_{i=1}^{n}i\sum_{x|gcd(i,n)}\mu (x)\\ =&\sum_{i=1}^n i\sum_{x|i,x|n}\mu(x)\\ =&\sum_{x|n}\mu(x)x^d\sum_{i=1}^{\frac{n}{x}}i^d \end{aligned} \]

以上都是一些套路，接下来才步入正题。
因为形如\(\sum_{i=1}^n i^d\)这种都是一个以\(n\)为自变量的，最高项次数为\(d+1\)的多项式。

到这步后，我们将后面的形式化，即将多项式表示出来，设\(a_i\)为相关系数，那么就有式子等于：

\[\sum_{x|n}\mu(x)x^d\sum_{i=0}^{d+1}a_i\lfloor\frac{n}{x}\rfloor^i \]

变化一下有：

\[\sum_{i=0}^{d+1}a_i\sum_{x|n}\mu(x)x^d\lfloor\frac{n}{x}\rfloor^i \]

这里面前面的\(a_i\)为多项式的系数，是未知的。
但其实因为我们知道多项式的形式为\(\sum_{i=0}^{d+1}a_ix^i\)，我们可以直接把\(x=1,x=2,\cdots,x=d+1\)的值求出来，然后高斯消元求解系数。
这里也可以利用拉格朗日插值来求解，代码是直接抄yyb(orz)的，思路应该是利用一下等式来求系数：

\[\sum_{i=0}^n y_i\prod_{i!=j}\frac{x-x_j}{x_i-x_j} \]

那么现在主要就是后面一部分的计算，我们令\(f(x)=\mu(x)x^d,g(x)=x^i\)，那么后面一部分可以写为：\(\sum_{x|n}f(x)g(\frac{n}{x})\)，这是狄利克雷卷积的的形式，因为\(f,g\)都为积性，那么\(h=f*g(n)\)也为积性。
所以\(h(n)=\sum_{x|n}\mu(x)x^d\lfloor\frac{n}{x}\rfloor^i\)也为积性函数，那么我们可以考虑单独素因子的贡献。显然每个素因子只会出现\(0\)次或者\(1\)次，否则贡献为\(0\)，那么有：

\[\begin{aligned} h(p^a)&=\sum_{j=0}^{a}\mu(p^j)p^{jd}(\frac{p^a}{p^j})^i\\ &=p^{ai}-p^{ai}p^{d-i} \end{aligned} \]

那么对于每个\(i,0\leq i\leq d+1\)，求出相应的\(h(n)\)，再与系数相乘最终结果就出来了。

感觉解法中将多项式形式化出来的想法很巧妙！没想到多项式还能这么用hhh，直接求解多项式的系数也是之前没想到的。之后对卷积的观察也很重要。
很好的一个题。细节参考代码：

/*
 * Author:  heyuhhh
 * Created Time:  2019/11/21 19:44:08
 */
#include <bits/stdc++.h>
#define MP make_pair
#define fi first
#define se second
#define sz(x) (int)(x).size()
#define all(x) (x).begin(), (x).end()
#define INF 0x3f3f3f3f
#define Local
#ifdef Local
  #define dbg(args...) do { cout << #args << " -> "; err(args); } while (0)
  void err() { std::cout << '\n'; }
  template<typename T, typename...Args>
  void err(T a, Args...args) { std::cout << a << ' '; err(args...); }
#else
  #define dbg(...)
#endif
void pt() {std::cout << '\n'; }
template<typename T, typename...Args>
void pt(T a, Args...args) {std::cout << a << ' '; pt(args...); }
using namespace std;
typedef long long ll;
typedef pair<int, int> pii;
//head
const int N = 1000 + 5, MOD = 1e9 + 7;

ll qpow(ll a, ll b) {
    ll ans = 1;
    while(b) {
        if(b & 1) ans = ans * a % MOD;
        a = a * a % MOD;
        b >>= 1;   
    }
    return ans;   
}
//求d次多项式系数
struct Lagrange {
    ll f[N], a[N], b[N];
    int d;
    void init(int _d) {
        d = _d;
        //y_i
        for(int i = 1; i <= d + 1; i++) f[i] = (f[i - 1] + qpow(i, d - 1)) % MOD; 
        b[0] = 1;
    }
    void work() {
        for(int i = 0; i <= d; i++) {
            for(int j = i + 1; j; j--) b[j] = (b[j - 1] + MOD - 1ll * b[j] * (i + 1) % MOD) % MOD;
            b[0] = 1ll * b[0] * (MOD - i - 1) % MOD;
        }
        for(int i = 0; i <= d; i++) {
            int s = f[i + 1], inv = qpow(i + 1, MOD - 2);
            for(int j = 0; j <= d; j++) if(i != j) s = 1ll * s * qpow((i - j + MOD) % MOD, MOD - 2) % MOD;
            b[0] = 1ll * b[0] * (MOD - inv) % MOD;
            for(int j = 1; j <= d + 1; j++) b[j] = (MOD - 1ll * (b[j] + MOD - b[j - 1]) * inv % MOD) % MOD;
            for(int j = 0; j <= d + 1; j++) a[j] = (a[j] + 1ll * s * b[j]) % MOD;
            for(int j = d + 1; j; j--) b[j] = (b[j - 1] + MOD - 1ll * b[j] * (i + 1) % MOD) % MOD;
            b[0] = 1ll * b[0] * (MOD - i - 1) % MOD;
        }   
    }
}A;

int d, w;
ll prod[N];

void run(){
    A.init(d + 1);
    A.work();
    for(int i = 0; i <= d + 1; i++) prod[i] = 1;
    while(w--) {
        int p, a; cin >> p >> a;
        for(int i = 0; i <= d + 1; i++) {
            ll res = (qpow(p, 1ll * a * i) - qpow(p, 1ll * a * i + d) * qpow(qpow(p, i), MOD - 2) % MOD + MOD) % MOD;
            prod[i] = prod[i] * res % MOD;
        }   
    }
    ll ans = 0;
    for(int i = 0; i <= d + 1; i++) ans = (ans + A.a[i] * prod[i] % MOD) % MOD;
    cout << ans << '\n';
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0); cout.tie(0);
    cout << fixed << setprecision(20);
    while(cin >> d >> w) run();
	return 0;
}