P4464 [国家集训队]JZPKIL

如果不知道伯努利数的点这里。

\[\sum_{i=1}^n (i,n)^x[i,n]^y \]

\[\sum_{i=1}^n (i,n)^{x-y} \times i^y \times n^y \]

\[n^y \sum_{d|n} d ^{x-y} \sum_{i=1}^n [(i,n)=d] i^y \]

\[n^y \sum_{d|n} d ^{x} \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} [(i,\frac{n}{d})=1] i^y \]

\[n^y \sum_{d|n} d ^{x} \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \sum_{s|(i,\frac{n}{d})}\mu(s)i^y \]

\[n^y \sum_{d|n} d ^{x} \sum_{s|\frac{n}{d}} \mu(s) \sum_{i=1}^{\frac{n}{ds}}(is)^y \]

\[n^y \sum_{d|n} d ^{x} \sum_{s|\frac{n}{d}} \mu(s) s^y \sum_{i=1}^{\frac{n}{ds}}i^y \]

设 \(k_i\) 为 \(y\) 次幂展开后 \(n^i\) 的系数

\[n^y \sum_{d|n} d ^{x} \sum_{s|\frac{n}{d}} \mu(s) s^y \sum_{i=1}^{y+1} k_i (\frac{n}{ds})^i \]

\[n^y \sum_{i=1}^{y+1} k_i\sum_{d|n} d ^{x} \sum_{s|\frac{n}{d}} \mu(s) s^y (\frac{n}{ds})^i \]

后面是卷积的形式，记为 \(f(i)\)。

记 \(f_1(n)=n^x,f_2(n)=\mu(n)n^y,f_3(n)=n^i\)。

那么 \(f=f_1*f_2*f_3\) ，也为一个积性函数,考虑求出 \(f(p^w)\) 。

\(s\) 只有等于 \(1\) 或 \(p\) 时 \(\mu(s)\not=0\),

当 \(s = 1\)

\[f(n)=\sum_{d|n} d ^{x} \sum_{s|\frac{n}{d}} \mu(s) s^y (\frac{n}{ds})^i \]

\[=n^i\sum_{d|n} d ^{x-i} \]

\[f(p^k)=p^{ki}\sum_{d|p^k}d^{x-i} \]

\[=p^{ki}\sum_{j=0}^k p^{j(x-i)} \]

当 \(s=p\)

\[f(n)=\sum_{d|n} d ^{x} \sum_{s|\frac{n}{d}} \mu(s) s^y (\frac{n}{ds})^i \]

\[=\sum_{d|n} d ^{x} \mu(p) p^y (\frac{n}{dp})^i \]

\[=-n^i\sum_{d|n} d ^{x-i} p ^{y-i} \]

\[f(p^k)=-p^{ki}\sum_{d|p^k}d^{x-i}p ^{y-i} \]

\[=-p^{ki}\sum_{j=0}^{k-1} p^{j(x-i)}p ^{y-i} \]

\[=-p^{ki+y-i}\sum_{j=0}^{k-1} p^{j(x-i)} \]

然后用 \(\text{Pollard~rho}\) 分解质因数后计算。

但计算等比数列和时不能用公式算，也不能暴力算（具体看代码注释）。

这份常数极大的代码不开o2会t。

#include <ctime>
#include <cstdio>
#include <cstdlib>
#include <algorithm>
using namespace std;
#define Lint __int128
#define LL long long

template<typename _T>
void Read( _T &x ) {
	x = 0; int f = 1;
	char s = getchar( );
	for( ; s < '0' || s > '9' ; s = getchar( ) ) f = s == '-' ? -f : f;
	for( ; s >= '0' && s <= '9' ; s = getchar( ) ) x = x * 10 + s - '0';
	x *= f;
}
template<typename _T>
void Write( _T x ) {
	if( x < 0 ) putchar( '-' ) , x = -x;
	if( x >= 10 ) Write( x / 10 );
	putchar( x % 10 + '0' );
}

const int MAXN = 3005 , Mod = 1e9 + 7;
LL n;
int t , x , y;

template<typename _T>
_T Quick_pow( _T x , _T po , _T p = Mod ) {
	_T Ans = 1; bool f = 0;
    if( po < 0 ) po = -po , f = 1;
	while( po ) {
		if( po % 2 )
			Ans = (Lint)Ans * x % p;
		x = (Lint)x * x % p;
		po /= 2;
	}
	return f ? Quick_pow( (int)Ans , Mod - 2 ) : Ans;
}
template<typename _T>
_T Gcd( _T x , _T y ) {
	return y == 0 ? x : Gcd( y , x % y );
}
template<typename _T>
_T Abs( _T x ) {
	return x >= 0 ? x : -x;
}
template<typename _T>
_T Inverse( _T x ) {
    return Quick_pow( x , Mod - 2 );
}

template<typename _T>
_T f( _T x , _T a , _T Mod ) {
	return ( (Lint)x * x + a ) % Mod;
}
LL Pollard_rho( LL n ) {
	LL r = rand( ) % ( n - 1 ) + 1;
	LL s = 0 , t = 0 , val , d;
	for( int k = 1 ; k ; k <<= 1 ) {
		val = 1; s = t;
		for( int i = 1 ; i <= k ; i ++ ) {
			t = f( t , r , n );
			val = (Lint)val * Abs( t - s ) % n;
			if( i % 127 == 0 && ( d = Gcd( val , n ) ) > 1 ) return d;
		}
		if( ( d = Gcd( val , n ) ) > 1 ) return d;
	}
}

