求解三元一次不定方程

#include "stdio.h"
#include "math.h"
#include "time.h"
#include "stdlib.h"

const int iAgentDim=3;
const int iRangL=1;
const int iRangR=10000;
const int a=5;
const int b=6;
const int c=7;
const int n=3000;

const int iPSONum=20;

int iStep=1000;
double w=0.9;
const double delta1=1;
const double delta2=1;

#define rnd(low,upper) ((rand()/(double)RAND_MAX)*((upper)-(low))+(low))
#define rdint(i) (rand()%(int)(i))
#define rndint(a,b) (rdint((int)(b)-(int)(a)+1)+(int)(a))
int gbest[iAgentDim];

class Agent
{
public:
int dpos[iAgentDim];
int dpbest[iAgentDim];
int dv[iAgentDim];
int m_dFitness;
int m_dBestfitness;

Agent()
{
srand((unsigned)(time(NULL)+rand()));
int i=0;
for(i=0;i<iAgentDim;i++)
{
dpos[i]=rndint(iRangL,iRangR);
dv[i]=dpbest[i]=dpos[i];
}
}

void UpdateFitness()
{

m_dFitness=abs(a*dpos[0]+b*dpos[1]-c*dpos[2]-n);
if(m_dFitness<m_dBestfitness)
{
m_dBestfitness=m_dFitness;
int i=0;
for(;i<iAgentDim;i++)
{
dpbest[i]=dpos[i];
}
}
}

void UpdatePos()
{
int i=0;
for(;i<iAgentDim;i++)
{
dv[i]=(int)(w*dv[i]+delta1*rnd(0,1)*(dpbest[i]-dpos[i])+delta2*rnd(0,1)*(gbest[i]-dpos[i]));
dpos[i]+=dv[i];
}

}
};

class PSO
{
private:
Agent agents[iPSONum];
int m_dBestFitness;
int m_iTempPos;
public:
void Init();
void Search();
int GetBestFitness(){return m_dBestFitness;}

};
void PSO::Search()
{
int k=0;
while(k<iStep)
{
m_iTempPos=9999;
int i;
for(i=0;i<iPSONum;i++)
{
if(m_dBestFitness>agents[i].m_dBestfitness)
{
m_dBestFitness=agents[i].m_dBestfitness;
m_iTempPos=i;
}
}
if(m_iTempPos!=9999)
{
int j;
for(j=0;j<iAgentDim;j++)
{
gbest[j]=agents[m_iTempPos].dpos[j];
}

}

for(i=0;i<iPSONum;i++)
{
agents[i].UpdatePos();
agents[i].UpdateFitness();
}
k++;
}

}

void PSO::Init()
{
int i=0;
m_dBestFitness=1000;
srand((unsigned)(time(NULL)+rand()));
for(;i<iPSONum;i++)
{
agents[i].m_dBestfitness=1000;
agents[i].UpdateFitness();
}
}

int main()
{
PSO pso;
pso.Init();
pso.Search();
for(int i=0;i<iAgentDim; i++)
{
if(gbest[i]<=0)main();
}
if(pso.GetBestFitness()==0)
{
for(int i=0;i<iAgentDim;i++)printf(" %d ",gbest[i]);
}
return 0;
}
很简单的PSO算法，拿来改了改，你们谁有好的方法也贴出来啊！！！！

Input
The input consists of several test cases. Each test case consists of four integers x, y, z, n in a single line. (0 ≤ x, y, z, n ≤ 200)
The input will terminated by the end of file.

Output
For each input test case, if nuanran can get n stones, output a single line with the times of getting x stones, the times of getting y stones and the times of throwing z stones, separated by spaces. Otherwise output -1 in a single line.
Please note your answer must fit into 32 signed bits.

Sample Input
2 4 3 5
2 4 2 5
Sample Output
2 1 1
-1

The Chinese remainder theorem is the key to solve

a_1 x_1 + a_2 x_2 + ... + a_r x_r = n

for integer solutions. You may want to be familiar with the theorem.

For the case r=2, i.e., a_1 x_1 + a_2 x_2 = n, the extended Euclid algorithm is your friend. I implemented it inside the code below.

For the case r=3, by using the extended Euclid algorithm twice, you arrive at the solution.

For the case r>3, you can still use the algorithm r-1 times to solve it. But it become so cubersome that you want to use the Chinese remainder theorem and its corresponding algorithm instead.

I am not using the theorem in the code.

=================================================

#include <iostream>
using namespace std;

int gcd(int a, int b, int& x, int& y)
{
int x2=1, y2=0, t, q, a2=a, b2=b;

x=0;
y=1;
while(b)
{
t = b;
q = a / b;
b = a% b;
a = t;

t = x;
x = x2-q*x;
x2 = t;

t = y;
y = y2-q*y;
y2 = t;
}

x=x2;
y=y2;

// make x>=0, y<=0
x=x+b2/a;
y=y-a2/a;

return a;
}

int main()
{
int a, b, c, x, y, z, n, v, w, x1, y1, x2, y2, t;

while(cin>>x>>y>>z>>n)
{
v=gcd(x, y, x1, y1);
w=gcd(v, z, x2, y2);

if(n%w)
cout<<-1<<endl;
else
{
t=n/w;

// calculate one solution: here a>=0, b<=0, c>=0
a=t*x2*x1;
b=t*x2*y1;
c=t*y2;

// adjust b and c to make b non-negative
c=-c;
t=-b/z+1;
b+=z*t;
c+=t*y;

cout<<a<<" "<<b<<" "<<c<<endl;
}
}

return 0;
}