BZOJ3601 一个人的数论 【数论 + 高斯消元】
题目链接
题解
挺神的
首先有
\[\begin{aligned}
f(n) &= \sum\limits_{x = 1}^{n} x^{d} [(x,n) = 1] \\
&= \sum\limits_{x = 1}^{n} x^{d} \sum\limits_{c|(x,n)}\mu(c) \\
&= \sum\limits_{c|n}\sum\limits_{x = 1}^{\frac{n}{c}} (cx)^{d} \mu(c) \\
&= \sum\limits_{c|n}\mu(c)c^{d}\sum\limits_{x = 1}^{\frac{n}{c}} x^{d} \\
\end{aligned}
\]
我们记
\[g(x) = \sum\limits_{i = 1}^{x}i^{d}
\]
然后就是最匪夷所思的地方，我们大力猜想这是关于\(x\)的一个\(d + 1\)次多项式
即
\[g(x) = \sum\limits_{i = 1}^{d + 1}a_ix^{i}
\]
只需高斯消元得出系数\(a_i\)
【upd:其实很显然，展开\(\sum\limits_{i = 0}^{x - 1}(x - i)^{d}\)，\(x^d\)有\(x\)项，合并后就是一个关于\(x\)的\(d + 1\)次多项式】
然后\(f(n)\)可以继续化简
\[\begin{aligned}
f(n) &= \sum\limits_{c|n}\mu(c)c^{d}g(\frac{n}{c}) \\
&= \sum\limits_{c|n}\mu(c)c^{d}\sum\limits_{i = 1}^{d + 1} a_i(\frac{n}{c})^{i} \\
&= \sum\limits_{i = 1}^{d + 1}a_i\sum\limits_{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}) = \sum\limits_{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) = \sum\limits_{i = 1}^{d + 1}a_i\prod_{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\n",ans);
}
int main(){
d = read(); w = read();
REP(i,w) p[i] = read(),k[i] = read();
cal();
work();
return 0;
}
