[LOJ6055]「from CommonAnts」一道数学题 加强版

题目

点这里看题目。

分析

首先可以娴熟地推倒一发式子:

\[\begin{split}&f_k(n)&\overset{\mathrm{def}}{=}f(k,n)\\\Rightarrow &f_k(n)&=\sum_{i=1}^{n-1}f_k(i)+n^k\\\end{split} \]

定义\(S_k(n)=\sum_{i=1}^{n}f_k(i)\),做差计算可以得到:

\[\begin{split}\Rightarrow &&S_k(n)&&=2S_k(n-1)+n^k\\&&&&=\dots\\&&&&=\sum_{i=1}^n2^{n-i}i^k\\&&&&=2^n\sum_{i=1}^n2^{-i}i^k\\&&&a&\overset{\mathrm {def}}{=}2^{-1}\\\Rightarrow &&S_k(n)&&=2^n\sum_{i=1}^na^ii^k\end{split} \]

自然我们可以得到:

\[f_k(n)=n^k+2^{n-1}\sum_{i=1}^{n-1}a^ii^k \]

为了方便,我们就可以定义出:

\[F_k(n)=\sum_{i=1}^{n-1}a^ii^k \]

问题于是转化为了:如何快速求出\(F_k(n)\)

60 pts ~ 80 pts

可以发现\(F_k(n)\)就是一个前缀和,我们尝试错位相减:

\[\begin{split}(1-a)F_k(n)&=\sum_{i=1}^na^ii^k-\sum_{i=1}^{n+1}a^i(i-1)^k\\&=\sum_{i=1}^na^i[i^k-(i-1)^k]-a^{n+1}n^k\\&=\sum_{i=1}^na^i\sum_{j=1}^k\binom k j(-1)^{j+1}i^{k-j}-a^{n+1}n^k\\&=\sum_{j=1}^k\binom k j(-1)^{j+1}\sum_{i=1}^na^ii^{k-j}-a^{n+1}n^k\\&=\sum_{j=1}^k\binom k j(-1)^{j+1}F_{k-j}(n)-a^{n+1}n^k\end{split} \]

并且当\(k=0\)时,\(F_0(n)=\sum_{i=1}^na^i=\frac{1-a^{n}}{1-a}\),可以快速计算。通过这个递推关系我们可以\(O(k^2)\)计算出\(F_k(n)\)。用一个任意模数 NTT 说不定可以优化到\(O(k\log_2k)\)

说不定就是说我也没试过

100 pts

考虑一个\(f_k(n)\)的不存在求和符的递推式:

\[\begin{split}f_k(n)-f_k(n-1)&=n^k+F_k(n-1)\\f_k(n)&=2f_k(n-1)+n^k-(n-1)^k\end{split} \]

尝试构造出一个\(g_k(n)\)满足:

\[f_k(n)+g_k(n)=2(f_k(n-1)+g_k(n-1)) \]

显然一个可行的解为:

\[g_k(n)=2g_k(n-1)-n^k+(n-1)^k \]

即可以解出:

\[f_k(n)=2^{n-1}(1+g_k(1))-g_k(n) \]

这也就意味着:

\[F_k(n)=-a^{n-1}(g_k(n)+n^k)+(g_k(1)+1) \]

定义\(G_k(n)=-2(g_k(n)+n^k)\),我们就有:

\[F_k(n)=a^nG_k(n)-aG_k(1) \]

又可以通过\(g_k(n)\)的通项公式推导得出:

\[G_k(1)=2G_k(0) \]

于是就可以得到

\[F_k(n)=a^nG_k(n)-G_k(0) \]

根据某些奇奇怪怪的东西我们可以知道\(G_k(n)\)是关于\(n\)的不超过\(k\)次的多项式。

i.e.作者自己也不会证明

这就意味着,我们只需要知道了\(G_k(n)\)\(k+1\)个点值,就可以为所欲为使用拉格朗日插值法计算出\(G_k\)在任何一个位置的点值。

另外,我们不难写出如下式子:

\[G_k(n)=2^nG_k(0)+2^nF_k(n) \]

这意味着,我们只需要计算出\(F_k(n)\),我们就可以得到\(G_k(n)\)。鉴于 \(F_k(n)\)有其简洁的递推形式,我们可以快速而方便地得到\(G_k(1)\dots G_k(k+1)\)\(G_k(0)\)的关系。

考虑到如下的高阶差分公式:

\[\Delta^kf(k)=\sum_{i=0}^k\binom{k}{i}(-1)^{k-i}f(k+i) \]

推导难度不大,自行定义算子\(\mathrm Rf(n)=f(n+1)\)并重定义差分,展开式子即可。

不过有一个很有用的性质——任何\(k\)次的多项式在经过\(k+1\)次的差分后,结果必然是\(0\)

这就是很像是求导了(差分本身就可以理解为离散中的导)。于是应有:

\[\Delta^{k+1}G_k(0)=\sum_{i=0}^{k+1}\binom{k+1}{i}(-1)^{k+1-i}G_k(i)=0 \]

根据这个式子我们可以列出\(G_k(1)\dots G_k(k+1)\)关于\(G_k(0)\)的关系。而\(G_k(1)\dots G_k(k+1)\)本身可以使用\(G_k(0)\)表示,因而我们就相当于得到了一个关于\(G_k(0)\)的一次方程,直接解出\(G_k(0)\)即可。

另外,还有一种列方程的方法。我们直接利用拉格朗日插值,可以得出:

\[G_k(0)=\sum_{i=1}^{k+1}G_k(i)\prod_{j\not=i}\frac{-j}{(i-j)} \]

