大意:求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;
}