HDU3292 佩尔pell方程、矩阵快速幂

Posted on | In , , ,

给出n，k，满足$x^2-ny^2=1$，若存在输出第$k$的$x\mod8191$。
$(n<29，k\le10^9)$

题目链接

题解

  • 佩尔pell方程裸题，至需要套套板子即可

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
#include<bits/stdc++.h>
using namespace std;
#define ll long long

const int maxn = 35;
const int mod = 8191;

int  n, k;
ll x1, y1;


int p;
struct martix
{
    int ma[4][4];
    martix friend operator *(const martix a, const martix b){
        martix ret;
        memset(ret.ma,0,sizeof(ret.ma));
        for(int i=0;i<p;i++){
            for(int j=0;j<p;j++){
                for(int k=0;k<p;k++)
                    ret.ma[i][j]=(ret.ma[i][j] + a.ma[i][k]*b.ma[k][j]%mod)%mod;
            }
        }
        return ret;
    }
}ans;

martix multipow(martix x, int k)
{
    memset(ans.ma, 0, sizeof(ans.ma));
    for(int i=0;i<p;i++)
        ans.ma[i][i]=1;
    while(k)
    {
        if(k&1) ans=ans*x;
        k>>=1;
        x=x*x;
    }
    return ans;
}

void seek(ll &x, ll &y, ll d)
{
    y=1;
    while(1){
        x=1LL*sqrt(d*y*y+1);
        if(x*x-d*y*y==1) break;
        y++;
     //   cout<<y<<endl;
    }
}

bool check(int x)
{
    int y=(int)sqrt(x);
    if(y*y==x) return true;
    return false;
}

int main()
{
    p=2;
    while(scanf("%d%d", &n, &k)!=EOF)
    {
        if(check(n)) {
            printf("No answers can meet such conditions\n");
            continue;
        }
        seek(x1, y1, n);
        //cout<<x1<<' '<<y1<<endl;
        martix now;
        now.ma[0][0]=x1;
        now.ma[0][1]=(n)*(y1);
        now.ma[1][0]=y1;
        now.ma[1][1]=x1;
        martix ans=multipow(now, k-1);
        int sum=0;
            sum = (sum + ans.ma[0][0]*x1%mod + ans.ma[0][1]*y1%mod) % mod;
        printf("%d\n", sum);
    }
}