int prinum = 12;
LL prime[ 20 ] = { 0 , 2 , 3 , 5 , 7 , 11 , 13 , 17 , 19 , 23 , 29 , 31 , 37 };
bool Miller_Rabin( LL x ) {
	if( x == 0 || x == 1 ) return 0;
	for( int i = 1 ; i <= prinum ; i ++ ) {
		if( x == prime[ i ] ) return 1;
		if( x < prime[ i ]  ) return 0;
		
		LL s = Quick_pow( prime[ i ] , x - 1 , x ) , tmp = x - 1;
		if( s != 1 ) return 0; 
		while( s == 1 && tmp % 2 == 0 ) {
			tmp /= 2;
			s = Quick_pow( prime[ i ] , tmp , x );
			if( s != 1 && s != x - 1 ) return 0;
		}
	}
	return 1;
}

int k; LL pi[ 100 ];
void Divide( LL x ) {
	if( x == 1 ) return;
	if( Miller_Rabin( x ) ) {
		pi[ ++ k ] = x; return;
	}
	if( x == 4 ) {
		pi[ ++ k ] = 2 , pi[ ++ k ] = 2;
		return;
	}
	LL p = x;
	while( p == x ) p = Pollard_rho( x );
	Divide( p ) , Divide( x / p );
} 

int Fac[ MAXN + 5 ] , Inv[ MAXN + 5 ];
void Init_Fac( ) {
    Fac[ 0 ] = 1;
    for( int i = 1 ; i <= MAXN ; i ++ ) Fac[ i ] = 1ll * Fac[ i - 1 ] * i % Mod;
    Inv[ MAXN ] = Inverse( Fac[ MAXN ] );
    for( int i = MAXN ; i >= 1 ; i -- ) Inv[ i - 1 ] = 1ll * Inv[ i ] * i % Mod;
}
int C( int n , int m ) {
    return 1ll * Fac[ n ] * Inv[ m ] % Mod * Inv[ n - m ] % Mod;
}

int Bernoulli[ MAXN + 5 ] , S[ MAXN + 5 ];
void Init_Bernoulli( ) {
    Bernoulli[ 0 ] = 1;
    for( int i = 1 ; i <= MAXN ; i ++ ) {
        int Sum = 0;
        for( int j = 0 ; j < i ; j ++ )
            Sum = ( Sum + 1ll * C( i + 1 , j ) * Bernoulli[ j ] % Mod ) % Mod;
        Bernoulli[ i ] = ( -1ll * Sum * Inverse( C( i + 1 , i ) ) % Mod + Mod ) % Mod;
    }Bernoulli[ 1 ] ++;
}

LL di[ MAXN + 5 ];
int wi[ MAXN + 5 ];

/*int Aq( int a , int q , int n ) {
    if( q == 1 ) return 1ll * a * n % Mod;
    return 1ll * a * ( 1 - Quick_pow( q , n ) + Mod ) % Mod * Inverse( 1 - q ) % Mod;
}*/
int F( int p , int w , int i ) {
//	int Ans = 0;
//	for( int j = 0 ; j <= w ; j ++ ) Ans = ( Ans + Quick_pow( p , w * i + j * ( x - i ) ) ) % Mod;
//	for( int j = 0 ; j <  w ; j ++ ) Ans = ( Ans - Quick_pow( p , w * i + y - i + j * ( x - i ) ) ) % Mod;
//	return ( Ans + Mod ) % Mod;

//	return ( ( Aq( Quick_pow( p , w * i ) , Quick_pow( p , x - i ) , w + 1 ) - Aq( Quick_pow( p , w * i + y - i ) , Quick_pow( p , x - i ) , w ) ) % Mod + Mod ) % Mod;

	int Ans = 0 , a , q = Quick_pow( p , x - i );
	a = Quick_pow( p , w * i );
	for( int j = 0 ; j <= w ; j ++ ) Ans = ( Ans + a ) % Mod , a = 1ll * a * q % Mod; 
	a = Quick_pow( p , w * i + y - i );
	for( int j = 0 ; j <  w ; j ++ ) Ans = ( Ans - a ) % Mod , a = 1ll * a * q % Mod; 
	return ( Ans + Mod ) % Mod;
}

/*
p^{wi}*p^{j(x-i)}
a1:p^{wi} q:p^(x-i) n:w+1

p^{wi+y-i}*p^{j(x-i)}
a1:p^{wi+y-i} q:p^(x-i) n:w
*/

int Solve( LL n , int x , int y ) {
	if( n == 0 ) return 0;
	k = 0; Divide( n );
	sort( pi + 1 , pi + k + 1 );
	
	S[ 0 ] = 0;
	for( int i = 0 ; i <= y ; i ++ )
		S[ y + 1 - i ] = 1ll * Inverse( y + 1 ) * C( y + 1 , i ) % Mod * Bernoulli[ i ] % Mod;

	int cnt = 0;
	for( int i = 1 ; i <= k ; i ++ ) {
		if( pi[ i ] != di[ cnt ] ) di[ ++ cnt ] = pi[ i ] , wi[ cnt ] = 1;
		else wi[ cnt ] ++;
	}

	int Ans = 0;
	for( int i = 0 ; i <= y + 1 ; i ++ ) {
		int tmp = 1;
		for( int j = 1 ; j <= cnt ; j ++ )
			tmp = 1ll * tmp * F( (int)( di[ j ] % Mod ) , wi[ j ] , i ) % Mod;
        Ans = ( Ans + 1ll * S[ i ] * tmp % Mod ) % Mod;
	}
	return 1ll * Quick_pow( int( n % Mod ) , y ) * Ans % Mod;
}

signed main() {
    srand( time( 0 ) );
    Init_Fac( );
    Init_Bernoulli( );

	Read( t );
    while( t -- ) {
        Read( n ) , Read( x ) , Read( y );
		Write( Solve( n , x , y ) ) , putchar('\n');
    }
    return 0;
}

------------恢复内容结束------------