#include "circuit.h"
#include <bits/stdc++.h>
using namespace std;
#define REP(i, j, k) for (int i = j; i < k; i++)
#define RREP(i, j, k) for (int i = j; i >= k; i--)
template <class T>
inline bool mnto(T &a, const T b) {return a > b ? a = b, 1 : 0;}
template <class T>
inline bool mxto(T &a, const T b) {return a < b ? a = b, 1 : 0;}
typedef long long ll;
#define FI first
#define SE second
typedef pair<int, int> ii;
typedef pair<ll, ll> pll;
#define ALL(x) x.begin(), x.end()
#define SZ(x) (int) x.size()
#define pb push_back
typedef vector<int> vi;
typedef vector<ll> vll;
typedef vector<ii> vii;
typedef tuple<int, int, int> iii;
#ifndef DEBUG
#define cerr if (0) cerr
#endif
#define submod(x) if (x >= MOD) x -= MOD
#define addmod(x) if (x < 0) x += MOD
const int MAXN = 200005;
const int INF = 1000000007;
const ll LINF = 1000000000000000007;
const int MOD = 1000002022;
namespace brute {
int n, m;
vi p, a;
vi adj[MAXN];
ll ways[MAXN];
ll dfs(int u) {
if (u >= n) {
ways[u] = 1;
return a[u - n];
}
vll dp(1, 1);
ways[u] = SZ(adj[u]);
for (int v : adj[u]) {
ll tmp = dfs(v);
vll ndp(SZ(dp) + 1, 0);
REP (i, 0, SZ(dp)) {
ndp[i + 1] += dp[i] * tmp % MOD;
submod(ndp[i + 1]);
}
REP (i, 0, SZ(dp)) {
ndp[i] += dp[i] * (ways[v] - tmp + MOD) % MOD;
submod(ndp[i]);
}
dp = ndp;
ways[u] = ways[u] * ways[v] % MOD;
}
ll res = 0;
REP (i, 1, SZ(dp)) {
res += dp[i] * i % MOD;
submod(res);
}
return res;
}
void init(int N, int M, vi P, vi A) {
n = N, m = M, p = P, a = A;
REP (i, 1, n + m) {
adj[p[i]].pb(i);
}
}
int count_ways(int L, int R) {
int l = L, r = R;
REP (i, l - n, r + 1 - n) {
a[i] ^= 1;
}
return dfs(0);
}
}
int n, m;
vi p, a;
vi adj[MAXN];
ll ways[MAXN];
void dfs(int u) {
ways[u] = max(1, SZ(adj[u]));
for (int v : adj[u]) {
dfs(v);
ways[u] = ways[u] * ways[v] % MOD;
}
}
ll pw[MAXN];
void dfs2(int u) {
if (adj[u].empty()) {
return;
}
vll pp(SZ(adj[u])), sp(SZ(adj[u]));
REP (i, 0, SZ(adj[u])) {
pp[i] = ways[adj[u][i]];
if (i) {
pp[i] = pp[i] * pp[i - 1] % MOD;
}
}
RREP (i, SZ(adj[u]) - 1, 0) {
sp[i] = ways[adj[u][i]];
if (i != SZ(adj[u]) - 1) {
sp[i] = sp[i] * sp[i + 1] % MOD;
}
}
REP (i, 0, SZ(adj[u])) {
ll tmp = (i == 0 ? 1 : pp[i - 1]) * (i + 1 >= SZ(adj[u]) ? 1 : sp[i + 1]) % MOD;
pw[adj[u][i]] = tmp * pw[u] % MOD;
dfs2(adj[u][i]);
}
}
#define MLR int mid = lo + hi >> 1, lc = u << 1, rc = u << 1 ^ 1
int sm[MAXN * 4], lz[MAXN * 4];
ll coef[MAXN * 4];
void propo(int u, int lo, int hi) {
if (lo == hi || lz[u] == 0) return;
MLR;
sm[lc] = coef[lc] - sm[lc];
addmod(sm[lc]);
sm[rc] = coef[rc] - sm[rc];
addmod(sm[rc]);
lz[lc] ^= 1;
lz[rc] ^= 1;
lz[u] = 0;
}
void init(int u = 1, int lo = 0, int hi = m - 1) {
if (lo == hi) {
coef[u] = pw[lo + n];
sm[u] = a[lo] * coef[u] % MOD;
lz[u] = 0;
return;
}
MLR;
init(lc, lo, mid);
init(rc, mid + 1, hi);
sm[u] = sm[lc] + sm[rc];
submod(sm[u]);
coef[u] = coef[lc] + coef[rc];
submod(coef[u]);
}
void upd(int s, int e, int u = 1, int lo = 0, int hi = m - 1) {
if (lo >= s && hi <= e) {
lz[u] ^= 1;
sm[u] = coef[u] - sm[u];
addmod(sm[u]);
return;
}
propo(u, lo, hi);
MLR;
if (s <= mid) {
upd(s, e, lc, lo, mid);
}
if (e > mid) {
upd(s, e, rc, mid + 1, hi);
}
sm[u] = sm[lc] + sm[rc];
submod(sm[u]);
}
void init(int N, int M, vi P, vi A) {
n = N, m = M, p = P, a = A;
if (N <= 1000 && M <= 1000) {
brute::init(N, M, P, A);
return;
}
REP (i, 1, n + m) {
adj[p[i]].pb(i);
}
dfs(0);
pw[0] = 1;
dfs2(0);
REP (i, n, n + m) {
cerr << i << ": " << pw[i] << '\n';
}
init();
brute::init(N, M, P, A);
}
int count_ways(int L, int R) {
if (n <= 1000 && m <= 1000) {
return brute::count_ways(L, R);
}
int l = L, r = R;
upd(l - n, r - n);
//assert(sm[1] == brute::count_ways(L, R));
return sm[1];
}
