5个不同的自然数， 它们中的任意二个之和都是平方数

[northwolves]

这是最大数最小的一组解： 7442, 28658, 148583, 177458, 763442

问题：

试列出10000000内的所有满足条件的解

6个自然数有没有解？

[northwolves]

非常感谢郭老师和mathe兄弟的关注。

平方数末尾二位：
00 01 04 09 16 21 24 25 29 36 41 44 49 56 61 64 69 76 81 84 89 96

我利用以上平方数的特点，做了一个很笨的循环，花了25分钟在 1000000内找到了http://www.iwr.uni-heidelberg.de/groups/ngg/People/winckler/PU/s051.html 提供的
两个解：

7442 28658 148583 177458 763442
32018 104882 188882 559343 956018

[northwolves]

a+b=x1*x1
a+c=x2*x2
b+c=x3*x3
则
b-a=x3^2-x2^2
c-b=x2^2-x1^2
c-a=x3^2-x1^2

20000000的平方根取整4472，10000000内的任意两平方数之差共有 3068382 个，用一个字节数组来表示，1 表示可这样表示，0表示非这种形式，做不下去了。

[northwolves]

然后我们在关系表中找出所有的前三个数字相同，但是最后一项不同的数据对，比如
<a,b,c;N>
<a,b,c;M>
最后计算
N+M-(a^2+b^2+c^2)
如果这个数是完全平方数，那么我们已经找到一组5个数的情况。
-------------
mathe兄，这个计算量应该很可观

[uzone]

关注...

[mathe]

还有一种特殊情况没有考虑，所有数可以乘上公因子,最后结果如下：
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>

#define abs(x) ((x)>0?(x):-(x))
#define HALFLIMIT 6325
#define LIMIT (HALFLIMIT*HALFLIMIT)
#define SUMBOUND (LIMIT*25)
//The limitation is 40000000*25
#ifdef WIN32
typedef __int64 longlong;
#else
typedef long long longlong;
#endif
#define MAXNUMS 73000000
int not_prime[LIMIT];
int pcount;
#define primelist not_prime

int data[MAXNUMS][4];
int dcount;
int cmpdata(const void *p, const void *q){
const int *pi=(const int *)p;
const int *qi=(const int *)q;
int i;
for(i=0;i<4;i++){
if(pi[i]<qi[i])return -1;
if(pi[i]>qi[i])return 1;
}
return 0;
}

int cmpsd(const void *p, const void *q){
const int *pi=(const int *)p;
const int *qi=(const int *)q;
int i;
for(i=0;i<3;i++){
if(pi[i]<qi[i])return -1;
if(pi[i]>qi[i])return 1;
}
return 0;
}

void pushdata(int X,int Y,int Z,int W)
{
int t;
data[dcount][0]=X;
data[dcount][1]=Y;
data[dcount][2]=Z;
data[dcount][3]=W;
dcount++;
}

class cn{
longlong r,i;
public:
cn(longlong r, longlong i){this->r=r;this->i=i;}
cn com()const{cn x(r,-i);return x;}
longlong sqrmod()const{return r*r+i*i;}
longlong rel()const{return r;}
longlong img()const{return i;}
cn operator+(const cn& x)const{cn y(r+x.r,i+x.i);return y;}
cn operator-(const cn& x)const{cn y(r-x.r,i-x.i);return y;}
cn operator*(const cn& x)const{cn y(r*x.r-i*x.i,r*x.i+i*x.r);return y;}
cn div(longlong x)const;
bool iszero()const{return r==0&&i==0;}
};


cn cn::div(longlong x)const{
longlong u,v;
u=r%x;
v=i%x;
if(u<0)u+=x;
if(v<0)v+=x;
if(u>x/2)u-=x;
if(v>x/2)v-=x;
u=(r-u)/x;
v=(i-v)/x;
return cn(u,v);
}

cn round(const cn& x, const cn& y){
cn z=x*y.com();
longlong r=y.sqrmod();
return z.div(r);
}

int *A, *B;
cn factor(const cn& x, const cn& y){
if(x.iszero())return y;
return factor(y-round(y,x)*x,x);
}

int powmod(longlong a, longlong n, int p){
longlong mul=a%p;
longlong base=1;
while(n){
if(n&1){
base*=mul;
base%=p;
}
mul*=mul;
mul%=p;
n>>=1;
}
return (int)base;
}

int rm1(int p){
int i;
int t=(p-1)/4;
for(i=2;i<p;i++){
int u=powmod(i,t,p);
if(((longlong)u*u+1)%p==0)return u;
}
}

cn defactor(int p){
int v=rm1(p);
cn x(v,1);
cn y(p,0);
return factor(x,y);
}

