题解：求解佩尔方程 (x^2 - D y^2 = 1) 的最小正整数解

一、题目核心分析

1. 问题定义

佩尔方程是形如 (x^2 - D y^2 = 1) 的二元二次不定方程，其中 (D) 是正整数且不是完全平方数（若 (D) 是完全平方数，方程无正整数解）。本题要求找到该方程的最小正整数解（即 (x) 和 (y) 均为正整数，且在所有解中 (x) 最小，对应 (y) 也最小）。

2. 关键算法：连分数展开法

佩尔方程的最小解可通过对 (\sqrt{D}) 进行循环连分数展开求得，这是数论中解决佩尔方程的标准方法，核心原理如下：

• 任何非完全平方数的平方根 (\sqrt{D}) 都可以展开为无限循环连分数：(\sqrt{D} = [a_0; (a_1, a_2, ..., a_k)])（括号内为循环节）；
• 佩尔方程的最小解对应循环连分数的周期结束时的渐进分数（即循环节完整遍历一次后的收敛项）；
• 渐进分数的分子和分母分别对应最小解 (x) 和 (y)。

3. 连分数展开核心公式

对 (\sqrt{D}) 展开连分数时，需迭代计算以下参数（初始值）：

• (a_0 = \lfloor \sqrt{D} \rfloor)（(\sqrt{D}) 的整数部分）；
• (m_0 = 0)，(d_0 = 1)；
• 迭代公式（第 (i) 步）： [ m_{i+1} = a_i \cdot d_i - m_i ] [ d_{i+1} = \frac{D - m_{i+1}^2}{d_i} ] [ a_{i+1} = \left\lfloor \frac{a_0 + m_{i+1}}{d_{i+1}} \right\rfloor ]
• 渐进分数计算：设第 (i) 个渐进分数为 (\frac{p_i}{q_i})，则： [ p_i = a_i \cdot p_{i-1} + p_{i-2} ] [ q_i = a_i \cdot q_{i-1} + q_{i-2} ] 初始值：(p_{-1}=1, p_0=a_0)；(q_{-1}=0, q_0=1)。
• 终止条件：当迭代到 (a_{i} = 2a_0) 时，循环节结束，此时的渐进分数 (\frac{p_{i-1}}{q_{i-1}}) 即为最小解（(x=p_{i-1}, y=q_{i-1})）。

4. 大整数处理

由于 (D) 最大为 (1000)，最小解的 (x) 和 (y) 可能非常大（如样例中 (D=13) 时，(x=649)，(y=180)；更大 (D) 可能超出 64 位整数范围），因此代码中使用vector模拟大整数，实现加法和乘法运算。

二、完整代码

#include <bits/stdc++.h>
using namespace std;

// 将整数转换为大整数（vector存储，低位在前）
vector<int> vti(int x) {
    vector<int> v;
    if (x == 0) { // 处理x=0的特殊情况
        v.push_back(0);
        return v;
    }
    while (x) {
        v.push_back(x % 10); // 低位存于vector前端
        x /= 10;
    }
    return v;
}

// 大整数加法（vector存储，低位在前）
vector<int> operator+(vector<int> a, vector<int> b) {
    if (a.size() < b.size()) swap(a, b); // 确保a更长
    int jw = 0; // 进位
    for (int i = 0; i < a.size(); i++) {
        jw += a[i];
        if (i < b.size()) jw += b[i];
        a[i] = jw % 10; // 当前位结果
        jw /= 10; // 更新进位
    }
    if (jw) a.push_back(jw); // 处理最后的进位
    return a;
}

