bzoj 3811: 玛里苟斯

3811: 玛里苟斯

Time Limit: 10 Sec  Memory Limit: 256 MB
Submit: 190  Solved: 95
[Submit][Status][Discuss]

Description

魔法之龙玛里苟斯最近在为加基森拍卖师的削弱而感到伤心，于是他想了一道数学题。

S 是一个可重集合，S={a1,a2,…,an}。

等概率随机取 S 的一个子集 A={ai1,…,aim}。

计算出 A 中所有元素异或 x, 求 xk 的期望。

 

 

Input

第一行两个正整数 n, k。

以下 n 行每行一个整数，表示 ai。

 

 

Output

如果结果是整数，直接输出。如果结果是小数（显然这个小数是有限的），输出精确值（末尾不加多余的 0）。

 

 

Sample Input

4 2
0
1
2
3

Sample Output

3.5

 

---

 

首先，对于题目的答案小于2⁶³，我们可以看出对于不同的k，a_i是不同的：

k=1,a_i<2⁶³
k=2,a_i<2³²
k=3,a_i<2²¹
k=4,a_i<2¹⁶
k=5,a_i<2¹³

首先，因此我们依据数据范围，对于不同的ai有不同的解法。

其次，题目要求输出精确的小数，由以下推导我们可以知道，小数最多只有.5,不可能存在.25,.125这样的：

对于求出来了其中一种情况x，他对答案的贡献是x^k*p，其中p是该情况出现的概率。
我们以x=2为例，其他k的情况类似：
我们把这个x按二进制分成了最多n位，也就是最大的数的最高位为n。那么他的k次方是(w0*2⁰+w1*2¹+w2*2²……+wn*2ⁿ)*(w0*2⁰+w1*2¹+w2*2²……+wn*2ⁿ)，展开就是w0*w0*2⁰*2⁰+w0*w1*2⁰*2¹……wn*wn*2ⁿ*2ⁿ。其中w_i是它每个二进制位的值。
对于每个项w_i*w_j*2^i+j，他存在的情况是w_i为1且w_j为1，也就是该位异或和为1。那么他的存在的概率为P/2^N(N为输入的数列数字数),P为该两位皆为1的情况数。
证明题外话：对于数列，我们可以先做个线性基，把线性相关的数字去掉，剩下的都是线性无关的数，这样方便处理。
我们可以构造一个异或（%2)的增广矩阵来求解这样的一个P的数量，x₁~x_n,即x_i代表这个数字是否在集合中，作为一个行向量X，把所有数字按按二进制位分解，每个数字占一列，构造一个矩阵MART.
那么我们求解的是X*MART=ANS，ANS为一个行向量，只有你要求的对应位为1，其余为0。
假如我们求解的是i位和j位,我们的矩阵对应方程差不多是这样的：
(2333 此处^为异或)
x₁*w(1,0)^x₂*w(2,0)^x₃*w(3,0) ……^x_N*w(N,0) =0;
x₁*w(1,1)^x₂*w(2,1)^x₃*w(3,1) ……^x_N*w(N,1) =0;
……
x₁*w(1,i)^x₂*w(2,i)^x₃*w(3,i) ……^x_N*w(N,i) =1;
……
x₁*w(1,j)^x₂*w(2,j)^x₃*w(3,j) ……^x_N*w(N,j) =1;
……
x₁*w(1,n)^x₂*w(2,n)^x₃*w(3,n) ……^x_N*w(N,n) =0;
下面还有N-n个方程，对应更高位。毕竟是个矩阵嘛行数等于列数，但它们的w全为0，等号右边也为0，对求解无影响就不列举了。