#define MAX_PRIME_COUNT 20
int pstack[MAX_PRIME_COUNT];
int pindex[MAX_PRIME_COUNT];
int prod[MAX_PRIME_COUNT];
int max_part;
longlong total_count;
struct pair{
int S,L;
}SL[50];
int cmpsl(const void *p, const void *q){
const struct pair *ps=(const struct pair *)p;
const struct pair *qs=(const struct pair *)q;
if(ps->S<qs->S)return -1;
if(ps->S>qs->S)return 1;
return 0;
}

int sc;
void insert_num(int lindex[MAX_PRIME_COUNT], int step){
cn num(1,0);
int i,j;
int u,v;
for(i=0;i<=step;i++){
cn p(A[pstack[i]],B[pstack[i]]);
cn q=p.com();
for(j=0;j<lindex[i];j++){
num=num*p;
}
for(;j<pindex[i];j++){
num=num*q;
}
}
u=abs((int)num.rel());
v=abs((int)num.img());
if(v<u){
i=v;
v=u;
u=i;
}
SL[sc].S=u*u;SL[sc].L=v*v;
sc++;
}

[mathe]

//product in M=prod[step+1] or M=2*prod[step+1];
//data in pstack^pindex
void trynum(int step){
int odd_index;
int count,i,j,k;
int lindex[MAX_PRIME_COUNT];
int need2=(longlong)2*prod[step+1]<SUMBOUND;
int need4=(longlong)4*prod[step+1]<SUMBOUND;
int need8=(longlong)8*prod[step+1]<SUMBOUND;
int M=prod[step+1];
if(step==0&&pindex[0]<5)return;//One prime number
if(step==1&&pindex[0]==1&&pindex[1]==1)return;//Two prime numbers
odd_index=0;
count=1;
for(i=0;i<=step;i++){
count*=pindex[i]+1;
if(pindex[i]&1)odd_index=i;
}
sc=0;
memset(lindex,0,sizeof(lindex));
do{
insert_num(lindex,step);
for(i=step;i>=0;i--){
if((i==odd_index&&lindex[i]==(pindex[i]-1)/2)||
lindex[i]==pindex[i]){
lindex[i]=0;
continue;
}
lindex[i]++;
break;
}
}while(i>=0);
if(count&1){
for(j=1;j<=step;j++){
lindex[j-1]=pindex[j]/2;
odd_index=j;
do{
insert_num(lindex,step);
for(i=step;i>=j;i--){
if((i==odd_index&&lindex[i]==(pindex[i]-1)/2)||
lindex[i]==pindex[i]){
lindex[i]=0;
continue;
}
lindex[i]++;
break;
}
}while(i>=j);
}
}
qsort(SL,sc,sizeof(SL[0]),cmpsl);
for(i=0;i<sc-2;i++){
for(j=i+1;j<sc-1;j++){
for(k=j+1;k<sc;k++){
int N,X,Y,Z,W;
if(SL[i].S+SL[j].S>SL[k].S){
N=SL[i].S+SL[j].S+SL[k].S;
if(need4&&(N&1)){
X=2*N-4*SL[k].S;
Y=2*N-4*SL[j].S;
Z=2*N-4*SL[i].S;
W=prod[step+1]*4-2*N;
pushdata(X,Y,Z,W);
}else if(!(N&1)){
N/=2;
X=N-SL[k].S;
Y=N-SL[j].S;
Z=N-SL[i].S;
W=prod[step+1]-N;
pushdata(X,Y,Z,W);
}
}
if(SL[i].S+SL[j].S>SL[k].L&&SL[k].L>SL[k].S){
N=SL[i].S+SL[j].S+SL[k].L;
if(need4&&(N&1)){
X=2*N-4*SL[k].L;
Y=2*N-4*SL[j].S;
Z=2*N-4*SL[i].S;
W=prod[step+1]*4-2*N;
if(W>Z)
pushdata(X,Y,Z,W);
}else if(!(N&1)){
N/=2;
X=N-SL[k].L;
Y=N-SL[j].S;
Z=N-SL[i].S;
W=prod[step+1]-N;
if(W>Z)
pushdata(X,Y,Z,W);
}
}
}
}
}
if(need2)for(i=0;i<sc-2;i++){
for(j=i+1;j<sc-1;j++){
for(k=j+1;k<sc;k++){
struct pair s1,s2,s3;
struct pair tmp;
int N,X,Y,Z,W;
s1.S=SL[i].L-SL[i].S,s2.S=SL[j].L-SL[j].S,s3.S=SL[k].L-SL[k].S;
s1.L=SL[i].L+SL[i].S,s2.L=SL[j].L+SL[j].S,s3.L=SL[k].L+SL[k].S;
if(s1.S>s2.S){
tmp=s1;s1=s2;s2=tmp;
}
if(s1.S>s3.S){
tmp=s1;s1=s3;s3=tmp;
}
if(s2.S>s3.S){
tmp=s2;s2=s3;s3=tmp;
}
if(s1.S+s2.S<=s3.S)
continue;
N=s1.S+s2.S+s3.S;
if(need8&&(N&1)){
X=2*N-4*s3.S;
Y=2*N-4*s2.S;
Z=2*N-4*s1.S;
W=prod[step+1]*8-2*N;
pushdata(X,Y,Z,W);
}else if(!(N&1)){
N/=2;
X=N-s3.S;
Y=N-s2.S;
Z=N-s1.S;
W=prod[step+1]*2-N;
pushdata(X,Y,Z,W);
}
if(s1.S+s2.S<=s3.L||s3.L==s3.S)
continue;
N=s1.S+s2.S+s3.L;
if(need8&&(N&1)){
X=2*N-4*s3.L;
Y=2*N-4*s2.S;
Z=2*N-4*s1.S;
W=prod[step+1]*8-2*N;
if(W>Z)
pushdata(X,Y,Z,W);
}else if(!(N&1)){
N/=2;
X=N-s3.L;
Y=N-s2.S;
Z=N-s1.S;
W=prod[step+1]*2-N;
if(W>Z)
pushdata(X,Y,Z,W);
}
}
}
}

}