这同样是\(G_k(1)\dots G_k(k+1)\)\(G_k(0)\)的关系,解方程就好。

进而我们可以得到\(G_k(1)\dots G_k(k+1)\),然后就可以插值得到\(G_k(n)\),最后计算出\(F_k(n)\)\(f_k(n)\)

时间复杂度\(O(k)\sim O(k\log_2m)\),其中\(m\)为模数。

感觉......本题的核心在其数学构造。 60 pts ~ 80 pts 中出现的错位相减是一个常见的推前缀和通项的方法。但是考场上好像把减的对象搞错了。 100 pts 中的构造很巧妙。对于\(G_k(n)\)的多项式性质似乎有点 “ 猜想 ” 的意味,直接证明好像有难度,如果可以证明可以评论,谢谢。

代码

#include <cstdio>

const int mod = 1e9 + 7, phi = mod - 1;
const int MAXK = 1e6 + 5;

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' );
}

int inv( const int );
int qkpow( int, int );

struct func
{
	int a, b;
	func() { a = b = 0; }
	func( const int C ) { a = 0, b = C; }
	func( const int A, const int B ) { a = A, b = B; }
	func operator + ( const func &g ) const { return func( ( a + g.a ) % mod, ( b + g.b ) % mod ); }
	func operator - ( const func &g ) const { return func( ( a - g.a + mod ) % mod, ( b - g.b + mod ) % mod ); }
	func operator * ( const func &g ) const { return func( ( 1ll * a * g.b % mod + 1ll * b * g.a % mod ) % mod, 1ll * b * g.b % mod ); }
	func operator / ( const int &g ) const { int t = inv( g ); return func( 1ll * a * t % mod, 1ll * b * t % mod ); }
	void operator += ( const func &g ) { *this = *this + g; }
	void operator -= ( const func &g ) { *this = *this - g; }
	void operator *= ( const func &g ) { *this = *this * g; }
	void operator /= ( const int &g ) { *this = *this / g; }
	
	int getX( const int y ) const { return 1ll * ( y - b + mod ) * inv( a ) % mod; }
	int getY( const int x ) const { return ( 1ll * a * x % mod + b ) % mod; }
};

func coe[MAXK];
int G[MAXK], neg[MAXK];
int fac[MAXK], finv[MAXK];
int prime[MAXK], pn, pwk[MAXK];
char S[MAXK];
int N, K;
bool isPrime[MAXK];

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

int C( const int n, const int m ) { return 1ll * fac[n] * finv[m] % mod * finv[n - m] % mod; }
int inv( const int n ) { return n <= K + 10 ? 1ll * finv[n] * fac[n - 1] % mod : qkpow( n, mod - 2 ); }

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

void init( const int siz )
{
	fac[0] = 1;
 	EulerSieve( siz );
	for( int i = 1 ; i <= siz ; i ++ ) fac[i] = 1ll * fac[i - 1] * i % mod;
	finv[siz] = qkpow( fac[siz], mod - 2 );
	for( int i = siz - 1 ; ~ i ; i -- ) finv[i] = 1ll * finv[i + 1] * ( i + 1 ) % mod;
}

int Lagrange( const int n )
{
	#define M K + 1
	if( n <= M ) return G[n];
	neg[0] = 1; int up = 1, ans = 0;
	for( int i = 2 ; i <= M ; i ++ ) up = 1ll * up * ( n - i ) % mod;
	for( int i = 1 ; i <= M ; i ++ ) neg[i] = 1ll * neg[i - 1] * inv( mod - i ) % mod;
	for( int i = 1 ; i <= M ; i ++ )
	{
		ans = ( ans + 1ll * up * finv[i - 1] % mod * neg[M - i] % mod * G[i] % mod ) % mod;
		if( i < M ) up = 1ll * up * ( n - i ) % mod * inv( n - i - 1 ) % mod;
	}
	return ans;
	#undef M
}

int main()
{
	int ind = 0;
	scanf( "%s %d", S, &K );
	for( int i = 0 ; S[i] ; i ++ )
		N = ( 10ll * N + S[i] - '0' ) % mod,
		ind = ( 10ll * ind + S[i] - '0' ) % phi;
	init( K + 10 ), coe[0] = func( 1, 0 );
	int F = 0, alph = inv( 2 ), pwa = 1, pw2 = 2;
	for( int i = 1 ; i <= K + 1 ; i ++ )
	{
		F = ( F + 1ll * pwa * pwk[i - 1] % mod ) % mod;
		coe[i] = func( pw2, 1ll * pw2 * F % mod );
		pwa = 1ll * pwa * alph % mod, pw2 = 2ll * pw2 % mod;
	}
	func tmp;
	for( int i = 0 ; i <= K + 1 ; i ++ )
	{
		if( ( K + 1 - i ) & 1 ) tmp -= coe[i] * C( K + 1, i ); 
		else tmp += coe[i] * C( K + 1, i );
	}
	G[0] = tmp.getX( 0 );
	for( int i = 1 ; i <= K + 1 ; i ++ )
		G[i] = coe[i].getY( G[0] );
	int GN = Lagrange( N );
	int FN = ( 1ll * qkpow( 2, phi - ind ) * GN % mod - G[0] + mod ) % mod;
	write( ( qkpow( N, K ) + 1ll * FN * qkpow( 2, ind - 1 + phi ) % mod ) % mod ), putchar( '\n' );
	return 0;
}
posted @ 2020-07-21 16:06  crashed  阅读(259)  评论(0编辑  收藏  举报