【算法竞赛学习笔记】佩尔方程-数学提升计划

title : 佩尔方程
date : 2021-10-31
tags : ACM,数学
author : Linno

佩尔方程

形如x2−dy2=1(d>1且d不为完全平方数)x^2-dy^2=1(d>1且d不为完全平方数)x2−dy2=1(d>1且d不为完全平方数)的不定方程称为佩尔方程
设两个解为(x1,y1)和(x2,y2)，那么有{x12−dy12=1x22−dy22=1得到(x12−dy12)(x22−dy22)=1展开得到x12x22−dx12y22−dy12x22+d2y12y22=1⇒(x12x22+d2y12y22)−d(dx12y22+y12x22)=1⇒(x12x22+d2y12y22+2dx1x2y1y2)−d(dx12y22+y12x22−2dx1x2y1y2)=1⇒(x1x2+dy1y2)−d(dx12y22+y12x22−2dx1x2y1y2)=1所以x3=x2x1+dy2y1,y3=x2y1+x1y2根据迭代公式，通解为：xn=xn−1x1+dyn−1y1,yn=xn−1y1+yn−1x1设两个解为(x_1,y_1)和(x_2,y_2)，那么有\\ \begin{cases} x_1^2-dy_1^2=1\\ x_2^2-dy_2^2=1\\ \end{cases} 得到(x_1^2-dy_1^2)(x_2^2-dy_2^2)=1\\ 展开得到x_1^2x_2^2-dx_1^2y_2^2-dy_1^2x_2^2+d^2y_1^2y_2^2=1\\ \Rightarrow (x_1^2x_2^2+d^2y_1^2y_2^2)-d(dx_1^2y_2^2+y_1^2x_2^2)=1\\ \Rightarrow (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-2dx_1x_2y_1y_2)=1\\ \Rightarrow (x_1x_2+dy_1y_2)-d(dx_1^2y_2^2+y_1^2x_2^2-2dx_1x_2y_1y_2)=1\\ 所以x_3=x_2x_1+dy_2y_1,y_3=x_2y_1+x_1y_2\\ 根据迭代公式，通解为：\\ x_n=x_{n-1}x_1+dy_{n-1}y_1,y_n=x_{n-1}y_1+y_{n-1}x_1\\ 设两个解为(x1,y1)和(x2,y2)，那么有{x12−dy12=1x22−dy22=1得到(x12−dy12)(x22−dy22)=1展开得到x12x22−dx12y22−dy12x22+d2y12y22=1⇒(x12x22+d2y12y22)−d(dx12y22+y12x22)=1⇒(x12x22+d2y12y22+2dx1x2y1y2)−d(dx12y22+y12x22−2dx1x2y1y2)=1⇒(x1x2+dy1y2)−d(dx12y22+y12x22−2dx1x2y1y2)=1所以x3=x2x1+dy2y1,y3=x2y1+x1y2根据迭代公式，通解为：xn=xn−1x1+dyn−1y1,yn=xn−1y1+yn−1x1

矩阵快速幂迭代法

[xnyn]=[x1dy1y1x1]n−1[x1y1]\begin{bmatrix}x_n\\y_n\end{bmatrix}=\begin{bmatrix}x_1&dy_1\\y_1&x_1\end{bmatrix}^{n-1}\begin{bmatrix}x_1\\y_1\end{bmatrix} [xnyn]=[x1y1dy1x1]n−1[x1y1]

这样我们就可以用矩阵快速幂来求取方程的第n个解，时间复杂度为O(logn)。

HDU3292 No more tricks, Mr Nanguo

#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=k的方程叫做第二类佩尔方程形如x^2-dy^2=k的方程叫做第二类佩尔方程形如x2−dy2=k的方程叫做第二类佩尔方程
推导：假设x2−dy2=k①的一个特解为p,q,设x2−dy2=1②整数解为x1,y1①×②得：(p2−dq2)∗(x12−dy12)=k所以通解为x=px1±dqy1,y=py1±qx1推导：\\ 假设x^2-dy^2=k①的一个特解为p,q,设x^2-dy^2=1②整数解为x_1,y_1\\ ①×②得：(p^2-dq^2)*(x_1^2-dy_1^2)=k\\ 所以通解为x=px_1\pm dqy_1,y=py_1\pm qx_1\\ 推导：假设x2−dy2=k①的一个特解为p,q,设x2−dy2=1②整数解为x1,y1①×②得：(p2−dq2)∗(x12−dy12)=k所以通解为x=px1±dqy1,y=py1±qx1

矩阵形式为：[xnyn]=[pdqqp]n−1[x1y1]矩阵形式为：\\\begin{bmatrix}x_n\\y_n\end{bmatrix}=\begin{bmatrix}p&dq\\q&p\end{bmatrix}^{n-1}\begin{bmatrix}x_1\\y_1\end{bmatrix} 矩阵形式为：[xnyn]=[pqdqp]n−1[x1y1]

参考资料：

https://blog.csdn.net/TheWayForDream/article/details/113772229