void enumerate(int step, int start){
int i,j;
for(i=start;i<pcount;i++){
int p=primelist[i];
pstack[step]=i;
prod[step+1]=prod[step];
for(j=1;prod[step+1]*(longlong)p<SUMBOUND;j++){
pindex[step]=j;
prod[step+1]*=p;
trynum(step);
enumerate(step+1,i+1);
}
}
}

int issqr(int x){
double d=sqrt(x);
int n=(int)d;
if(n*n==x)return 1;
n=n++;
if(n*n==x)return 1;
return 0;
}

int comfac2(int a, int b){
if(a==0)return b;
return comfac2(b%a,a);
}

int comfac3(int d[3]){
int c=comfac2(d[0],d[1]);
return comfac2(c,d[2]);
}

int main(){
int i,j;
not_prime[0]=not_prime[1]=1;
for(i=2;i<HALFLIMIT;i++){
if(!not_prime[i]){
if(i%4==1){
primelist[pcount++]=i;
}
for(j=i*i;j<LIMIT;j+=i)
not_prime[j]=1;
}
}
A = new int[pcount];
B = new int[pcount];
for(i=0;i<pcount;i++){
cn factors=defactor(primelist[i]);
A[i]=abs((int)factors.rel());
B[i]=abs((int)factors.img());
}
prod[0]=1;
enumerate(0,0);
fprintf(stderr,"4 pairs count=%d\n",dcount);
qsort(data,dcount,sizeof(data[0]),cmpdata);
for(j=0,i=1;i<dcount;i++){
int k;
int h;
for(k=0;k<3;k++){
if(data[i][k]!=data[j][k])break;
}
if(k<3){
j=i;
continue;
}
for(k=j;k<i;k++){
if(issqr(data[k][3]+data[i][3])){
int t;
for(t=0;t<4;t++){
printf("%d,",data[k][t]);
}
printf("%d\n",data[i][3]);
}
}
h=comfac3(data[i]);
if(h>1){
int u;
for(u=2;u<=h/2;u++){
if(h%u==0){
int d[4];
int *r;
d[0]=data[i][0]/u;
d[1]=data[i][1]/u;
d[2]=data[i][2]/u;
r=(int *)bsearch(d,data,dcount,sizeof(data[0]),cmpsd);
if(r){
while(cmpsd(d,r-4)==0)r-=4;
while(cmpsd(d,r)==0){
int t;
for(t=0;t<4;t++){
printf("%d,",data[i][t]);
}
printf("%d\n",r[3]*u);
r+=4;
}
}
}
}
}
}
}

[mathe]

程序大概运行1分半钟。
但是不知道为什么会找到一些重复解（怀疑还有些bug)
让程序过滤掉重复解后，可以找到总共6916个解。

[mathe]

