题目链接
题解
挺神的
首先有
[egin{aligned}
f(n) &= sumlimits_{x = 1}^{n} x^{d} [(x,n) = 1] \
&= sumlimits_{x = 1}^{n} x^{d} sumlimits_{c|(x,n)}mu(c) \
&= sumlimits_{c|n}sumlimits_{x = 1}^{frac{n}{c}} (cx)^{d} mu(c) \
&= sumlimits_{c|n}mu(c)c^{d}sumlimits_{x = 1}^{frac{n}{c}} x^{d} \
end{aligned}
]
我们记
[g(x) = sumlimits_{i = 1}^{x}i^{d}
]
然后就是最匪夷所思的地方，我们大力猜想这是关于(x)的一个(d + 1)次多项式
即
[g(x) = sumlimits_{i = 1}^{d + 1}a_ix^{i}
]
只需高斯消元得出系数(a_i)
【upd:其实很显然，展开(sumlimits_{i = 0}^{x - 1}(x - i)^{d})，(x^d)有(x)项，合并后就是一个关于(x)的(d + 1)次多项式】
然后(f(n))可以继续化简
[egin{aligned}
f(n) &= sumlimits_{c|n}mu(c)c^{d}g(frac{n}{c}) \
&= sumlimits_{c|n}mu(c)c^{d}sumlimits_{i = 1}^{d + 1} a_i(frac{n}{c})^{i} \
&= sumlimits_{i = 1}^{d + 1}a_isumlimits_{c|n}mu(c)c^{d}(frac{n}{c})^{i}
end{aligned}
]
后面是一个狄利克雷卷积
(F(x) = mu(x)x^{d})是一个积性函数，(F(x) = x^{i})显然也是一个积性函数
两个积性函数的狄利克雷卷积依旧是一个积性函数
所以我们只需计算(n)的所有质因子的函数值乘起来
所以我们记
[h(p^{k}) = sumlimits_{c|p^{k}}mu(c)c^{d}(frac{p^{k}}{c})^{i}
]
显然只有(mu(1))和(mu(p))两项
化简得
[h(p^{k}) = p^{ki}(1 - p^{d - i})
]
可以(O(1))计算
所以式子就化为
[f(n) = sumlimits_{i = 1}^{d + 1}a_iprod_{i=1}^{w}h(p_i^{k_i})
]
(O(dw))计算即可
总复杂度(O(d^3 + dw))
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<map>
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define mp(a,b) make_pair<int,int>(a,b)
#define cls(s) memset(s,0,sizeof(s))
#define cp pair<int,int>
#define LL long long int
using namespace std;
const int maxn = 105,maxm = 1005,INF = 1000000000,P = 1000000007;
inline int read(){
int out = 0,flag = 1; char c = getchar();
while (c < 48 || c > 57){if (c == '-') flag = -1; c = getchar();}
while (c >= 48 && c <= 57){out = (out << 3) + (out << 1) + c - 48; c = getchar();}
return out * flag;
}
int w,d,p[maxm],k[maxm],a[maxn];
int A[maxn][maxn],N;
inline int qpow(int a,LL b){
if (b < 0) b += P - 1;
int re = 1;
for (; b; b >>= 1,a = 1ll * a * a % P)
if (b & 1) re = 1ll * re * a % P;
return re;
}
void gause(){
for (int i = 1; i <= N; i++){
int j = i;
/*for (int k = i + 1; k <= N; k++)
if (A[k][i] > A[j][i]) j = k;
if (j != i) for (int k = i; k <= N + 1; k++) swap(A[j][k],A[i][k]);*/
for (j = i + 1; j <= N; j++){
int t = 1ll * A[j][i] * qpow(A[i][i],P - 2) % P;
for (int k = i; k <= N + 1; k++)
A[j][k] = ((A[j][k] - 1ll * A[i][k] * t % P) % P + P) % P;
}
}
for (int i = N; i; i--){
for (int j = i + 1; j <= N; j++)
A[i][N + 1] = ((A[i][N + 1] - 1ll * a[j] * A[i][j] % P) % P + P) % P;
a[i] = 1ll * A[i][N + 1] * qpow(A[i][i],P - 2) % P;
}
}
void cal(){
N = d + 1;
for (int x = 1; x <= N; x++){
A[x][N + 1] = (A[x - 1][N + 1] + qpow(x,d)) % P;
for (int j = 1; j <= N; j++) A[x][j] = qpow(x,j);
}
gause();
int s1 = 0,s2 = 0;
for (int i = 1; i <= N; i++) s1 = (s1 + 1ll * a[i] * qpow(5,i) % P) % P;
for (int i = 1; i <= 5; i++) s2 = (s2 + qpow(i,d)) % P;
}
void work(){
int ans = 0;
for (int i = 1; i <= N; i++){
int tmp = a[i];
for (int j = 1; j <= w; j++)
tmp = 1ll * tmp * qpow(p[j],1ll * k[j] * i) % P * (((1 - qpow(p[j],d - i)) % P + P) % P) % P;
ans = (ans + tmp) % P;
}
printf("%d
",ans);
}
int main(){
d = read(); w = read();
REP(i,w) p[i] = read(),k[i] = read();
cal();
work();
return 0;
}