[LOJ6055]「from CommonAnts」一道数学题 加强版
题目
点这里看题目。
分析
首先可以娴熟地推倒一发式子:
定义\(S_k(n)=\sum_{i=1}^{n}f_k(i)\),做差计算可以得到:
自然我们可以得到:
为了方便,我们就可以定义出:
问题于是转化为了:如何快速求出\(F_k(n)\)?
60 pts ~ 80 pts
可以发现\(F_k(n)\)就是一个前缀和,我们尝试错位相减:
并且当\(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)\)的不存在求和符的递推式:
尝试构造出一个\(g_k(n)\)满足:
显然一个可行的解为:
即可以解出:
这也就意味着:
定义\(G_k(n)=-2(g_k(n)+n^k)\),我们就有:
又可以通过\(g_k(n)\)的通项公式推导得出:
于是就可以得到
根据某些奇奇怪怪的东西我们可以知道\(G_k(n)\)是关于\(n\)的不超过\(k\)次的多项式。
i.e.作者自己也不会证明
这就意味着,我们只需要知道了\(G_k(n)\)的\(k+1\)个点值,就可以为所欲为使用拉格朗日插值法计算出\(G_k\)在任何一个位置的点值。
另外,我们不难写出如下式子:
这意味着,我们只需要计算出\(F_k(n)\),我们就可以得到\(G_k(n)\)。鉴于 \(F_k(n)\)有其简洁的递推形式,我们可以快速而方便地得到\(G_k(1)\dots G_k(k+1)\)与\(G_k(0)\)的关系。
考虑到如下的高阶差分公式:
推导难度不大,自行定义算子\(\mathrm Rf(n)=f(n+1)\)并重定义差分,展开式子即可。
不过有一个很有用的性质——任何\(k\)次的多项式在经过\(k+1\)次的差分后,结果必然是\(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(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;
}