HDU 3292 (佩尔方程矩阵快速幂)

最新推荐文章于 2021-10-31 12:43:38 发布

morejarphone

最新推荐文章于 2021-10-31 12:43:38 发布

阅读量390

点赞数

CC 4.0 BY-SA版权

分类专栏：数论/推导

本文链接：https://blog.youkuaiyun.com/morejarphone/article/details/52945734

数论/推导专栏收录该内容

26 篇文章

订阅专栏

本文介绍了一种求解佩尔方程x²-ny²=1第k大解的方法。通过找到最小特解，利用矩阵快速幂进行迭代计算。文章提供了一个C++实现示例。

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

题目链接：点击这里

题意：求解佩尔方程 $x^2-ny^2=1$ 的第k大解。

首先暴力求出佩尔方程的最小特解，然后根据迭代式子

[x k y k] = [x 1 y 1 d \times y 1 x 1] k - 1 [x 1 y 1]

$\begin{bmatrix} x_k \\ y_k\\ \end{bmatrix} =\begin{bmatrix} x_1 & d\times y_1\\ y_1 & x_1 \\ \end{bmatrix} ^{k-1} \begin{bmatrix} x_1\\ y_1\\ \end{bmatrix}$
直接用矩阵快速幂求出结果即可。

#include <bits/stdc++.h>
using namespace std;
#define maxn 500005
#define INF 1e15
#define mod 8191

struct m {
    long long a[2][2];
    m operator * (const m &gg) const {
        m ans;
        memset (ans.a, 0, sizeof ans.a);
        for (int i = 0; i < 2; i++) {
            for (int j = 0; j < 2; j++) {
                for (int l = 0; l < 2; l++) {
                    ans.a[i][j] += a[i][l]*gg.a[l][j];
                    ans.a[i][j] %= mod;
                }
            }
        }
        return ans;
    }
    void show () {
        for (int i = 0; i < 2; i++) {
            for (int j = 0; j < 2; j++)
                cout << a[i][j] << " ";
            cout << endl;
        }
    }
};

m qpow (m a, long long kk) {
    m ans;
    int i, j;
    for(i = 0; i < 2; ++i)
        for(j = 0; j < 2; ++j)
            ans.a[i][j] = (i == j ? 1 : 0);
    for(; kk; kk >>= 1) {
        if(kk&1) ans = ans*a;
        a = a*a;
    }
    return ans;
}

int n, k;

void solve (long long x, long long y, long long k) {
    m ans;
    ans.a[0][0] = x, ans.a[0][1] = n*y;
    ans.a[1][0] = y, ans.a[1][1] = x;
    ans = qpow (ans, k-1);
    cout << (x*ans.a[0][0]+y*ans.a[0][1]) % mod << endl;
}

int main () {
    while (cin >> n >> k) {
        if (n == 1 || n == 4 || n == 9 || n == 16 || n == 25) {
            cout << "No answers can meet such conditions" << endl;
            continue;
        }
        for (long long y = 1; ; y++) {
            long long x = sqrt (1LL*n*y*y+1);
            if (x*x-y*y*n == 1) { //cout << x << " " << y << endl;
                solve (x, y, k);
                goto out;
            }
        }
        //cout << "No answers can meet such conditions" << endl;
        out: ;
    }
    return 0;
}