最后一步弄错了，修改一下main函数：
int main(){
int i,j;
not_prime[0]=not_prime[1]=1;
for(i=2;i<HALFLIMIT;i++){
if(!not_prime[i]){
if(i%4==1){
primelist[pcount++]=i;
}
for(j=i*i;j<LIMIT;j+=i)
not_prime[j]=1;
}
}
A = new int[pcount];
B = new int[pcount];
for(i=0;i<pcount;i++){
cn factors=defactor(primelist[i]);
A[i]=abs((int)factors.rel());
B[i]=abs((int)factors.img());
}
prod[0]=1;
enumerate(0,0);
fprintf(stderr,"4 pairs count=%d\n",dcount);
qsort(data,dcount,sizeof(data[0]),cmpdata);
for(j=0,i=1;i<dcount;i++){
int k;
int h;
for(k=0;k<3;k++){
if(data[i][k]!=data[j][k])break;
}
if(k<3){
j=i;
continue;
}
if(data[i][3]==data[j][3]){
data[i][3]=-1;
continue;
}
for(k=j;k<i;k++){
if(data[k][3]==-1)continue;
if(issqr(data[k][3]+data[i][3])){
int t;
for(t=0;t<4;t++){
printf("%d,",data[k][t]);
}
printf("%d\n",data[i][3]);
}
}
h=comfac3(data[i]);
if(h>1){
int u;
for(u=2;u*u<h/2;u++){
if(h%(u*u)==0){
int d[4];
int *r;
int v=u*u;
d[0]=data[i][0]/v;
d[1]=data[i][1]/v;
d[2]=data[i][2]/v;
r=(int *)bsearch(d,data,dcount,sizeof(data[0]),cmpsd);
if(r){
while(cmpsd(d,r-4)==0)r-=4;
while(cmpsd(d,r)==0){
int t;
if(r[3]!=-1&&issqr(data[i][3]+r[3]*v)){
for(t=0;t<4;t++){
printf("%d,",data[i][t]);
}
printf("%d\n",r[3]*v);
}
r+=4;
}
}
}
}
}
}
}

[mathe]

找到一个很严重的BUG了，代码：
s1.S=SL[i].L-SL[i].S,s2.S=SL[j].L-SL[j].S,s3.S=SL[k].L-SL[k].S;
s1.L=SL[i].L+SL[i].S,s2.L=SL[j].L+SL[j].S,s3.L=SL[k].L+SL[k].S;
错了
其中，每个数应该分别弄成
sqr(sqrt(*.L)-sqrt(*.S)) (sqr(x):=x*x)
修改以后，程序运行53秒，共得出1483个不同结果。

[mathe]

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>

#define abs(x) ((x)>0?(x):-(x))
#define HALFLIMIT 6325
#define LIMIT (HALFLIMIT*HALFLIMIT)
#define SUMBOUND (LIMIT*25)
//The limitation is 40000000*25
#ifdef WIN32
typedef __int64 longlong;
#else
typedef long long longlong;
#endif
#define MAXNUMS 73000000
int not_prime[LIMIT];
int pcount;
#define primelist not_prime

int data[MAXNUMS][4];
int dcount;
int cmpdata(const void *p, const void *q){
const int *pi=(const int *)p;
const int *qi=(const int *)q;
int i;
for(i=0;i<4;i++){
if(pi[i]<qi[i])return -1;
if(pi[i]>qi[i])return 1;
}
return 0;
}

int cmpsd(const void *p, const void *q){
const int *pi=(const int *)p;
const int *qi=(const int *)q;
int i;
for(i=0;i<3;i++){
if(pi[i]<qi[i])return -1;
if(pi[i]>qi[i])return 1;
}
return 0;
}

void pushdata(int X,int Y,int Z,int W)
{
int t;
data[dcount][0]=X;
data[dcount][1]=Y;
data[dcount][2]=Z;
data[dcount][3]=W;
dcount++;
}

class cn{
longlong r,i;
public:
cn(longlong r, longlong i){this->r=r;this->i=i;}
cn com()const{cn x(r,-i);return x;}
longlong sqrmod()const{return r*r+i*i;}
longlong rel()const{return r;}
longlong img()const{return i;}
cn operator+(const cn& x)const{cn y(r+x.r,i+x.i);return y;}
cn operator-(const cn& x)const{cn y(r-x.r,i-x.i);return y;}
cn operator*(const cn& x)const{cn y(r*x.r-i*x.i,r*x.i+i*x.r);return y;}
cn div(longlong x)const;
bool iszero()const{return r==0&&i==0;}
};


cn cn::div(longlong x)const{
longlong u,v;
u=r%x;
v=i%x;
if(u<0)u+=x;
if(v<0)v+=x;
if(u>x/2)u-=x;
if(v>x/2)v-=x;
u=(r-u)/x;
v=(i-v)/x;
return cn(u,v);
}

cn round(const cn& x, const cn& y){
cn z=x*y.com();
longlong r=y.sqrmod();
return z.div(r);
}

int *A, *B;
cn factor(const cn& x, const cn& y){
if(x.iszero())return y;
return factor(y-round(y,x)*x,x);
}

int powmod(longlong a, longlong n, int p){
longlong mul=a%p;
longlong base=1;
while(n){
if(n&1){
base*=mul;
base%=p;
}
mul*=mul;
mul%=p;
n>>=1;
}
return (int)base;
}

int rm1(int p){
int i;
int t=(p-1)/4;
for(i=2;i<p;i++){
int u=powmod(i,t,p);
if(((longlong)u*u+1)%p==0)return u;
}
}

cn defactor(int p){
int v=rm1(p);
cn x(v,1);
cn y(p,0);
return factor(x,y);
}

