佩尔方程

关于佩尔方程

佩尔方程是具有\(x^2-ny^2=1\)形式的丢番图方程（不定方程）

当\(n\)为完全平方数的时候，这个方程只有平凡解\((\pm1,0)\)，对于其他情况拉格朗日证明了佩尔方程总有非平凡解。而这些解都可以由\(\sqrt{n}\)的连分数求出

关于连分数怎么求，比如\(\sqrt{7}=2+\frac{1}{r_1},a_0=2\)

\(r_1=\frac{\sqrt{7}+2}{3}=1+\frac{1}{r_2},a_1=1\)

\(r_2=\frac{\sqrt{7}+1}{2}=1+\frac{1}{r_2},a_2=1\)

\(r_3=\frac{\sqrt{7}+1}{3}=1+\frac{1}{r_3},a_3=1\)

\(r_4=\sqrt{7}+2=4+\frac{1}{r_4},a_4=4\)

从这开始就是\([4,1,1,1]\)这样的周期了

我们取第一个周期

我们得到了\([2:1,1,1,1]\)

这时我们的连分数为\(\frac{8}{3}\),\(h_3=8,k_3=3\)，发现\(8^2-7*3^2=1\)

这就是佩尔方程的最小解了

其实对于\(\frac{h_n}{k_n}\)的\(h_n,k_n\)也有递推式

\(h_n=a_nh_{n-1}+h_{n-2}\)

\(k_n=a_nk_{n-1}+k_{n-2}\)

我们由最小解可以推得通解

此时的最小解为\(x_1=8,y_1=3\)

有

\(x_{i+1}=x_1x_i+ny_1y_i\)

\(y_{i+1}=x_1y_i+y_1x_i\)

可以推出

\(\begin{bmatrix}x_k\\y_k\end{bmatrix}=\begin{bmatrix}x_1~~~~ny_1\\y_1~~~~x_1\end{bmatrix}^{k-1}\begin{bmatrix}x_1\\y_1\end{bmatrix}\)

可以由矩阵快速幂得到第k项通解

例题 hdu2281-Square Number

传送门

解题思路：

\(x^2=\frac{n(2n+1)(n+1)}{6n}=\frac{(2n+1)(n+1)}{6}=\frac{1}{3}(n^2+\frac{3}{2}n+\frac{9}{16}-\frac{1}{16})=\frac{1}{3}(n+\frac{3}{4})^2-\frac{1}{48}\)

所以\((4n+3)^2-48x^2=1\)

这是显然的佩尔方程形式，而这个最小解为\(4n+3=7,x=1\)即\(n=1,x=1\)

可以用最小解去递推通解，还要判一下是否\(n\)满足为整数

打了个表看了下最多跑到第十六项 暴力来就好了

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/hash_policy.hpp>
#include <ext/pb_ds/tree_policy.hpp>
#include <ext/pb_ds/trie_policy.hpp>
using namespace __gnu_pbds;
using namespace std;
// freopen("k.in", "r", stdin);
// freopen("k.out", "w", stdout);
// clock_t c1 = clock();
// std::cerr << "Time:" << clock() - c1 <<"ms" << std::endl;
//#pragma comment(linker, "/STACK:1024000000,1024000000")
mt19937 rnd(time(NULL));
#define de(a) cout << #a << " = " << a << endl
#define rep(i, a, n) for (int i = a; i <= n; i++)
#define per(i, a, n) for (int i = n; i >= a; i--)
#define ls ((x) << 1)
#define rs ((x) << 1 | 1)
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> PII;
typedef pair<double, double> PDD;
typedef pair<char, char> PCC;
typedef pair<ll, ll> PLL;
typedef vector<int> VI;
#define inf 0x3f3f3f3f
const ll INF = 0x3f3f3f3f3f3f3f3f;
const ll MAXN = 1e5 + 7;
const ll MAXM = 4e5 + 7;
const ll MOD = 1e9 + 7;
const double eps = 1e-7;
ll x[105], y[105];
ll ansn[105];
int main()
{
    ll n = 48;
    x[1] = 7, y[1] = 1;
    ansn[1] = 1;
    for (int i = 2; i <= 20; i++)
    {
        x[i] = x[i - 1] * x[1] + n * y[1] * y[i - 1];
        if ((x[i] - 3) % 4 == 0)
            ansn[i] = (x[i] - 3) / 4;
        else
            ansn[i] = -1;
        y[i] = x[1] * y[i - 1] + y[1] * x[i - 1];
    }
    ll N;
    while (~scanf("%lld", &N) && N)
    {
        ll ans = -inf;
        for (int i = 1; i <= 16; i++)
        {
            if (ansn[i] == -1)
                continue;
            if (N < ansn[i])
                break;
            ans = i;
        }
        printf("%lld %lld\n", ansn[ans], y[ans]);
    }
    return 0;
}

例题 Street Numbers

传送门

知道题意之后就跟上面一模一样的过程了

代码