其中w(i,j)，j代表数列第i个数字的第j位。
我们经过高斯消元以后弄成上三角,可以得出有t个f[i][i]非零的行,那么只有这t个xi是有唯一解的，其他的x_i有多解。因为%2，所以其他x_i有2解，那么总共就有2^N-t个解，即P=2^N-t;
那么w_i*w_j存在的概率为p/2^N=1/(2^t);
那么每项对答案的贡献为2^i+j-t。
对于i≠j,i+j≥1,t≤2，因此最多小数点后1位.5。
对于i==j,上述方程只有一行初始为1，所以t≤1,i+j≥0，也是最多.5。
上述结论扩展成k=n的情况，无非是方程等于1的行增加为n。
如果我们求解的位(w_i*w_j*w_t.....）中有k位不同，那么方程等于1的行就缩小为k行，因此i+j+t....≥0+1+……k-1=k*(k-1)/2>=k-1=t-1,k≥2。
∴Σp-t>=-1,p为所求解对应位的位号。对应k=2的贡献2^i+j-t，K=n贡献为2^Σp-t>=2^-1。所以最多有一位小数.5。
k=1显而易见最多/2所以k=1也是最多也是.5。
至此证明完毕。

接下来我们对应数据范围，提出三种不同范围下的解法：

k=1时，就求x的期望，它是线性的。对于所有数的所有子集的异或和，我们从他们某个二进制位看。如果有数字该位为1,那么该位有奇数个1和偶数个1的概率是相等的，皆为1/2。为奇数个为该位为1的情况，为1/2。如果没有则不可能为1，该位始终为0。
那么我们只要把每个数字或一下，然后乘1/2就是答案了。

k=2时，就要按位做了。先做个线性基剔除线性无关的向量。然后把在线性基每位唯一化(即消消元弄成每位上只有一个数字有1)，这样做方便处理独立性。按照上面证明的展开，那么对于两位i,j，由k=1可得，如果这两位独立的话（即不在同一个数字中相应位为1），那么等概率出现(0,0),(0,1),(1,0),(1,1)，(1,1)即w_i，w_j存在的概率为1/4，不独立的话只可能等概率出现(0,0)和(1,1)，(1,1)概率为1/2。
那么枚举32位*32位,算算他们均为1的概率P*2^i+j,作为答案贡献加入答案中。

k≥3时，由于a_i<2²¹ ,即使枚举0~a_i也最多百万级别的时间复杂度。我们可以类似线性基求第k小枚举出所有可能出现的数字。于是我们可以考虑做完线性基后，唯一化，构造新数组存不为0的元素中。
然后枚举每个数是否加入异或贡献中，这样可以求出所有数异或的所有值了。借用之前题目的结论，每个不同的异或值有P=2^n-|μ|选择元素的方法，n为原本数组个数，|μ|为新构造的数组个数。那么每个值出现的概率就是P/总选择方法数=P/(2ⁿ)=1/(2^|μ|)的概率，值 乘 概率即为对答案的贡献。
这里需要做个两位的高精来保证数字不出现错误。

#include<bits/stdc++.h>
#define clr(x) memset(x,0,sizeof(x))
#define clr_1(x) memset(x,-1,sizeof(x))
#define LL unsigned long long
#define mod 1000000007
using namespace std;
const int N=2e5+10;
LL a[N],liner[100],per[100];
bool f[40][40];
bool bits[40];
int n,m,k,cnt;
LL ans,res,keepbit;
void calc(int bit)
{
    clr(liner);
    LL p;
    for(int i=1;i<=n;i++)
    {
        for(int j=bit;j>=0;j--)
        {
            p=a[i];
            if(p>>j)
            {
                if(liner[j]) p^=liner[j];
                else
                {
                    liner[j]=p;
                    break;
                }
            }
        }
    }
    for(int i=bit;i>=0;i--)
    {
        for(int j=i-1;j>=0;j--)
            if((liner[i]>>j)&1) liner[j]^=liner[i];
    }
    cnt=0;
    for(int i=bit;i>=0;i--)
        if(liner[i])
            per[++cnt]=liner[i];
    return ;
}
void solve1(int n)
{
    ans=0,res=0;
    for(int i=1;i<=n;i++)
        ans|=a[i];
    if(ans&1)
        res=1;
    ans>>=1;
    return ;
}
void solve2(int n)
{
    clr(bits);
    ans=res=0;
    bool zero;
    for(int i=32;i>=0;i--)
    {
        zero=0;
        for(int j=1;j<=cnt;j++)
        {
            f[i][j]=(per[j]>>i)&1;
            zero|=f[i][j];
        }
        bits[i]=zero;
    }
    int p;
    for(int i=32;i>=0;i--)
        for(int j=32;j>=0;j--)
            if(bits[i]>0 && bits[j]>0)
            {
                p=1;
                for(int l=1;l<=cnt;l++)
                    if(f[i][l]^f[j][l])
                    {
                        p++;
                        break;
                    }
                if(i+j-p>=0)
                    ans+=(1LL<<(i+j-p));
                else
                {
                    res++;
                }
                ans+=res>>1;
                res&=1;
            }
    return ;
}

void dfs(int pos,LL num)
{

    if(pos>cnt)
    {
        //两位高精度计算u高位，v低位，因为这东西太大了。k>=3的情况下ans和res也算是一个高精度的高位和低位关系。
        LL u=0,v=1;
        for(int i=1;i<=k;i++)
        {
            u*=num;
            v*=num;
            u+=v>>cnt;
            v&=keepbit;
        }
        ans+=u;
        res+=v;
        ans+=res>>cnt;
        res&=keepbit;
        return ;
    }
    dfs(pos+1,num^per[pos]);
    dfs(pos+1,num);
    return ;
}
void solve3(int n,int k)
{
    ans=0;
    res=0;
    keepbit=(1LL<<cnt)-1;
    dfs(1,0);
    return ;
}
int main()
{
    scanf("%d%d",&n,&k);
    for(int i=1;i<=n;i++)
        scanf("%llu",&a[i]);
    if(k==1)
        solve1(n);
    else if(k==2)
    {
        calc(32);
        solve2(n);
    }
    else
    {
        calc(61/k+1);
        solve3(n,k);
    }
    printf("%llu",ans);
    if(res) printf(".5");
    printf("\n");
    return 0;
}