「WC2021」斐波那契

题目

点这里看题目。

分析

有点厉害啊......

首先注意到，由于递推关系和询问参数 \(a,b\) 没有关系，我们可以使用斐波那契数列 \(f\) 来给出 \(F\) 的“通项”。换言之，也就是：

\[F_n=af_n+(b-a)f_{n-1} \]

Note.

简单的理解方法：考虑转移矩阵的幂次 \(\begin{bmatrix}1&1\\1&0\end{bmatrix}^n=\begin{bmatrix}f_n&f_{n-1}\\ f_{n-1}&f_{n-2}\end{bmatrix}\)。

这里，我们允许 \(f\) 出现负下标。负下标的情况可以由递推式给出：

\[f_{n-2}=f_{n}-f_{n-1},n\le 1 \]

比如，有 \(f_{-1}=0,f_{-2}=1\)。

去掉一些容易处理的边界，我们现在考虑 \(a,b,b-a\not\equiv 0\pmod m\) 的情况。所以，我们只需要解决：

\[uf_n\equiv vf_{n-1}\pmod m \]

先做一点显然的处理：对于 \(u,v,m\) 同时除掉 \(\gcd(u,v,m)\)，得到 \(u',v',m'\)。

如果有着良好的互质性，比如 \(m'\) 是一个质数，那么接下来我们可以直接将 \(f\) 和 \(u,v\) 分开，从而做一些预处理。不过，虽然 \(\gcd(u',v',m')=1\)，我们却还无法保证 \(u'\perp m'\) 或者 \(v'\perp m'\)，因此还需要再做尝试。

很容易想到把 \(\gcd(u',m')\) 和 \(\gcd(v',m')\) 给除掉。但是怎么保证除了之后不会出现分数？结合我们的目的和斐波那契自身的性质 \(f_{n-1}\perp f_n\)，我们不难猜到这样一个结论：

\[\begin{aligned} \gcd(u',m')&=\gcd(f_{n-1},m')=d_2\\ \gcd(v',m')&=\gcd(f_n,m')=d_1 \end{aligned} \]

为了节约篇幅，证明缩在了这里：

Proof.

刚拿到这个条件的时候，我们注意到：同余条件可以产生 \(\gcd\)，而反过来比较困难。因此，我们选择将仅有的同余式转化成 \(\gcd\) 的关系：

\[\gcd(u'f_n,m')=\gcd(v'f_{n-1},m')=a \]

把 \(a\) 除到 \(\gcd\) 里面去，并考虑 \(a\) 是如何“分配”而来的。也就是说：

\[\begin{aligned} &\gcd\left(\frac{u'd_1}{a}\cdot \frac{f_n}{d_1},\frac{m'}{a}\right)\\ =&\gcd\left(\frac{v'd_2}{a}\cdot \frac{f_{n-1}}{d_2},\frac{m'}{a}\right)\\ =&1 \end{aligned} \]

这将会给出：

\[\gcd(d_1u',d_2v',m')=a \]

由于 \(f_{n-1}\perp f_n\)，显然有 \(d_1\perp d_2\)。结合 \(\gcd(u',v',m')=1\)，我们可以知道 \(d_1d_2=a\)。此时我们已经完成了 \(d_2|\gcd(u',m')\) 和 \(d_1|\gcd(v',m')\)。

接下来证明 \(\gcd(u',m')=d_2\)。这一点只需要结合 \(\gcd(u'f_n,m')=d_1d_2,d_1\perp d_2\) 即可得到。证明 \(\gcd(v',m')=d_1\) 类似。

所以，这样处理过后我们又将 \(u,v\) 和 \(f\) 分开了。对于每个 \(m'|m\) 我们可以处理出所有可能的三元组 \((d_1,d_2,(f_nd_1^{-1})(f_{n-1}d_2^{-1})^{-1}\bmod (m'd_1^{-1}d_2^{-1}))\)。由于斐波那契的循环节长度是 \(O(m')\) 的，因此对于每个质因子暴力枚举即可。询问直接查一下就结束了。

时间复杂度为 \(O(\sigma (m)\log m+n\log m)\)。

代码

由于我偷懒没有精细实现，这份代码在洛谷上会 TLE，在 LOJ 上才能过

#include <map>
#include <set>
#include <cstdio>
#include <vector>
#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 -- )

const int MAXN = 1e6 + 5, MAXD = 2e3 + 5;

template<typename _T>
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>
void write( _T x ) {
	if( x < 0 ) putchar( '-' ), x = -x;
	if( 9 < x ) write( x / 10 );
	putchar( x % 10 + '0' );
}

struct Triple {
	int a, b, c;

	Triple(): a( -1 ), b( -1 ), c( -1 ) {}
	Triple( int A, int B, int C ): a( A ), b( B ), c( C ) {}

	inline bool operator < ( const Triple &q ) const {
		return a == q.a ? ( b == q.b ? c < q.c : b < q.b ) : a < q.a;
	}
};

std :: set<std :: pair<int, int> > ever;
std :: map<Triple, int> occ[MAXD];

int dvs[MAXD], tot = 0;

int N, M;

inline int Gcd( int x, int y ) { for( int z ; y ; z = x, x = y, y = z % y ); return x; }

inline int Exgcd( const int &a, const int &b, int &x, int &y ) {
	if( ! b ) return x = 1, y = 0, a;
	int d = Exgcd( b, a % b, y, x );
	y -= x * ( a / b ); return d;
}

inline int Inv( const int &a, const int &mod ) {
	static int x, y;
	Exgcd( a, mod, x, y );
	return ( x % mod + mod ) % mod;
}

void Enumerate( const int &mod ) {
	ever.clear();
	dvs[++ tot] = mod;
	int x = 1, y = 0, p, q;
	for( int n = 0 ; ; n ++ ) {
		if( ever.find( { x, y } ) != ever.end() ) break;
		ever.insert( { x, y } );
		p = Gcd( x, mod ), q = Gcd( y, mod );
		Triple nxt( p, q, 1ll * ( y / q ) * Inv( x / p, mod / p / q ) % ( mod / p / q ) );
		if( occ[tot].find( nxt ) == occ[tot].end() )
			occ[tot][nxt] = n;
		y = ( x + y ) % mod, std :: swap( x, y );
	}
}

int main() {
	read( N ), read( M );
	rep( i, 1, M )
		if( M % i == 0 )
			Enumerate( i );
	while( N -- ) {
		int a, b, m = M;
		read( a ), read( b );
		if( a == 0 ) { 
			puts( "0" ); 
			continue; 
		}
		if( b == 0 ) {
			puts( "1" );
			continue;
		}
		b = ( a - b + m ) % m;
		int d = Gcd( Gcd( a, b ), m );
		a /= d, b /= d, m /= d;
		int d1 = Gcd( a, m ), d2 = Gcd( b, m );
		int idx = std :: lower_bound( dvs + 1, dvs + 1 + tot, m ) - dvs;
		Triple qry( d2, d1, 1ll * ( a / d1 ) * Inv( b / d2, m / d1 / d2 ) % ( m / d1 / d2 ) );
		if( occ[idx].find( qry ) == occ[idx].end() ) 
			puts( "-1" );
		else
			write( occ[idx][qry] ), putchar( '\n' );
	}
	return 0;
}