No more tricks, Mr Nanguo HDU - 3292 【佩尔方程矩阵快速幂求第K大解】

最新推荐文章于 2020-02-11 23:54:39 发布

Conchpeng

最新推荐文章于 2020-02-11 23:54:39 发布

阅读量395

点赞数

CC 4.0 BY-SA版权

分类专栏：数论之佩尔方程

本文链接：https://blog.youkuaiyun.com/Haipai1998/article/details/81610851

数论之佩尔方程专栏收录该内容

7 篇文章

订阅专栏

本文介绍了一种利用矩阵快速幂解决佩尔方程第K大解的方法，并通过一个具体的HDU-3292题目进行实现。文章详细展示了矩阵乘法、矩阵快速幂等关键步骤的代码实现。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

No more tricks, Mr Nanguo HDU - 3292

题意：输入D,K ，求 $x^2 -D*y^2==1$ 的第K大X解

用佩尔方程的矩阵表示，马上会总结到（佩尔方程专题总结）

用JAVA写失败了。不知道问题出在哪里

#include<cstdio>
#include<vector>
#include<cmath>
#include<math.h>
#include<time.h>
#include<string>
#include<string.h>
#include<iostream>
#include<algorithm>
#include<map>
#define PI acos(-1.0)
#define pb push_back
#define F first
#define S second
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
const int N=1005;
const int MOD=8191;
const double eps=1e-10;

typedef struct{
    ll mat[20][20];
}matrix;

// ¾ØÕó³Ë·¨
matrix matrix_mul(matrix a,matrix b,ll y){
    matrix c;
    memset(c.mat,0,sizeof(c.mat));
    for(ll i=0;i<=y;i++){
        for(ll j=0;j<=y;j++){
            for(ll k=0;k<=y;k++){
                c.mat[i][j]+=(a.mat[i][k]%MOD*b.mat[k][j]%MOD)%MOD;
                c.mat[i][j]%=MOD;
            }
        }
    }
    return c;
}

// ¿ìËÙÃÝ
matrix matrix_quick_power(matrix a,ll k ,ll y){
    matrix b;
    memset(b.mat,0,sizeof(b.mat));
    for(int i=0;i<=y;i++)
        b.mat[i][i]=1ll;
    while(k){
        if(k%2==1)  b=matrix_mul(b,a,y);
        k=k/2;
        a=matrix_mul(a,a,y);
    }
    return b;

}
ll d,k;
void work(ll &x1,ll &y1){
    y1=1;
    while(true){
        ll x_2=1+d*y1*y1;
        if(x_2>0){
            x1=sqrt(x_2);
            if(x1*x1-d*y1*y1==1)    return ;
        }
        ++y1;
    }
}
bool is_square(ll x){
    ll t=sqrt(x);
    for(ll i=t-1;i<=t+1;i++)
        if(i*i==x)  return true;
    return false;
}
int main(void){
    while(scanf("%lld%lld",&d,&k)==2){
        matrix a,res;
        ll x1,y1;
        if(is_square(d)){
            printf("No answers can meet such conditions\n");
            continue;
        }
        work(x1,y1);
        a.mat[1][1]=x1;
        a.mat[2][1]=y1;
        a.mat[1][2]=d*y1;
        a.mat[2][2]=x1;
        res=matrix_quick_power(a,k-1,5);
        ll ansx;
        ansx=res.mat[1][1]*x1%MOD+res.mat[1][2]*y1%MOD;
        ansx%=MOD;
        printf("%lld\n",ansx);
    }

    return 0;
}