c语言求佩尔方程的解设计思路,c语言版 佩尔方程求最小正整数解及第k解(矩阵快速幂)...
佩爾方程講解連接:
若一個丟番圖方程具有以下的形式:
且
為正整數,則稱此方程為佩爾方程(英文:Pell's equation 德文:Pellsche Gleichung)
若
是完全平方數,則這個方程式只有解
(實際上對任意的
,
都是解)。對于其余情況,拉格朗日證明了佩爾方程總有解。而這些解可由
的連分數求出。
設
?是
的連分數表示:
的漸近分數列,由連分數理論知存在?
?使得(pi,qi) 為佩爾方程的解。取其中最小的?
,將對應的 (pi,qi) 稱為佩爾方程的基本解,或最小解,記作(x1,y1) ,則所有的解(xi,yi) 可表示成如下形式:
或者由以下遞推公式得到:
?
——————————————————分割線————————————————————
求得佩爾方程最小正整數解后,由公式
及
可求得第k解(X1,Y1為最小正整數解)。
到這里你可能會想用遞歸的方法求解Xk及Yk。可是事實上如果k的值很大的話,就會花費好多時間。所以在這里求解的時候,用矩陣快速冪便可節約很多時間。
現在構造矩陣,如下圖
swun oj 里的一題,請參考,以便理解
#include
#include
#include
using namespace std;
typedef __int64 ll;
#define Mod 1000000007
ll x,y,n,k;
struct PellAns
{
ll p,q;
};
struct Node
{
ll g,h;
};
struct Matrix
{
ll a[2][2];
void init()
{
a[0][0]=x%Mod;a[0][1]=y%Mod;
a[1][0]=(n%Mod*y%Mod%Mod)%Mod;a[1][1]=x%Mod;
}
};
//矩陣乘法
Matrix matrix_mul(Matrix a,Matrix b)
{
ll i,j,k;
Matrix ans;
for(i=0;i<2;i++)
{
for(j=0;j<2;j++)
{
ans.a[i][j]=0;
for(k=0;k<2;k++)
ans.a[i][j]=(ans.a[i][j]%Mod+(a.a[i][k]%Mod*b.a[k][j]%Mod)%Mod)%Mod;
}
}
return ans;
}
//矩陣快速冪
Matrix mult(Matrix a,ll b)
{
Matrix ans;
ans.a[0][0]=1;ans.a[0][1]=0;
ans.a[1][0]=0;ans.a[1][1]=1;
while(b)
{
if(b & 1)
ans=matrix_mul(ans,a);
b>>=1;
//cout<
a=matrix_mul(a,a);
}
return ans;
}
//求佩爾方程最小正整數解...模板
PellAns Solve( ll n1)
{
PellAns s[4];
Node w[4];
int a[4];
s[0].p=0; s[0].q=1;
s[1].p=1; s[1].q=0;
a[0]=(ll)floor(sqrt( (double)n ));
a[2]=a[0];
w[1].g=0;w[1].h=1;
while( 1 )
{
w[2].g = -w[1].g+a[2]*w[1].h;
w[2].h = (n1-w[2].g*w[2].g)/w[1].h;
a[3] = (ll)floor( (double)(w[2].g+a[0])/w[2].h );
s[2].p = a[2]*s[1].p+s[0].p;
s[2].q = a[2]*s[1].q+s[0].q;
if( (s[2].p*s[2].p-n1*s[2].q*s[2].q) == 1 &&s[2].p>0&&s[2].q>0 )
return s[2];
w[0]=w[1];w[1]=w[2];
a[2]=a[3];
s[0]=s[1];s[1]=s[2];
}
}
int main()
{
PellAns ans;
// freopen("a.in","r",stdin);
//freopen("1.out","w",stdout);
while( ~scanf("%I64d%I64d",&n,&k) )
{
if(sqrt(double(n))*sqrt(double(n))==n) {
printf("No solution\n");continue;
}
ans = Solve(n);//求得佩爾方程最小正整數解
x=ans.p%Mod,y=ans.q%Mod;
Matrix tmp,ans1;
tmp.init(); //初始化
ans1=mult(tmp,(k-1)%Mod);
ll x1=x%Mod;
x=((ans1.a[0][0]%Mod*x%Mod)%Mod+(ans1.a[1][0]%Mod*y%Mod)%Mod)%Mod;
y=((ans1.a[0][1]%Mod*x1%Mod)+(ans1.a[1][1]%Mod*y%Mod)%Mod)%Mod;
printf("%I64d,%I64d %I64d,%I64d\n",ans.p,ans.q,x%Mod,y%Mod);
}
return 0;
}
總結
以上是生活随笔為你收集整理的c语言求佩尔方程的解设计思路,c语言版 佩尔方程求最小正整数解及第k解(矩阵快速幂)...的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: c语言中有哪些函数关系,C语言中有哪些常
- 下一篇: C语言输出最后一个空格去掉,新人提问:如