[51nod 1847]奇怪的数学题

题目

点这里看题目。

分析

是挺奇怪的......

以下定义质数集合为\(P\)，\(p_i\)为第\(i\)个质数。

定义\(mp(x)\)为\(x\)的最小质因子，则可以得到：

\[sgcd(a,b)=\frac{\gcd(a,b)}{mp(\gcd(a,b))} \]

这个比较显然。然后可以娴熟地变换式子得到：

\[\begin{aligned}\sum_{i=1}^n\sum_{j=1}^n sgcd(i,j)^k&=\sum_{d=2}^n\left(\frac d{mp(d)}\right)^k \sum_{i=1}^{\lfloor\frac n d\rfloor}\sum_{j=1}^{\lfloor\frac n d\rfloor}[\gcd(i,j)=1]\\&=\sum_{d=2}^n\left(\frac d{mp(d)}\right)^k \left(2\sum_{i=1}^{\lfloor\frac n d\rfloor}\varphi(i)-1\right)\end{aligned} \]

后面的求\(\varphi\)的前缀和的部分可以用杜教筛等筛法高速求出。这里推荐杜教筛，因为它可以方便地记忆化。

并且可以发现，\(\varphi\)的前缀和的上限实际上可以整除分块。所以，如果可以告诉求出\(\left(\frac d{mp(d)}\right)^k\)，我们就可以整除分块了。这也是杜教筛比较占优势的原因——一次杜教筛就可以筛出所有需要的值，并且存下来。

问题变成了如何解决前一部分的求和。考虑使用 min_25。

min_25 首先需要算出在质数位置的贡献，不难发现这些位置的贡献为\(1^k=1\)。因此只需要筛出质数数量即可。

其它的位置需要筛。首先有常见操作：

设\(f(x)=x^k\)，\(g(a,b)\)为前\(a\)个数进行了\(b\)轮埃筛之后的\(f\)的总贡献。

转移略。这是常见操作。

不过仔细想想，我们会发现转移出现的\(g(\lfloor\frac a{p_b}\rfloor, b-1)-g(p_{b-1},b-1)\)实际上就是 " 最小质因数为\(p_b\)的数的贡献 " （最小质因数少乘上一次）。因此我们实际操作的时候就不用写 min_25 的第二步，而是直接每轮累加起来，就可以得到合数的贡献。

而质数的贡献已经知道是质数个数了。所以对于每次询问，我们可以直接将合数和质数的贡献加起来回答。

但是，\(g(a,0)=\sum_{i=2}^a i^k\)，这个该怎么计算呢？

可以用第二类斯特林数来处理：

\[\begin{aligned} \sum_{i=1}^n i^k &=\sum_{i=1}^n \sum_{j=1}^k {k\brace j} i^{\underline j}\\ &=\sum_{j=1}^k {k\brace j}\sum_{i=1}^n i^{\underline j}\\ &=\sum_{j=1}^k{k\brace j}\frac{(n+1)^{\underline{j+1}}}{j+1} \end{aligned}\]

其中一步推导用到了 " 离散微积分 " 的东西，可以参考[第二类斯特林数]自然数幂求和。

需要注意的是，由于这道题取模方法是自然溢出，因此不能求逆元。但由于求幂和的时候结果一定是整数，所以在\((n+1)^{\underline{j+1}}\)里面一定有一个\(j+1\)的倍数。我们可以先用余数把它找出来，除掉\(j+1\)之后再将剩下的乘起来即可。

代码

#include <map>
#include <cmath>
#include <cstdio>
using namespace std;

typedef long long LL;
typedef unsigned int ui;

const int MAXN = 1e5 + 5, MAXK = 55;

template<typename _T>
void read( _T &x )
{
	x = 0;char s = getchar();int f = 1;
	while( s > '9' || s < '0' ){if( s == '-' ) f = -1; s = getchar();}
	while( s >= '0' && s <= '9' ){x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar();}
	x *= f;
}

template<typename _T>
void write( _T x )
{
	if( x < 0 ){ putchar( '-' ); x = ( ~ x ) + 1; }
	if( 9 < x ){ write( x / 10 ); }
	putchar( x % 10 + '0' );
}

map<int, ui> mp;

