2016CCPC网预 HDU 5833 Zhu and 772002

2016CCPC网预02

HDU 5833 Zhu and 772002

高斯消元，数论，异或方程组

传送门：HDU

---

题意

给你一些数，从中选出至少1个，每个数只能用一次，使得他们的乘积是完全平方数。问选法数。

---

思路

训练指南原题改编，P160。下为书上讲解：

“不含大于2000的素因子”提示我们考虑每个数的唯一分解式，用01向量表示—个数， 再用n个01变量xi来表示我们的选择，其中xi=1表示要选第i个数，xi=0表示不选它，则可对每个素数的幂列出一个模2的方程。
这话听起来比较抽象，让我们分析一下题目中的例子。4个整数4,6,10,15的素因子只 有2, 3, 5这3种，首先把这些整数写成01向量的的形式，即4=22∗30∗50,即(2,0,0)；6=21∗31∗50,即(1,1,0)；10=21∗30∗51,即(1,0,1)；15=20∗31∗51,即(0,1,1)。
选出来的数乘积为22x1+x2+x3∗3x2+x4∗5x3+x4。如果要让这个数是完全平方数，每个幂都应该是偶数，即

⎧⎩⎨⎪⎪x2+x3≡0(mod2)x2+x4≡0(mod2)x3+x3≡0(mod2)

注意，这也是一个线性方程组，只是代数系统变成了Z2（模2的剩余系）。可是第一 个方程里的2x1不见了。这是因为2x1总是偶数，所以没必要写在方程里。同理，3x1会变成x1，任意变量的系数非0即1。还可以把这个方程组看成是如下xor方程组

⎧⎩⎨⎪⎪x2xorx3≡0x2xorx4≡0x3xorx3≡0

需要求解方程组解的组数。可以求秩，即求到自由变元的个数。xor方程组是很好消元的，因为不需要做乘法和除法，只需要做xor;每次也不需要找 绝对值最大的系数(每个系数不是0就是1)，任意一个系数为1即可实现消元。
最后，假设自由变量有f个，则线性方程组的解共有2f个，因为每个自由变量可以取0和1。比如，刚才的方程组对应的增广矩阵消元后为

⎡⎣⎢0001101010|01|01|0⎤⎦⎥⇒⎡⎣⎢0001001100|01|00|0⎤⎦⎥

有两个自由变量x1和x4，有两个有界变量x2和x3，因此—共有22=4种选法。注意，本题不允许一个整数都不选，因此最终答案需要减1。

换句话说，我们需要做的就是对给的数做素数分解，因为素数大小不超过2000，2000以内素数有303个，所以开个二维数组，记录每个题给数的素数因子个数，个数偶数记为0，奇数记为1(素数分解时不断异或1就行了)。得到二维数组实际就是系数矩阵。求矩阵的秩就可以。高斯消元的模板。

代码

#include <iostream>#include <cstdio>#include <algorithm>#include <cstdlib>#include <cstring>#include <cmath>#include <vector>#include <queue>#include <stack>#include <iomanip>#include <string>using namespace std;const int MAXN=305;const int oo=1000000007;typedef int Matrix[2107][2107];typedef long long LL;Matrix A;int ra_nk(int m,int n)//A是异或方程组系数矩阵 返回秩{    int i=0,j=0,k,r,u;    while(i<m&&j<n)    {        r=i;        for(int k=i;k<m;k++)        {            if(A[k][j]) { r=k;break; }        }        if(A[r][j])        {            if(r!=i) for(k=0;k<=n;k++) swap(A[r][k],A[i][k]);            //消元后第i行的第一个非0列是第j列，且第u>i行的第j列均为0            for(u=i+1;u<m;u++) if(A[u][j])                for(k=i;k<=n;k++) A[u][k]^=A[i][k];            i++;        }        j++;    }    return i;}int vis[2107];int prime[2107];int gen_primes(int m){    memset(vis,0,sizeof(vis));    int cnt=0;    for(int i=2; i < m; i++)    {        if(!vis[i])        {            prime[cnt++]=i;            for(int j=i*i; j<m; j+=i)                vis[j]=1;        }    }    return cnt;}LL powmod(LL a,LL n){    LL mod=1e9+7;    LL res=1;    while(n){        if(n&1) res=res*a%mod;        a=a*a%mod;        n>>=1;    }    return res;}int main(){    int prime_n=gen_primes(2107);    int T;    scanf("%d",&T);    for(int t=1;t<=T;t++)    {        memset(A,0,sizeof(A));        int n;        scanf("%d",&n);        int maxp=0;        for(int i=0;i<n;i++)        {            LL x;            scanf("%lld",&x);            for(int j=0;j<prime_n;j++)            {                while(x%prime[j]==0)                {                    maxp=max(maxp,j);                    x/=prime[j];                    A[j][i]^=1;                }            }        }        int r=ra_nk(maxp+1,n);        printf("Case #%d:\n",t);        printf("%lld\n",powmod(2,(LL)n-r)-1);    }}