I am trying to solve http://www.spoj.com/problems/GS . Here is my approach.

Let u, v denote the start and end vertices. [0 based index]

Let E[i] = expected number of roads to cross while going from vertex i to vertex v [0 <= i <= N-1]

Let preferred value of an edge (i, j) be P(i, j).

Now, from vertex i (other than vertex v), there are degree(i) possibilities.

Let the neighbors of i be i1, i2, ..., ik.

Thus, E[i] = [(1+E[i1])*P(i, i1) + ... + (1+E[ik])*P(i, ik)] / [P(i, i1) + ... + P(i, ik)]

For vertex v, E[v] = 0

So, we get a linear equation for each vertex. Overall, we get N equations with N unknowns, which can be solved by elimination.

However, I am getting Wrong Answer. I think there is a bug in my implementation of Gauss Jordan Elimination (this is the first time I am coding that). Can somebody help me?

Here is my code.

typedef long double db;
typedef vector<db> vd;
typedef vector<vd> vvd;
#define EPS 1e-9

db GaussJordan(vvd &A, vvd &b) {
	int N = A.size();
	assert(A[0].size() == N && b.size() == N);
	int M = b[0].size();
	db det = 1.0;
	for (int i = 0; i < N; i++) {
		if (fabs(A[i][i]) < EPS) {
			int j = i+1;
			for (; j < N; j++)
				if (fabs(A[j][i]) >= EPS)
					break;
			if (j == N)
			{
				assert(1 == 0);
				exit(-1);
			}
			for (int k = 0; k < N; k++)
				swap(A[i][k], A[j][k]);
			for (int k = 0; k < M; k++)
				swap(b[i][k], b[j][k]);
			det *= -1.0;
		}
		assert(fabs(A[i][i]) >= EPS);
		db pivot = A[i][i];
		det *= pivot;
		for (int k = 0; k < N; k++)
			A[i][k] /= pivot;
		for (int k = 0; k < M; k++)
			b[i][k] /= pivot;
		for (int j = i+1; j < N; j++) {
			db mul = A[j][i];
			for (int k = 0; k < N; k++)
				A[j][k] -= A[i][k]*mul;
			for (int k = 0; k < M; k++)
				b[j][k] -= b[i][k]*mul;
		}			
	}
	for (int i = N-1; i >= 0; i--) {
		for (int j = i-1; j >= 0; j--) {
			db mul = A[j][i];
			for (int k = 0; k < N; k++)
				A[j][k] -= A[i][k]*mul;
			for (int k = 0; k < M; k++)
				b[j][k] -= b[i][k]*mul;
		}
	}
	return det;
}

vector<pii> adj[20];

int main() {
	int T;
	scanf("%d", &T);
	for (int t = 0; t < T; t++) {
		int N, u, v;
		scanf("%d %d %d", &N, &u, &v);
		for (int i = 0; i < N-1; i++) {
			int x, y, p;
			scanf("%d %d %d", &x, &y, &p);
			adj[x-1].push_back(pii(y-1, p));
			adj[y-1].push_back(pii(x-1, p));
		}
		vvd A;
		vvd b;
		for (int i = 0; i < N; i++) {
			db tmp[N];
			for (int j = 0; j < N; j++)
				tmp[j] = 0.0;
			if (i == v-1)
			{
				tmp[i] = 1.0;
				A.push_back(vd(tmp, tmp+N));
				vd t;
				t.push_back(0.0);
				b.push_back(t);
				continue;
			}
			int sum = 0;
			for (int j = 0; j < adj[i].size(); j++)
				sum += adj[i][j].second;
			for (int j = 0; j < adj[i].size(); j++)
				tmp[adj[i][j].first] = -(db)adj[i][j].second/(db)sum;
			tmp[i] = 1.0;
			A.push_back(vd(tmp, tmp+N));
			vd t;
			t.push_back(1.0);
			b.push_back(t);
		}
		GaussJordan(A, b);
		printf("%.5lLf\n", b[u-1][0]);
	}
	return 0;
}