题面图片真是大到离谱……
题目要求的是
\(\begin{align*}\sum\limits_{i=1}^N i^d[gcd(i,n) == 1] &= \sum\limits_{i=1}^N i^d \sum\limits_{p \mid gcd(i,n)} \varphi(p) \\ &= \sum\limits_{p|n} \varphi(p) p^d \sum\limits_{i=1}^{\frac{n}{p}} i^d\end{align*}\)
然后就不会做了qwq，后面的自然数次幂和似乎和前面的\(\varphi(p)p^d\)没什么关系
注意到\(f(x) = \sum\limits_{i=1}^x i^d\)可以写成一个\(d+1\)次多项式，即\(f(x) = \sum\limits_{i=0}^{d+1} a_ix^i\)，将这个代入上面的式子
\(\begin{align*}\sum\limits_{p|n} \varphi(p) p^d \sum\limits_{i=1}^{\frac{n}{p}} i^d &= \sum\limits_{p|n} \varphi(p) p^d \sum\limits_{i=0}^{d+1} a_i(\frac{n}{p})^i \\ &= \sum\limits_{i=0}^{d+1} a_i n^i \sum\limits_{p \mid n} \varphi(p)p^{d-i} \end{align*}\)
\(\varphi(p)p^{d-i}\)是一个积性函数，\(\varphi(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;
}