#include <bits/stdc++.h>
// #include <ext/pb_ds/assoc_container.hpp>
#define asc(x) x.begin(), x.end()
#define des(x) x.rbegin(), x.rend()
#define pb push_back
#define eb emplace_back
#define ll long long
#define db long double
#define fori(a, b) for (int i = (a); i < (b); ++i)
#define forj(a, b) for (int j = (a); j < (b); ++j)

using namespace std;
// using namespace __gnu_pbds;
// using ordered_set = tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update>;

mt19937 gen(228);

const int N = 1e6 + 5, MOD = 998244353;

int pw(int x, int y) {
    if (x == 1) return 1;
    if (y == 1) return x;

    int _x = pw(x, y / 2);

    _x = (1ll * _x * _x) % MOD;
    if (y % 2) {
        _x = (1ll * _x * x) % MOD;
    }

    return _x;
}

int sz[N], f[N], frv[N];
bool need[N];
vector<int> g[N], dp[N];

int c(int n, int k) {
    if (k > n || k < 0) return 0;

    return 1ll * f[n] * frv[k] % MOD * frv[n - k] % MOD;
}

void dfs(int v) {
    sz[v] = 1;

    for (int u : g[v]) {
        dfs(u);

        sz[v] += sz[u];
    }


    if (!need[v]) {
        int mx = 1;

        dp[v].assign(mx + 1, 0);
        dp[v][1] = 1;

        for (int u : g[v]) {
            vector<int> ndp(sz[u] + mx + 1, 0);

            fori(0, mx + 1) {
                if (!dp[v][i]) continue;

                forj(0, sz[u] + 1) {
                    ndp[i + j] += (1ll * dp[v][i] * dp[u][j]) % MOD;

                    if (ndp[i + j] >= MOD) {
                        ndp[i + j] -= MOD;
                    }
                }
            }

            swap(dp[v], ndp);

            mx += sz[u];
        }

        return;
    }

    vector<vector<int>> dpp(2, vector<int>(3, 0));
    dpp[1][0] = 1, dpp[0][2] = 1;

    int mx = 1;

    for(int u : g[v]) {
        vector<vector<int>> ndp(mx + sz[u] + 1, vector<int>(3, 0));

        for (int cnt = 0; cnt < 3; ++cnt) {
            fori(0, mx + 1) {
                if (!dpp[i][cnt]) continue;

                forj(0, sz[u] + 1) {
                    for (int add = 0; add <= j; ++add) {
                        if (add) {
                            ndp[i + j - add][min(cnt + 1, 2)] += (1ll * dpp[i][cnt] * dp[u][j] % MOD * c(j, add)) % MOD;

                            if (ndp[i + j - add][min(cnt + 1, 2)] >= MOD) {
                                ndp[i + j - add][min(cnt + 1, 2)] -= MOD;
                            }
                        } else {
                            ndp[i + j - add][cnt] += (1ll * dpp[i][cnt] * dp[u][j] % MOD * c(j, add)) % MOD;

                            if (ndp[i + j - add][cnt] >= MOD) {
                                ndp[i + j - add][cnt] -= MOD;
                            }
                        }
                    }
                }
            }
        }

        mx += sz[u];

        swap(ndp, dpp);
    }

    dp[v].assign(mx + 1, 0);
    fori(0, sz[v] + 1) {
        dp[v][i] = dpp[i][2];
    }
}

void solve() {
    int n, k; cin >> n >> k;

    fori(0, k) {
        int x; cin >> x; --x;

        need[x] = 1;
    }

    fori(1, n) {
        int p; cin >> p; --p;

        g[p].pb(i);
    }

    dfs(0);

    ll ans = 0;
    fori(0, n + 1) {
        ans += dp[0][i];
    }

    cout << ans % MOD;
}

signed main() {
    ios_base::sync_with_stdio(0), cin.tie(0), cout.tie(0);

    f[0] = frv[0] = 1;
    fori(1, N) {
        f[i] = 1ll * f[i - 1] * i % MOD;
    }

    frv[N - 1] = pw(f[N - 1], MOD - 2);

    for (int i = N - 2; i > 0; --i) {
        frv[i] = 1ll * frv[i + 1] * (i + 1) % MOD;
    }

    int tt = 1;
    while (tt--) {
        solve();
        cout << "\n";
    }
}