自动搬运

查看原文

来自洛谷，原作者为

	Soulist 再见了

搬运于2025-08-24 22:00:28，当前版本为作者最后更新于2020-06-29 10:32:29，作者可能在搬运后再次修改，您可在原文处查看最新版

自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

以下是正文

---

多次查询，每次给定 $n,x,y$，求：

$$\sum_{i=1}^n \gcd(i,n)^x\textrm{lcm}(i,n)^y \pmod {10^9+7} $$

• $T\le 100,1\le n\le 10^{18},x,y\le 3000$

---

$\rm Sol:$

答案是：

$$\sum_{i=1}^n i^y\cdot n^y\gcd(i,n)^{x-y}$$

枚举 $\gcd=d$，忽略 $n^y$ 则有：

$$\begin{aligned} &\sum_{d|n} d^x\sum_{j=1}^{\frac{n}{d}}[\gcd(j,\frac{n}{d})=1]j^{y} \\&=\sum_{d|n} d^x\sum_{j=1}^{n/d}\sum_{k|j,k|\frac{n}{d}} \mu(k)j^y \\&=\sum_{d|n} d^x\sum_{k|\frac{n}{d}}\mu(k) \sum_{k|j} j^y \\&=\sum_{d|n} d^x\sum_{k|\frac{n}{d}}\mu(k)k^y \sum_{j=1}^{n/(kd)} j^y \\&=\sum_{d|n} d^x \sum_{k|\frac{n}{d}}\mu(k)k^y S_y(\frac{n}{kd}) \end{aligned}$$

我们将 $S_y(\frac{n}{kd})$ 视为关于其的 $y+1$ 次多项式，那么此处有：

$$\begin{aligned} &\sum_{d|n} d^x \sum_{k|\frac{n}{d}}\mu(k)k^y \sum_{i=0}^{y+1}f_i (\frac{n}{kd})^i \\&=\sum_{i=0}^{y+1}f_i\sum_{d|n} d^x \sum_{k|\frac{n}{d}}\mu(k)k^y(\frac{n}{kd})^i \end{aligned}$$

注意到后者可以视为 $f,g,h$ 的迪利克雷卷积，令 $f(d)=d^x,g(d)=\mu(d)d^y,h^c(d)=d^c$，这个卷积满足分配律和交换律以及结合律，同时注意到三者均为积性函数，所以卷积结果也是积性函数。

设结果为 $\Z(x)$，于是我们所求为 $\Z(n)$，那么显然有：

1. $\Z(1)=1$
2. $\Z(p)=p^x+p^c-p^y$
3. $\Z(p^k)=f\cdot h^c(p^k)-f\cdot h^c(p^{k-1})\times p^y$

于是只需要考虑 $f\cdot h^c(p^k)$，其为：

$$\sum_{i+j=k} f(p^i)\cdot h^c(p^j)$$

所以通过 MR 进行暴力质因数分解，然后我们 $\mathcal O(\log n)$ 的枚举每个质因子（include 次幂）然后暴力卷积并计算答案即可。

还有一个问题，如何计算自然数幂的和的多项式？

自然数幂和多项式

构造多项式 $G_n(x)=\frac{e^{nx}-1}{1-e^{-x}}$，一堆神仙分析可以得到 $G_n(x)[x^k]$ 就是 $\sum_{i=1}^n i^k$（其实应该是根据 $G_n(x)[x^k]$ 是 $\sum_{i=1}^n i^k$ 推出 $G_n(x)$ 的表达式的...）

将 $G_n(x)$ 裂项成 $\frac{xe^x}{e^x-1}\times \frac{e^{nx}-1}{x}$，前者的 EGF 序列记为 $\{B_0,B_1,B_2...\}$，后者显然是 EGF 序列 $\{\frac{n}{1},\frac{n^2}{2},\frac{n^3}{3}...\}$

所以我们会发现其实有 $G_n(x)[x^k]$ 是这两个数列的二项卷积的结果，所以有：

$$\begin{aligned} &G_n[x^k]=\sum_{i+j=k} B_i \frac{n^{j+1}}{j+1}\binom{k}{i} \\&= \frac{1}{k+1}\sum_{i=0}^{k} B_i \binom{k+1}{i}n^{k-i+1} \end{aligned}$$

换而言之如果能够得到伯努利数那么就可以知道自然数幂和多项式

然后只需要递推伯努利数即可。

注意到 $\frac{xe^x}{e^x-1}\times \frac{e^x-1}{x}=e^x$ 所以 $\mathbf{EGF}\{B_0,B_1...\}\times \mathbf{EGF}\{\frac{1}{1},\frac{1}{2}...\}=\mathbf{EGF}\{1,1,1...\}$

所以有：

$$\begin{aligned} &\sum_{i+j=k} B_i\times \frac{1}{j+1}\times \binom{k}{i}=1 \\& \sum_{i}^k B_i\times\binom{k+1}{i}=k+1 \end{aligned}$$

于是得到 $B_k=1-\sum_{i<k} \frac{B_i}{k+1}\binom{k+1}{i}$

综上，本题复杂度为 $\mathcal O(T\cdot (n^{\frac{1}{4}}\log n+y\log^2 n)+y^2)$

---

然后这道题卡常，你需要尽可能的省略快速幂，如果你是 $\log^3$ 然后还是我这种大常数玩家就会死得非常惨。

然后千万不要写龟速乘，然后挑一个夜深人静的时候交一发，多洗一下脸，没准就过了。

不过我在 darkbzoj 这道题过得非常稳...但是咕咕就死得非常惨

$Code:$

