LGPP4067题解

直接计算太困难了，考虑转化。

可以转化为原储能表的和减去原储能表中不大于 \(k\) 的部分，然后减去数量乘上 \(k\) 即可。零次和与一次和可以同时统计。

原储能表的元素和非常好算啊，直接拆位即可，复杂度 \(O(\log n)\)。

我们假设存在一个 \(t\) 满足 \(2^{t-1}<k\leq 2^t\)。

将行列分为长度为 \(2^t\) 的若干段，发现每一行只有第 \(i\) 行对应第 \(i\) 列才有可能小于 \(k\)。（因为只有此时高位异或才为 \(0\)）

而让一个数去异或 \([0,2^t-1]\) 中的所有数，得到的结果一定还是 \([0,2^t-1]\) 中的所有数（对于 \(a\ne b\) 有 \(a\texttt{ xor } x\ne b\texttt{ xor } x\)）。于是只要两端中有一段是整的那么贡献都非常好算。

我们还需要考虑的是不完整的段对不完整的段的贡献。

长度不是 \(2\) 的幂好像就没有优秀的性质了，于是我们直接把这两段都拆成 \(\log k\) 段长度为 \(2\) 的幂的段。

考虑两个区间 \([a\times 2^x,(a+1)\times 2^x-1]\) 和 \([b\times 2^y,(b+1)\times 2^y-1]\) 的答案。（方便计算有 \(x>y\)）

如果高位异或起来大于 \(k\) 的对应的高位，那么贡献显然为 \(0\)。

如果高位异或起来小于 \(k\) 的对应的高位那么贡献显然是 \(2^{x+y}\)。

考虑到 \(b\) 的 \([y,x-1]\) 位异或到第一个区间上其实是没有任何影响的（因为原本是哪些数那还是哪些数），剩下的部分的贡献就是 \((k\bmod 2^{x})\times 2^y\)。

#include<cstdio>
#include<cctype>
typedef unsigned long long ull;
const int M=64;
int T;ull n,m,k,mod;int l1,l2,k1[M],k2[M];ull a[M],b[M];
inline ull mul(ull n,ull m){
	return 1ull*(n%mod)*(m%mod)%mod;
}
inline ull GetSum(const ull&n){
	return n&1?mul(n,n+1>>1):mul(n>>1,n+1);
}
inline ull read(){
	ull n(0);char s;while(!isdigit(s=getchar()));while(n=n*10+(s&15),isdigit(s=getchar()));return n;
}
inline ull calc(const ull&n,const ull&m){
	ull ans(0);
	for(int i=0;i<64;++i)if((1ull<<i)<n||(1ull<<i)<m){
		ull c1=n&(1ull<<i)-1,c2=m&(1ull<<i)-1;
		if(n>>i&1)c1=(n+(1ull<<i)-c1>>1)-(1ull<<i)+c1;else c1=n-c1>>1;
		if(m>>i&1)c2=(m+(1ull<<i)-c2>>1)-(1ull<<i)+c2;else c2=m-c2>>1;
		ans=(ans+mul(mul(c1,m-c2)+mul(n-c1,c2),1ull<<i))%mod;
	}
	return ans;
}
inline void f(const ull&a,const int&k1,const ull&b,const int&k2,const ull&k,ull&ans,ull&cnt){
	if(((a^b)>>k1)>(k>>k1))return;
	if(((a^b)>>k1)<(k>>k1)){
		ans=(ans+mul(mod-(GetSum((1ull<<k1)-1)+mul(((a>>k1)^(b>>k1))<<k1,1ull<<k1))%mod,1ull<<k2))%mod;
		cnt=(cnt+mul(mod-(1ull<<k1)%mod,1ull<<k2))%mod;
	}
	else{
		ans=(ans+mul(mod-(GetSum(k&(1ull<<k1)-1)+mul(((a>>k1)^(b>>k1))<<k1,(k&(1ull<<k1)-1)+1))%mod,1ull<<k2))%mod;
		cnt=(cnt+mul(mod-((k&(1ull<<k1)-1)+1)%mod,1ull<<k2))%mod;
	}
}
signed main(){
	T=read();
	while(T--){
		n=read();m=read();k=read();mod=read();if(n>m)n^=m^=n^=m;
		ull ans=calc(n,m),cnt=mul(n,m),len(1);while(len<k)len<<=1;
		if((n+len-1)/len!=(m+len-1)/len){
			ans=(ans+mul(n,mod-GetSum(k)))%mod;cnt=(cnt+mul(n,mod-(k+1)%mod))%mod;
		}
		else{
			ull x1=((n+len-1)/len-1)*len,x2=((m+len-1)/len-1)*len;l1=l2=0;
			ans=(ans+mul(x1,mod-GetSum(k)))%mod;cnt=(cnt+mul(x1,mod-(k+1)%mod))%mod;
			for(int i=63;i>=0;--i){
				if(x1+(1ull<<i)<=n)a[++l1]=x1&len-1,k1[l1]=i,x1+=1ull<<i;
				if(x2+(1ull<<i)<=m)b[++l2]=x2&len-1,k2[l2]=i,x2+=1ull<<i;
			}
			for(int i=1;i<=l1;++i)for(int j=1;j<=l2;++j){
				if(k1[i]>k2[j])f(a[i],k1[i],b[j],k2[j],k,ans,cnt);
				else f(b[j],k2[j],a[i],k1[i],k,ans,cnt);
			}
		}
		printf("%llu\n",(ans+mul(mod-k%mod,cnt))%mod);
	}
}