#define MAX_PRIME_COUNT 20
int pstack[MAX_PRIME_COUNT];
int pindex[MAX_PRIME_COUNT];
int prod[MAX_PRIME_COUNT];
int max_part;
longlong total_count;
struct pair{
int S,L,s,l;
}SL[50];
int cmpsl(const void *p, const void *q){
const struct pair *ps=(const struct pair *)p;
const struct pair *qs=(const struct pair *)q;
if(ps->S<qs->S)return -1;
if(ps->S>qs->S)return 1;
return 0;
}

int sc;
void insert_num(int lindex[MAX_PRIME_COUNT], int step){
cn num(1,0);
int i,j;
int u,v;
for(i=0;i<=step;i++){
cn p(A[pstack[i]],B[pstack[i]]);
cn q=p.com();
for(j=0;j<lindex[i];j++){
num=num*p;
}
for(;j<pindex[i];j++){
num=num*q;
}
}
u=abs((int)num.rel());
v=abs((int)num.img());
if(v<u){
i=v;
v=u;
u=i;
}
SL[sc].S=u*u;SL[sc].L=v*v;
SL[sc].s=u;SL[sc].l=v;
sc++;
}

//product in M=prod[step+1] or M=2*prod[step+1];
//data in pstack^pindex

[mathe]

void trynum(int step){
int odd_index;
int count,i,j,k;
int lindex[MAX_PRIME_COUNT];
int need2=(longlong)2*prod[step+1]<SUMBOUND;
int need4=(longlong)4*prod[step+1]<SUMBOUND;
int need8=(longlong)8*prod[step+1]<SUMBOUND;
int M=prod[step+1];
if(step==0&&pindex[0]<5)return;//One prime number
if(step==1&&pindex[0]==1&&pindex[1]==1)return;//Two prime numbers
odd_index=0;
count=1;
for(i=0;i<=step;i++){
count*=pindex[i]+1;
if(pindex[i]&1)odd_index=i;
}
sc=0;
memset(lindex,0,sizeof(lindex));
do{
insert_num(lindex,step);
for(i=step;i>=0;i--){
if((i==odd_index&&lindex[i]==(pindex[i]-1)/2)||
lindex[i]==pindex[i]){
lindex[i]=0;
continue;
}
lindex[i]++;
break;
}
}while(i>=0);
if(count&1){
for(j=1;j<=step;j++){
lindex[j-1]=pindex[j]/2;
odd_index=j;
do{
insert_num(lindex,step);
for(i=step;i>=j;i--){
if((i==odd_index&&lindex[i]==(pindex[i]-1)/2)||
lindex[i]==pindex[i]){
lindex[i]=0;
continue;
}
lindex[i]++;
break;
}
}while(i>=j);
}
}
qsort(SL,sc,sizeof(SL[0]),cmpsl);
for(i=0;i<sc-2;i++){
for(j=i+1;j<sc-1;j++){
for(k=j+1;k<sc;k++){
int N,X,Y,Z,W;
if(SL[i].S+SL[j].S>SL[k].S){
N=SL[i].S+SL[j].S+SL[k].S;
if(need4&&(N&1)){
X=2*N-4*SL[k].S;
Y=2*N-4*SL[j].S;
Z=2*N-4*SL[i].S;
W=prod[step+1]*4-2*N;
pushdata(X,Y,Z,W);
}else if(!(N&1)){
N/=2;
X=N-SL[k].S;
Y=N-SL[j].S;
Z=N-SL[i].S;
W=prod[step+1]-N;
pushdata(X,Y,Z,W);
}
}
if(SL[i].S+SL[j].S>SL[k].L&&SL[k].L>SL[k].S){
N=SL[i].S+SL[j].S+SL[k].L;
if(need4&&(N&1)){
X=2*N-4*SL[k].L;
Y=2*N-4*SL[j].S;
Z=2*N-4*SL[i].S;
W=prod[step+1]*4-2*N;
if(W>Z)
pushdata(X,Y,Z,W);
}else if(!(N&1)){
N/=2;
X=N-SL[k].L;
Y=N-SL[j].S;
Z=N-SL[i].S;
W=prod[step+1]-N;
if(W>Z)
pushdata(X,Y,Z,W);
}
}
}
}
}
#define SQR(x) ((x)*(x))
if(need2)for(i=0;i<sc-2;i++){
for(j=i+1;j<sc-1;j++){
for(k=j+1;k<sc;k++){
struct pair s1,s2,s3;
struct pair tmp;
int N,X,Y,Z,W;
s1.S=SQR(SL[i].l-SL[i].s),s2.S=SQR(SL[j].l-SL[j].s),s3.S=SQR(SL[k].l-SL[k].s);
s1.L=SQR(SL[i].l+SL[i].s),s2.L=SQR(SL[j].l+SL[j].s),s3.L=SQR(SL[k].l+SL[k].s);
if(s1.S>s2.S){
tmp=s1;s1=s2;s2=tmp;
}
if(s1.S>s3.S){
tmp=s1;s1=s3;s3=tmp;
}
if(s2.S>s3.S){
tmp=s2;s2=s3;s3=tmp;
}
if(s1.S+s2.S<=s3.S)
continue;
N=s1.S+s2.S+s3.S;
if(need8&&(N&1)){
X=2*N-4*s3.S;
Y=2*N-4*s2.S;
Z=2*N-4*s1.S;
W=prod[step+1]*8-2*N;
pushdata(X,Y,Z,W);
}else if(!(N&1)){
N/=2;
X=N-s3.S;
Y=N-s2.S;
Z=N-s1.S;
W=prod[step+1]*2-N;
pushdata(X,Y,Z,W);
}
if(s1.S+s2.S<=s3.L||s3.L==s3.S)
continue;
N=s1.S+s2.S+s3.L;
if(need8&&(N&1)){
X=2*N-4*s3.L;
Y=2*N-4*s2.S;
Z=2*N-4*s1.S;
W=prod[step+1]*8-2*N;
if(W>Z)
pushdata(X,Y,Z,W);
}else if(!(N&1)){
N/=2;
X=N-s3.L;
Y=N-s2.S;
Z=N-s1.S;
W=prod[step+1]*2-N;
if(W>Z)
pushdata(X,Y,Z,W);
}
}
}
}

}

