PROBLEM LINK:
Author: me
Tester: radoslav192
Editorialist: also me
DIFFICULTY:
Super-hard
PREREQUISITES:
- sweepline Mo
- Manipulations of trees
PROBLEM:
Given an n-node DAG composed of an outwards arborescence with m additional edges, solve queries of the form: given [l,r], count the number of i,j such that l\le i<j\le r and that node i can reach node j.
QUICK EXPLANATION:
After the Sweepline Mo reduction, we can relabel the nodes in such a way so that the nodes any node can reach can be represented by O(m) intervals, and the nodes that can reach any node can also be represented by O(m) intervals. We then use a generalized segment tree to solve O(nm) interval add, O(n\sqrt n) point value query operations in O(\frac{n\log n\sqrt n}{W(\sqrt n/m)}).
EXPLANATION:
The first step for these interval pair counting problems is sweepline Mo’s algorithm. If you don’t know about sweepline Mo you can refer to one of my blog posts about it.
Sweepline Mo effectively reduces this problem into allowing the rapid computation of two things:
- The number of nodes in [1,i] that can reach x.
- The number of nodes in [1,i] that x can reach.
(presumed for some monotonically increasing i.) We will deal with these separately; we call the former Prefix-to-Node, and the latter Node-to-Prefix.
Prefix-to-Node
Consider the DFS order of the given tree. A prime attribute of the DFS order is that for a tree, the nodes that an arbitrary node can reach can be represented by exactly one interval of the DFS order. This interval is the representative interval of a node. We can extend this to our DAG, too. We first calculate the DFS order of the underlying tree.
Theorem 1: The nodes a node can reach on the DAG can be represented by O(m) disjoint intervals of the DFS order.
We can prove this constructively. A node can reach an edge if and only if it can reach the starting point of this edge. Define the critical nodes of a node u as the set \{u\}\cup\{y:u\text{ can reach extra edge }(x,y)\}.
Lemma 1: I claim that the union of the representative intervals of the critical nodes of u are exactly the nodes that u can reach.
Consider the path from u to a node it can reach, v. The last critical node on this path represents which representative interval v is contained in. In reality, we don’t need to use critical nodes to implement Prefix-to-Node; we can calculate the intervals in reverse topological order and merge all the intervals of nodes adjacent to a node with sweepline.
Node-to-Prefix
We now slightly change what intervals are. They no longer are intervals of the DFS order, but instead a path from a node to one of its descendants. Like the DFS order interval, the ancestors of a node on a tree is exactly one of these intervals.
Theorem 2: The nodes that can reach a node on the DAG can be represented by O(m) disjoint node-descendant paths.
We prove this through contribution: each extra link will contribute at most one such path to all the nodes its destination reaches, and there are O(m) such extra links. The actual construction is a bit messier. We calculate the paths for a node in topological order, ensuring that when we calculate node u we already know the paths of its father fa_u. First, take the path in fa_u that ends with fa_u and extend it to u. Then, add the paths the extra edges’ contribution, like this.
Define the critical ancestors of node u as the set \{x:\text{extra edge }(x,y)\text{ can reach }u\}; we can calculate this by iterating over the DAG, once again, in topological order. To add the extra contribution, we iterate over all x such that (x,u) is an extra edge. Then we take the critical ancestors of all such x. Finally, for all critical ancestors, we use brute force and keep taking the father of that node until we reach a node already in the paths, and we add a new path from the last node not in the paths to said critical ancestor.
The time complexity of this is correct: a node can be a critical ancestor of the start of an extra edge at most m^2 times, so the time complexity of this initialization is O(nm^2) with a low constant. With some merging tactics we can do O(nm), but it is not required for this problem.
The initial problem is now effectively reduced to doing the following things:
- Querying the weight of a node.
- Adding 1 to all nodes in a subtree.
- Adding 1 to all nodes in a node-descendant path.
The first operation is done O(n\sqrt n) times, and the second/third operation is done O(nm) times. We can reduce both to operations on a linear array, and what we need to do now is O(nm) range add operations and O(n\sqrt n) point query operations. Directly employing a segment tree gives a O(n(m+\sqrt n)\log n) solution; directly employing square-root decomposition gives a O(nm\sqrt n) solution. Consider balancing between a segment tree and square-root decomposition by doing B layers of decomposition. This way, modification is O(Bn^{1/B}) and queries are O(B). Analyzing the best B gives a time complexity of O(\frac{n\log n\sqrt n}{W(\sqrt n/m)}). Keeping in mind that W(n) is approximately \ln n-\ln\ln n, this value is almost O(n\log m\sqrt n), with quite a low constant factor because we can’t simultaneously maximize the Prefix-to-Node intervals and the Node-to-Prefix intervals.
In practice, B=3 works well.
SOLUTIONS:
Setter's Solution
// writer: w33z8kqrqk8zzzx33
#include <bits/stdc++.h>
using namespace std;
// begin fast read template by CYJian (source: https://www.luogu.com.cn/paste/i11c3ppx)
namespace io {
const int __SIZE = (1 << 21) + 1;
char ibuf[__SIZE], *iS, *iT, obuf[__SIZE], *oS = obuf, *oT = oS + __SIZE - 1, __c, qu[55]; int __f, qr, _eof;
#define Gc() (iS == iT ? (iT = (iS = ibuf) + fread (ibuf, 1, __SIZE, stdin), (iS == iT ? EOF : *iS ++)) : *iS ++)
inline void flush () { fwrite (obuf, 1, oS - obuf, stdout), oS = obuf; }
inline void gc (char &x) { x = Gc(); }
inline void pc (char x) { *oS ++ = x; if (oS == oT) flush (); }
inline void pstr (const char *s) { int __len = strlen(s); for (__f = 0; __f < __len; ++__f) pc (s[__f]); }
inline void gstr (char *s) { for(__c = Gc(); __c < 32 || __c > 126 || __c == ' ';) __c = Gc();
for(; __c > 31 && __c < 127 && __c != ' '; ++s, __c = Gc()) *s = __c; *s = 0; }
template <class I> inline bool gi (I &x) { _eof = 0;
for (__f = 1, __c = Gc(); (__c < '0' || __c > '9') && !_eof; __c = Gc()) { if (__c == '-') __f = -1; _eof |= __c == EOF; }
for (x = 0; __c <= '9' && __c >= '0' && !_eof; __c = Gc()) x = x * 10 + (__c & 15), _eof |= __c == EOF; x *= __f; return !_eof; }
template <class I> inline void print (I x) { if (!x) pc ('0'); if (x < 0) pc ('-'), x = -x;
while (x) qu[++ qr] = x % 10 + '0', x /= 10; while (qr) pc (qu[qr --]); }
struct Flusher_ {~Flusher_(){flush();}}io_flusher_;
} using io::pc; using io::gc; using io::pstr; using io::gstr; using io::gi; using io::print;
// end fast read template by CYJian
#define iter(i, a, b) for(int i=(a); i<(b); i++)
#define rep(i, a) iter(i, 0, a)
#define rep1(i, a) iter(i, 1, (a)+1)
#define log(a) cerr<<"\033[32m[DEBUG] "<<#a<<'='<<(a)<<" @ line "<<__LINE__<<"\033[0m"<<endl
#define all(a) a.begin(), a.end()
#define fi first
#define se second
#define pb push_back
#define mp make_pair
using ll=long long;
using pii=pair<int, int>;
//#define int ll
const int MOD = 1000000007;
const int MAX_N = 300005;
const int MAX_Q = 300005;
const int MAX_M = 25;
const int B_MO=547;
int n, m;
struct gpp { int a, b; };
// Mo stuff
ll ans[MAX_Q], fans[MAX_Q];
struct q {
int l, r, i;
const bool operator<(const q&o)const{
return(l/B_MO==o.l/B_MO)?((l/B_MO%2)?(r<o.r):(r>o.r)):(l<o.l);
}
} qs[MAX_Q];
struct q2{ int l, r, i, c; };
vector<q2> iv[MAX_N];
// Graph
int elist_node_begin[MAX_N];
int elist[MAX_N];
int dfnL[MAX_N], dfnR[MAX_N], father[MAX_N];
// toplogical sort
namespace dagify {
pii more_eges[MAX_N+MAX_M];
int elist2_node_begin[MAX_N];
int elist2[MAX_N+MAX_M];
int idx, topseq[MAX_N];
int ind[MAX_N];
void topify() {
rep1(u, n)
iter(i, elist2_node_begin[u], elist2_node_begin[u+1]) {
// cerr << u << ' ' << v << endl;
ind[elist2[i]]++;
}
queue<int> q;
rep1(u, n)
if(!ind[u])
q.push(u);
while(!q.empty()) {
int u = q.front();
q.pop();
topseq[idx++] = u;
// cerr << u << endl;
iter(i, elist2_node_begin[u], elist2_node_begin[u+1])
if(!(--ind[elist2[i]]))
q.push(elist2[i]);
}
// cerr << idx << endl;
assert(idx == n);
}
}
// data
namespace data_graph {
int reverse_graph_extra_node_begin[MAX_N];
int reverse_graph_extra[MAX_N];
pii extra1[MAX_M];
pii extra2[MAX_M];
pii eges[MAX_N], egegeg[MAX_N];
int seq[MAX_N];
int clo;
}
// prefix to node data
namespace data_PTN {
pii igtmp[(MAX_N+MAX_M*MAX_M)*2];
int ivStart[MAX_N], ivEnd[MAX_N];
int tivs;
gpp ivs[MAX_N*MAX_M];
}
// node to prefix critical data
namespace data_crit {
int critStart[MAX_N], critEnd[MAX_N];
int crit[MAX_N*MAX_M];
int crits3;
int cstck[MAX_M];
int td2;
pii *dfs_init_crit_L_ptr = data_graph::extra2;
}
// node to prefix data
namespace data_NTP {
int ivStart2[300005], ivEnd2[MAX_N];
gpp ivs2[MAX_N*MAX_M];
int tivs2;
gpp tmp[MAX_M];
int ctmpl;
bool ocd[MAX_N];
pii *init_NTP_L_ptr = data_graph::extra2;
}
// generic dfs
void dfs1(int u, int f) {
using namespace data_graph;
dfnL[u] = clo++;
father[u] = f;
seq[dfnL[u]] = u;
iter(i, elist_node_begin[u], elist_node_begin[u+1])
dfs1(elist[i], u);
dfnR[u] = clo;
}
// dfs to init node to prefix
void dfs_init_NTP(int u) {
using namespace data_graph;
using namespace data_NTP;
ivStart2[u] = tivs2;
rep(i, ctmpl)
ivs2[tivs2++] = tmp[i];
int octmpl = ctmpl;
assert(!ocd[u]);
ocd[u] = 1;
iter(i, data_crit::critStart[u], data_crit::critEnd[u]) {
int Lse = data_crit::crit[i];
// cerr << Lse << endl;
int pos = Lse;
while(pos && !ocd[pos]) {
ocd[pos] = 1;
pos = father[pos];
}
assert(pos);
if(Lse != pos) {
tmp[ctmpl] = {Lse, pos};
ivs2[tivs2++] = tmp[ctmpl];
ctmpl++;
}
}
ivs2[tivs2++] = {u, 0};
ivEnd2[u] = tivs2;
iter(i, elist_node_begin[u], elist_node_begin[u+1])
dfs_init_NTP(elist[i]);
ocd[u] = 0;
iter(i, octmpl, ctmpl) {
int lo = tmp[i].a, hi = tmp[i].b;
while(lo != hi) {
ocd[lo] = 0;
lo = father[lo];
}
}
ctmpl = octmpl;
}
// node to prefix crit init
// updated to some extent
void init_NTP() {
static bool us[MAX_N];
using namespace data_crit;
critStart[1] = critEnd[1] = crits3;
rep1(_, n-1) {
int u = dagify::topseq[_];
critStart[u] = crits3;
iter(j, critStart[father[u]], critEnd[father[u]]) {
if(!us[crit[j]]) {
us[crit[j]]=1;
crit[crits3++]=crit[j];
}
}
// crit[crits3++] = crit[j];
iter(i, data_graph::reverse_graph_extra_node_begin[u], data_graph::reverse_graph_extra_node_begin[u+1]) {
int v = data_graph::reverse_graph_extra[i];
if(!us[v]){
us[v]=1;
crit[crits3++]=v;
}
// crit[crits3++] = v;
iter(j, critStart[v], critEnd[v])
if(!us[crit[j]]){
us[crit[j]]=1;
crit[crits3++]=crit[j];
}
// crit[crits3++] = crit[j];
}
iter(i, critStart[u], crits3)
us[crit[i]]=0;
critEnd[u] = crits3;
}
dfs_init_NTP(1);
}
// prefix to node topological init
void init_PTN() {
using namespace data_graph;
using namespace data_PTN;
for(int i=n-1; i>=0; i--) {
int u = dagify::topseq[i];
int bg = 0;
igtmp[bg++] = {dfnL[u], -1};
igtmp[bg++] = {dfnL[u]+1, 1};
iter(ed, elist_node_begin[u], elist_node_begin[u+1])
iter(j, ivStart[elist[ed]], ivEnd[elist[ed]]) {
// cerr << bg << endl;
igtmp[bg++] = {ivs[j].a, -1};
igtmp[bg++] = {ivs[j].b, 1};
}
pii *L = lower_bound(extra1, extra1+m, pii(dfnL[u], 0));
// cout << i << ' ' << (R-2)->fi << ' ' << dfnR[i] << endl;
while(L->fi == dfnL[u]) {
// cout << i << ' ' << L->fi << endl;
if(!(dfnL[u] <= dfnL[L->se] && dfnL[L->se] < dfnR[u]))
iter(j, ivStart[L->se], ivEnd[L->se]) {
// cerr << bg << endl;
igtmp[bg++] = {ivs[j].a, -1};
igtmp[bg++] = {ivs[j].b, 1};
}
L++;
}
if(bg != 2) sort(igtmp, igtmp+bg);
ivStart[u] = tivs;
int las = -1, sm = 0;
rep(i, bg) {
if(igtmp[i].se == -1) {
if(sm == 0)
las = igtmp[i].fi;
sm++;
} else {
sm--;
if(!sm)
ivs[tivs++] = {las, igtmp[i].fi};
}
}
ivEnd[u] = tivs;
assert(!sm);
}
// if(n == 29)
// rep1(i, n) {
// cout << i << ": " << ivEnd[i] - ivStart[i] << endl;
// iter(j, ivStart[i], ivEnd[i]) cout << " " << ivs[j].a << ' ' << ivs[j].b << endl;
// }
}
// graph read
void initGraph() {
using namespace data_graph;
gi(n), gi(m);
assert(n<=300000);
int ecec = 0;
int ecec_dagify = 0;
int ecec_egegeg = 0;
iter(i, 2, n+1) {
int f; gi(f);
eges[ecec++] = {f, i};
// elist[f].push_back(i);
dagify::more_eges[ecec_dagify++] = {f, i};
}
{
sort(eges, eges+ecec);
int pt = 0;
rep1(i, n) {
elist_node_begin[i] = pt;
while(pt != ecec && eges[pt].fi == i) {
elist[pt] = eges[pt].se;
pt++;
}
}
elist_node_begin[n+1] = pt;
}
dfs1(1, 0);
rep(i, m) {
gi(extra1[i].fi), gi(extra1[i].se);
// cout << extra[i].fi << ' ' << extra[i].se << endl;
extra2[i] = {extra1[i].se, extra1[i].fi};
dagify::more_eges[ecec_dagify++] = extra1[i];
egegeg[ecec_egegeg++] = {extra1[i].se, extra1[i].fi};
// reverse_graph_extra[extra1[i].se].push_back(extra1[i].fi);
extra1[i].fi = dfnL[extra1[i].fi];
extra2[i].fi = dfnL[extra2[i].fi];
}
{
sort(dagify::more_eges, dagify::more_eges+ecec_dagify);
int pt = 0;
rep1(i, n) {
dagify::elist2_node_begin[i] = pt;
while(pt != ecec_dagify && dagify::more_eges[pt].fi == i) {
dagify::elist2[pt] = dagify::more_eges[pt].se;
pt++;
}
}
dagify::elist2_node_begin[n+1] = pt;
}
{
sort(egegeg, egegeg+ecec_egegeg);
int pt = 0;
rep1(i, n) {
reverse_graph_extra_node_begin[i] = pt;
while(pt != ecec_egegeg && egegeg[pt].fi == i) {
reverse_graph_extra[pt] = egegeg[pt].se;
pt++;
}
}
reverse_graph_extra_node_begin[n+1] = pt;
}
// rep1(i, n) cout << dfnL[i] << " \n"[i==n];
// rep1(i, n) cout << dfnR[i] << " \n"[i==n];
sort(extra1, extra1+m);
sort(extra2, extra2+m);
extra1[m] = extra2[m] = {1e9, 1e9};
dagify::topify();
init_PTN();
init_NTP();
// if(n == 29)
// rep1(i, n) {
// cout << i << ": " << data_NTP::ivEnd2[i] - data_NTP::ivStart2[i] << endl;
// iter(j, data_NTP::ivStart2[i], data_NTP::ivEnd2[i])
// cout << " " << data_NTP::ivs2[j].a << ' ' << data_NTP::ivs2[j].b << endl;
// }
}
// Mo structure
constexpr int BFDS=64;
constexpr int BFDS2 = BFDS*BFDS;
// general segtree
class FastDS {
public:
int b0[MAX_M*4], b1[MAX_N/BFDS+5], b2[MAX_N];
void _add(int l, int r) {
if(l/BFDS == r/BFDS)
iter(i, l, r)
b1[i]++;
else {
int nl = l+((l%BFDS)?(BFDS-l%BFDS):0);
int nr = r-(r%BFDS);
iter(i, l, nl)
b1[i]++;
iter(i, nr, r)
b1[i]++;
nl /= BFDS; nr /= BFDS;
iter(i, nl, nr)
b0[i]++;
}
}
void add(int l, int r) {
if(l/BFDS == r/BFDS)
iter(i, l, r)
b2[i]++;
else {
int nl = l+((l%BFDS)?(BFDS-l%BFDS):0);
int nr = r-(r%BFDS);
iter(i, l, nl)
b2[i]++;
iter(i, nr, r)
b2[i]++;
nl /= BFDS; nr /= BFDS;
_add(nl, nr);
}
}
int qry(int i) {
return b0[i/BFDS2] + b1[i/BFDS] + b2[i];
}
};
// modularization
namespace PrefixToNode {
using namespace data_PTN;
// TODO: multi-layer array decomposition
// how many nodes in the current prefix can reach a node?
FastDS inner;
// int query(int i) {
// return ;
//}
void add(int i) {
iter(j, ivStart[i], ivEnd[i])
inner.add(ivs[j].a, ivs[j].b);
}
}
namespace NodeToPrefix {
using namespace data_NTP;
// TODO: multi-layer tree decomp
// how many nodes in the current prefix can a node reach?
FastDS inner;
// int tmp[MAX_N];
//int query(int i) {
// return 0;
// return tmp[i];
// return ;
//}
void add(int i) {
iter(j, ivStart2[i], ivEnd2[i]) {
// int d = ivs2[j].a;
// while(d != ivs2[j].b) tmp[d]++, d = father[d];
if(ivs2[j].b)
inner.add(dfnL[ivs2[j].b]+1, dfnL[ivs2[j].a]+1);
else
inner.add(0, dfnL[ivs2[j].a]+1);
}
}
}
// Sweepline Mo ----------------------
//int globalquery(int i) {
// cout << "qry " << i << " res " << PrefixToNode::query(i) << " " << NodeToPrefix::query(i) << endl;
// return PrefixToNode::query(i) + NodeToPrefix::query(i);
//}
#define globalquery(i) (\
PrefixToNode::inner.b0[dfnL[i]/BFDS2]+PrefixToNode::inner.b1[dfnL[i]/BFDS]+PrefixToNode::inner.b2[dfnL[i]]\
+NodeToPrefix::inner.b0[dfnL[i]/BFDS2]+NodeToPrefix::inner.b1[dfnL[i]/BFDS]+NodeToPrefix::inner.b2[dfnL[i]]\
-NodeToPrefix::inner.b0[dfnR[i]/BFDS2]-NodeToPrefix::inner.b1[dfnR[i]/BFDS]-NodeToPrefix::inner.b2[dfnR[i]])
#define queryPTN(i) (\
PrefixToNode::inner.b0[dfnL[i]/BFDS2]+PrefixToNode::inner.b1[dfnL[i]/BFDS]+PrefixToNode::inner.b2[dfnL[i]])
#define queryNTP(i) (\
NodeToPrefix::inner.b0[dfnL[i]/BFDS2]+NodeToPrefix::inner.b1[dfnL[i]/BFDS]+NodeToPrefix::inner.b2[dfnL[i]]\
-NodeToPrefix::inner.b0[dfnR[i]/BFDS2]-NodeToPrefix::inner.b1[dfnR[i]/BFDS]-NodeToPrefix::inner.b2[dfnR[i]])
//int globalquery(int i) {
// int ans = -NodeToPrefix::inner.qry(dfnR[i]);
// int x = dfnL[i], y = dfnL[i]/BFDS, z = dfnL[i]/BFDS2;
// ans += PrefixToNode::inner.b2[x] + PrefixToNode::inner.b1[y] + PrefixToNode::inner.b0[z];
// ans += NodeToPrefix::inner.b2[x] + NodeToPrefix::inner.b1[y] + NodeToPrefix::inner.b0[z];
// return ans;
//}
void globaladd(int i) {
// cout << "add " << i << endl;
PrefixToNode::add(i);
NodeToPrefix::add(i);
}
ll lptrmoveL[MAX_N], lptrmoveR[MAX_N];
ll rptrmoveL[MAX_N], rptrmoveR[MAX_N];
signed main() {
initGraph();
// return 0;
// if(n == 19) {
// rep1(i, n) cout << i << " " << dfnL[i] << " " << dfnR[i] << endl;
// }
int q; gi(q);
// q = min(q, 100000);
// if(q == 1) {
// int l, r; gi(l), gi(r);
// int ans = 0;
// iter(i, l, r+1) {
// ans += globalquery(i);
// globaladd(i);
// }
// cout << ans + r - l + 1 << endl;
// return 0;
// }
rep(i,q){
gi(qs[i].l),gi(qs[i].r),qs[i].i=i;
// fans[i] = qs[i].r - qs[i].l + 1;
}
sort(qs, qs + q);
int cl = 1, cr = 0;
rep(i, q) {
// cout << qs[i].l << ' ' << qs[i].r << endl;
if(cr < qs[i].r) {
iv[cl-1].push_back({cr + 1, qs[i].r, i, 0}); // -1
cr = qs[i].r;
}
if(qs[i].l < cl) {
iv[cr].push_back({qs[i].l, cl - 1, i, 1}); // 1
cl = qs[i].l;
}
if(qs[i].r < cr) {
iv[cl-1].push_back({qs[i].r + 1, cr, i, 2}); // 1
cr = qs[i].r;
}
if(cl < qs[i].l) {
iv[cr].push_back({cl, qs[i].l - 1, i, 3}); // -1
cl = qs[i].l;
}
}
for(int i = 1; i <= n; i++) {
rptrmoveR[i] = rptrmoveR[i-1] + queryPTN(i);
rptrmoveL[i] = rptrmoveL[i-1] + queryPTN(i);
// ansPrev[i] = ansPrev[i-1] + globalquery(i);
globaladd(i);
lptrmoveR[i] = lptrmoveR[i-1] + queryNTP(i);
lptrmoveL[i] = lptrmoveL[i-1] + queryNTP(i);
for(q2& c14compliance:iv[i]) {
const int& l=c14compliance.l;
const int& r=c14compliance.r;
const int& x=c14compliance.i;
const int& c=c14compliance.c;
ll g = 0;
if(c == 0)
for(int p = l; p <= r; p++)
g -= queryPTN(p);
if(c == 1)
for(int p = l; p <= r; p++)
g += queryNTP(p);
if(c == 2)
for(int p = l; p <= r; p++)
g += queryPTN(p);
if(c == 3)
for(int p = l; p <= r; p++)
g -= queryNTP(p);
ans[x] += g;
}
// ansThere[i] = ansThere[i - 1] + globalquery(i);
}
// cout << rptrmoveR[n]-(rptrmoveL[n]-rptrmoveL[n-1]) << ' ' << rptrmoveR[n-1] << endl;
// return 0;
ll curans=0;
cl = 1, cr = 0;
for(int i = 0; i < q; i++) {
curans += ans[i];
if(cr < qs[i].r) {
curans += rptrmoveR[qs[i].r] - rptrmoveR[cr];
cr = qs[i].r;
}
if(qs[i].l < cl) {
curans -= lptrmoveL[cl-1] - lptrmoveL[qs[i].l-1];
cl = qs[i].l;
}
if(qs[i].r < cr) {
curans -= rptrmoveL[cr] - rptrmoveL[qs[i].r];
cr = qs[i].r;
}
if(cl < qs[i].l) {
curans += lptrmoveR[qs[i].l-1] - lptrmoveR[cl-1];
cl = qs[i].l;
}
fans[qs[i].i] += curans;
}
for(int i = 0; i < q; i++) {
print(fans[i]);
pc('\n');
}
}
Tester's Solution
//Radoslav's solution
#include <bits/stdc++.h>
#define endl '\n'
#define SZ(x) ((int)x.size())
#define ALL(V) V.begin(), V.end()
#define L_B lower_bound
#define U_B upper_bound
#define pb push_back
using namespace std;
template<class T, class T1> int chkmin(T &x, const T1 &y) { return x > y ? x = y, 1 : 0; }
template<class T, class T1> int chkmax(T &x, const T1 &y) { return x < y ? x = y, 1 : 0; }
const int MAXN = 1 << 19;
const int MAXM = 22;
inline int64_t hilbertOrder(int x, int y, int pw, int rot) {
if (pw == 0) {
return 0;
}
int hpow = 1 << (pw-1);
int seg = (x < hpow) ? (
(y < hpow) ? 0 : 3
) : (
(y < hpow) ? 1 : 2
);
seg = (seg + rot) & 3;
const int rotateDelta[4] = {3, 0, 0, 1};
int nx = x & (x ^ hpow), ny = y & (y ^ hpow);
int nrot = (rot + rotateDelta[seg]) & 3;
int64_t subSquareSize = int64_t(1) << (2*pw - 2);
int64_t ans = seg * subSquareSize;
int64_t add = hilbertOrder(nx, ny, pw-1, nrot);
ans += (seg == 1 || seg == 2) ? add : (subSquareSize - add - 1);
return ans;
}
struct query {
int l, r, idx;
int64_t order;
query() { }
query(int _l, int _r, int i) {
l = _l;
r = _r;
idx = i;
order = hilbertOrder(l, r, 21, 0);
}
};
int read_int();
int n, q_cnt, m;
vector<int> adj[MAXN];
vector<pair<int, int> > bonus[MAXN];
pair<int, int> e[MAXM];
vector<pair<int, int>> q[MAXN];
query qs[MAXN];
void read() {
n = read_int();
m = read_int();
for(int i = 2; i <= n; i++) {
int tmp;
tmp = n - read_int() + 1;
adj[tmp].pb(n - i + 1);
}
for(int i = 0; i < m; i++) {
int u, v;
u = n - read_int() + 1;
v = n - read_int() + 1;
e[i] = {u, v};
bonus[u].pb({v, i});
}
q_cnt = read_int();
for(int i = 0; i < q_cnt; i++) {
int l, r;
r = n - read_int() + 1;
l = n - read_int() + 1;
q[r].pb({l, i});
qs[i] = query(l, r, i);
}
}
int ver[MAXN];
int st[MAXN], en[MAXN], dfs_time;
int reachable[MAXN];
void pre_dfs(int u) {
st[u] = ++dfs_time;
ver[st[u]] = u;
reachable[u] = 0;
for(auto ed: bonus[u]) {
reachable[u] |= (1 << ed.second);
}
for(int v: adj[u]) {
pre_dfs(v);
reachable[u] |= reachable[v];
}
en[u] = dfs_time;
}
int64_t sum[MAXM][MAXN];
int cnt[MAXM][MAXN];
int lb_pntr[MAXM][MAXN], sz[MAXN];
bool upper(int u, int v) {
return st[u] <= st[v] && st[v] <= en[u];
}
int64_t answer[MAXN];
int64_t query_important(int i, int l) {
int pos = lb_pntr[i][l], csz = sz[i];
if(csz - 1 == pos) return 0;
pos++;
int c = 0, p = pos - 1;
int64_t s = 0;
while(pos <= csz) {
c += cnt[i][pos];
s += sum[i][pos];
pos += pos & -pos;
}
return s - c * 1ll * p;
}
void update_important(int i, int r) {
int pos = lb_pntr[i][r];
if(pos != 0) {
int add = pos;
while(pos) {
sum[i][pos] += add;
cnt[i][pos]++;
pos -= pos & -pos;
}
}
}
int64_t delta[MAXN];
vector<tuple<int, int, int, int>> lq[MAXN], rq[MAXN];
int64_t valid[MAXN], invalid[MAXN];
struct sqrt_decompostion {
static const int B = 200;
int block_add[1601];
int single_add[MAXN];
void add(int l, int r) {
if(l / B == r / B) {
while(l <= r) single_add[l++]++;
return;
}
while(l % B != 0) single_add[l++]++;
while(r % B != B - 1) single_add[r--]++;
for(int from = l / B, to = r / B; from <= to; from++) {
block_add[from]++;
}
}
int query(int pos) {
return block_add[pos / B] + single_add[pos];
}
} ds1, ds2;
int query_mo_valid(int p) {
return ds2.query(en[p]) - ds2.query(st[p]);
}
int query_mo_invalid(int p) {
return ds1.query(st[p]);
}
void update_mo(int p) {
ds1.add(st[p], en[p]);
ds2.add(st[p], n);
}
void solve() {
dfs_time = 0;
pre_dfs(n);
for(int i = 0; i < m; i++) {
for(int _ = n; _ >= 1; _--) {
int u = ver[_];
for(auto ed: bonus[u]) {
reachable[u] |= reachable[ed.first];
}
for(int v: adj[u]) {
reachable[u] |= reachable[v];
}
}
}
vector<int> perm(m);
iota(ALL(perm), 0);
sort(ALL(perm), [](int i, int j) { return make_pair(st[e[i].second], i) < make_pair(st[e[j].second], j); });
for(int i = 0; i < m; i++) {
vector<int> subtree(en[e[i].second] - st[e[i].second] + 1);
for(int _ = st[e[i].second]; _ <= en[e[i].second]; _++) {
subtree[_ - st[e[i].second]] = ver[_];
}
sort(ALL(subtree));
int p = 0;
for(int q = 1; q <= n; q++) {
while(p < SZ(subtree) && subtree[p] < q) p++;
lb_pntr[i][q] = p;
}
sz[i] = SZ(subtree) + 1;
}
// M = 0 reduction
for(int r = 1; r <= n; r++) {
int last = -1;
for(int i: perm) {
if((reachable[r] & (1 << i)) && !upper(r, e[i].second) && st[e[i].second] > last) {
update_important(i, r);
last = en[e[i].second];
}
}
for(auto it: q[r]) {
int idx = it.second, l = it.first;
for(int i: perm) {
answer[idx] += query_important(i, l);
}
}
}
sort(qs, qs + q_cnt, [](query a, query b) { return a.order < b.order; });
int cl = 1, cr = 0;
for(int i = 0; i < q_cnt; i++) {
if (cr < qs[i].r) {
rq[cl - 1].pb({cr + 1, qs[i].r, i, -1});
cr = qs[i].r;
}
if (qs[i].l < cl) {
lq[cr].pb({qs[i].l, cl - 1, i, 1});
cl = qs[i].l;
}
if (qs[i].r < cr) {
rq[cl - 1].pb({qs[i].r + 1, cr, i, 1});
cr = qs[i].r;
}
if (cl < qs[i].l) {
lq[cr].pb({cl, qs[i].l - 1, i, -1});
cl = qs[i].l;
}
}
for(int r = 1; r <= n; r++) {
valid[r] = valid[r - 1] + query_mo_valid(r);
update_mo(r);
invalid[r] = invalid[r - 1] + query_mo_invalid(r);
for(auto it: lq[r]) {
int ql = get<0>(it), qr = get<1>(it), i = get<2>(it), c = get<3>(it);
int64_t contribution = 0;
for(int p = ql; p <= qr; p++) {
contribution += query_mo_invalid(p);
}
delta[i] += c * contribution;
}
for(auto it: rq[r]) {
int ql = get<0>(it), qr = get<1>(it), i = get<2>(it), c = get<3>(it);
int64_t contribution = 0;
for(int p = ql; p <= qr; p++) {
contribution += query_mo_valid(p);
}
delta[i] += c * contribution;
}
}
int64_t ans = 0;
cl = 1, cr = 0;
for(int i = 0; i < q_cnt; i++) {
if(cr < qs[i].r) {
ans += valid[qs[i].r] - valid[cr];
cr = qs[i].r;
}
if(qs[i].l < cl) {
ans -= invalid[cl - 1] - invalid[qs[i].l - 1];
cl = qs[i].l;
}
if(qs[i].r < cr) {
ans -= valid[cr] - valid[qs[i].r];
cr = qs[i].r;
}
if(cl < qs[i].l) {
ans += invalid[qs[i].l - 1] - invalid[cl - 1];
cl = qs[i].l;
}
ans += delta[i];
answer[qs[i].idx] += ans;
}
for(int i = 0; i < q_cnt; i++) {
cout << answer[i] << endl;
}
}
int main() {
ios_base::sync_with_stdio(false);
cin.tie(nullptr);
read();
solve();
return 0;
}
const int maxl = 100000;
char buff[maxl];
int ret_int, pos_buff = 0;
void next_char() { if(++pos_buff == maxl) fread(buff, 1, maxl, stdin), pos_buff = 0; }
int read_int()
{
ret_int = 0;
for(; buff[pos_buff] < '0' || buff[pos_buff] > '9'; next_char());
for(; buff[pos_buff] >= '0' && buff[pos_buff] <= '9'; next_char())
ret_int = ret_int * 10 + buff[pos_buff] - '0';
return ret_int;
}