大家好,我是小典,我来为大家解答以上问题。佩尔方程通解的证明,佩尔方程很多人还不知道,现在让我们一起来看看吧!
方程没有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;
}
本文到此讲解完毕了,希望对大家有帮助。