#include <bits/stdc++.h>
using namespace std;
#define FOR(i,k,n) for(int i = (int)(k); i < (int)(n); i++)
#define REP(i,n) FOR(i,0,n)
#define ALL(a) a.begin(), a.end()
#define MS(m,v) memset(m,v,sizeof(m))
typedef long long ll;
typedef long double ld;
typedef vector<int> vi;
typedef vector<string> vs;
typedef pair<int, int> pii;
const int MOD = 1000000007;
const int INF = MOD + 1;
const ld EPS = 1e-12;
template<class T> T &chmin(T &a, const T &b) { return a = min(a, b); }
template<class T> T &chmax(T &a, const T &b) { return a = max(a, b); }

/*--------------------template--------------------*/


typedef ld Data;
typedef vector<Data> Array;
typedef vector<Array> Matrix;

bool is_zero(Data dat) { return (abs(dat) < EPS); }

Matrix operator-(Matrix mat) {
	REP(i, mat.size()) REP(j, mat[0].size()) mat[i][j] = -mat[i][j];
	return mat;
}

Matrix operator+(Matrix lhs, const Matrix &rhs) {
	REP(i, lhs.size()) REP(j, lhs[0].size()) lhs[i][j] = lhs[i][j] + rhs[i][j];
	return lhs;
}

Matrix operator-(Matrix lhs, const Matrix &rhs) {
	REP(i, lhs.size()) REP(j, lhs[0].size()) lhs[i][j] = lhs[i][j] - rhs[i][j];
	return lhs;
}

Matrix operator*(const Matrix &lhs, const Matrix &rhs) {
	Matrix res(lhs.size(), Array(rhs[0].size(), 0));
	REP(i, lhs.size()) REP(j, rhs[0].size()) REP(k, rhs.size())
		res[i][j] = fma(lhs[i][k], rhs[k][j], res[i][j]);
	return res;
}

Array operator+(const Array &lhs, const Array &rhs) {
	Array res = lhs;
	REP(i, rhs.size()) res[i] = res[i] + rhs[i];
	return res;
}

Array operator*(const Matrix &lhs, const Array &rhs) {
	Array res(lhs.size());
	REP(i, lhs.size())
		res[i] = inner_product(begin(lhs[i]), end(lhs[i]), begin(rhs), Data(0));
	return res;
}

Array operator*(Data scalar, const Array &rhs) {
	Array res(rhs.size());
	REP(i, rhs.size())
		res[i] = scalar * rhs[i];
	return res;
}

Matrix transpose(const Matrix &A) {
	const int n = A[0].size(), m = A.size();
	Matrix res(n, Array(m, 0));
	REP(i, n)REP(j, m) res[i][j] = A[j][i];
	return res;
}

Data inner_product_ex(const Array& a, const Array& b, const Matrix& m) {
	Array b_prime = m*b;
	return inner_product(begin(a), end(a), begin(b_prime), Data(0));
}

Data norm(const Array &a) {
	return inner_product(begin(a), end(a), begin(a), Data(0));
}

Data norm(const Array &a, const Matrix &M) {
	return inner_product_ex(a, a, M);
}

int rankMat(Matrix A) {
	const int n = A.size(), m = A[0].size();
	int r = 0;
	for (int i = 0; r < n && i < m; ++i) {
		int pivot = r;
		for (int j = r + 1; j < n; ++j)
			if (abs(A[j][i]) > abs(A[pivot][i])) pivot = j;
		swap(A[pivot], A[r]);
		if (is_zero(A[r][i])) continue;
		for (int k = m - 1; k >= i; --k)
			A[r][k] = A[r][k] / A[r][i];
		for (int j = r + 1; j < n; ++j)
			for (int k = m - 1; k >= i; --k)
				A[j][k] = fma(-A[r][k], A[j][i], A[j][k]);
		++r;
	}
	return r;
}

Data det(Matrix A) {
	const int n = A.size();
	Data D = Data(1);
	for (int i = 0; i < n; ++i) {
		int pivot = i;
		for (int j = i + 1; j < n; ++j)
			if (abs(A[j][i]) > abs(A[pivot][i])) pivot = j;
		swap(A[pivot], A[i]);
		D = D * A[i][i] * Data(i != pivot ? -1 : 1);
		if (is_zero(A[i][i])) break;
		for (int j = i + 1; j < n; ++j)
			for (int k = n - 1; k >= i; --k)
				A[j][k] = A[j][k] - A[i][k] * A[j][i] / A[i][i];
	}
	return D;
}

