思路
转化题意
注意到这个“没有出现过的最小正整数”很不舒服。
改定义：
\[b_{x,y}=a_{x+1,y+1}-1 \]则我们有
\[b_{x,y}=\operatorname{mex}(\{b_{x,k}|0\le k<y\}\cup\{b_{k,y}|0\le k<x\}) \]而我们所求的即
\[\sum_{i=x_1}^{x_2}\sum_{j=y_1}^{y_2}[a_{i,j}\le k]a_{i,j} \\=\sum_{i=x_1-1}^{x_2-1}\sum_{j=y_1-1}^{y_2-1}[b_{i,j}<k](b_{i,j}+1)\]设
\[f_k(x,y)=\sum_{i<x}\sum_{j<y}[b_{i,j}<k] \]\[g_k(x,y)=\sum_{i<x}\sum_{j<y}[b_{i,j}<k]b_{i,j} \]\[h_k=f_k+g_k \]于是即求
\[h_k(x_2,y_2)-h_k(x_2,y_1-1)-h_k(x_1-1,y_2)+h_k(x_1-1,y_1-1) \]因此现在的任务就是快速计算 \(f,g\)。
我们指出如下事实：
\[b_{x,y}=x\oplus y \]其中 \(\oplus\) 表示异或。
证明：
构造 \(2\) 堆石子的 nim 游戏 \(a_1,a_2\)，设其 SG 值为 \(c_{a_1,a_2}\)。
由 SG 函数的定义式，我们得到 \(c_{a_1,a_2}=\operatorname{mex}(\{c_{a_1,k}|0\le k<a_2\}\cup\{c_{k,a_2}|0\le k<a_1\})\)。
于是使用数学归纳法容易证明，\(b_{x,y}=c_{x,y}\)。
由 SG 定理，我们得到 \(c_{a_1,a_2}=a_1\oplus a_2\)。
因此，\(b_{x,y}=x\oplus y\)。
于是
\[f_k(x,y)=\sum_{i<x}\sum_{j<y}[i\oplus j<k] \]\[g_k(x,y)=\sum_{i<x}\sum_{j<y}[i\oplus j<k](i\oplus j) \]这不就是储能表吗！
因此数位 dp 是能做的。当然，这并不优美。
考虑有没有更加优美的做法。
异或卷积
我们发现和式内部仅和 \(i\oplus j\) 有关，而枚举的却是 \(i,j\)。
因此考虑能否对每个 \(z=x\oplus y\)，反推出有多少个合法的 \((x,y)\) 有序二元组满足此点。
定义数列 \(f,g\) 的异或卷积为数列
\[h_z=\sum_{x\oplus y=z}f_xg_y \]不妨记为 \(h=f\times_\oplus g\)。
显然这个柿子的组合意义为枚举 \(f,g\) 对应下标异或值为 \(z\) 的方案数，其中 \(f,g\) 分别表示其作为子式时变量异或和为某值的方案数。
考虑计算 \(f_k(x,y),g_k(x,y)\) 的过程。
\[f_k(x,y)=\sum_{i<x}\sum_{j<y}[i\oplus j<k]=\sum_{z<k}\sum_{i<x}\sum_{j<y}[i\oplus j=z] \]\[g_k(x,y)=\sum_{i<x}\sum_{j<y}[i\oplus j<k](i\oplus j)=\sum_{z<k}z\sum_{i<x}\sum_{j<y}[i\oplus j=z] \]于是记
\[f_i=[i<x] \]\[g_i=[i<y] \]\[h=f\times_\oplus g \]则
\[h_z=\sum_{i<x}\sum_{j<y}[i\oplus j=z] \]因此
\[f_k(x,y)=\sum_{z<k}h_z \]\[g_k(x,y)=\sum_{z<k}zh_z \]特殊形式的异或卷积
接下来的任务就是观察如何快速计算与记录这个异或卷积的结果了。
因为值域较大，所以异或卷积显然不能暴力使用基于快速沃尔什变换的集合对称差卷积实现，我们考虑观察性质。
不妨设在讨论范围内，一个数列非零值的下标总是小于 \(2^m\)。
我们发现，对于数列 \(f,g\)，其前一半或后一半下标对应的值总是相等，另一半下标对应地递归满足相同性质。
简单讨论容易验证，异或卷积保持这个性质不变，因此 \(h=f\times_\oplus g\) 也满足这个性质。
于是我们可以使用 \(m+1\) 对数，来表示这样的一个数列（记录是前一半还是后一半为定值；记录具体是什么值）；对于这么一对数列的异或卷积，我们总可以在 \(O(m)\) 的时间内实现，并返回一个仍保持此性质的数列。
这个结论有一个经典的习题：loj3073「2019 集训队互测 Day 2」序列。感兴趣可以去做一做。
因此，我们就可以优美地在 \(O(\log v)\) 的时间内解决每一问了。
Code
翻译没翻到位，原来是对 \(10^9+7\) 取模的……
代码跑得很快。
下面是核心代码。
const ullt Mod=1e9+7;
typedef ConstMod::mod_ullt<Mod>modint;
typedef std::vector<std::pair<bol,modint> >vec;
typedef std::vector<modint>modvec;
const uint M=32;ullt Len[M+1];
modvec User(vec a){
modvec Ans(M+1);
for(uint i=0;i<=M;i++)Ans[i]=a[i].second*Len[i];
for(uint i=M;i;i--)Ans[i-1]+=Ans[i];
return Ans;
}
vec mul(vec A,vec B){
vec Ans(M+1);modvec UserA=User(A),UserB=User(B);
for(uint i=0;i<M;i++)
{
Ans[i].first=A[i].first==B[i].first;
Ans[i+1].second+=Ans[i].second+A[i].second*B[i].second*Len[i];
Ans[i].second+=A[i].second*UserB[i+1]+B[i].second*UserA[i+1];
}
Ans[M].second+=A[M].second+B[M].second;
return Ans;
}
vec get(ullt r){
vec Ans(M+1);
for(uint i=0;i<M;i++)if(r>=Len[i])Ans[i]={0,1},r-=Len[i];else Ans[i]={1,0};
Ans[M].second=r;return Ans;
}
modint sum(ullt l,ullt r){return(l+r-1)*(r-l)/2;}
modint find(vec A,ullt k){
modint ans;ullt from=0;
for(uint i=0;i<M;i++)if(A[i].first)
{
if(k>Len[i])ans+=sum(from+Len[i]+1,from+k+1)*A[i].second,k=Len[i];
}else{
if(k<=Len[i]){
ans+=sum(from+1,from+k+1)*A[i].second;
return ans;
}
ans+=sum(from+1,from+Len[i]+1)*A[i].second;
k-=Len[i];
from+=Len[i];
}
if(k)ans+=(from+1)*A[M].second;
return ans;
}