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; }
- BZOJ 3601 一个人的数论 莫比乌斯反演+高斯消元
- [莫比乌斯反演 高斯消元 数学技巧] BZOJ 3601 一个人的数论
- [bzoj3601] 一个人的数论 [莫比乌斯反演+高斯消元]
- 【bzoj3601】一个人的数论 莫比乌斯反演+高斯消元
- bzoj 3601: 一个人的数论 高斯消元&莫比乌斯反演
- 【BZOJ3601】一个人的数论 高斯消元+莫比乌斯反演
- 【BZOJ3601】一个人的数论,莫比乌斯反演+高斯消元
- bzoj3601 一个人的数论 莫比乌斯反演+高斯消元
- BZOJ 3601 一个人的数论 ——莫比乌斯反演 高斯消元
- BZOJ3601 一个人的数论
- BZOJ3601:一个人的数论(莫比乌斯反演+伯努利数)
- BZOJ_P1013 [JSOI2008]球形空间产生器sphere(数论+高斯消元)
- BZOJ3601 一个人的数论
- BZOJ-1013 球形空间产生器sphere 高斯消元+数论推公式
- BZOJ-1013 球形空间产生器sphere 高斯消元+数论推公式
- 【数论】bzoj3601一个人的数论
- BZOJ 3601: 一个人的数论
- 【bzoj 3601】一个人的数论 (莫比乌斯反演+伯努利数)
- [BZOJ3601] 一个人的数论 - 拉格朗日插值/伯努利数,狄利克雷卷积
- 高斯消元裸题 记录模板 BZOJ 1013