BZOJ3601 一个人的数论 莫比乌斯反演、高斯消元

传送门

---

题面图片真是大到离谱……

题目要求的是

(egin{align*}sumlimits_{i=1}^N i^d[gcd(i,n) == 1] &= sumlimits_{i=1}^N i^d sumlimits_{p mid gcd(i,n)} mu(p) \ &= sumlimits_{p|n} mu(p) p^d sumlimits_{i=1}^{frac{n}{p}} i^dend{align*})

然后就不会做了qwq，后面的自然数次幂和似乎和前面的(mu(p)p^d)没什么关系

注意到(f(x) = sumlimits_{i=1}^x i^d)可以写成一个(d+1)次多项式，即(f(x) = sumlimits_{i=0}^{d+1} a_ix^i)，将这个代入上面的式子

(egin{align*}sumlimits_{p|n} mu(p) p^d sumlimits_{i=1}^{frac{n}{p}} i^d &= sumlimits_{p|n} mu(p) p^d sumlimits_{i=0}^{d+1} a_i(frac{n}{p})^i \ &= sumlimits_{i=0}^{d+1} a_i n^i sumlimits_{p mid n} mu(p)p^{d-i} end{align*})

(mu(p)p^{d-i})是积性函数，(mu(p))又需要保证(p)中每一个质因数的指数只能为(1)，给出(n)的方式又是质因数分解之后的方式，所以可以比较方便地计算所有质数的贡献。

所以唯一的问题就是如何求出(a_i)。这个高斯消元和拉格朗日插值都可以。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
//This code is written by Itst
using namespace std;

const int MAXN = 1007 , MOD = 1e9 + 7;
int N , D , W , prm[MAXN] , e[MAXN] , gauss[107][107];

inline int poww(long long a , int b){
	int times = 1;
	if(b < 0) b += MOD - 1;
	while(b){
		if(b & 1) times = times * a % MOD;
		a = a * a % MOD;
		b >>= 1;
	}
	return times;
}

void init(){
	int sum = 0;
	for(int i = 0 ; i <= D + 1 ; ++i){
		int cur = 1;
		sum = (sum + poww(i , D)) % MOD;
		for(int j = 0 ; j <= D + 1 ; ++j){
			gauss[i][j] = cur;
			cur = 1ll * cur * i % MOD;
		}
		gauss[i][D + 2] = sum;
	}
	for(int i = 0 ; i <= D + 1 ; ++i){
		int j = i;
		while(j <= D + 1 && !gauss[j][i])
			++j;
		swap(gauss[i] , gauss[j]);
		int inv = poww(gauss[i][i] , MOD - 2);
		for(int k = i ; k <= D + 2 ; ++k)
			gauss[i][k] = 1ll * gauss[i][k] * inv % MOD;
		while(++j <= D + 1)
			if(gauss[j][i])
				for(int k = D + 2 ; k >= i ; --k)
					gauss[j][k] = (gauss[j][k] - 1ll * gauss[i][k] * gauss[j][i] % MOD + MOD) % MOD;
	}
	for(int i = D + 1 ; i >= 0 ; --i)
		for(int j = i - 1 ; j >= 0 ; --j)
			gauss[j][D + 2] = (gauss[j][D + 2] - 1ll * gauss[j][i] * gauss[i][D + 2] % MOD + MOD) % MOD;
}

int main(){
#ifndef ONLINE_JUDGE
	//freopen("in","r",stdin);
	//freopen("out","w",stdout);
#endif
	cin >> D >> W;
	N = 1;
	for(int i = 1 ; i <= W ; ++i){
		cin >> prm[i] >> e[i];
		N = 1ll * N * poww(prm[i] , e[i]) % MOD;
	}
	init();
	int ans = 0;
	for(int i = 0 ; i <= D + 1 ; ++i){
		int sum = 1;
		for(int j = 1 ; j <= W ; ++j)
			sum = 1ll * sum * (1 - poww(prm[j] , D - i) + MOD) % MOD;
		ans = (ans + 1ll * sum * poww(N , i) % MOD * gauss[i][D + 2]) % MOD;
	}
	cout << ans;
	return 0;
}