佩尔方程

第一类佩尔方程

定义：形如x2−dy2=1x^2 - dy^2 = 1x2−dy2=1(d>1,且d不是完全平方数)
要求第一类佩尔方程的解都是正整数解,也即(x,y),x>0,y>0(x,y),x>0,y>0(x,y),x>0,y>0

为什么要求d不是完全平方数呢？
我们假设d是完全平方数
则x2−dy2=1x^2 - dy^2 = 1x2−dy2=1等价于(x+d∗y)(x−d∗y)=1(x+\sqrt d*y)(x-\sqrt d*y) = 1(x+d∗y)(x−d∗y)=1
所以它的解只有可能是(1,0),(-1,0)，不是正整数解。

为什么d>1呢?因为如果d<=1，解也不会是正整数解，自己证明试试看。

显然佩尔方程x2−dy2=1x^2 - dy^2 = 1x2−dy2=1有无数多个解，那么如何求最小的正整数解呢

我们把式子化简一下就变成了，x=1+dy2x=\sqrt{1+dy{2}}x=1+dy2,我们从y=1开始枚举y，如果1+dy2\sqrt{1+dy{2}}1+dy2是整数的话，此时的x,y就是最小的正整数，也就得到了最小正整数解

那么如何得到所有解呢，也就是说我们最小正整数解是(x1,y1)，那么如何得到(x2,y2),(x3,y3),…？

这里给出通解的迭代公式
xn=xn−1x1+dyn−1y1x_n = x_{n-1}x_1 + dy_{n-1}y_1xn=xn−1x1+dyn−1y1
yn=xn−1y1+yn−1x1y_n = x_{n-1}y_1 + y_{n-1}x_1yn=xn−1y1+yn−1x1
这个通解公式怎么推导的呢？我们要知其然，而知其所以然

推导：我们令(x1,y1),(x2,y2)为方程的两个特解
x12−dy12=1x_1^2-dy_1^2=1x12−dy12=1----①
x22−dy22=1x_2^2-dy_2^2=1x22−dy22=1----②
两式想乘得到
(x12−dy12)∗(x22−dy22)=1(x_1^2-dy_1^2)*(x_2^2-dy_2^2)=1(x12−dy12)∗(x22−dy22)=1
展开得到
x12x22−dx12y22−dy12x22+d2y12y22=1x_1^2x_2^2-dx_1^2y_2^2-dy_1^2x_2^2+d^2y_1^2y_2^2=1x12x22−dx12y22−dy12x22+d2y12y22=1
-> [x12x22+d2y12y22]−d[dx12y22+y12x22]=1[x_1^2x_2^2+d^2y_1^2y_2^2]-d[dx_1^2y_2^2+y_1^2x_2^2]=1[x12x22+d2y12y22]−d[dx12y22+y12x22]=1
-> [x12x22+d2y12y22+2dx1x2y1y2]−d[dx12y22+y12x22+2x1x2y1y2]=1[x_1^2x_2^2+d^2y_1^2y_2^2+2dx_1x_2y_1y_2]-d[dx_1^2y_2^2+y_1^2x_2^2+2x_1x_2y_1y_2]=1[x12x22+d2y12y22+2dx1x2y1y2]−d[dx12y22+y12x22+2x1x2y1y2]=1
-> (x1x2−dy1y2)2−d(x1y2+x2y1)2=1(x_1x_2-dy_1y_2)^2-d(x_1y_2+x_2y_1)^2=1(x1x2−dy1y2)2−d(x1y2+x2y1)2=1
所以x3=x2x1+dy2y1,y3=x2y1+x1y2x_3=x_2x_1+dy_2y_1,y_3=x_2y_1+x_1y_2x3=x2x1+dy2y1,y3=x2y1+x1y2
我们把下标换为n就得到
xn=xn−1x1+dyn−1y1x_n = x_{n-1}x_1 + dy_{n-1}y_1xn=xn−1x1+dyn−1y1
yn=xn−1y1+yn−1x1y_n = x_{n-1}y_1 + y_{n-1}x_1yn=xn−1y1+yn−1x1

