题目链接

题解

挺神的

首先有

\[\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;

}

巴特西

BZOJ3601 一个人的数论【数论 + 高斯消元】

题目链接

题解

最新文章

热门文章

巴特西

BZOJ3601 一个人的数论 【数论 + 高斯消元】

题目链接

题解

最新文章

热门文章

BZOJ3601 一个人的数论【数论 + 高斯消元】