考虑将点分为 $A,B$ 两类。其中 $[x\in A]\iff\exists_{y\neq x},y|x$。
那么我们删去的点只可能在 $B$ 类中，且当前 $x\in B$ 可删当且仅当存在 $y\in B,z\in A$ 使得 $z|x\land z|y$。
那么对于 $z\in A,x\in B,z|x$，连边 $(z,x)$。这样得到一个二分图，每次删除操作相当于选择某个 $z\in A$ 且 $z$ 的度数 $\geq 2$，然后删去某个和 $z$ 相邻的点。
对二分图中每个连通块单独考虑，假设该连通块 $B$ 类有 $m$ 个点，那么发现答案能达到上界 $m-1$：以某个 $B$ 类点开始跑一棵生成树，然后从叶子往根删即可。
现在考虑统计方案数。只需统计单个连通块内的方案数，然后连通块间乘个组合数即可。
发现连通块内能删到上界的等价条件是：每次删完某个 $B$ 类点后，剩下的 $B$ 类点仍然是连通的。进一步地，若我们枚举最后剩下的是哪个 $B$ 类点 $lst$，我们只需保证当前所有未删的 $B$ 类点都与 $lst$ 连通即可。
数据范围提示我们可能是指数级做法。那么先挖掘到一条可能有用的性质：$A$ 类点至多 $20$ 个，因为 $>20$ 的 $A$ 类点至多有 $1$ 条邻边，那么它是没用的。
进一步地，由于 $x,2x$ 不可能同时成为 $A$ 类点，所以现在有用的 $A$ 类点至多 $10$ 个。同时这也是紧界：考虑 $A$ 类点为 $11,12,\cdots,20$。
把删点变成加点，虽然是等价的，但是会直观很多，且助于思考。
考虑 DP。发现我们可以只关注 $A$ 类点的添加过程（具体来说，是 $A$ 类点在 $B$ 类点添加序列的什么位置变化了，且变化是什么），而通过组合数算出 $B$ 类点的添加方案数。可以做到时间复杂度 $O(n2^{n/6})$。
#include<bits/stdc++.h>
#define N 65
using namespace std;
namespace modular
{
const int mod=1000000007;
inline int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
inline int dec(int x,int y){return x-y<0?x-y+mod:x-y;}
inline int mul(int x,int y){return 1ll*x*y%mod;}
inline void Add(int &x,int y){x=x+y>=mod?x+y-mod:x+y;}
inline void Dec(int &x,int y){x=x-y<0?x-y+mod:x-y;}
inline void Mul(int &x,int y){x=1ll*x*y%mod;}
inline int poww(int a,int b){int ans=1;for(;b;Mul(a,a),b>>=1)if(b&1)Mul(ans,a);return ans;}
}using namespace modular;
int n,a[N],fa[N];
int fac[N],ifac[N];
bool AB[N],vis[N];
int find(int x)
{
return x==fa[x]?x:(fa[x]=find(fa[x]));
}
pair<int,int> solve(int rt)
{
static int f[1050],w[1050],from[N];
vector<int> A,B;
for(int i=1;i<=n;i++)
if(find(i)==rt) (AB[i]?B:A).push_back(a[i]),vis[i]=1;
while(!A.empty()&&A.back()>20) A.pop_back();
int maxn=(1<<A.size())-1;
memset(f,0,sizeof(int)*(maxn+1));
for(int i=0;i<B.size();i++)
{
from[i]=0;
for(int j=0;j<A.size();j++)
if(!(B[i]%A[j])) from[i]|=(1<<j);
}
for(int S=0;S<=maxn;S++)
{
w[S]=0;
for(int j=0;j<B.size();j++)
if((from[j]&S)==from[j]) w[S]++;
}
for(int i=0;i<B.size();i++)
Add(f[from[i]],mul(fac[B.size()-1],ifac[B.size()-w[from[i]]]));//C(B-1,w[from[i]]-1)*fac[w[from[i]-1]]
for(int S=0;S<maxn;S++)
{
if(!f[S]) continue;
for(int j=0;j<B.size();j++)
{
if((S&from[j])&&((S&from[j])!=from[j]))
{
int T=S|from[j];
Add(f[T],mul(mul(fac[B.size()-w[S]-1],ifac[B.size()-w[T]]),f[S]));
//C(x,y)*fac[y]
}
}
}
return make_pair(B.size()-1,f[maxn]);
}
int main()
{
cin>>n;
for(int i=1;i<=n;i++) cin>>a[i];
sort(a+1,a+n+1);
for(int i=1;i<=n;i++)
for(int j=1;j<i;j++)
if(!(a[i]%a[j])){AB[i]=1;continue;}
for(int i=1;i<=n;i++) fa[i]=i;
for(int i=1;i<=n;i++)
if(!AB[i]) for(int j=i+1;j<=n;j++)
if(!(a[j]%a[i])) fa[find(i)]=find(j);
fac[0]=1;
for(int i=1;i<=n;i++) fac[i]=mul(fac[i-1],i);
ifac[n]=poww(fac[n],mod-2);
for(int i=n;i>=1;i--) ifac[i-1]=mul(ifac[i],i);
int ans=1,s=0;
for(int i=1;i<=n;i++)
{
if(AB[i]&&find(i)==i&&!vis[i])
{
auto pr=solve(i);
Mul(ans,mul(pr.second,ifac[pr.first]));
s+=pr.first;
}
}
Mul(ans,fac[s]);
printf("%d\n",ans);
return 0;
}