void enumerate(int step, int start){
int i,j;
for(i=start;i<pcount;i++){
int p=primelist[i];
pstack[step]=i;
prod[step+1]=prod[step];
for(j=1;prod[step+1]*(longlong)p<SUMBOUND;j++){
pindex[step]=j;
prod[step+1]*=p;
trynum(step);
enumerate(step+1,i+1);
}
}
}

[mathe]

int issqr(int x){
double d=sqrt(x);
int n=(int)d;
if(n*n==x)return 1;
n=n++;
if(n*n==x)return 1;
return 0;
}

int comfac2(int a, int b){
if(a==0)return b;
return comfac2(b%a,a);
}

int comfac3(int d[3]){
int c=comfac2(d[0],d[1]);
return comfac2(c,d[2]);
}

int main(){
int i,j;
not_prime[0]=not_prime[1]=1;
for(i=2;i<HALFLIMIT;i++){
if(!not_prime[i]){
if(i%4==1){
primelist[pcount++]=i;
}
for(j=i*i;j<LIMIT;j+=i)
not_prime[j]=1;
}
}
A = new int[pcount];
B = new int[pcount];
for(i=0;i<pcount;i++){
cn factors=defactor(primelist[i]);
A[i]=abs((int)factors.rel());
B[i]=abs((int)factors.img());
}
prod[0]=1;
enumerate(0,0);
fprintf(stderr,"4 pairs count=%d\n",dcount);
qsort(data,dcount,sizeof(data[0]),cmpdata);
for(j=0,i=1;i<dcount;i++){
int k;
int h;
for(k=0;k<3;k++){
if(data[i][k]!=data[j][k])break;
}
if(k<3){
j=i;
continue;
}
/*
if(data[i][3]==data[j][3]){
data[i][3]=-1;
continue;
}*/
for(k=j;k<i;k++){
// if(data[k][3]==-1)continue;
if(issqr(data[k][3]+data[i][3])){
int t;
for(t=0;t<4;t++){
printf("%d,",data[k][t]);
}
printf("%d\n",data[i][3]);
}
}
h=comfac3(data[i]);
if(h>1){
int u;
for(u=2;u*u<h/2;u++){
if(h%(u*u)==0){
int d[4];
int *r;
int v=u*u;
d[0]=data[i][0]/v;
d[1]=data[i][1]/v;
d[2]=data[i][2]/v;
r=(int *)bsearch(d,data,dcount,sizeof(data[0]),cmpsd);
if(r){
while(cmpsd(d,r-4)==0)r-=4;
while(cmpsd(d,r)==0){
int t;
if(issqr(data[i][3]+r[3]*v)){
for(t=0;t<4;t++){
printf("%d,",data[i][t]);
}
printf("%d\n",r[3]*v);
}
r+=4;
}
}
}
}
}
}
}

[mathe]

上面的计算量不算很大，如果你通过hash表来保存数据。
其实我后来重新想了一下。没有必要保存这个关系表。
首先找到所有4个数的情况，将它们完全保存下来。
然后对它们进行排序。最后对所有前面3个数相同的数对进行检测最后两个数字之后是否是平方数就可了。
时间复杂度同4个数的情况的数目相比基本上线性。

[mathe]

