Hdu 3292 No more tricks, Mr Nanguo

来源：互联网发布：功夫软件手机版编辑：程序博客网时间：2024/06/03 16:44

大意：求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-10const 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;}