这样的话我们只要知道第一个特解(x1,y1)(x1,y1)(x1,y1)，和第n-1个特解(xn−1,yn−1)(x_{n-1},y_{n-1})(xn−1,yn−1),就能够求出第n个解，但是很明现如果暴力求第n个解的话，时间复杂度是O(n),如果n很大的话就会超时，我们这时候可以用矩阵快速幂logn的复杂度来就快速求取，第n个解。

步骤：

step1:先通过暴力找到最小正整数解，也就是第一组解(x1,y1)，也有一种快速找最小正整数解的方法，叫做连分数法，感兴趣的小伙伴自行百度。

代码：

#include<bits/stdc++.h>
typedef long long ll;
void solve(ll &x,ll &y,ll d)
{y=1;for(y=1;;y++){tempc = ceil(sqrt(d*y*y+1));tempf = floor(sqrt(d*y*y+1));if(tempc==tempf) break;//只有整数的ceil和floor相同，也就是向上取整和向下取整}x=tempc;return ;
}

step2:根据迭代公式找到所有我们要找的解，也就是
xn=xn−1x1+dyn−1y1x_n = x_{n-1}x_1 + dy_{n-1}y_1xn=xn−1x1+dyn−1y1
yn=xn−1y1+yn−1x1y_n = x_{n-1}y_1 + y_{n-1}x_1yn=xn−1y1+yn−1x1

方法1:暴力迭代法
代码：

#include<bits/stdc++.h>
typedef long long ll;
using namespace std;
void solve(ll &x,ll &y,ll d)
{ll tempc,tempf; for(y=1;;y++){tempc = ceil(sqrt(d*y*y+1));tempf = floor(sqrt(d*y*y+1));if(tempc==tempf) break;}x=tempc;return ;
}
int main()
{ll x1,y1,d;cin>>d;//输入佩恩方程的d solve(x1,y1,d);ll xxi=x1,yyi=y1;ll xi,yi;int n; for(int i=2;i<=n;i++)//这里可以用于暴力求第n个解输入n就行 {xi = xxi*x1+d*yyi*y1;yi = xxi*y1+yyi*x1;xxi = xi;yyi = yi; }cout<<xi<<' '<<yi<<endl;return 0;
}

方法2:矩阵快速幂迭代法
如果我们需要用矩阵快速幂来求的话，我们就需要把通解转化为矩阵的表达形式。

[xnyn]\begin{bmatrix} x_n\\ y_n \end{bmatrix}[xnyn]=[x1dy1y1x1]n−1\begin{bmatrix} x_1&dy_1 \\ y_1&x_1 \end{bmatrix}^{n-1}[x1y1dy1x1]n−1[x1y1]\begin{bmatrix} x_1\\ y_1 \end{bmatrix}[x1y1]
这样我们就可以用矩阵快速幂来求取方程的第n个解，时间复杂度是O(logn),矩阵快速幂就类似于普通的快速幂，只不过我们的底数是一个矩阵，而不是一个常数，如果不了解矩阵快速幂的，请出门右转(矩阵快速幂)。

代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll Maxn=10;
/*矩阵乘法*/
void Matrix_ksm(ll a[][Maxn],ll b[][Maxn],ll size)
{ll temp[Maxn][Maxn]={0};//注意一定要初始化为0 for(ll i=1;i<=size;i++)for(ll j=1;j<=size;j++)for(ll k=1;k<=size;k++){temp[i][j]+=a[i][k]*b[k][j];}for(ll i=1;i<=size;i++)for(ll j=1;j<=size;j++)a[i][j]=temp[i][j];return ;
}
void solve(ll &x,ll &y,ll d)
{ll tempc,tempf; for(y=1;;y++){tempc = ceil(sqrt(d*y*y+1));tempf = floor(sqrt(d*y*y+1));if(tempc==tempf) break;}x=tempc;return ;
}
int main()
{ll x1,y1;ll n,d;printf("请输入d:\n");cin>>d;solve(x1,y1,d);//我们先求第一个解，也就是最小正整数解
//  cout<<x1<<' '<<y1<<'\n';printf("请输入要求的第n个矩阵,n=:\n");cin>>n;ll ans[Maxn][Maxn]={0};for(ll i=1;i<=Maxn;i++) ans[i][i]=1;ll dishu[Maxn][Maxn];dishu[1][1]=x1;dishu[1][2]=d*y1;dishu[2][1]=y1;dishu[2][2]=x1;ll zhishu=n-1;/*矩阵快速幂*/ while(zhishu){if(zhishu&1){zhishu--;Matrix_ksm(ans,dishu,2);}else {zhishu>>=1;Matrix_ksm(dishu,dishu,2);}}ll xn = ans[1][1]*x1+ans[1][2]*y1;ll yn = ans[2][1]*x1+ans[2][2]*y1;cout<<xn<<' '<<yn<<endl;return 0;
}

第一类佩尔方程例题讲解

有句熟话说得好“光说不练假把式”，下面我们就来分析几道例题

Street Numbers
题目大意:找出n，m使得1+2+3+....+(n−1)=(n+1)+(n+2)....+m1+2+3+....+(n-1)=(n+1)+(n+2)....+m1+2+3+....+(n−1)=(n+1)+(n+2)....+m,输出前十个解(n,m)。每个数字占是个字符，向右对齐，且每一排升序排列。

分析：我们不难发现，等式的左右都是一个等差数列，那么我们处理一下就可以得到这样一个等式，(2m+1)2−8n2=1(2m+1)^2-8n^2=1(2m+1)2−8n2=1,我们令x=(2m+1),y=n,d=8x=(2m+1),y=n,d=8x=(2m+1),y=n,d=8,这不就是一个裸的佩尔方程吗？我们只需要找到x，y就能找到m，n。

AC代码：

#include<iostream>
#include<cmath>
#include<cstdio>
using namespace std;
typedef long long ll;
/*找佩尔方程的第一个解，也就是正整数解*/
void solve(ll &x,ll &y,ll d)
{ll tempc,tempf; for(y=1;;y++){tempc = ceil(sqrt(d*y*y+1));tempf = floor(sqrt(d*y*y+1));if(tempc==tempf) break;}x=tempc;return ;
}
int main()
{ll x1,y1,d=8;solve(x1,y1,d);ll xxi=x1,yyi=y1;ll xi,yi;/*通过迭代公式依次找到前十个解，注意这道题答案，第一个解不算在里面，所以前十个解是从第2到第11个解*/for(int i=2;i<=11;i++)//这里可以用于暴力求第n个解输入n就行 {xi = xxi*x1+d*yyi*y1;yi = xxi*y1+yyi*x1;xxi = xi;yyi = yi; printf("%10lld%10lld\n",yyi,(xxi-1)/2);/*注意题目要求每一个数字占十个字符，每一排升序排列，所以我们先输出n，再输出m，因为(n<m),*/}return 0;
}

我们再来一道板子题No more tricks, Mr Nanguo

题目大意：就是一开始乐队有n个人，可以站成一个边长为x的正方形，现在走了一个人，把剩余的n-1个人分成d(这个d题目会给出)个边长为y的正方形站列，显然n有无穷多个解，对应的x，y也有无穷多组解，现在让你求第k(k题目会给出)小的正整数解x

这道题自己尝试着写一下，我就不写题解了，不过提醒一点，运算的时候注意取模，不然数据太大long long都可能超,由于k会比较大，所以不建议用暴力迭代法会超时，用矩阵快速幂来解决。

AC代码：

#include<iostream>
#include<cmath>
#include<cstdio>
using namespace std;
typedef long long ll;
const int mod=8191;
const ll Maxn=10;
/*矩阵乘法*/
void Matrix_ksm(ll a[][Maxn],ll b[][Maxn],ll size)
{ll temp[Maxn][Maxn]={0};//注意一定要初始化为0 for(ll i=1;i<=size;i++)for(ll j=1;j<=size;j++)for(ll k=1;k<=size;k++){/*在矩阵运算的时候，我们把每一个可能超范围的地方都取模*/temp[i][j]=(temp[i][j]%mod+a[i][k]%mod*b[k][j]%mod)%mod;}for(ll i=1;i<=size;i++)for(ll j=1;j<=size;j++)a[i][j]=temp[i][j];return ;
}
void solve(ll &x,ll &y,ll d)
{ll tempc,tempf; for(y=1;;y++){tempc = ceil(sqrt(d*y*y+1));tempf = floor(sqrt(d*y*y+1));if(tempc==tempf) break;}x=tempc;return ;
}
int main()
{ll x1,y1;ll n,d;while(scanf("%lld %lld",&d,&n)!=EOF){if(d<=1||ceil(sqrt(d))==floor(sqrt(d))){cout<<"No answers can meet such conditions"<<'\n';continue;} solve(x1,y1,d);//我们先求第一个解，也就是最小正整数解 ll ans[Maxn][Maxn]={0};for(ll i=1;i<=Maxn;i++) ans[i][i]=1;ll dishu[Maxn][Maxn];dishu[1][1]=x1;dishu[1][2]=d*y1;dishu[2][1]=y1;dishu[2][2]=x1;ll zhishu=n-1;/*矩阵快速幂*/ while(zhishu){if(zhishu&1){zhishu--;Matrix_ksm(ans,dishu,2);}else {zhishu>>=1;Matrix_ksm(dishu,dishu,2);}}ll xn = ans[1][1]%mod*x1%mod+ans[1][2]%mod*y1%mod;ll yn = ans[2][1]%mod*x1%mod+ans[2][2]%mod*y1%mod;cout<<xn%mod<<'\n';}return 0;
}

第二类佩尔方程

定义:形如x2−dy2=kx^2 - dy^2 = kx2−dy2=k的方程叫做第二类佩恩方程
第二类佩尔方程可以看成是第一类佩尔的拓展，也相对复杂一点

解第二类佩恩方程，需要用到第一类方程的解。

推导：
我们假设x2−dy2=k−−−−−①x^2 - dy^2 = k-----①x2−dy2=k−−−−−①的一个特解为p，qp，qp，q
设x2−dy2=1−−−−−②x^2 - dy^2 = 1-----②x2−dy2=1−−−−−②的正整数解为x1,y1x1,y1x1,y1，①×②式就可以得到下面式子。
->(p2−dq2)∗(x12−dy12)=k(p^2-dq^2)*(x_1^2-dy_1^2)=k(p2−dq2)∗(x12−dy12)=k
->(px1±dqy1)2−d(py1±qx1)2=k(px_1\pm dqy_1)^2-d(py_1\pm qx_1)^2 =k(px1±dqy1)2−d(py1±qx1)2=k
所以通解为,x=px1±dqy1x=px_1\pm dqy_1x=px1±dqy1,y=py1±qx1y=py_1\pm qx_1y=py1±qx1
我们假设②式中d=2，那么x2−2y2=1x^2-2y^2=1x2−2y2=1，解得最小正整数解是x1=3,y1=2x_1=3,y_1=2x1=3,y1=2,那么①式的解就是xn=3xn−1±4yn−1x_n=3x_{n-1}\pm 4y_{n-1}xn=3xn−1±4yn−1,yn=2xn−1±3yn−1y_n=2x_{n-1}\pm 3y_{n-1}yn=2xn−1±3yn−1。

遇到第二类佩尔方程也不用害怕，只需要和第一类佩尔方程联系起来求出通解即可。我们更具第二类佩尔方程的通解，
x=px1+dqy1x=px_1+dqy_1x=px1+dqy1,y=py1+qx1y=py_1+ qx_1y=py1+qx1，可以转化为矩阵乘法的形式
[xnyn]\begin{bmatrix} x_n\\ y_n \end{bmatrix}[xnyn]=[pdqqp]n−1\begin{bmatrix} p &dq \\ q&p \end{bmatrix}^{n-1}[pqdqp]n−1[x1y1]\begin{bmatrix} x_{1}\\ y_{1} \end{bmatrix}[x1y1],这样就方便用矩阵快速幂求第k个解

请关注我看更多数论分支知识体系讲解，也别忘点个赞额！

佩尔方程(超详细推导+例题讲解) 每日一遍，算法再见！相关推荐

传输线方程—超详细推导！
自己学的时候前面的知识都忘了,相当于0基础重学,推导如果有误请告知!
【转】流体动力学控制方程（详细推导）
[转]流体动力学控制方程(详细推导) diyhoos 原文链接:https://blog.csdn.net/qq_42020563/article/details/80940387 CFD建立在流体力 ...
傅里叶级数和傅里叶变换超详细推导（DR_CAN）
傅里叶级数和傅里叶变换超详细推导(DR_CAN) Part I 三角函数的正交性 Part Ⅱ周期为2π\piπ的 f(x)的傅里叶展开 Part Ⅲ 周期为"2L"的函数展开为傅 ...
一文让你彻底搞懂最小二乘法（超详细推导）
要解决的问题在工程应用中,我们经常会用一组观测数据去估计模型的参数,模型是我们根据先验知识定下的.比如我们有一组观测数据(xi,yi)(x_i,y_i)(xi,yi)(一维),通过一些数据分析我 ...
图像处理——几种简单的旋转变换的超详细推导过程（点在同一坐标系的变换）（一）
图像处理--几种简单的旋转变换的超详细推导过程(同一坐标系)(一) 本文主要推导了二维和三维坐标系中的绕点和绕轴的旋转变换,推导过程比较详细,希望可以给大家提供一些帮助. 一.绕原点的旋转(二维) 二 ...
详细推导HMM模型之：EM算法
详细推导HMM模型+直观例子+项目应用(三):EM算法进行参数估计 HMM模型介绍 HMM参数估计两种使用情景需要EM的情景 EM算法求解HMM中的参数回顾前向和后向算法 A.B.π\piπ的估 ...
支持向量机(SVM)----超详细原理分析讲解
文章目录支持向量机(SVM) 直观的本质理解几个基础概念决策超平面的求解(SVM模型的推导) 最大硬间隔的寻找与公式构建拉格朗日乘数法的应用使用对偶问题求解一个小例子(求解决策超平面与决策 ...
卡尔曼滤波(kalman filter)超详细推导
1. 概率论相关知识这一节主要回忆概率论的一些相关基础知识,包括全概率公式.贝叶斯公式.协方差矩阵.多维高斯分布等等,对这些熟悉的可以直接跳到第2节看贝叶斯滤波 1.1 条件概率 P(x,y)=P( ...
xgboost实例_XGBoost超详细推导，终于有人讲明白了！
- XGB中树结点分裂的依据是什么? - 如何计算树节点的权值? - 为防止过拟合,XGB做了哪些改进? 相信看到这篇文章的各位对XGBoost都不陌生,的确,XGBoost不仅是各大数据科学比赛的必 ...

佩尔方程(超详细推导+例题讲解) 每日一遍，算法再见！

这里写目录标题

佩尔方程

第一类佩尔方程

第一类佩尔方程例题讲解

第二类佩尔方程

佩尔方程(超详细推导+例题讲解) 每日一遍，算法再见！相关推荐

最新文章

热门文章