void Gauss(int n, double a[][4], int ell[])
{
	double s[4];
	int i, j, k, tmp;
	double r, rmax, smax, xmult;
	for (i=0;i<n;i++)
	{
		ell[i] = i;
		smax = 0;
		for (j=0;j<n;j++)
		{
			smax = max(smax, a[i][j]);
		}
		s[i] = smax;
	}
	for (k=0;k<n-1;k++)
	{
		rmax = 0;
		for (i=k;i<n;i++)
		{
			r = a[ell[i]][k]/s[ell[i]];
			r = fabs(r);
			if (r> rmax)
			{
				rmax = r;
				j = i;
			}
		}
		tmp = ell[j];
		ell[j] = ell[k];
		ell[k] = tmp;
		for(i=k+1;i<n;i++)
		{	
			xmult = a[ell[i]][k]/a[ell[k]][k];
			a[ell[i]][k] = xmult;
			for (j=k+1;j<n;j++)
			{
				a[ell[i]][j] = a[ell[i]][j] = xmult * a[ell[k]][j];
			}
		}
	}
}




void Solve(int n, double a[][4], int ell[], double b[], double x[])
{
	int i, k, j;
	double sum;
	for (k=0; k<n-1; k++)
	{
		for (i=k+1;i<n;i++)
		{
			b[ell[i]] = (b[ell[i]] - (a[ell[i]][k] * b[ell[k]]));
	
		}
	}
	x[n] = b[ell[n]]/a[ell[n]][n];
	for (i=n-2; i>=0; i--)
	{
		sum = b[ell[i]];
		for(j = i+1; j<n; j++)
		{
			sum = sum - a[ell[i]][j] * x[j];
		}
		x[i] = sum/a[ell[i]][i];
	}
	

}