佩尔方程通解的证明(佩尔方程)

导读 大家好,我是小典,我来为大家解答以上问题。佩尔方程通解的证明,佩尔方程很多人还不知道,现在让我们一起来看看吧!方程没有10亿以内的解...

大家好,我是小典,我来为大家解答以上问题。佩尔方程通解的证明,佩尔方程很多人还不知道,现在让我们一起来看看吧!

方程没有10亿以内的解,手动计算应该是行不通的

我写了一个程序,计算出方程的一组特解:

x=379516400906811930638014896080

y=12055735790331359447442538767

算法是求√991的渐进连分数h0/k0,h1/k1...h(l-1)/k(l-1),l是连分数的最小正周期,若l是偶数,则h(l-1)和k(l-1)就是pell方程的一组特解,否则h(2l - 1)和k(2l- 1)是方程的一组特解。证明可以在任何一本数论书上找到,不再赘述,也可以参考wiki百科的pell方程条目

程序附在后面,C++写的,任何一个windows平台下的C++编译器都可以编译(VC++,DEVCPP等等),涉及到大整数的运算比较麻烦,用matlab写起来会更简单

程序清单:

#include <stdio.h>

#include <stdlib.h>

#include <memory.h>

#include <math.h>

#define LEN 200

#define BASE 10000000 //大整数用10^7进制来计算

typedef __int64 INT;

struct BIG //大整数结构定义,整数每7位拼成一个int64,这样加减乘计算耗时较少

{

INT a[LEN];//存放大整数的数组

int len; //标记大整数的长度,不能超过200

}h[3],k[3];

INT a[1 << 14],c[1 << 14],q[1 << 14];

/*实际是计算:

h[i] = n*h[i - 1] + h[i - 2];

k[i] = n*k[i - 1] + k[i - 2];

这是连分数的递推形式

*/

void go(INT n)

{

int i;

for(i = 0;i < h[1].len;i++)

h[2].a[i] = h[1].a[i] * n + h[0].a[i];

for(i = 0;i < h[1].len;i++)

h[2].a[i + 1] += h[2].a[i] / BASE,h[2].a[i] %= BASE;

if(h[2].a[h[2].len])

h[2].len++;

for(i = 0;i < k[1].len;i++)

k[2].a[i] = k[1].a[i] * n + k[0].a[i];

for(i = 0;i < k[1].len;i++)

k[2].a[i + 1] += k[2].a[i] / BASE,k[2].a[i] %= BASE;

if(k[2].a[k[2].len])

k[2].len++;

for(i = 0;i < 2;i++)

h[i] = h[i + 1],k[i] = k[i + 1];

}

void print(BIG &a) //把大整数输出

{

int i;

printf("%I64d",a.a[a.len - 1]);

for(i = a.len - 2;i >= 0;i--)

printf("%07I64d",a.a[i]);

}

int main()

{

int i,j,len;

bool use = false;

c[0] = 0,q[0] = 1,a[0] = 31;

for(i = 1;;i++) //寻找连分数循环节,j-i是最小循环周期

{

c[i] = a[i - 1] * q[i - 1] - c[i - 1];

q[i] = (991 - c[i] * c[i])/q[i - 1];

a[i] = (a[0] + c[i])/q[i];

if(!use)

{

for(j = 0;j < i;j++)

{

if(c[i] == c[j] && a[i] == a[j] && q[i] == q[j])

{

len = i - j;

break;

}

}

if(j < i)

break;

}

}

if(len & 1) //循环节长度是奇数,长度要乘以2

len <<= 1;

for(i = 0;i < 3;i++)

{

memset(h[i].a,0,sizeof(h[i].a));

memset(k[i].a,0,sizeof(k[i].a));

h[i].len = k[i].len = 0;

}

h[0].a[0] = a[0],h[0].len = 1;

h[1].a[0] = a[0] * a[1] + 1;

h[1].len = 1;

k[0].a[0] = 1,k[0].len = 1;

k[1].a[0] = a[1],k[1].len = 1;

for(i = 2;i < len;i++)

go(a[i]);

printf("x = ");

print(h[2]);

printf(" ");

printf("y = ");

print(k[2]);

printf(" ");

system("PAUSE");

return 0;

}

本文到此讲解完毕了,希望对大家有帮助。

最新文章