深圳推荐企业网站制作维护,论坛网站开发开题报告,邢台学校网站建设报价,域名备案查询网站备案信息查询佩尔方程讲解连接#xff1a;若一个丢番图方程具有以下的形式#xff1a;且为正整数#xff0c;则称此方程为佩尔方程(英文#xff1a;Pells equation 德文#xff1a;Pellsche Gleichung) 若是完全平方数#xff0c;则这个方程式只有解(实际上对任意的#xff0c;都是解…佩尔方程讲解连接若一个丢番图方程具有以下的形式且为正整数则称此方程为佩尔方程(英文Pells equation 德文Pellsche Gleichung) 若是完全平方数则这个方程式只有解(实际上对任意的都是解)。对于其余情况拉格朗日证明了佩尔方程总有解。而这些解可由的连分数求出。 设 是的连分数表示:的渐近分数列由连分数理论知存在 使得(pi,qi) 为佩尔方程的解。取其中最小的 将对应的 (pi,qi) 称为佩尔方程的基本解或最小解记作(x1,y1) 则所有的解(xi,yi) 可表示成如下形式 或者由以下递推公式得到 ——————————————————分割线————————————————————求得佩尔方程最小正整数解后由公式及可求得第k解(X1Y1为最小正整数解)。 到这里你可能会想用递归的方法求解Xk及Yk。可是事实上如果k的值很大的话就会花费好多时间。所以在这里求解的时候用矩阵快速幂便可节约很多时间。现在构造矩阵如下图swun oj 里的一题请参考以便理解#include #include #include using namespace std;typedef __int64 ll;#define Mod 1000000007ll 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(i0;i2;i){for(j0;j2;j){ans.a[i][j]0;for(k0;k2;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)ansmatrix_mul(ans,a);b1;//coutamatrix_mul(a,a);}return ans;}//求佩尔方程最小正整数解...模板PellAns Solve( ll n1){PellAns s[4];Node w[4];int a[4];s[0].p0; s[0].q1;s[1].p1; s[1].q0;a[0](ll)floor(sqrt( (double)n ));a[2]a[0];w[1].g0;w[1].h1;while( 1 ){w[2].g -w[1].ga[2]*w[1].h;w[2].h (n1-w[2].g*w[2].g)/w[1].h;a[3] (ll)floor( (double)(w[2].ga[0])/w[2].h );s[2].p a[2]*s[1].ps[0].p;s[2].q a[2]*s[1].qs[0].q;if( (s[2].p*s[2].p-n1*s[2].q*s[2].q) 1 s[2].p0s[2].q0 )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);//求得佩尔方程最小正整数解xans.p%Mod,yans.q%Mod;Matrix tmp,ans1;tmp.init(); //初始化ans1mult(tmp,(k-1)%Mod);ll x1x%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;}
