Description

给定一个正整数 \(n,k\)，求不定方程 \(x^2 - n \cdot y^2 = 1\) 的解，其中 \(x,y\) 为正整数。

由于方程解很多，所以你的解需要使 \(x\) 在所有解中为第 \(k\) 小的。

由于这个解可能很大，所以请输出你得到的解中的 \(x \bmod 8191\) 的值。

如果无解，输出 No answers can meet such conditions 。

Hint

\(2\le n \le 29\)

\(1\le k \le 10^9\)

Solution

首先，这个方程其实有一个名字，叫做 “佩尔（Pell）方程” 。

很容易判断无解的情况，即当 \(\sqrt{n} \in \text{N}^+\)（\(n\) 为完全平方数时）是无解的。

然后确定这个方程是有解的之后，我们先找出其最小解 \((x_1,y_1)\)。

怎么找？数据范围不大，暴力即可。

之后的解，存在一个迭代公式：

下面给出证明。

证明：

\(x_i^2 - ny_i^2\)

\(= (x_{i-1} x_1 + ny_{i-1} y_1)^2 - n(x_{i-1} y_1 + y_{i-1} x_1)^2\)

\(= (x_{i-1}^2 x_1^2 + n^2 y_{i-1}^2 y_1^2 + 2n x_{i-1} x_1 y_{i-1} y_1) - (nx_{i-1}^2 y_1^2 + ny_{i-1}^2 x_1^2 +2n x_{i-1} x_1 y_{i-1} y_1)\)

\(= x_{i-1}^2 (x_1^2 - ny_1^2) + ny_{i-1}^2 (ny_1^2 - x_1^2)\)

\(\because (x_1,y_1)\) 是方程的一组解，

\(\therefore x_1^2 - ny_1^2 = 1\)

\(\therefore x_i^2 - ny_i^2 = x_{i-1}^2 (x_1^2 - ny_1^2) + ny_{i-1}^2 (ny_1^2 - x_1^2) =x_{i-1}^2 +ny_{i-1}^2\)

证毕。

那么接下来，我们就可以由一个最小的正整数解，递推推出第 \(k\) 个解。

然而本题中， \(\Theta(k)\) 复杂度的算法是过不了的。

但是前面的辛苦并没有白费，我们惊喜地发现，这个递推式子可以写成矩阵的形式：

使用快速幂优化，递推复杂度直降至 \(\Theta(\log k)\) ，可以通过此题。

Code

#include<iostream>
#include<cmath>
#include<cstring>

using std::cin;
using std::cout;
using std::endl;
const int mod=8191;

struct matrix {
	int element[2][2];
	int r,c;
	inline matrix(int _r,int _c):r(_r),c(_c) {
		memset(element,0,sizeof element);
	}
	inline int* operator [](const int &p) {
		return element[p];
	}
};
inline matrix unit_matrix(int s) {
	matrix ret(s,s);
	for(register int i=0;i<s;i++)
		ret[i][i]=1;
	return ret;
}
inline matrix operator *(matrix A,matrix B) {
	matrix C(A.r,B.c);
	for(register int i=0;i<A.r;i++)	
		for(register int j=0;j<B.c;j++)
			for(register int k=0;k<A.c;k++)
				(C[i][j]+=A[i][k]*B[k][j]%mod)%=mod;
	return C;
}
matrix fastpow(matrix a,int b) {
	if(!b) return unit_matrix(a.r);
	matrix t=fastpow(a,b>>1);
	if(b&1) return t*t*a;
	else return t*t;
}

signed main() {
	int N,K;
	while(cin>>N>>K) {
		int SQRTN=floor(sqrt(N*1.0));
		if(SQRTN*SQRTN==N) {
			cout<<"No answers can meet such conditions"<<endl;
			continue;
		}
		int y=1,x;
		for(;;y++) {
			x=floor(sqrt(N*y*y+1));
			if(x*x-N*y*y==1) break;
		}
		matrix base(2,2);
		base[0][0]=x;
		base[0][1]=N*y;
		base[1][0]=y;
		base[1][1]=x;
		
		matrix sol(2,1);
		sol[0][0]=x;
		sol[1][0]=y;
		
		sol=fastpow(base,K-1)*sol;
		cout<<sol[0][0]<<endl;
	}
}