高斯消元法——HDU 5833

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

• 分析：
  给出n个数，从中取一个数或者取多个数相乘得到完全平方数，问有多少种取法。

• 题解：
  因为每个数的大小为$[0, 10^{18}]$，每个数的质因子大小不超过2000，又2000以内的质数只有303个，所以我们可以先预处理一遍：把每个数质因子分解存进数组，若含有同一个质因子的个数为偶数则存为0，为奇数则存为1。如果某个数质因子分解得到的数组里存的全是0，那么这个数本身就是完全平方数；如果多个数的质因子数组互相异或得到的新数组里全是0，那么这几个数相乘能得到一个完全平方数。所以我们可以把题意转换成求一个线性方程组的解：
  $a_{11}x_1+a_{12}x_2+...+a_{1n}x_n=0$
  $a_{21}x_1+a_{22}x_2+...+a_{2n}x_n=0$
  …
  $a_{n1}x_1+a_{n2}x_2+...+a_{nn}x_n=0$
  其中$a_{ij}$表示第i个质数在第j个数里含有的个数是奇数(为1)还是偶数(为0)，$X_i$表示第i个数有没有被选。
  求出这个方程系数矩阵的秩rank ，则基础解系中含有n-rank个解向量，那么总共有$2^{n-rank}$种取法(每一个解向量都有取和不取两种情况，任意多个解向量组合都是该方程的解)，再减去全部不取，即全为0的情况，则为题目的答案。

• 参考代码：

#include <cstdio>
#include <cstring>
#include <algorithm>
#define ll long long
#define mod 1000000007
using namespace std;
const int N=2000;
const int M=310;
int prime[N+1],cnt;
int n,t,mat[M][M];
ll a[M];
void getPrime()      //预处理找出前几个质数
{
    for(int i=2;i<=N;i++)
    {
        if(!prime[i])prime[++cnt]=i;    //第cnt个质数
        for(int j=1;j<=cnt&&prime[j]<=N/i;j++)
        {
            prime[prime[j]*i]=1;        //用当前的数乘以之前找到的所有质数来排除后面的合数。
            if(i%prime[j]==0)break;     //如果当前的数与之前的质数不互质则跳出
        }
    }
}
int Rank(int c[][M])
{                            //异或版的高斯消元求秩
    int i=0,j=0,k,r,u;
    while(i<=cnt&&j<=n)
    {
        r=i;
        while(c[r][j]==0&&r<=cnt)r++;
        if(c[r][j])
        {
            swap(c[i],c[r]);
            for(u=i+1;u<=cnt;u++)
            {    
                if(c[u][j])
                {    
                    for(k=i;k<=n;k++)
                        c[u][k]^=c[i][k];
                }
            }
            i++;
        }   
        j++;
    }
    return i;
}
int solve()
{
    memset(mat,0,sizeof mat);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=cnt;j++)
        {
            ll tmp=a[i];
            while(tmp%prime[j]==0)
            {
                tmp/=prime[j];
                mat[j][i]^=1;
            }
        }
    int b=n-Rank(mat);  //b个解向量
    ll ans=1;
    ll k=2;
    while(b)           //快速幂
    {
        if(b&1)
        {
            ans=ans*k%mod;
        }
        k=k*k%mod;
        b>>=1;
    }
    return ans-1;//减去全为0的解
}
int main() 
{
    getPrime();
    scanf("%d",&t);
    for(int cas=1;cas<=t;cas++)
    {
        scanf("%d",&n);
        for(int i=1;i<=n;i++)
            scanf("%lld",&a[i]);
        printf("Case #%d:\n%d\n",cas,solve());
    }
    return 0;
}