ui G[MAXN << 1], g1[MAXN << 1], gk[MAXN << 1];
ui ps[MAXN], pks[MAXN], pk[MAXN];
ui S[MAXK][MAXK];
int val[MAXN << 1];
int id1[MAXN], id2[MAXN];
int phi[MAXN], prime[MAXN], pn;
int N, K, s, cnt;
bool isPrime[MAXN];

int& ID( const int x ) { return x <= s ? id1[x] : id2[N / x]; }

ui qkpow( ui base, int indx )
{
	ui ret = 1;
	while( indx )
	{
		if( indx & 1 ) ret *= base;
		base *= base, indx >>= 1;
	}
	return ret;
}

void EulerSieve( const int siz )
{
	phi[1] = 1, isPrime[1] = true;
	for( int i = 2 ; i <= siz ; i ++ )
	{
		if( ! isPrime[i] ) prime[++ pn] = i, phi[i] = i - 1;
		for( int j = 1 ; j <= pn && 1ll * i * prime[j] <= siz ; j ++ )
		{
			isPrime[i * prime[j]] = true;
			if( ! ( i % prime[j] ) ) { phi[i * prime[j]] = phi[i] * prime[j]; break; }
			phi[i * prime[j]] = phi[i] * ( prime[j] - 1 );
		}
	}
	for( int i = 1 ; i <= pn ; i ++ ) pks[i] = pks[i - 1] + ( pk[i] = qkpow( prime[i], K ) );
	for( int i = 1 ; i <= siz ; i ++ ) ps[i] = ps[i - 1] + phi[i];
}

ui getS( const int a )
{
	ui ret = 0, t;
	for( int i = 1, tmp ; i <= K ; i ++ )
	{
		tmp = ( a + 1 ) % ( i + 1 ), t = S[K][i];
		for( int j = 1 ; j <= tmp ; j ++ ) t *= ( a - j + 2 );
		t *= ( a - tmp + 1 ) / ( i + 1 );
		for( int j = tmp + 2 ; j <= i + 1 ; j ++ ) t *= ( a - j + 2 );
		ret += t;
	}
	return ret;
}

ui SPhi( const int n )
{
	if( n <= s ) return ps[n];
	if( mp[n] ) return mp[n];
	ui ret;
	if( n & 1 ) ret = ( ui ) ( n + 1 ) / 2 * n;
	else ret = ( ui ) n / 2 * ( n + 1 );
	for( int l = 2, r ; l <= n ; l = r + 1 )
	{
		r = n / ( n / l );
		ret -= ( r - l + 1 ) * SPhi( n / l );
	}
	return mp[n] = ret;
}

int main()
{
	read( N ), read( K );
	s = sqrt( N );
	EulerSieve( s );
	for( int i = 1 ; i <= K ; i ++ ) S[i][i] = 1, S[i][0] = 0;
	for( int i = 2 ; i <= K ; i ++ )
		for( int j = 1 ; j <= K ; j ++ )
			S[i][j] = S[i - 1][j - 1] + ( ui ) j * S[i - 1][j];
	for( int l = 1, r, v ; l <= N ; l = r + 1 )
	{
		r = N / ( v = N / l ), val[++ cnt] = v; 
		g1[ID( v ) = cnt] = v - 1, gk[cnt] = getS( v ) - 1;
	}
	for( int j = 1, k ; j <= pn ; j ++ )
		for( int i = 1 ; i <= cnt && 1ll * prime[j] * prime[j] <= val[i] ; i ++ )
		{
			k = ID( val[i] / prime[j] );
			g1[i] -= g1[k] - ( j - 1 ); 
			gk[i] -= ( ui ) pk[j] * ( gk[k] - pks[j - 1] );
			G[i] += gk[k] - pks[j - 1];
		}
	ui ans = 0, pre = 0, nxt;
	for( int l = 2, r, v ; l <= N ; l = r + 1 )
	{
		r = N / ( v = N / l );
		nxt = G[ID( l )] + g1[ID( l )];
		ans += ( ui ) ( 2u * SPhi( v ) - 1 ) * ( nxt - pre ), pre = nxt;
	}
	write( ans ), putchar( '\n' );
	return 0;
}