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;
}