佩尔方程+矩阵快速幂整理

佩尔方程最基本的形式(我目前了解到的)：

x² - d × y² = 1

首先找到一组最小正整数解：(x₁,y₁)

解的递推式为：

X_n  = X_n-1  × X₁ + d × Y_n-1 ×Y₁

Y_n  = X_n-1  × Y₁ + Y_n-1  × X₁

矩阵快速幂递推：

例题：佩尔方程+矩阵快速幂(HDU3292)

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
typedef long long ll;
const int mod = 8191;

struct matrix
{
    int m[2][2];
    matrix()
    {
        for(int i=0;i<2;i++)
            for(int j=0;j<2;j++)
                m[i][j] = 0;
        m[1][1] = 1;
        m[0][0] = 1;
    }
    matrix operator *(const matrix& a) const
    {
        matrix temp;
        for(int i=0;i<2;i++)
        {
            for(int j=0;j<2;j++)
            {
                temp.m[i][j] = 0;
                for(int k=0;k<2;k++)
                {
                    temp.m[i][j]+=m[i][k]*a.m[k][j]%mod;
                    temp.m[i][j]%=mod;
                }
            }
        }
        return temp;
    }
};

matrix ks(matrix a,int b)
{
    matrix ans = matrix();
    while(b)
    {
        if(b&1)
        {
            ans = ans*a;
        }
        a = a*a;
        b>>=1;
    }
    return ans;
}

//传入d后暴力求解
pair<int,int> find(int d)
{
    int x,y = 1;
    while(1)
    {
        ll tmp = (ll)y*y*d+1;
        int temp = sqrt(tmp);
        if((ll)temp*temp==tmp)
        {
            x = temp;
            break;
        }
        y++;
    }
    return make_pair(x,y);
}


int main()
{
    int n,k;
    while(scanf("%d %d",&n,&k)!=EOF)
    {
        int temp = sqrt(n*1.0);
        if(temp*temp == n)
        {
            printf("No answers can meet such conditions
");
            continue;
        }
        pair<int,int> base = find(n);
        int x1 = base.first,y1 = base.second;
        //printf("%d %d
",base.first,base.second);
        matrix d;
        d.m[0][0] = x1;
        d.m[0][1] = n * y1;
        d.m[1][0] = y1;
        d.m[1][1] = x1;
        matrix ans = ks(d,k-1);
        int ansx = (ans.m[0][0]*x1%mod+ans.m[0][1]*y1%mod)%mod;
        int ansy = (ans.m[1][0]*x1%mod+ans.m[1][1]*y1%mod)%mod;
        printf("%d
",ansx);
    }
}