据说是NP问题。
很明显我们要先预处理出来勾股数对。
但由于数过于大，所以常规的枚举是解决不了问题的。
但也貌似没有什么很好的办法可以立马找到一个数的勾股数对。
所以只能缩减枚举范围。
已知：
\[\begin{aligned} (x-y)^2+4xy&=(x+y)^2\\ (x^2-y^2)^2+4x^2y^2&=(x^2+y^2)^2\\ (x^2-y^2)^2+(2xy)^2&=(x^2+y^2)^2 \end{aligned}\]所以可以把 \(A\) 看成 \((x-y)\)，\(B\) 看成 \(2xy\)，然后枚举 \(i,j\)，代表式子中的 \(x,y\)，即可缩小枚举范围。因为要满足 \(x^2-y^2\le10^6\hspace{0.3cm}2xy\le10^6\)。枚举复杂度约为 \(O(n)\)。
然后对于每一对勾股数对，连接一条边，即可得到一个森林。
这样就可以进行树形dp了。
设 \(dp_{i,(0,1)}\) 代表选与不选数 \(i\) 的方案数。
如果不选，方案即为自己儿子所有方案的乘积
\[dp_{i,0}=\prod_{(i,j)} dp_{j,0}+dp_{j,1} \]若选，则自己儿子不能选，乘上自己儿子不选的方案数与选这个数的子集数（不算空集）即可。
\[dp_{i,1}=2^{num_i}\prod_{(i,j)}dp_{j,0} \]然后对每棵树树形dp即可。
复杂度 \(O(n)\)。
但是交上去不对，那是因为在这之前的所有假设都是基于建出来的是棵树。
那接下来就麻烦了。
但是通过楼上大佬的模拟，发现在本题的数据范围内建出来的仅仅是棵仙人掌，这样NP问题就可以转换为普通问题了，那么我们就可以进行仙人掌dp。
首先利用Tarjan思想，对于每个整棵仙人掌遍历一遍，建立时间戳，寻找环。
找到环以后，我们考虑断掉环上的一条边，并记录下来边两边的点。
然后对于每一棵仙人掌，对记录下来的点的状态进行搜索，然后在搜索终端判断两边的点是否同为 \(1\)。如果满足都不同为 \(1\)，进行一次树形dp。
那么这棵仙人掌所带来的方案数即为 \(dp_{i,0}+dp_{i,1}\)。
最后把所有仙人掌所得来的答案相乘，减去全部为空集情况，即减去一，就可得到答案。
#include<iostream>
#include<cstdio>
#include<vector>
#define int long long
using namespace std;
const int N=1e6+5;
const int mod=1e9+7;
int n,num[N],p[N],vis[N],ins[N],dfn[N],cnt,dp[N][2],chose[N],col[N],color;
vector<int>a[N],t;
int gcd(int x,int y)
{
if(y==0)return x;
return gcd(y,x%y);
}
void insert(int x)//记录断掉的边两边的点。
{
if(!vis[x])
{
vis[x]=1;
t.push_back(x);
}
}
void build(int x,int fa)//找环。
{
dfn[x]=++cnt;
int len=a[x].size();
for(int i=0;i<len;i++)
{
if(a[x][i]==fa)continue;
if(!dfn[a[x][i]])build(a[x][i],x);
else if(dfn[a[x][i]]<dfn[x])insert(a[x][i]),insert(x);
}
}
bool check(int x,int fa)//判断断掉的边的两点是否符合条件
{
col[x]=color;
int len=a[x].size();
for(int i=0;i<len;i++)
{
if(a[x][i]==fa)continue;
if(chose[x]&&chose[a[x][i]])return false;
if(col[a[x][i]]!=color&&!check(a[x][i],x))return false;
}
return true;
}
void get_ans(int x,int fa)//树形dp。
{
col[x]=color;
int len=a[x].size();
dp[x][0]=1;
dp[x][1]=p[num[x]]-1;
if(vis[x])
{
if(chose[x])dp[x][0]=0;
else dp[x][1]=0;
}
for(int i=0;i<len;i++)
{
if(a[x][i]==fa||col[a[x][i]]==color)continue;
if(col[a[x][i]]!=color)get_ans(a[x][i],x);
dp[x][1]*=(dp[a[x][i]][0]);
dp[x][1]%=mod;
dp[x][0]*=(dp[a[x][i]][0]+dp[a[x][i]][1]);
dp[x][0]%=mod;
}
}
int dfs(int x,int step)//搜索。
{
if(step==t.size())
{
color++;
if(check(x,0))
{
color++;
get_ans(x,0);
return dp[x][0]+dp[x][1];
}
return 0;
}
chose[t[step]]=0;
int sum=0;
sum+=dfs(x,step+1);
chose[t[step]]=1;
sum+=dfs(x,step+1);
return sum%mod;
}
int solve(int x)
{
t.clear();
build(x,0);
int sum=dfs(x,0);
return sum;
}
signed main()
{
scanf("%lld",&n);
for(int i=1,x;i<=n;i++)scanf("%lld",&x),num[x]++;
p[0]=1;
for(int i=1;i<=(N-5);i++)p[i]=p[i-1]*2,p[i]%=mod;
//以下为预处理。
for(int i=1;i<=(N-5)/2;i++)
{
for(int j=i+1;2*i*j<=(N-5)&&j*j-i*i<=(N-5);j++)
{
if(((i&1)^(j&1))&&gcd(i,j)==1)
{
int u=j*j-i*i,v=2*i*j;
if(!num[u]||!num[v])continue;
a[u].push_back(v);
a[v].push_back(u);
}
}
}
int ans=1;
for(int i=1;i<=N-5;i++)if(num[i]&&!col[i])ans=ans*solve(i)%mod;
printf("%lld",ans-1);
return 0;
}