一个多变量丢番图方程

mathe · 发表于 2008-11-4 14:10:01

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？欢迎注册

×

http://topic.csdn.net/u/20081102 ... d-971e4ff1e0ea.html 中,northwolves问了一个问题，
i)对于四个不同自然数$a<=b<=c<=d$,满足$a^2+b^2+c^2+d^2=abcd$,求出一亿以内所有的(a,b,c,d)组合
ii)对于n个不同自然数$a_1<=a_2<=...<=a_n$,满足$a_1^2+a_2^2+...+a_n^2=a_1a_2...a_n$,求出一亿以内所有组合$(a_1,a_2,...,a_n)$

开始northwolves和medie2005试用了很多不错的方法求得了（特别是问题一的）部分解，然后medie2005发现一个很有趣的规律，这些解中很多数字会大量重复出现。在帖子中50楼，litaoye发现了一个重要的规律：
假设我们已经找到$a^2+b^2+c^2+d^2=abcd$的一组解(a,b,cd)
我们知道d是方程$X^2-abcX+a^2+b^2+c^2=0$的一个解，而这个方程还有另外一个解$e=abc-d$.
如果e<d,那么我们就可以将这个解降到另外一个更小一点的解(a,b,c,e)(这里它们顺序可能需要重派）。
所有反过来，如果对于某个解(a,b,c,d),我们分别假设其中a是某个解(b,c,d,x)对应的e;b是某个解(a,c,d,x)对应的e等等就可以通过一个基础解得出很多其他派生解。
所以对于这个问题，我们可以求出所有那些最基本的解，然后其他解都可以通过上面方法派生出来就可以了。
litaoye还分析出来对于4个变量的情况，所有其他解都可以通过基础解(2,2,2,2)直接或间接派生出来。

我们可以进一步分析更加一般的情况，所以如果有某个基础解$(a_1,a_2,...,a_n)$无法通过其他解派生出来（就是一个最基本的解），由于$a_n$是方程
$X^2-a_1a_2...a_{n-1}X+a_1^2+a_2^2+...+a_{n-1}^2=0$的一个解，这个方程有两个解
$X={a_1a_2...a_{n-1}+-sqrt(a_1^2a_2^2...a_{n-1}^2-4(a_1^2+a_2^2+...+a_{n-1}^2))}/2$
由于这个解是最基本的，所以$a_n$必须是其中较小的解，由此我们得到
$a_n={a_1a_2...a_{n-1}-sqrt(a_1^2a_2^2...a_{n-1}^2-4(a_1^2+a_2^2+...+a_{n-1}^2))}/2>=a_{n-1}$ (1)
于是我们可以得到
$(a_1a_2...a_{n-1}-2a_{n-1})^2>=a_1^2a_2^2...a_{n-1}^2-4(a_1^2+a_2^2+...+a_{n-1}^2)$而且$a_1a_2...a_{n-2}>=2$
即
$2a_{n-1}^2+a_1^2+...+a_{n-2}^2>=a_1a_2...a_{n-2}a_{n-1}^2$
所以
$a_{n-1}<=sqrt({a_1^2+a_2^2+...+a_{n-2}^2}/{a_1a_2....a_{n-2}-2})$ (2)
由于$2a_{n-1}^2+a_1^2+...+a_{n-2}^2<=na_{n-1}^2$，所以我们还可以得到
$2<=a_1a_2...a_{n-2}<=n$ (3)
而穷举满足(3)的自然数$a_1,a_2,...,a_{n-2}$非常容易，此后可以根据(2)式穷举对应的$a_{n-1}$,再根据(1)式判断是否存在对应的$a_n$

mathe · 发表于 2008-11-4 14:17:24

再贴一个用gmp程序实现的对应代码，可以用来解出对于给定n的所有基本解或前面若干个解（比如前10000个解）。

#include <stdio.h>
#include <gmp.h>
#include <math.h>
#include <set>
using namespace std;
#ifndef N
#define N 4
#endif
#define MAX_COUNT 10000
typedef struct sDataEntry{
mpz_t datas[N];
mpz_t nc;
sDataEntry(){
int i;
for(i=0;i<N;i++)mpz_init(datas[i]);
mpz_init(nc);
}
sDataEntry(const sDataEntry& d){
int i;
for(i=0;i<N;i++)mpz_init(datas[i]);
mpz_init(nc);
for(i=0;i<N;i++){
mpz_set(datas[i],d.datas[i]);
}
mpz_set(nc,d.nc);
}
~sDataEntry(){
int i;
for(i=0;i<N;i++)mpz_clear(datas[i]);
mpz_clear(nc);
}
void set_nc(){
int i;
mpz_t tmp;
mpz_init(tmp);
mpz_set_ui(nc,0);
for(i=0;i<N;i++){
mpz_mul(tmp,datas[i],datas[i]);
mpz_add(nc,nc,tmp);
}
mpz_clear(tmp);
}
bool operator<(const struct sDataEntry& s)const{
int i;
int c=mpz_cmp(nc,s.nc);
if(c!=0)return c<0;
for(i=N-1;i>=0;i--){
int c=mpz_cmp(datas[i],s.datas[i]);
if(c!=0)return c<0;
}
return false;
}
bool operator==(const struct sDataEntry& s)const{
int i;
for(i=N-1;i>=0;i--){
int c=mpz_cmp(datas[i],s.datas[i]);
if(c!=0)return false;
}
return true;
}
}DataEntry;
set<DataEntry> results;
DataEntry t;
int td[N];
void search_last()
{
long long mul=1LL;
long long sum=0LL;
int i;
if(td[N-3]==1)return;///No solution available
for(i=0;i<N-2;i++){
mul*=td[i];
sum+=(long long)td[i]*td[i];
}
if(mul<=2LL)return;///No solution available;
long long d=sum/(mul-2LL);
int up=(int)sqrt(d+0.01);
for(i=td[N-3];i<=up;i++){
td[N-2]=i;
long long a=sum+(long long)i*i;
long long b=mul*i;
long long c=b*b-4*a;
if(c<0)continue;
int u=(int)sqrt((double)c);
if(u*(long long)u<c)u=u+1;
if(u*(long long)u==c){///Found a solution now
if((b-u)&1)continue;///Invalid
td[N-1]=(int)((b-u)/2);
int j;
for(j=0;j<N;j++){
mpz_set_ui(t.datas[j],td[j]);
}
t.set_nc();
results.insert(t);
}
}
}
void search(int last)
{
int i;
if(last==N-2){
search_last();
return;
}
long long mul=1LL;
long long mul2;
double ub;
for(i=0;i<last;i++){
mul*=td[i];
}
if(mul>N)return;
ub=pow((double)N/(double)mul,1.0/(N-2-last))+0.01;
i=1;if(last>0)i=td[last-1];
for(;i<=(int)ub;i++){
td[last]=i;
search(last+1);
}
}
int tcount=0;
void dump(const DataEntry& d)
{
int i;
for(i=0;i<N;i++){
mpz_out_str(stdout,10,d.datas[i]);
printf("\t");
}
printf("\n");
}
void dumpall()
{
set<DataEntry>::iterator it;
for(it=results.begin();it!=results.end();++it){
dump(*it);
}
}
bool gens(const DataEntry& d)///return false when tcount reaches MAX_COUNT
{
mpz_t sum;
mpz_t tmp;
mpz_init(sum);mpz_init(tmp);
int i,j;
mpz_set_ui(sum,0);
for(i=0;i<N;i++){
mpz_mul(tmp,d.datas[i],d.datas[i]);
mpz_add(sum,sum,tmp);
}
for(i=0;i<N;i++){
if(i>0&&mpz_cmp(d.datas[i],d.datas[i-1])==0)
continue;
mpz_mul(tmp,d.datas[i],d.datas[i]);
mpz_sub(tmp,sum,tmp);
mpz_div(tmp,tmp,d.datas[i]);
if(mpz_cmp(tmp,d.datas[N-1])>0){
for(j=0;j<i;j++){
mpz_set(t.datas[j],d.datas[j]);
}
for(j=i;j<N-1;j++){
mpz_set(t.datas[j],d.datas[j+1]);
}
mpz_set(t.datas[j],tmp);
t.set_nc();
set<DataEntry>::iterator lit;
lit=results.find(t);
if(lit==results.end()){
results.insert(t);
tcount++;
}
}
}
mpz_clear(sum);mpz_clear(tmp);
return tcount<MAX_COUNT;
}
void gen_more()
{
set<DataEntry>::iterator it;
for(it=results.begin();it!=results.end();++it){
if(!gens(*it))
break;
}
mpz_t maxt;
mpz_init(maxt);
mpz_set(maxt,it->nc);
set<DataEntry>::iterator nit;
for(nit=results.begin();nit!=results.end();++nit){
if(mpz_cmp(nit->datas[N-1],maxt)<=0)
dump(*nit);
}
mpz_clear(maxt);
}
int main()
{
int i;
search(0);
#ifdef IONLY
dumpall();
#else
gen_more();
#endif
}

复制代码

我们知道n=2时方程没有非平凡解（也就是只有全0的解）
而利用这个程序还可以知道在100以内，在n分别为
2,6,9,11,12,15,16,18,20,21,24,29,32,33,36,41,42,45,48,50,51,56,57,60,66,72,76,77,81,82,84,90,96,99时没有非平凡解。而对于100以内（包括100）的其他n都有非平凡解。

mathe · 发表于 2008-11-5 08:04:30

附件中文件给出了10000以内所有对应方程没有非平凡解的n,总共1360个,毫无规律,是否同素数的分布规律有得一拼

,而且可能频率也差不多

sqm.tar.gz (3.14 KB, 下载次数: 5)

无心人 · 发表于 2008-11-5 08:34:05

比如n = 2
应该证明确实没非平凡解
但是n = 6是否在整个正整数区域均无非平凡解呢？

mathe · 发表于 2008-11-5 08:45:16

原帖由 无心人 于 2008-11-5 08:34 发表
比如n = 2
应该证明确实没非平凡解
但是n = 6是否在整个正整数区域均无非平凡解呢？

是的,已经用计算机证明了

无心人 · 发表于 2008-11-5 10:10:55

我的意思是严格证明
毕竟某些方程的最小整数解是
很大的

mathe · 发表于 2008-11-5 11:44:51

是严格证明。前面已经通过数学方法将最小解规约到一个很小的范围。然后用计算机穷举那个小范围就可以了。

无心人 · 发表于 2008-11-5 14:08:59

知道了

hujunhua · 发表于 2010-1-30 13:29:40

又看到一棵3叉树，也许会于我有用，收下了。

wjy · 发表于 2021-10-3 22:42:38

oeis.org/A146974的第一个链接被标记为无法打开，楼主能修复它吗

账号		自动登录	找回密码
密码			欢迎注册

[转载] 一个多变量丢番图方程

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。