佩尔方程
关于佩尔方程
佩尔方程是具有\(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}\)
用这个去算就好了
