HDU 3292 Pell方程第k小解

HDU 3292 Pell方程第k小解

链接：HDU3292

方法：Pell方程递推式 + 矩阵快速幂

AC code：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod = 8191;
const int NUM = 2;
int N = 2;
struct mat{
    int a[NUM][NUM];
    void init(){
        memset(a, 0, sizeof(a));
        for(int i = 0; i < N; ++i) a[i][i] = 1;
    }
    mat operator * (const mat &b) const{
        mat c;
        for(int i = 0; i < N; ++i){
            for(int j = 0; j < N; ++j){
                c.a[i][j] = 0;
                for(int k = 0; k < N; ++k){
                    c.a[i][j] = (c.a[i][j] + (long long)a[i][k] * b.a[k][j] % mod) % mod;
                }
            }
        }
        return c;
    }
};
mat fast_pow(mat a, int b){
    mat res;
    res.init();
    while(b > 0){
        if(b & 1){
            res = res * a;
        }
        a = a * a;
        b >>= 1;
    }
    return res;
}

/*

[ xn+1 ] = [x0  d * y0] [ xn ]
[ yn+1 ] = [y0  x0]     [ yn ]
*/
inline bool issq(int d){
    int sq = floor(sqrt(d + 0.5));
    return (sq * sq == d);
}
void Find(int &x, int &y, int d){
    y = 1;
    while(1){
        x = floor(sqrt((long long)d * y * y + 1.5));
        if((ll)x * x == (long long)d * y * y + 1LL) break;
        ++y;
    }
}
int main(){
    int d, k;
    while(cin >> d >> k){
        if(issq(d)){
            puts("No answers can meet such conditions");
            continue ;
        }
        int x, y;
        Find(x, y, d);
        if(k == 1){
            printf("%d\n", x);
            continue ;
        }
        mat A;
        A.a[0][0] = x % mod, A.a[0][1] = (ll)d * y % mod;
        A.a[1][0] = y % mod, A.a[1][1] = x % mod;
        A = fast_pow(A, k - 1);
        int ansx = ((ll)A.a[0][0] * x  % mod + (ll)A.a[0][1]  * y % mod) % mod;
        printf("%d\n", ansx);
    }
    return 0;
}