上面分解平方数的代码有个BUG,我说怎么这么慢呢
int rm1(int p){
int i;
int t=(p-1)/4;
for(i=2;i<p;i++){
int u=powmod(i,t,p);
>> if((u*u+1)%p==0)return u;
改成：
if(((longlong)u*u+1)%p==0)return u;
}
}

[mathe]

我写了个程序，找出所有任何4个数之和不超过1000140625的5个数组合，
共有764组，在我的计算机上花费55秒。
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>

#define abs(x) ((x)>0?(x):-(x))
#define HALFLIMIT 6325
#define LIMIT (HALFLIMIT*HALFLIMIT)
#define SUMBOUND (LIMIT*25)
//The limitation is 40000000*25
#ifdef WIN32
typedef __int64 longlong;
#else
typedef long long longlong;
#endif
#define MAXNUMS 51000000
int not_prime[LIMIT];
int pcount;
#define primelist not_prime

int data[MAXNUMS][4];
int dcount;
int cmpdata(const void *p, const void *q){
const int *pi=(const int *)p;
const int *qi=(const int *)q;
int i;
for(i=0;i<4;i++){
if(pi[i]<qi[i])return -1;
if(pi[i]>qi[i])return 1;
}
return 0;
}

void pushdata(int X,int Y,int Z,int W)
{
int t;
if(X>Y){
t=X;X=Y;Y=t;
}
if(X>Z){
t=X;X=Z;Z=t;
}
if(Y>Z){
t=Y;Y=Z;Z=t;
}
data[dcount][0]=X;
data[dcount][1]=Y;
data[dcount][2]=Z;
data[dcount][3]=W;
dcount++;
}

class cn{
longlong r,i;
public:
cn(longlong r, longlong i){this->r=r;this->i=i;}
cn com()const{cn x(r,-i);return x;}
longlong sqrmod()const{return r*r+i*i;}
longlong rel()const{return r;}
longlong img()const{return i;}
cn operator+(const cn& x)const{cn y(r+x.r,i+x.i);return y;}
cn operator-(const cn& x)const{cn y(r-x.r,i-x.i);return y;}
cn operator*(const cn& x)const{cn y(r*x.r-i*x.i,r*x.i+i*x.r);return y;}
cn div(longlong x)const;
bool iszero()const{return r==0&&i==0;}
};


cn cn::div(longlong x)const{
longlong u,v;
u=r%x;
v=i%x;
if(u<0)u+=x;
if(v<0)v+=x;
if(u>x/2)u-=x;
if(v>x/2)v-=x;
u=(r-u)/x;
v=(i-v)/x;
return cn(u,v);
}

cn round(const cn& x, const cn& y){
cn z=x*y.com();
longlong r=y.sqrmod();
return z.div(r);
}

int *A, *B;
cn factor(const cn& x, const cn& y){
if(x.iszero())return y;
return factor(y-round(y,x)*x,x);
}

int powmod(longlong a, longlong n, int p){
longlong mul=a%p;
longlong base=1;
while(n){
if(n&1){
base*=mul;
base%=p;
}
mul*=mul;
mul%=p;
n>>=1;
}
return (int)base;
}

int rm1(int p){
int i;
int t=(p-1)/4;
for(i=2;i<p;i++){
int u=powmod(i,t,p);
if(((longlong)u*u+1)%p==0)return u;
}
}

cn defactor(int p){
int v=rm1(p);
cn x(v,1);
cn y(p,0);
return factor(x,y);
}

#define MAX_PRIME_COUNT 20
int pstack[MAX_PRIME_COUNT];
int pindex[MAX_PRIME_COUNT];
int prod[MAX_PRIME_COUNT];
int max_part;
longlong total_count;

int S[50],L[50];
int sc;
void insert_num(int lindex[MAX_PRIME_COUNT], int step){
cn num(1,0);
int i,j;
int u,v;
for(i=0;i<=step;i++){
cn p(A[pstack[i]],B[pstack[i]]);
cn q=p.com();
for(j=0;j<lindex[i];j++){
num=num*p;
}
for(;j<pindex[i];j++){
num=num*q;
}
}
u=abs((int)num.rel());
v=abs((int)num.img());
if(v<u){
i=v;
v=u;
u=i;
}
S[sc]=u*u;L[sc]=v*v;
sc++;
}

