#ifdef LOCAL
#define _GLIBCXX_DEBUG
#endif

#include <vector>

// #pragma GCC optimize("O3, Ofast")
// #pragma GCC target("avx2")

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>

using namespace std;


template<class T>
istream &operator>>(istream &in, vector<T> &a) {
    for (auto &i: a) {
        cin >> i;
    }
    return in;
}

template<class T1, class T2>
istream &operator>>(istream &in, pair<T1, T2> &a) {
    cin >> a.first >> a.second;
    return in;
}

template<class T>
ostream &operator<<(ostream &out, vector<T> &a) {
    for (auto &i: a) {
        cout << i << ' ';
    }
    cout << '\n';
    return out;
}

void solve();

template<class T>
using oset = __gnu_pbds::tree<T, __gnu_pbds::null_type, less<T>, __gnu_pbds::rb_tree_tag,
    __gnu_pbds::tree_order_statistics_node_update>;

signed main() {
#ifdef LOCAL
    freopen(R"(C:\work\untitled\input.txt)", "r", stdin);
    freopen(R"(C:\work\untitled\output.txt)", "w", stdout);

    auto begin_ = chrono::steady_clock::now();
#else
    cin.tie(0);
    cout.tie(0);
    std::ios_base::sync_with_stdio(0);
#endif
    int t_ = 1;
    // cin >> t_;
    while (t_--) {
        solve();
    }
#ifdef LOCAL
    auto end_ = chrono::steady_clock::now();

    cout << "\n\n===========\n";
    cout << chrono::duration_cast<chrono::milliseconds>(end_ - begin_).count() << " ms";
#endif
}

using ll = long long;
#define int ll

using pii = pair<int, int>;
using vi = vector<int>;
using vp = vector<pii>;
using vv = vector<vi>;

const int maxn = 1e6 + 1;
const int mod = 998244353;

struct md {
    int x = 0;

    md(ll p = 0) {
        x = p;
    }

    md operator+(md o) {
        if (x + o.x >= mod) {
            return x + o.x - mod;
        }
        return x + o.x;
    }

    md operator *(md o) {
        return {1ll * x * o.x % mod};
    }

    md operator^(ll n) {
        if (n == 0) {
            return 1;
        }
        if (n % 2 == 0) {
            md g = *this ^ (n / 2);
            return g * g;
        }
        return *this * (*this ^ (n - 1));
    }
};

vi g[maxn];

vector<vector<md> > dp(maxn);
vi hot(maxn);
vector<md> fact(2 * maxn);
vi sm(maxn);
vi sz(maxn);

md c(int n, int k) {
    return fact[n] * (fact[k] ^ (mod - 2)) * (fact[n - k] ^ (mod - 2));
}

void dfs(int v) {
    sm[v] = hot[v];
    sz[v] = 1;
    for (int u: g[v]) {
        dfs(u);
        sm[v] += sm[u];
        sz[v] += sz[u];
    }
    dp[v].resize(sz[v] + 1);
    if (!hot[v]) {
        vector<vector<md> > d(g[v].size() + 1, vector<md>(sz[v] + 1));
        d[0][0] = 1;
        int i = 0;
        for (int u: g[v]) {
            i++;
            for (int now = sz[v]; now >= 0; --now) {
                for (int k = 0; k <= min(now, sz[u]); ++k) {
                    d[i][now] = d[i][now] + d[i - 1][now - k] * dp[u][k];
                }
            }
        }
        dp[v] = d[g[v].size()];
    } else {
        md d[g[v].size() + 1][3][sz[v] + 1];
        d[0][2][1] = 1;
        d[0][0][0] = 1;
        int i = 0;
        for (int u: g[v]) {
            i++;
            d[i][0][0] = 1;
            for (int now = 1; now <= sz[v]; ++now) {
                for (int k = 0; k <= min(sz[u], now); ++k) {
                    for (int use = 1; use <= sz[u] - k && k + use <= now; ++use) {
                        d[i][2][now] = d[i][2][now] + d[i - 1][2][now - k - use] * dp[u][k] * c(sz[u] - k, use);
                        d[i][2][now] = d[i][2][now] + d[i - 1][1][now - k - use] * dp[u][k] * c(sz[u] - k, use);
                        d[i][1][now] = d[i][1][now] + d[i - 1][0][now - k - use] * dp[u][k] * c(sz[u] - k, use);
                    }
                    d[i][0][now] = d[i][0][now] + d[i - 1][0][now - k] * dp[u][k] * md(1);
                    d[i][1][now] = d[i][1][now] + d[i - 1][1][now - k] * dp[u][k] * md(1);
                    d[i][2][now] = d[i][2][now] + d[i - 1][2][now - k] * dp[u][k] * md(1);
                }
            }
        }
        for (int i = 0; i <= sz[v]; ++i) {
            dp[v][i] = d[g[v].size()][2][i];
        }
    }
}


void solve() {
    fact[0] = 1;
    for (int i = 1; i < 2 * maxn; ++i) {
        fact[i] = fact[i - 1] * md(i);
    }
    int n, k;
    cin >> n >> k;
    for (int i = 0; i < k; ++i) {
        int x;
        cin >> x;
        hot[x] = 1;
    }
    // cout << ((md(3) ^ (mod - 2)) * md(3)).x << '\n';
    for (int i = 2; i <= n; ++i) {
        int p;
        cin >> p;
        g[p].push_back(i);
    }
    dfs(1);
    md res;
    for (int i = 0; i < dp[1].size(); ++i) {
        res = res + dp[1][i];
    }
    cout << res.x << '\n';
}
