Hdu 3292 No more tricks, Mr Nanguo

大意：求x^2-d*y^2 = 1的第K解。

思路：

佩尔方程的求解，我们知道佩尔方程有解的情况是d不是完全平方数，根据的佩尔方程的递推公式那么第K解就是：

xn     =  [x1 dy1] ^n-1     x1

yn    =  [y1 x1   ]            y1

我们知道了d、k，用矩阵快速幂算出xn-1,d*yn-1，然后分别于与x1和y1相乘即可。

#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;

#define eps 1e-10

const int MAXN = 31;
const int BASE = 8191;
const int n = 2;

typedef struct
{
	int A[MAXN][MAXN];
}M;

M a, per;

int d, k;

void init(int x0, int y0)
{
	a.A[0][0] = x0%BASE, a.A[0][1] = d*y0%BASE;
	a.A[1][0] = y0%BASE, a.A[1][1] = x0%BASE;
	for(int i = 0; i < n; i++)
	{
		for(int j = 0; j < n; j++)
		{
			per.A[i][j] = (i == j);
		}
	}
}

void search(int &x0, int &y0, int d)
{
	int x, y = 1;
	for(;;)
	{
		x = (int)sqrt(d*y*y+1.0);
		if(x*x == d*y*y+1) break;
		y++;
	}
	x0 = x, y0 = y;
}

M Mul(M a, M b)
{
	M c;
	memset(c.A, 0, sizeof(c.A));
	for(int i = 0; i < n; i++)
	{
		for(int j = 0; j < n; j++)
		{
			for(int k = 0; k < n; k++)
			{
				c.A[i][j] += a.A[i][k] * b.A[k][j];
			}
			c.A[i][j] %= BASE;
		}
	}
	return c;
}

M power(int k)
{
	M p, ans = per;
	p = a;
	while(k)
	{
		if(k & 1) ans = Mul(ans, p);
		p = Mul(p, p);
		k /= 2;
	}
	return ans;
}

int check(int n)
{
	for(int i = 1; n > 0; i += 2) n -= i; // 1 + 3 + 5 + 7 + .... + (2*n-1) = n^2;
	return n == 0;
}

void solve()
{
	int x0, y0;
	if(check(d)) { printf("No answers can meet such conditions\n"); return ;}
	search(x0, y0, d);
	init(x0, y0);
	M ans = power(k-1);
	int q = (ans.A[0][0]*x0%BASE + ans.A[0][1]*y0%BASE)%BASE;
	printf("%d\n", q);
}

int main()
{
	while(~scanf("%d%d", &d, &k))
	{
		solve();
	}
	return 0;
}