佩尔方程

关于佩尔方程

佩尔方程是具有(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)

可以推出

(egin{bmatrix}x_k\y_kend{bmatrix}=egin{bmatrix}x_1~~~~ny_1\y_1~~~~x_1end{bmatrix}^{k-1}egin{bmatrix}x_1\y_1end{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
", 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
", 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
",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
",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	", x, n, y);
            printf("%lld-%lld=1
", 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
", 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项

(egin{bmatrix}x_k\y_kend{bmatrix}=egin{bmatrix}x_1~~~~ny_1\y_1~~~~x_1end{bmatrix}^{k-1}egin{bmatrix}x_1\y_1end{bmatrix})

用这个去算就好了