Problem: ButCoder Inc. is a startup company whose main business is the development and operation of the programming competition site "ButCoder".

There are N members of the company including the president, and each member except the president has exactly one direct boss. Each member has a unique ID number from 1 to N, and the member with the ID i is called Member i. The president is Member 1, and the direct boss of Member i (2 ≤ i ≤ N) is Member b_i (1 ≤ b_i < i).

All the members in ButCoder have gathered in the great hall in the main office to take a group photo. The hall is very large, and all N people can stand in one line. However, they have been unable to decide the order in which they stand. For some reason, all of them except the president want to avoid standing next to their direct bosses.

How many ways are there for them to stand in a line that satisfy their desires? Find the count modulo 10^9+7.

Constraints

* 2 ≤ N ≤ 2000
* 1 ≤ b_i < i (2 ≤ i ≤ N)

Input

Input is given from Standard Input in the following format:


N
b_2
b_3
:
b_N


Output

Print the number of ways for the N members to stand in a line, such that no member except the president is next to his/her direct boss, modulo 10^9+7.

Examples

Input

4
1
2
3


Output

2


Input

3
1
2


Output

0


Input

5
1
1
3
3


Output

8


Input

15
1
2
3
1
4
2
7
1
8
2
8
1
8
2


Output

97193524

Solution (in CPP):

Solution:
#include <bits/stdc++.h>
#define ADD(a, b) a = (a + ll(b)) % mod
#define MUL(a, b) a = (a * ll(b)) % mod
#define MAX(a, b) a = max(a, b)
#define MIN(a, b) a = min(a, b)
#define rep(i, a, b) for(int i = int(a); i < int(b); i++)
#define rer(i, a, b) for(int i = int(a) - 1; i >= int(b); i--)
#define all(a) (a).begin(), (a).end()
#define sz(v) (int)(v).size()
#define pb push_back
#define sec second
#define fst first
#define debug(fmt, ...) Debug(__LINE__, ":", fmt, ##__VA_ARGS__)
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> pi;
typedef pair<ll, ll> pl;
typedef pair<int, pi> ppi;
typedef vector<int> vi;
typedef vector<ll> vl;
typedef vector<vl> mat;
typedef complex<double> comp;
void Debug() {cout << '\n'; }
template<class FIRST, class... REST>void Debug(FIRST arg, REST... rest){
	cout<<arg<<" ";Debug(rest...);}
template<class T>ostream& operator<<(ostream& out,const vector<T>& v) {
	out<<"[";if(!v.empty()){rep(i,0,sz(v)-1)out<<v[i]<<", ";out<<v.back();}out<<"]";return out;}
template<class S, class T>ostream& operator<<(ostream& out,const pair<S, T>& v){
	out<<"("<<v.first<<", "<<v.second<<")";return out;}
const int MAX_N = 200010;
const int MAX_V = 100010;
const double eps = 1e-6;
const ll mod = 1000000007;
const int inf = 1 << 29;
const ll linf = 1LL << 60;
const double PI = 3.14159265358979323846;
///////////////////////////////////////////////////////////////////////////////////////////////////

ll fac[MAX_N], inv[MAX_N], fiv[MAX_N]; //fiv:inv(fac(i))

ll mod_pow(ll a, ll n) {
	if(n == 0) return 1;
	ll res = mod_pow(a, n / 2);
	if(n % 2 == 0) return res * res % mod;
	else return a * res % mod * res % mod;
}

ll invf(ll a) {
	return mod_pow(a, mod - 2);
}

void C_init(int n) {
	fac[0] = fac[1] = 1; inv[1] = 1;
	fiv[0] = fiv[1] = 1;
	rep(i, 2, n + 1) {
		fac[i] = fac[i - 1] * i % mod;
		inv[i] = mod - inv[mod % i] * (mod / i) % mod;
		fiv[i] = fiv[i - 1] * inv[i] % mod;
	}
}

ll C(int a, int b) { //assume a >= b
	if(a < b || a < 0 || b < 0) return 0;
	return fac[a] * fiv[b] % mod * fiv[a - b] % mod;
}

int N;
vi G[MAX_N];
ll dp[2010][2010][3];
int ts[2010];

void loop(int v) {

	int n = sz(G[v]);
	ts[v] = 1;
	rep(i, 0, n) {
		int to = G[v][i]; 
		loop(to);
		ts[v] += ts[to];
	}

	int m = ts[v];

	vector<vl> tdp[2] = {vector<vl>(m + 1, vl(3, 0)), vector<vl>(m + 1, vl(3, 0))};
	tdp[0][1][0] = 1;
	int now = 0, nex = 1;
	int cnt = 1;
	rep(i, 0, n) {
		int to = G[v][i];

		rep(j, 0, cnt + ts[to] + 1) {
			rep(k, 0, 3) {
				tdp[nex][j][k] = 0;
			}
		}

		rep(j, 0, cnt + 1) {
			rep(k, 0, ts[to] + 1) {
				rep(l, 0, 3) ADD(tdp[nex][j + k][0], tdp[now][j][0] * dp[to][k][l] % mod);
				rep(l, 0, 3) ADD(tdp[nex][j + k][1], tdp[now][j][1] * dp[to][k][l] % mod);
				rep(l, 0, 3) ADD(tdp[nex][j + k][2], tdp[now][j][2] * dp[to][k][l] % mod);

				if(k < ts[to]) {
					ADD(tdp[nex][j + k][1], tdp[now][j][0] * dp[to][k + 1][0] % mod * 2 % mod);
					ADD(tdp[nex][j + k][1], tdp[now][j][0] * dp[to][k + 1][1] % mod);
					ADD(tdp[nex][j + k][2], tdp[now][j][1] * dp[to][k + 1][0] % mod);
					ADD(tdp[nex][j + k][2], tdp[now][j][1] * dp[to][k + 1][1] % mod * fiv[2] % mod);
				}
			}
		}
		cnt += ts[to];
		swap(now, nex);
	}
	rep(j, 0, m + 1) 
		rep(k, 0, 3) {
			dp[v][j][k] = tdp[now][j][k];
		}
}

void solve() {
	cin >> N;
	C_init(N);
	rep(i, 0, N - 1) {
		int a; cin >> a; a--;
		G[a].pb(i + 1);
	}
	loop(0);

	ll res = 0;
	rep(j, 0, N + 1) {
		rep(k, 0, 3) {
			if((N - j) % 2 == 0) ADD(res, dp[0][j][k] * fac[j] % mod);
			else ADD(res, mod - dp[0][j][k] * fac[j] % mod);
		}
	}
	cout << res << "\n";
}

int main() {
#ifndef LOCAL
	ios::sync_with_stdio(false);
    cin.tie(0);
#endif
    cout << fixed;
	cout.precision(20);
	srand((unsigned int)time(NULL));
#ifdef LOCAL
	//freopen("in.txt", "wt", stdout); //for tester
    freopen("in.txt", "rt", stdin);
#endif	
	solve();
#ifdef LOCAL
    cerr << "Time elapsed: " << 1.0 * clock() / CLOCKS_PER_SEC << " s.\n";
#endif
	return 0;
}