// 大整数乘法（vector存储，低位在前）
vector<int> operator*(vector<int> a, vector<int> b) {
    if (a.empty() || b.empty()) return {};
    vector<int> c(a.size() + b.size() - 1, 0); // 乘积初始长度：a.len + b.len - 1
    // 逐位相乘，累加结果
    for (int i = 0; i < a.size(); i++) {
        for (int j = 0; j < b.size(); j++) {
            c[i + j] += a[i] * b[j];
        }
    }
    // 处理进位
    int jw = 0;
    for (int i = 0; i < c.size(); i++) {
        jw += c[i];
        c[i] = jw % 10;
        jw /= 10;
    }
    if (jw) c.push_back(jw); // 最后的进位
    return c;
}

// 大整数 + 普通整数
vector<int> operator+(vector<int> a, int b) {
    return a + vti(b);
}

// 大整数 * 普通整数
vector<int> operator*(vector<int> a, int b) {
    return a * vti(b);
}

// 输出大整数（vector存储，低位在前，需反向输出）
void wrt(vector<int> v) {
    for (int i = v.size() - 1; i >= 0; i--) {
        putchar(v[i] + '0');
    }
}

// 比较两个大整数是否相等（vector存储，低位在前）
bool eq(vector<int> a, vector<int> b) {
    if (a.size() != b.size()) return false;
    for (int i = 0; i < a.size(); i++) {
        if (a[i] != b[i]) return false;
    }
    return true;
}

// 求解佩尔方程x² - D y² = 1的最小正整数解
void sol() {
    int D;
    scanf("%d", &D);
    int a0 = sqrt(D); // a0 = floor(sqrt(D))
    int m = 0, d = 1, a = a0;

    // 初始化渐进分数的分子p和分母q：p0=a0, p1=a0*a1 + p0；q同理
    vector<int> p_prev_prev = vti(1);  // p_{-1} = 1
    vector<int> p_prev = vti(a0);      // p_0 = a0
    vector<int> q_prev_prev = vti(0);  // q_{-1} = 0
    vector<int> q_prev = vti(1);       // q_0 = 1

    // 检查初始项是否为解（特殊情况，如D=2时a0=1，p0=1, q0=1：1-2*1=-1≠1，不满足）
    vector<int> x_sq = p_prev * p_prev;
    vector<int> y_sq_D = (q_prev * q_prev) * D;
    vector<int> right = y_sq_D + 1; // 右边：D*y² + 1
    if (eq(x_sq, right)) {
        wrt(p_prev);
        putchar(' ');
        wrt(q_prev);
        putchar('\n');
        return;
    }

    // 迭代展开连分数，计算渐进分数
    for (;;) {
        // 迭代更新m, d, a
        m = a * d - m;
        d = (D - m * m) / d;
        a = (a0 + m) / d;

        // 计算当前渐进分数的p和q：p_i = a*i * p_{i-1} + p_{i-2}
        vector<int> p_curr = (p_prev * a) + p_prev_prev;
        vector<int> q_curr = (q_prev * a) + q_prev_prev;

        // 检查是否满足佩尔方程：p_curr² - D*q_curr² == 1
        x_sq = p_curr * p_curr;
        y_sq_D = (q_curr * q_curr) * D;
        right = y_sq_D + 1;
        if (eq(x_sq, right)) {
            wrt(p_curr);
            putchar(' ');
            wrt(q_curr);
            putchar('\n');
            return;
        }

        // 更新前两项，准备下一轮迭代
        p_prev_prev = p_prev;
        p_prev = p_curr;
        q_prev_prev = q_prev;
        q_prev = q_curr;
    }
}

int main() {
    int tc;
    scanf("%d", &tc);
    while (tc--) {
        sol();
    }
    return 0;
}

三、关键逻辑详解

1. 大整数处理模块

由于佩尔方程的最小解可能超出普通整数（甚至64位整数）范围，代码使用vector<int>存储大整数，低位在前（如数字123存储为[3,2,1]），方便加法和乘法的进位处理：

• vti(int x)：将普通整数转换为大整数（vector存储）；
• operator+：实现两个大整数的加法，处理进位；
• operator*：实现两个大整数的乘法，先逐位相乘累加，再统一处理进位；
• wrt(vector<int> v)：反向输出vector中的元素，得到正确的数字顺序。

