题意:求x^2+n*y^2=1 按x排序第k大的解。
题目就是求佩尔方程,用矩阵连乘的方法,当n为完全平方数的时候无解。
如果我们求出Pell方程的最小正整数解后,就可以根据递推式求出所有的解。
则根据上式我们可以构造矩阵,然后就可以快速幂了。
这样就可以求出第k大的解。
HDU3292题就要用到上面的矩阵方法求第k大的解。
#include <iostream> #include<cstdio> #include<cstring> #include<algorithm> #include<cmath> using namespace std; #define M 8191 const int MAX=2; typedef struct { long long m[MAX][MAX]; } Matrix; Matrix I= {1,0, 0,1 }; Matrix P; Matrix matrixmul(Matrix a,Matrix b) //矩阵乘法 { int i,j,k; Matrix c; for (i = 0 ; i < MAX; i++) for (j = 0; j < MAX; j++) { c.m[i][j] = 0; for (k=0; k<MAX; k++) c.m[i][j]+=((a.m[i][k]%(M))*(b.m[k][j]%(M)))%(M); c.m[i][j] %=M; } return c; } Matrix quickpow(long long n) { Matrix m = P, b = I; while (n >= 1) { if (n & 1) b = matrixmul(b,m); n = n >> 1; m = matrixmul(m,m); } return b; } int main() { long long n,k,x,y; while(~scanf("%I64d%I64d",&n,&k)) { if(n==4||n==9||n==16||n==25) { puts("No answers can meet such conditions"); continue; } for(y=1;; y++) { x=(long long)sqrt(double(y*y*n+1)); if(x*x==y*y*n+1) break; } P.m[0][0]=P.m[1][1]=x%M; P.m[0][1]=(n*y)%M,P.m[1][0]=y%M; Matrix ans=quickpow(k-1); printf("%I64d\n",((ans.m[0][0]*x)%M+(ans.m[0][1]*y)%M)%M); } return 0; }