#include<bits/stdc++.h>
using namespace std ;
#define Next( i, x ) for( register int i = head[x]; i; i = e[i].next )
#define rep( i, s, t ) for( register int i = (s); i <= (t); ++ i )
#define drep( i, s, t ) for( register int i = (t); i >= (s); -- i )
#define re register
#define int __int128
int gi() {
	char cc = getchar() ; int cn = 0, flus = 1 ;
	while( cc < '0' || cc > '9' ) {  if( cc == '-' ) flus = - flus ; cc = getchar() ; }
	while( cc >= '0' && cc <= '9' )  cn = cn * 10 + cc - '0', cc = getchar() ;
	return cn * flus ;
}
const int M = 3000 + 5 ; 
const int P = 1e9 + 7 ;
int n, fx, fy, fc, maxn, top, cnt, num, tp[M], st[M], w[M], fac[M], inv[M], B[M], f[M] ; 
int mul( int a, int b, int p ) {
	return a * b % p ; 
}
int fpow( int x, int k, int p ) {
	int ans = 1, base = x ; 
	while(k) {
		if(k & 1) ans = mul(ans, base, p) ;
		base = mul(base, base, p), k >>= 1 ; 
	} return ans ; 
}
int sw( int x ) { return ( ( 1ll * rand() << 16ll ) | rand() ) % x + 1 ; }
bool M_R( int p ) {
	if( p == 2 || p == 3 ) return 1 ; 
	if( p == 1 || ( p % 2 == 0 ) ) return 0 ; 
	re int d = p - 1, s = 0 ; while( !( d & 1 ) ) ++ s, d >>= 1 ; 
	rep( i, 0, 5 ) {
		re int a = sw( p - 3 ) + 2 ; 
		re int x = fpow( a, d, p ), y = 0 ;
		for( register int j = 0; j < s; ++ j ) {
			y = mul( x, x, p ) ; if( y == 1 && ( x != 1 ) && x != ( p - 1 ) ) return 0 ;
			x = y ;
		} if( y != 1 ) return 0 ;
	} return 1 ; 
}
int gcd( int x, int y ) {
	return ( x == 0 ) ? y : gcd( y % x, x ) ;
}
int work( int p ) {
	re int k = 2, x = sw(p - 1), y = x, d = 1, c = sw(99997) ;
	for( re int i = 1; d == 1; ++ i ) {
		x = ( mul( x, x, p ) + c ) % p ;
		d = gcd( (x > y) ? x - y : y - x, p ) ;
		if( i == k ) y = x, k <<= 1 ; 
	} return d ; 
}
void Pollard_Rho(int p) {
	if( p == 1 ) return ; 
	if( M_R(p) ) { st[++ top] = p; return ; }
	int x = p ; while( x == p ) x = work(x) ;
	Pollard_Rho(x), Pollard_Rho(p / x) ;
}
int g[M], gg[M], fg[M], fgg[M], fff[M] ; 
int c( int x, int y ) {
	if( x < y ) return 0 ;
	return fac[x] * inv[y] % P * inv[x - y] % P ;  
}
signed main()
{
	srand(time(NULL)) ;
	int T = gi() ; maxn = 3002 ; B[0] = 1, fac[0] = inv[0] = 1 ; 
	rep( i, 1, maxn ) fac[i] = fac[i - 1] * i % P, inv[i] = fpow( fac[i], P - 2, P ) ;
	rep( i, 1, maxn ) {
		for( re int j = 0; j < i; ++ j ) B[i] = ( B[i] - c(i + 1, j) * B[j] % P + P ) % P ; 
		B[i] = B[i] * fpow( i + 1, P - 2, P ) % P, ++ B[i] ;  
	}
	while( T-- ) {
		n = gi(), fx = gi(), fy = gi() ; int k = fy ; top = 0 ; 
		for( re int i = 0; i <= k; ++ i ) f[k - i + 1] = B[i] * c(k + 1, i) % P * fpow( k + 1, P - 2, P ) % P ; 
		Pollard_Rho(n), sort(st + 1, st + top + 1) ; cnt = 0 ; 
		rep( i, 1, top ) {
			if( st[i] != st[i - 1] ) tp[++ cnt] = st[i] ; 
			++ w[cnt] ; 
		} int Ans = 0 ; 
		rep( j, 1, cnt ) fg[j] = fpow( tp[j], fx, P ), fff[j] = fpow(tp[j], fy, P) ; 
		for( re int ic = 0; ic <= k + 1; ++ ic ) {
			int ans = 1 ;
			for( re int j = 1; j <= cnt; ++ j ) {
				int tk = w[j], qwq = 0, de = 0 ; 
				g[0] = 1, gg[0] = 1 ; fgg[j] = ( !ic ) ? 1 : fgg[j] * tp[j] % P ;  
				int dg = fg[j], dgg = fgg[j] ; 
				for( re int o = 1; o <= w[j]; ++ o ) 
					g[o] = g[o - 1] * dg % P, gg[o] = gg[o - 1] * dgg % P ; 
				for( re int o = 0; o <= w[j]; ++ o )
					qwq = ( qwq + g[o] * gg[(w[j] - o)] % P ) % P ; 
				for( re int o = 0; o < w[j]; ++ o ) 
					de = ( de + g[o] * gg[(w[j] - o - 1)] % P ) % P ; 
				ans = ans * ( qwq - de * fff[j] % P + P ) % P ; 
			}
			Ans = ( Ans + ans * f[ic] % P + P ) % P ; 
		} rep( i, 1, cnt ) w[i] = tp[i] = 0 ;  
		cout << (long long)(Ans * fpow( n, fy, P ) % P) << endl ; 
	}
	return 0 ;
}