【 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;

}

巴特西

[51nod 1847]奇怪的数学题

【 51nod 1847 】奇怪的数学题

题目

分析

代码

最新文章

热门文章