如果n、m、k都是2的冪次方，答案非常好統計。於是容易想到數位dp，考慮每一位是否卡限制即可，即設f[i][0/1][0/1][0/1]為第i位是/否卡n、m、k的限制時，之前的位的總貢獻；g[i][0/1][0/1][0/1]為第i位是/否卡n、m、k的限制時，之前的位的方案數。為了方便可以改為統計小於k的貢獻再減去。

　　莫名其妙的搞錯了很多地方，簡直調一年，不知道在幹啥。

　　人醜常數大，根本沒辦法。

　　（突然發現以前大部分數位dp都是直接按位計數就搞出來了……這個題應該也行。

#include<iostream> 
#include<cstdio>
#include<cmath>
#include
 
<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;
#define ll long long
#define N 64
ll read()
{
    ll x=0,f=1;char c=getchar();
    while (c<‘0‘||c>‘9‘) {if (c==‘-‘) f=-1;c=getchar();}
    while (c>=‘0‘&&c<=‘9‘) x=(x<<1)+(x<<3)+(c^48
 
),c=getchar();
    return x*f;
}
int T,p,f[N][2][2][2],g[N][2][2][2],a[N],b[N],c[N],q[N];
ll n,m,k;
void inc(int &x,int y,int p){x+=y;if (x>=p) x-=p;}
void calc(ll n,ll m,ll k,int p)
{
    memset(f,0,sizeof(f));memset(g,0,sizeof(g));
    int t=-1;ll x=max(max(n,m),k);
    while (x) t++,x>>=1
 
;
    for (int i=0;i<=t;i++) a[i]=(n&(1ll<<i))>0;
    for (int i=0;i<=t;i++) b[i]=(m&(1ll<<i))>0;
    for (int i=0;i<=t;i++) c[i]=(k&(1ll<<i))>0;
    g[t+1][1][1][1]=1;
    for (register int i=t;~i;i--)
        for (register int x=0;x<=1;x++)
            for (register int y=0;y<=1;y++)
                for (register int z=0;z<=1;z++)
                    for (register int u=x;u<=1;u++)
                        for (register int v=y;v<=1;v++)
                            for (register int w=z;w<=1;w++)
                            {
                                int t=(!w|c[i]&z)?(((!u|a[i]&x)&(!v|b[i]^y))+((!v|b[i]&y)&(!u|a[i]^x))):0;
                                inc(f[i][x][y][z],1ll*q[i]*g[i+1][u][v][w]%p*t%p,p);
                                if (!w|c[i]^z) t+=(((!u|a[i]&x)&(!v|b[i]&y))+((!u|a[i]^x)&(!v|b[i]^y)));
                                inc(f[i][x][y][z],1ll*f[i+1][u][v][w]*t%p,p);
                                inc(g[i][x][y][z],1ll*g[i+1][u][v][w]*t%p,p);
                            }
}
int main()
{
#ifndef ONLINE_JUDGE
    freopen("bzoj4513.in","r",stdin);
    freopen("bzoj4513.out","w",stdout);
    const char LL[]="%I64d\n";
#else
    const char LL[]="%lld\n";
#endif
    T=read();
    while (T--)
    {
        n=read()-1,m=read()-1,k=read(),p=read();
        q[0]=1;for (int i=1;i<=63;i++) q[i]=(q[i-1]<<1)%p;
        int ans=0;
        calc(n,m,max(n,m)<<1,p);
        for (int x=0;x<=1;x++)
            for (int y=0;y<=1;y++)
                for (int z=0;z<=1;z++)
                inc(ans,f[0][x][y][z],p),inc(ans,p-1ll*k%p*g[0][x][y][z]%p,p);
        calc(n,m,k,p);
        for (int x=0;x<=1;x++)
            for (int y=0;y<=1;y++)
                for (int z=0;z<=1;z++)
                inc(ans,p-f[0][x][y][z],p),inc(ans,1ll*k%p*g[0][x][y][z]%p,p);
        cout<<ans<<endl;
    }
    return 0;
}

BZOJ4513 SDOI2016儲能表（數位dp）