习题之高斯消元 + 分解质因数 + 快速幂

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

题意：n个数中任选m（m >= 1）个数相乘，使得乘积是一个平方数的方案数

解析：首先可以对每个数分解质因数，显然，想要乘积是一个平方数，那么所有的质因子的个数都应该是偶数。（注：此题不能用dp，因为如果是这样的一组数 2 ，3， 6 那么当dp到六的时候，是需要乘以前两个的，最后还是要枚举全部的方案）

如何保证质因数是偶数呢？

矩阵a[j][i]，表示第i个数中第j个数出现的次数，如果是奇数就为1，如果是偶数就为0

这样就有了方程组

a11 * k1 + a12 * k2 + a13 * k3 …… + a1n * kn = 0

a21 * k1 + a22 * k2 + a23 * k3 …… + a2n * kn = 0

……

am1 * k1 + am2 * k2 + am3 * k3 …… + amn * kn = 0

得到行列式a[m]
只要解出行列式的秩，其次n -行列式的秩 == 自由元的个数。

自由元的个数就是有些点取不取都行，那么方案数就是2^(n - r) - 1,减去一定会存在的那个全零解

关于最后一列，对于其次方程有或者没有都是可行的。

求方案数是个大的幂，要用快速幂来求

代码：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <sstream>
#include <string>
#include <algorithm>
#include <list>
#include <map>
#include <vector>
#include <queue>
#include <stack>
#include <cmath>
#include <cstdlib>
using namespace std;
typedef long long LL;
int const maxn = 2005;
int const mod = 1000000007;
bool isPrime[maxn];
LL primeList[maxn],primeCount = 0;
void primeInit(LL n)
{

memset(isPrime,true,sizeof(isPrime));//初始化认为全部是素数
int m = sqrt(n + 0.5);
for(int i = 2; i <= m;  i ++)
{
if(isPrime[i])//判断是素数
{
for(int j = i * i; j <= n; j += i)
{

isPrime[j] = false;
}
}
}
for(int i = 2; i <= n; i ++)
{
if(isPrime[i])
{

primeList[primeCount] = i;
primeCount ++;
}

}
}
long long get_inv(LL x,LL m)
{
LL r,y;
for(r = 1,y = m; y; x = x * x % mod,y >>= 1)
{
if(y % 2 == 1)
r = r *  x % mod;
}
return r;
}
int A[maxn][maxn];
int rank2(int n,int m)
{
int i = 0,j = 0,r;
while(i < n && j <m )
{
for(r = i; r < n; r ++)
{
if(A[r][j])
{
break;
}
}
if(A[r][j])
{
if(i != r)
{
for(int k = 0; k < m; k ++)
{
swap(A[i][k],A[r][k]);
}

}
for(r = i+ 1; r < n; r  ++)
{
if(A[r][j])
{
for(int k = 0; k < m; k ++)
{
A[r][k] ^=A[i][k];
}
}
}
i ++;//只有当A[r][j] != 0 时，才有i++。才用到了当前行
}
j ++;
}
return i;
}
int main()
{
// freopen("in.txt","r",stdin);
int T,n;
primeInit(2005);
scanf("%d",&T);
for(int t = 1; t <= T; t ++)
{
memset(A,0,sizeof(A));
int maxp = 0;
long long x;
scanf("%d",&n);
for(int i = 0; i < n; i ++)
{
scanf("%I64d",&x);
for(int j = 0; j < primeCount ; j ++)
{

while(x % primeList[j] == 0)
{
maxp = max(maxp,j);
x = x / primeList[j] ;
A[j][i] ^= 1;
}
}
}
int r = rank2(maxp + 1,n + 1);
printf("Case #%d:\n",t);
printf("%I64d\n",get_inv(2,n - r) - 1);

}
return 0;
}