#include <iostream>
#include <stdio.h>
#include <cmath>
#include <algorithm>
#include <string>
#include <cstring>
#include <string.h>
#include <map>
#include <queue>
#include <stack>
#include <list>
#include <cctype>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <set>
#include <vector>
#include <cstdlib>
#include <time.h>
using namespace std;
/* freopen("k.in", "r", stdin);
freopen("k.out", "w", stdout); */
// clock_t c1 = clock();
// std::cerr << "Time:" << clock() - c1 <<"ms" << std::endl;
//#pragma comment(linker, "/STACK:1024000000,1024000000")
#define de(a) cout << #a << " = " << a << endl
#define rep(i, a, n) for (int i = a; i <= n; i++)
#define per(i, a, n) for (int i = n; i >= a; i--)
#define ls ((x) << 1)
#define rs ((x) << 1 | 1)
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> PII;
typedef pair<double, double> PDD;
typedef pair<ll, ll> PLL;
typedef vector<int, int> VII;
#define inf 0x3f3f3f3f
const ll INF = 0x3f3f3f3f3f3f3f3f;
const ll MAXN = 2e6 + 7;
const ll MAXM = 4e5 + 7;
const ll MOD = 1e9 + 7;
const double eps = 1e-7;
const double pi = acos(-1.0);
ll x[20], y[20];
int main()
{
    vector<pair<ll, ll> > ans;
    x[1] = 3, y[1] = 1;
    ll n = 8;
    for (int i = 2; i <= 20; i++)
    {
        x[i] = x[i - 1] * x[1] + n * y[1] * y[i - 1];
        y[i] = x[1] * y[i - 1] + y[1] * x[i - 1];
        if ((x[i] - 1) % 2 == 0)
            ans.push_back({(x[i] - 1) / 2, y[i]});
        if (ans.size() == 10)
            break;
    }
    for (int i = 0; i < 10; i++)
        printf("%10lld%10lld\n", ans[i].second, ans[i].first);
    return 0;
}

推佩尔方程的板子

连分数版

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/hash_policy.hpp>
#include <ext/pb_ds/tree_policy.hpp>
#include <ext/pb_ds/trie_policy.hpp>
using namespace __gnu_pbds;
using namespace std;
// freopen("k.in", "r", stdin);
// freopen("k.out", "w", stdout);
// clock_t c1 = clock();
// std::cerr << "Time:" << clock() - c1 <<"ms" << std::endl;
//#pragma comment(linker, "/STACK:1024000000,1024000000")
mt19937 rnd(time(NULL));
#define de(a) cout << #a << " = " << a << endl
#define rep(i, a, n) for (int i = a; i <= n; i++)
#define per(i, a, n) for (int i = n; i >= a; i--)
#define ls ((x) << 1)
#define rs ((x) << 1 | 1)
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> PII;
typedef pair<double, double> PDD;
typedef pair<char, char> PCC;
typedef pair<ll, ll> PLL;
typedef vector<int> VI;
#define inf 0x3f3f3f3f
const ll INF = 0x3f3f3f3f3f3f3f3f;
const ll MAXN = 1e5 + 7;
const ll MAXM = 4e5 + 7;
const ll MOD = 1e9 + 7;
const double eps = 1e-7;
ll a[20000];
bool pell_minimum_solution(ll n, ll &x0, ll &y0)
{
    ll m = (ll)sqrt((double)n);
    double sq = sqrt(n);
    int i = 0;
    if (m * m == n)
        return false; //当n是完全平方数则佩尔方程无解
    a[i++] = m;
    ll b = m, c = 1;
    double tmp;
    do
    {
        c = (n - b * b) / c;
        tmp = (sq + b) / c;
        a[i++] = (ll)(floor(tmp));
        b = a[i - 1] * c - b;
        //printf("%lld %lld %lld\n",a[i-1],b,c);
    } while (a[i - 1] != 2 * a[0]);
    ll p = 1, q = 0;
    for (int j = i - 2; j >= 0; j--)
    {
        ll t = p;
        p = q + p * a[j];
        q = t;
        //printf("a[%d]=%lld %lld %lld\n",j,a[j],p,q);
    }
    if ((i - 1) % 2 == 0)
    {
        x0 = p;
        y0 = q;
    }
    else
    {
        x0 = 2 * p * p + 1;
        y0 = 2 * p * q;
    }
    return true;
}
int main()
{
    ll n, x, y;
    //x^2-ny^2=1
    while (~scanf("%lld", &n))
    {
        if (pell_minimum_solution(n, x, y))
        {
            printf("%lld^2-%lld*%lld^2=1\t", x, n, y);
            printf("%lld-%lld=1\n", x * x, n * y * y);
        }
    }
}

暴力

int ax[MAXN], ay[MAXN];
void pell(int &x, int &y, int n)
{ //暴力寻找pell方程最小解
    y = 1;
    while (true)
    {
        x = (ll)sqrt(n * y * y + 1);
        if (x * x - n * y * y == 1)
            break;
        y++;
    }
}
int main()
{
    int n;
    while (scanf("%d", &n) != EOF)
    {
        int m = (int)sqrt((double)n);
        if (m * m == n)
        { //d不能为完全平方数
            cout << "No Solution" << endl;
            continue;
        }
        int x = 0, y = 0;
        pell(x, y, n); //暴力找到最小解
        printf("%d %d\n", x, y);
    }
    return 0;
}

递推式

for (int i = 2; i <= 20; i++)
{
    x[i] = x[i - 1] * x[1] + n * y[1] * y[i - 1];
    y[i] = x[1] * y[i - 1] + y[1] * x[i - 1];
}

矩阵递推求第k项

\(\begin{bmatrix}x_k\\y_k\end{bmatrix}=\begin{bmatrix}x_1~~~~ny_1\\y_1~~~~x_1\end{bmatrix}^{k-1}\begin{bmatrix}x_1\\y_1\end{bmatrix}\)

用这个去算就好了