#ifndef LOCAL
#pragma GCC Optimize("O3, Ofast")
#pragma GCC Target("avx,avx2,unroll-loops")
#endif
#include <bits/stdc++.h>

using namespace std;
using ll = long long;
using ld = long double;

#define f first
#define s second
#define all(x) (x).begin(), (x).end()
#define mp make_pair
#define pb push_back
#define pii pair<int, int>
#ifndef LOCAL
#define endl "\n"
#endif

struct F {
    vector<int> t;
    F(int n) {
        t.assign(n, 0);
    }
    int get(int r) {
        int ans = 0;
        for (; r >= 0; r = (r&(r + 1)) - 1) {
            ans += t[r];
        }
        return ans;
    }
    void upd(int i, int x) {
        for (; i < t.size(); i = (i|(i + 1))) {
            t[i] += x;
        }
    }
};

mt19937 rnd(time(0));

const int MOD = 998244353;

int power(int x, int n) {
    int ans = 1;
    while (n) {
        if (n&1) {
            ans = (ans * 1ll * x) % MOD;
        }
        n /= 2;
        x = (x * 1ll * x) % MOD;
    }
    return ans;
}

int add(int a, int b) {
    a += b;
    if (a >= MOD) {
        a -= MOD;
    }
    return a;
}

int mul(int a, int b) {
    return (a * 1ll * b) % MOD;
}
int sub(int a, int b) {
    a -= b;
    if (a < 0) {
        a += MOD;
    }
    return a;
}

void solve() {
    int n, k;
    cin >> n >> k;
    // int pw[k + 1][k + 1];
    // for (int x = 0; x <= k; ++x) {
    //     pw[x][0] = 1;
    //     for (int y = 1; y <= k; ++y) {
    //         pw[x][y] = mul(pw[x][y - 1], x);
    //     }
    // }
    // auto power = [&](int x, int n) {
    //     if (n <= k) {
    //         return pw[x][n];
    //     }
    //     int cnt = n % k;

    // };
    vector<int> vertex(k);
    for (auto &x : vertex) {
        cin >> x;
        // x = (rnd() % n) + 1;
        --x;
    }
    vector<int> act = vertex;
    vector<int> a(n), isA(n);
    for (auto &x : act) {
        a[x] = 1;
        isA[x] = 1;
    }
    vector<vector<int>> graph(n);
    for (int i = 1; i < n; ++i) {
        int p;
        cin >> p;
        // p = max(rnd() % i, rnd() % i) + 1;
        --p;
        graph[p].pb(i);
    }
    vector<int> tin(n), tout(n), sz(n, 1);
    int Tm = 0;
    auto dfs = [&](int v, int prev, auto &&self) -> void {
        tin[v] = Tm++;
        for (auto &u : graph[v]) {
            self(u, v, self);
            sz[v] += sz[u];
            a[v] += a[u];
        }
        tout[v] = Tm;
    };
    dfs(0, 0, dfs);
    vector<vector<pii>> dp(4010, vector<pii>(k + 2));
    int H = 1;
    for (int i = 0; i < dp[0].size(); ++i) {
        dp[0][i].f = 1;
        dp[0][i].s = 0;
    }
    vector<vector<int>> M(2, vector<int>(2)), M2 = M;
    auto Sdfs = [&](int v, int way, auto &&self) -> pii { // ind + k
        if (a[v] == 0) {
            return mp(0, sz[v]);
        }
        vector<pii> kids;
        int cntActive = 0;
        int hisSz = 0;
        int activeInd = -1;
        int tot = 0;
        if (!isA[v]) {
            for (auto &u : graph[v]) {
                if (a[u]) {
                    kids.pb(self(u, way + isA[v], self));
                    activeInd = kids.size() - 1;
                    ++cntActive;
                    hisSz += sz[u];
                } else {
                    tot += sz[u];
                }
            }
            kids.pb(mp(0, tot));
        } else {
            for (auto &u : graph[v]) {
                kids.pb(self(u, way + isA[v], self));
                if (a[u]) {
                    activeInd = kids.size() - 1;
                    ++cntActive;
                    hisSz += sz[u];
                }
            }
        }
        if (cntActive == 1 && !isA[v]) {
            return mp(kids[activeInd].f, kids[activeInd].s + sz[v] - hisSz);
        }
        assert(H < dp.size());
        int ret = H++;
        for (int S = 0; S <= way; ++S) {
            M[0][0] = 1;
            M[0][1] = 0;
            M[1][0] = 0;
            M[1][1] = 0;
            for (auto &u : kids) {
                M2[0][0] = 0;
                M2[0][1] = 0;
                M2[1][0] = 0;
                M2[1][1] = 0;
                for (int z = 0; z <= 1; ++z) {
                    M2[0][z] = add(M2[0][z], mul(M[0][z], dp[u.f][S].f * 1ll * power(S + 1, u.s) % MOD));
                    M2[1][z] = add(M2[1][z], mul(M[1][z], dp[u.f][S].f * 1ll * power(S + 1, u.s) % MOD));
                    M2[1][z^1] = add(M2[1][z^1], mul(M[0][z], sub(dp[u.f][S + 1].f * 1ll * power(S + 2, u.s) % MOD, dp[u.f][S].f * 1ll * power(S + 1, u.s) % MOD)));
                    M2[0][z^1] = add(M2[0][z^1], mul(M[0][z], dp[u.f][S].s * 1ll * power(S + 1, u.s) % MOD));
                    M2[1][z^1] = add(M2[1][z^1], mul(M[1][z], dp[u.f][S].s * 1ll * power(S + 1, u.s) % MOD));
                    M2[1][z] = add(M2[1][z], mul(M[0][z], sub(dp[u.f][S + 1].s * 1ll * power(S + 2, u.s) % MOD, dp[u.f][S].s * 1ll * power(S + 1, u.s) % MOD)));
                }
                M = M2;
            }
            if (isA[v]) {
                if (S) {
                    dp[ret][S - 1].f = add(dp[ret][S - 1].f, M[0][0] * 1ll * (S + 1) % MOD);
                    dp[ret][S - 1].s = add(dp[ret][S - 1].s, M[0][1] * 1ll * (S + 1) % MOD);
                }
                dp[ret][S].s = add(dp[ret][S].s, (M[0][0] + M[1][1]) * 1ll * (S + 1) % MOD);
                dp[ret][S].f = add(dp[ret][S].f, (M[0][1] + M[1][0]) * 1ll * (S + 1) % MOD);
            } else {
                dp[ret][S].f = add(dp[ret][S].f, M[0][0] * 1ll * (S + 1) % MOD);
                dp[ret][S].s = add(dp[ret][S].s, M[0][1] * 1ll * (S + 1) % MOD);
            }
        }
        return mp(ret, 0);
    };
    int ret = Sdfs(0, 1, Sdfs).f;
    cout << (dp[ret][0].f - dp[ret][0].s + MOD) % MOD << endl;
}

signed main() {
    #ifdef LOCAL
        freopen("in.txt", "r", stdin);
        freopen("out.txt", "w", stdout);
    #else
        ios::sync_with_stdio(false);
        cin.tie(0);
    #endif
    int t = 1;
    // cin >> t;
    while (t--) {
        solve();
    }
}