hdu 5648 DZY Loves Math 组合数+深搜(子集法)

题目链接:http://acm.hdu.edu.cn/showproblem.php?pid=5648

题意:给定n,m(1<= n,m <= 15,000),求Σgcd(i|j,i&j);(1 <= i <= n,1<=j<=m);

至多三组数据，至多两组数据max(n,m) > 2000.至多一组数据max(n,m) > 8000;

很多题解是用递推打表，将数据压缩250倍，即[i][j]:代表[1...250*i][1...250*j]，之后零头就直接暴力求解；这样时间复杂度为O(T*250*max(n,m));

但是这并不完美。。代码长度接近5W b..如果卡代码长度的话，就呵呵了；

事先说明代码是在BC看的AC代码，来源于NanoApe。特在此orz..

思路:题解讲用枚举子集法，但是当时并没有看懂。。没有去细想在这个gcd里面是要确定什么，枚举什么。这就是深搜的思路了；

原本看到gcd想的就是莫队反演。。莫队反演是容易求出同一个gcd里面的对数；但是这道题gcd里面并不是原始的a,b;而是a,b的位或和位与运算以后的值；

并且数据规模不大，(至少从组数及n,m最大值的限制上来看)；

现在我们就来枚举i|j,因为i&j能够在i|j枚举1，这就是题解上面的子集枚举的对象；

下面分情况来看如何深搜；i表示两个数位或运算后的值；并且默认n >= m

<1>当i<= m时，这时易知对于我们实际的值x,y一定是<=m的；直接在范围内枚举,即可dfs1();

dfs1里面a 代表位或的值，b代表位与的值；容易看出b其实位元1是在a中选的；至于对ans的贡献是每一个a位元为1，b为0的位可以是两种，所以1<<(c[a]-c[b]);对于b为0的情况，会出现x = 0 || y = 0所以减去这二者；(d在dfs2中解释)；

注意并不是i <= n作为循环结束的标志，位或最大为2*n-1,这也是开始运算出mx的另一个用处；

<2>当i > m时，没有a里面的每一个1，一定有一个数(x||y)的位元是1的，这就可以分成x,y,x&&y;对应三个if;
对于这个位或值，如果前面有一个位1，原本x可以承担(即x += 1<<p 仍然<=n)但是并有没有取，而是被y承担了,这时之后任意位元为1的x均可以承担，这时对于x来说，进入dfs1就可以随便取a里面剩下的位了；这样就保证了所取得的值<=n;这就是为什么在判断是否能进入dfs1对后面的位元1进行组合数取的条件；即加的是(1<<(p+1))-1;

<3>对于dfs2中的p == -1的情况，是否需要求解?答案是要的，因为能到p == -1表示里面的x,y是符合情况的；如果不符合情况，就会直接在前面就无法dfs2下去了；直接计算即可；至于d的含义，因为当a > m时，前面大的位元1一定是x承担的，并且这承担的在后面组合数枚举的时候，虽然没有在高位上枚举，但是c[a]还是记录了这高位的，d即表示这没有枚举的高位a比b多出来的个数；由于枚举到的位可能在m范围内，但是之后却没有枚举，所以出现b = x&y;

很神奇~~~hhhhh

1326ms 1824k 自己手写__builtin_popcount()用时1263ms 因为库函数是查表的，并没有真正计算；

#include<bits/stdc++.h>
using namespace std;
#define rep0(i,l,r) for(int i = (l);i < (r);i++)
#define rep1(i,l,r) for(int i = (l);i <= (r);i++)
#define rep_0(i,r,l) for(int i = (r);i > (l);i--)
#define rep_1(i,r,l) for(int i = (r);i >= (l);i--)
#define MS0(a) memset(a,0,sizeof(a))
#define MS1(a) memset(a,-1,sizeof(a))
#define MSi(a) memset(a,0x3f,sizeof(a))
#define inf 0x3f3f3f3f
#define lson l, m, rt << 1
#define rson m+1, r, rt << 1|1
typedef pair<int,int> PII;
#define A first
#define B second
#define MK make_pair
typedef __int64 ll;
template<typename T>
void read1(T &m)
{
    T x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
    m = x*f;
}
template<typename T>
void read2(T &a,T &b){read1(a);read1(b);}
template<typename T>
void read3(T &a,T &b,T &c){read1(a);read1(b);read1(c);}
template<typename T>
void out(T a)
{
    if(a>9) out(a/10);
    putchar(a%10+'0');
}
int T,kase = 1,i,j,k,n,m,x,y,a,b,d;
#define N 1<<16
int c[N];
ll ans;
inline ll gcd(ll a,ll b){return b == 0?a:gcd(b,a%b);}
void dfs1(int p)
{
    while(p >= 0) {
        if(!(a&(1<<p))){p--;continue;}
        b += (1<<p);
        dfs1(p-1);
        b -= (1<<p);
        p--;
    }
    if(b) ans += 1LL*gcd(a,b)*(1<<(c[a]-c[b]-d));
    else ans += 1LL*gcd(a,b)*((1<<(c[a]-c[b]-d))-(x==0?1:0)-(y==0?1:0));//x,y分别是各自的实际值;
}
void dfs2(int p)
{
    if(p == -1){
        if(x == 0 || y == 0) return ;
        ans += gcd(x|y,x&y);
        return ;
    }
    while(p >= 0 && !(a&(1<<p))) p--;
    // **+(1<<p+1)-1是为了说明之后任意的枚举a中的位元1度不会超出n,m的范围
    if(x+(1<<p+1)-1 <= n && y+(1<<p+1)-1 <= m){ 
        b = x&y;
　　　　 dfs1(p);
        return ;
    }
    d++;//d表示x中1的个数比y中1的个数多d个
    if(x+(1<<p) <= n) x += 1<<p,dfs2(p-1), x -= 1<<p;// x承担
    if(y+(1<<p) <= m) y += 1<<p,dfs2(p-1), y -= 1<<p;// y承担
    d--;
    if(x+(1<<p) <= n && y+(1<<p) <= m)// x,y均承担
        x += 1<<p, y += 1<<p, dfs2(p-1), x -= 1<<p,y -= 1<<p;
}
int pop(unsigned x) //使用的是分治的思想
{
    x = x - ((x>>1) & 0x55555555);
    x = (x & 0x33333333) + ((x >> 2) & 0x33333333);
    x = (x + (x >> 4)) & 0x0F0F0F0F;
    x = x + (x >> 8);
    x = x + (x >> 16);
    return x & 0x0000003F;
}
int main()
{
    rep0(i,0,N) c[i] = pop(i);//__builtin_popcount
    read1(T);
    while(T--){
        ans = 0;
        read2(n,m);
        if(n < m) swap(n,m);
        int mx = 0;
        while((1<<mx) <= n) mx++;
        rep1(i,1,(1<<mx)-1){
            if(i <= m) a = i,b = 0,x = 0,dfs1(mx);
            else a = i,b = 0,x = 0,y = 0,d = 0,dfs2(mx);
        }
        out(ans);
        puts("");
    }
    return 0;
}