Array gauss_jordan(const Matrix& A, const Array& b)
{
	int n = A.size();
	Matrix B(n, Array(n + 1));
	REP(i, n)REP(j, n) B[i][j] = A[i][j];
	REP(i, n) B[i][n] = b[i];
	REP(i, n)
	{
		int pivot = i;
		FOR(j, i, n)
		{
			if (abs(B[j][i]) > abs(B[pivot][i])) pivot = j;
		}
		swap(B[i], B[pivot]);

		if (abs(B[i][i]) < EPS) return Array();
		FOR(j, i + 1, n + 1) B[i][j] /= B[i][i];
		REP(j, n)
		{
			if (i != j)
			{
				FOR(k, i + 1, n + 1) B[j][k] -= B[j][i] * B[i][k];
			}
		}
	}
	Array x(n);
	REP(i, n) x[i] = B[i][n];
	return x;
}

typedef int Weight;
struct Edge
{
	int from, to; Weight cost;
	bool operator < (const Edge& e) const { return cost < e.cost; }
	bool operator > (const Edge& e) const { return cost > e.cost; }
};
typedef vector<Edge> Edges;
typedef vector<Edges> Graph;
void add_edge(Graph &g, int from, int to, Weight cost)
{
	g[from].push_back(Edge{ from, to, cost });
}

void dijkstra(Graph &g, int s, vector<Weight> &d)
{
	d.assign(g.size(), INF);
	d[s] = 0;
	typedef pair<Weight, int> P;
	priority_queue<P, vector<P>, greater<P>> que;
	que.push(P(0, s));
	while (!que.empty())
	{
		Weight dist = que.top().first;
		int v = que.top().second;
		que.pop();
		if (d[v] < dist) continue;
		REP(i, g[v].size())
		{
			Edge e = g[v][i];
			if (d[e.to] > d[v] + e.cost)
			{
				d[e.to] = d[v] + e.cost;
				que.push(P(d[e.to], e.to));
			}
		}
	}
}

int main()
{
	cin.sync_with_stdio(false); cout << fixed << setprecision(10);
	int n, s, t;
	while (cin >> n >> s >> t, n)
	{
		s--; t--;
		vi sign(n);
		REP(i, n) cin >> sign[i];
		Graph g(n);
		REP(i, n)REP(j, n)
		{
			int c; cin >> c;
			if (!c) continue;
			add_edge(g, i, j, c);
		}
		vector<Weight> dist;
		dijkstra(g, t, dist);
		if (dist[s] == INF)
		{
			cout << "impossible" << endl;
			continue;
		}
		Matrix A(n, Array(n));
		Array b(n);
		REP(i, n)
		{
			if (i == t)
			{
				REP(j, n)
				{
					A[i][j] = (i == j ? 1 : 0);
				}
				b[i] = 0;
			}
			else if (sign[i])
			{
				vi v;
				int d = INF;
				int sum = 0;
				for (auto e : g[i])
				{
					int dst = e.to;
					if (dist[dst] + e.cost > d) continue;
					if (dist[dst] + e.cost == d)
					{
						sum += e.cost;
						v.push_back(dst);
					}
					else if (dist[dst] + e.cost < d)
					{
						d = dist[dst] + e.cost;
						v.clear();
						sum = e.cost;
						v.push_back(dst);
					}
				}
				if (v.size() == 0) continue;
				for (auto j : v)
				{
					A[i][j] = -1;
				}
				A[i][i] += v.size();
				b[i] = sum;
			}
			else
			{
				vi v;
				int sum = 0;
				for (auto e : g[i])
				{
					sum += e.cost;
					v.push_back(e.to);
				}
				if (v.size() == 0) continue;
				for (auto j : v)
				{
					A[i][j] = -1;
				}
				A[i][i] += v.size();
				b[i] = sum;
			}
		}
		Array x = gauss_jordan(A, b);
		cout << x[s] << endl;
	}
	return 0;
}