//product in N=prod[step+1] or N=2*prod[step+1];
//data in pstack^pindex
void trynum(int step){
int odd_index;
int count,i,j,k;
int lindex[MAX_PRIME_COUNT];
int need2=(longlong)2*prod[step+1]<SUMBOUND;
int need4=(longlong)4*prod[step+1]<SUMBOUND;
int need8=(longlong)8*prod[step+1]<SUMBOUND;
if(step==0&&pindex[0]<5)return;//One prime number
if(step==1&&pindex[0]==1&&pindex[1]==1)return;//Two prime numbers
odd_index=0;
count=1;
for(i=0;i<=step;i++){
count*=pindex[i]+1;
if(pindex[i]&1)odd_index=i;
}
sc=0;
memset(lindex,0,sizeof(lindex));
do{
insert_num(lindex,step);
for(i=step;i>=0;i--){
if((i==odd_index&&lindex[i]==(pindex[i]-1)/2)||
lindex[i]==pindex[i]){
lindex[i]=0;
continue;
}
lindex[i]++;
break;
}
}while(i>=0);
if(count&1){
for(j=1;j<=step;j++){
lindex[j-1]=pindex[j]/2;
odd_index=j;
do{
insert_num(lindex,step);
for(i=step;i>=j;i--){
if((i==odd_index&&lindex[i]==(pindex[i]-1)/2)||
lindex[i]==pindex[i]){
lindex[i]=0;
continue;
}
lindex[i]++;
break;
}
}while(i>=j);
}
}
for(i=0;i<sc-2;i++){
for(j=i+1;j<sc-1;j++){
for(k=j+1;k<sc;k++){
int N,X,Y,Z,W;
if(S[i]+S[j]<=S[k]||S[i]+S[k]<=S[j]||S[j]+S[k]<=S[i])
continue;
N=S[i]+S[j]+S[k];
if(need4&&(N&1)){
X=2*N-4*S[i];
Y=2*N-4*S[j];
Z=2*N-4*S[k];
W=prod[step+1]*4-2*N;
pushdata(X,Y,Z,W);
}else if(!(N&1)){
N/=2;
X=N-S[i];
Y=N-S[j];
Z=N-S[k];
W=prod[step+1]-N;
pushdata(X,Y,Z,W);
}
}
}
}
if(need2)for(i=0;i<sc-2;i++){
for(j=i+1;j<sc-1;j++){
for(k=j+1;k<sc;k++){
int i1,i2,i3;
int N,X,Y,Z,W;
i1=L[i]-S[i],i2=L[j]-S[j],i3=L[k]-S[k];
if(i1+i2<=i3||i1+i3<=i2||i2+i3<=i1)
continue;
N=i1+i2+i3;
if(need8&&(N&1)){
X=2*N-4*i1;
Y=2*N-4*i2;
Z=2*N-4*i3;
W=prod[step+1]*8-2*N;
pushdata(X,Y,Z,W);
}else if(!(N&1)){
N/=2;
X=N-i1;
Y=N-i2;
Z=N-i3;
W=prod[step+1]*2-N;
pushdata(X,Y,Z,W);
}
}
}
}

}

void enumerate(int step, int start){
int i,j;
for(i=start;i<pcount;i++){
int p=primelist[i];
pstack[step]=i;
prod[step+1]=prod[step];
for(j=1;prod[step+1]*(longlong)p<SUMBOUND;j++){
pindex[step]=j;
prod[step+1]*=p;
trynum(step);
enumerate(step+1,i+1);
}
}
}

int issqr(int x){
double d=sqrt(x);
int n=(int)d;
if(n*n==x)return 1;
n=n++;
if(n*n==x)return 1;
return 0;
}

int main(){
int i,j;
not_prime[0]=not_prime[1]=1;
for(i=2;i<HALFLIMIT;i++){
if(!not_prime[i]){
if(i%4==1){
primelist[pcount++]=i;
}
for(j=i*i;j<LIMIT;j+=i)
not_prime[j]=1;
}
}
A = new int[pcount];
B = new int[pcount];
for(i=0;i<pcount;i++){
cn factors=defactor(primelist[i]);
A[i]=abs((int)factors.rel());
B[i]=abs((int)factors.img());
}
prod[0]=1;
enumerate(0,0);
qsort(data,dcount,sizeof(data[0]),cmpdata);
for(j=0,i=1;i<dcount;i++){
int k;
for(k=0;k<3;k++){
if(data[i][k]!=data[j][k])break;
}
if(k<3){
j=i;
continue;
}
for(k=j;k<i;k++){
if(issqr(data[k][3]+data[i][3])){
int t;
for(t=0;t<4;t++){
printf("%d,",data[k][t]);
}
printf("%d\n",data[i][3]);
}
}
}
}

[mathe]

发现还漏了一部分解。
对于N=a^2+f^2=b^2+e^2=c^2+d^2
a<b<c<d<e<f
除了使用
S=(a^2+b^2+c^2)/2
外，我们还应该测试
S=(a^2+b^2+d^2)/2这组结果，
只是这组结果计算出来的W可能比X,Y,Z中某个数小，这种情况简单抛弃就可以了。
经过这样的修改后就可以找到所有结果了。

[heixia108]

mark

[softwarewander]

牛！

[ahjoe]

先学好数学再来...

[jinwei1984]

mathe 真是大哥级人物啊