「HDU 3292」 No more tricks, Mr Nanguo

Description

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

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

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

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

Hint

(2le n le 29)

(1le k le 10^9)

Solution

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

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

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

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

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

[x_i = x_{i-1} cdot x_1 + ny_{i-1} cdot y_1 ]

[y_i = x_{i-1} cdot y_1 + y_{i-1} cdot x_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))

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

( herefore x_1^2 - ny_1^2 = 1)

( herefore 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)) 复杂度的算法是过不了的。

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

[egin{bmatrix} x_k \ y_k end{bmatrix} = egin{bmatrix} x_1 & ny_1 \ y_1 & x_1 end{bmatrix}^{k-1} imes egin{bmatrix} x_1 \ y_1 end{bmatrix} ]

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

Code

/*
 * Author : _Wallace_
 * Source : https://www.cnblogs.com/-Wallace-/
 * Problem : HDU 3292 No more tricks, Mr Nanguo
 */
#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;
	}
}