393
HDU 3292 pell方程
題意：求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;
}
最後更新：2017-04-03 16:48:46