Template for Mo's algorithm, uses two (or four if different behaviour on each side) functors for addition and deletion. also two kinds of ordering options are available, one is the classical Mo's algorithm ((n+q)sqrt n time complexity), the other is Hilbert-Mo (n sqrt q time complexity). You can define other ordering options as you wish

mosalg.cpp

#include<iostream>
#include<algorithm>
#include<vector>
#include<cmath>
using namespace std;

struct Query
{
	size_t idx,s,e;
	long long ord;
};

struct default_ordering
{
	int sq;
	default_ordering(int sz):sq(sqrt(sz)){}
	void operator()(Query& in)
	{
		in.ord=in.s/sq;
	}
	bool operator()(Query a,Query b)
	{
		if(a.ord!=b.ord)
			return a.ord<b.ord;
		return a.e<b.e;
	}
};

struct hilbert_ordering
{
	int pw;
	hilbert_ordering(int sz):
	pw(int(log2l(sz)+1)){}
	int64_t hilbertOrder(int x, int y, int pow, int rotate);
	void operator()(Query& in)
	{
		in.ord=hilbertOrder(in.s,in.e,pw,0);
	}
	bool operator()(Query a,Query b)
	{
		if(a.ord!=b.ord)
			return a.ord<b.ord;
		return a.e<b.e;
	}
};

int64_t
hilbert_ordering::hilbertOrder(int x, int y, int pow, int rotate) {
	if (pow == 0) {
		return 0;
	}
	int hpow = 1 << (pow-1);
	int seg = (x < hpow) ? (
		(y < hpow) ? 0 : 3
	) : (
		(y < hpow) ? 1 : 2
	);
	seg = (seg + rotate) & 3;
	const int rotateDelta[4] = {3, 0, 0, 1};
	int nx = x & (x ^ hpow), ny = y & (y ^ hpow);
	int nrot = (rotate + rotateDelta[seg]) & 3;
	int64_t subSquareSize = int64_t(1) << (2*pow - 2);
	int64_t ans = seg * subSquareSize;
	int64_t add = hilbertOrder(nx, ny, pow-1, nrot);
	ans += (seg == 1 || seg == 2) ? add : (subSquareSize - add - 1);
	return ans;
}

template<class T, class R,
class AddL, class SubL,
class AddR=AddL, class SubR=SubL>
struct Mo
{
	size_t n,curr;
	vector<T> values;
	vector<Query> queries;
	vector<R> ans;
	AddL addl; AddR addr;
	SubL subl; SubR subr;
	Mo(size_t nn):n(nn),curr(0){}
	Mo(vector<T>& vec)
	:n(size(vec)),values(vec),curr(0){}
	void addquery(size_t s, size_t e)
	{
		queries.push_back({curr++,s,e});
		ans.push_back(R{});
	}
	template<class Ord=default_ordering>
	vector<R> process()
	{
		Ord order(n);
		for(auto&q:queries)order(q);
		sort(begin(queries),end(queries),
		order);
		R res{};
		size_t s=queries[0].s, e = queries[0].e;
		//cout<<s<<" "<<e<<endl;
		for(int i=s; i<=e; i++)res=addr(res,values[i]);
		ans[queries[0].idx]=res;
		for(int i=1;i<queries.size(); i++)
		{
			while(s < queries[i].s)
			res=subl(res,values[s++]);
			while(s > queries[i].s)
			res=addl(res,values[--s]);
			while(e < queries[i].e)
			res=addr(res,values[++e]);
			while(e > queries[i].e)
			res=subr(res,values[e--]);
			//cout<<s<<" "<<e<<endl;
			ans[queries[i].idx] = res;
		}
		return ans;
	}
};

// hilbert-mo from https://codeforces.com/blog/entry/61203