「BZOJ4833」最小公倍佩尔数
题目
对于非负整数 \(n\)，令 \((1+\sqrt 2)^n=e_n+f_n\sqrt 2\)。
设 \(g_n=\operatorname{lcm}_{1\le k\le n}f_k\)。给定正整数 \(n\) 和质数 \(p\)，求 \(\sum_{k=1}^n kg_k\bmod p\)。
多组测试，所有测试点满足 \(1\le n\le 10^6,1\le \sum n\le 3\times 10^6\)，且 \(\forall 1\le k\le n,p\nmid f_k\)。
分析
容易得到递推关系：
\[\begin{aligned}
\begin{pmatrix}
e_0\\
f_0
\end{pmatrix}
&=
\begin{bmatrix}
1\\
0
\end{bmatrix}\\
\begin{pmatrix}
e_n\\
f_n
\end{pmatrix}
&=
\begin{bmatrix}
1&2\\
1&1
\end{bmatrix}
\begin{pmatrix}
e_{n-1}\\
f_{n-1}
\end{pmatrix},\forall n\ge 1
\end{aligned}
\]
注意到答案最后和 \(e\) 没关系，所以我们尝试消掉 \(e\)，得到只关于 \(f\) 的式子：
\[\begin{aligned}
f_0&=0\\
f_{n}&=1+f_{n-1}+2\sum_{k=0}^{n-2}f_{k},\forall n\ge 1
\end{aligned}
\]
进而得到一个递推关系：
\[\begin{aligned}
f_0&=0,f_1=1\\
f_n&=2f_{n-1}+f_{n-2},\forall n\ge 2
\end{aligned}
\]
注意到这个 \(f\) 有两个很好的性质：
-
\(\gcd(f_n,f_{n-1})=1\)，考察 \(\gcd\) 时必看。
-
\(f_{n+m}=f_{m+1}f_n+f_{m}f_{n-1}\)，考察转移矩阵可以得到这一点。
由这两点，将其和 Fibonacci 数列进行类比，我们可以得到：
\[\gcd(f_{n},f_{m})=f_{\gcd(n,m)},\forall n,m\in \mathbb Z_+
\]
Proof.
不妨设 \(m=n+k,k\in \mathbb N\)。
则：
\[\begin{aligned} \gcd(f_n,f_m) &=\gcd(f_n,f_{n+k})\\ &=\gcd(f_n,f_{k+1}f_n+f_{k}f_{n-1})\\ &=\gcd(f_n,f_{n-1}f_k)\\ &=\gcd(f_n,f_k)=\gcd(f_n,f_{m-n}) \end{aligned} \]进而可以得到 \(\gcd(f_n,f_m)=\gcd(f_n,f_{m-\lfloor\frac m n\rfloor n})=\gcd(f_n,f_{m\bmod n})\)。反复利用此关系即可证明。
最后，\(g_n\) 可以直接根据 min-max 反演通过 \(f\) 的 \(\gcd\) 计算得到：
\[g_n=\prod_{d=1}^nf_d^{\sum_{k=1}^{\lfloor\frac n d\rfloor}\mu(k)}
\]
\(g_n\) 可以在 \(O(n\log n)\) 的时间内处理出来，于是问题就解决了。
代码
#include <cstdio>
#include <iostream>
#define rep( i, a, b ) for( int i = (a) ; i <= (b) ; i ++ )
#define per( i, a, b ) for( int i = (a) ; i >= (b) ; i -- )
typedef long long LL;
const int MAXN = 1e6 + 5;
template<typename _T>
inline void Read( _T &x ) {
x = 0; char s = getchar(); bool f = false;
while( ! ( '0' <= s && s <= '9' ) ) { f = s == '-', s = getchar(); }
while( '0' <= s && s <= '9' ) { x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar(); }
if( f ) x = -x;
}
template<typename _T>
inline void Write( _T x ) {
if( x < 0 ) putchar( '-' ), x = -x;
if( 9 < x ) Write( x / 10 );
putchar( x % 10 + '0' );
}
int tag[MAXN];
int f[MAXN], mu[MAXN];
int N, mod;
inline int Qkpow( int, LL );
inline int Inv( const int &a ) { return Qkpow( a, mod - 2 ); }
inline int Mul( int x, const int &v ) { return 1ll * x * v % mod; }
inline int Sub( int x, const int &v ) { return ( x -= v ) < 0 ? x + mod : x; }
inline int Add( int x, const int &v ) { return ( x += v ) >= mod ? x - mod : x; }
inline int& MulEq( int &x, const int &v ) { return x = 1ll * x * v % mod; }
inline int& SubEq( int &x, const int &v ) { return ( x -= v ) < 0 ? ( x += mod ) : x; }
inline int& AddEq( int &x, const int &v ) { return ( x += v ) >= mod ? ( x -= mod ) : x; }
inline int Qkpow( int base, LL indx ) {
int ret = 1;
while( indx ) {
if( indx & 1 ) MulEq( ret, base );
MulEq( base, base ), indx >>= 1;
}
return ret;
}
inline void Init( const int &n ) {
mu[1] = 1;
for( int i = 1 ; ( i << 1 ) <= n ; i ++ )
for( int j = i << 1 ; j <= n ; j += i )
mu[j] -= mu[i];
}
int main() {
Init( 1e6 );
int T; Read( T );
while( T -- ) {
Read( N ), Read( mod );
f[0] = 0, f[1] = 1;
rep( i, 1, N ) tag[i] = 1;
rep( i, 2, N ) f[i] = Add( Mul( 2, f[i - 1] ), f[i - 2] );
for( int d = 1 ; d <= N ; d ++ ) {
int fInv = Inv( f[d] );
for( int j = 1 ; d * j <= N ; j ++ )
if( mu[j] == 1 )
MulEq( tag[j * d], f[d] );
else if( mu[j] == -1 )
MulEq( tag[j * d], fInv );
}
int ans = 0;
rep( i, 2, N ) MulEq( tag[i], tag[i - 1] );
rep( i, 1, N ) AddEq( ans, Mul( i, tag[i] ) );
Write( ans ), putchar( '\n' );
}
return 0;
}