2. 连分数展开与渐进分数计算

（1）初始参数初始化

• (a0 = \lfloor \sqrt{D} \rfloor)：(\sqrt{D}) 的整数部分（如 (D=13) 时，(a0=3)）；
• 渐进分数初始值：(p_{-1}=1, p_0=a0)；(q_{-1}=0, q_0=1)（对应第一个渐进分数 (\frac{p_0}{q_0} = a0)）。

（2）迭代过程

每次迭代执行以下步骤：

1. 按公式更新 (m, d, a)（连分数的循环项参数）；
2. 按渐进分数公式计算当前的 (p_curr)（分子）和 (q_curr)（分母）；
3. 验证是否满足佩尔方程：计算 (p_curr^2) 和 (D \cdot q_curr^2 + 1)，若相等则为最小解；
4. 若不满足，更新前两项渐进分数，进入下一轮迭代。

（3）终止条件

当计算出的 (p_curr) 和 (q_curr) 满足 (p_curr^2 - D \cdot q_curr^2 = 1) 时，停止迭代，此时的 (p_curr) 和 (q_curr) 即为最小正整数解。这是因为连分数展开的周期性保证了第一个满足方程的渐进分数就是最小解。

3. 样例验证（以 (D=13) 为例）

• 初始：(a0=3)，(p0=3)，(q0=1)（(3^2 -13 \cdot 1^2 =9-13=-4≠1)，不满足）；
• 第一次迭代：(m=31-0=3)，(d=(13-9)/1=4)，(a=(3+3)/4=1) → (p1=13+1=4)，(q1=1*1+0=1)（(16-13=3≠1)）；
• 第二次迭代：(m=14-3=1)，(d=(13-1)/4=3)，(a=(3+1)/3=1) → (p2=14+3=7)，(q2=11+1=2)（(49-134=49-52=-3≠1)）；
• 第三次迭代：(m=13-1=2)，(d=(13-4)/3=3)，(a=(3+2)/3=1) → (p3=17+4=11)，(q3=12+1=3)（(121-139=121-117=4≠1)）；
• 第四次迭代：(m=13-2=1)，(d=(13-1)/3=4)，(a=(3+1)/4=1) → (p4=111+7=18)，(q4=13+2=5)（(324-1325=324-325=-1≠1)）；
• 第五次迭代：(m=14-1=3)，(d=(13-9)/4=1)，(a=(3+3)/1=6) → (p5=618+11=119)，(q5=65+3=33)（(119²-1333²=14161-13*1089=14161-14157=4≠1)）；
• 第六次迭代：(m=61-3=3)，(d=(13-9)/1=4)，(a=(3+3)/4=1) → (p6=1119+18=137)，(q6=1*33+5=38)（不满足）；
• ... 持续迭代直到循环节结束，最终得到 (p=649)，(q=180)（(649²-13180²=421201-1332400=421201-421200=1)），即为最小解。

四、复杂度分析

• 时间复杂度：对于每个 (D)，连分数的循环节长度 (k) 满足 (k = O(\sqrt{D} \log D))，每次迭代的大整数运算（加法、乘法）复杂度为 (O(L^2))（(L) 为大整数的位数，最大约为 (O(\log x))，(x) 为最小解的 (x) 值）。由于 (D \leq 1000)，(T \leq 100)，整体复杂度完全可控；
• 空间复杂度：大整数存储的空间复杂度为 (O(L))（(L) 为最大数字的位数），仅需常数级额外空间，效率极高。

五、注意事项

1. (D) 不是完全平方数：题目已保证，无需额外判断；
2. 大整数处理的细节：如低位在前的存储方式、进位处理、相等比较（需自定义eq函数）；
3. 初始项验证：部分特殊 (D) 可能在初始项 (p0, q0) 就满足方程（虽实际少见，但代码需覆盖）。