$O(n)$的做法：

显而易见的结论：

质因子最多的一个数一定能写成$2^x\times3^y$的形式$(y\leq 1)$。

因为如果有$>3$的质因子，那么可以用$2^2$替代。

如果$y\leq2$，那么可以用$2^3$替代。

所以在搜索的某一步：

使得$b_i$不变的数是一样的；

使得$b_i/2$的数是一样的；

使得$b_i/3$的数是一样的。

搜索的时候放在一起就可以了。

此时可以写出$O(C^n_{2n})$的代码

第一步优化：

设$dp_{i,x,y}$表示选了$i$个数，当前$gcd$是$2^x\times3^y$的方案数；

设$f_{x,y}$表示$1\to n$中有多少个数的因子包含$2^x\times3^y$。

若填完这个位置$gcd$不变：

$dp_{i+1,x,y}=dp_{i+1,x,y}+dp_{i,x,y}\times(f_{x,y}-i)$

若填完这个位置$gcd/2$：

$dp_{i+1,x-1,y}=dp_{i+1,x-1,y}+dp_{i,x,y}\times(f_{x-1,y}-f_{x,y})$

若填完这个位置$gcd/3$：

$dp_{i+1,x,y-1}=dp_{i+1,x,y-1}+dp_{i,x,y}\times(f_{x,y-1}-f_{x,y})$

使用记忆化搜索可以使得时间复杂度为$O(n\log n)$。

但是我们还可以优化：

第二步优化：

考虑优化记忆化搜索：

可以发现当前还没填的数，如果填了$gcd$不变的话，那么这个数填在之后任意一个位置都是可以的，直接用组合数计算即可，然后进入下一层。

最终复杂度：

阶乘预处理$O(n)$。

逆元预处理$O(n)$。

记忆化搜索$O(\log n)$。

总时间复杂度$O(n)$。

可以过$n\leq 10^7$的数据。

代码：

​x
#include <bits/stdc++.h>
using namespace std;
#define LL long long
#define G 11111111
#define mod 1000000007
inline int read(){
    int x=0,f=1;
    char c=getchar();
    while (c<'0' || c>'9'){
        if (c=='-') f=-1;
        c=getchar();
    }
    while (c>='0' && c<='9'){
        x=x*10+c-48;
        c=getchar();
    }
    return x*f;
}
LL f[31][2];
LL fac[G],inv[G];
int n,M;
LL dfs(int left,int x,int y){
    LL tmp[31][2]; //tmp是寄存数组，所以必须是局部变量。用来暂时保存f数组的值
    for (int i=0;i<=M;++i){ //保存f数组的值
        for (int j=0;j<=1;++j){
            tmp[i][j]=f[i][j];
        }
    }
    LL tot=0; //tot记录f数组内的和
    for (int i=x;i<=M;++i){ 
        for (int j=y;j<=1;++j){
            tot+=f[i][j]; 
            f[i][j]=0;
        }
    }
    if (!x && !y){ 
        for (int i=0;i<=M;++i){ //还原f数组的值
            for (int j=0;j<=1;++j){
                f[i][j]=tmp[i][j];
            }
        }
        return fac[tot]; //直接返回阶乘
    }
    LL ans=0;
    if (x) ans+=dfs(left-tot,x-1,y); //如果可以使gcd/=2，递归
    if (y) ans+=dfs(left-tot,x,y-1); //如果可以使gcd/=3，递归
    for (int i=0;i<=M;++i){ //还原f数组的值
        for (int j=0;j<=1;++j){
            f[i][j]=tmp[i][j];
        }
    }
    return fac[left-1]*inv[left-tot]%mod*tot%mod*ans%mod; //返回结果
}
int main(){
    n=read();
    fac[0]=1;
    for (int i=1;i<=n;++i) fac[i]=fac[i-1]*i%mod; //预处理阶乘
    inv[0]=inv[1]=1;
    for (int i=2;i<=n;++i) inv[i]=(mod-mod/i)*inv[mod%i]%mod;
    for (int i=1;i<=n;++i) inv[i]=inv[i]*inv[i-1]%mod; //预处理逆元的阶乘
    //预处理f数组，具体见上面解析
    f[0][0]=n;
    for (M=1;f[M-1][0];++M){
        f[M][0]=f[M-1][0]>>1;
    } 
    M-=2;
    for (int i=0;i<=M;++i){
        f[i][0]-=f[i+1][0];
    }
    f[0][1]=n/3;
    for (int i=1;i<M;++i){
        f[i][1]=f[i-1][1]>>1;
    }
    for (int i=0;i<M;++i){
        f[i][1]-=f[i+1][1];
    }
    for (int i=0;i<=M;++i){
        f[i][0]-=f[i][1];
    }
    //预处理结束
    LL ans=dfs(n,M,0); //不考虑3的情况
    if ((1<<(M-1))*3<=n) ans+=dfs(n,M-1,1); //考虑3的情况
    cout<<ans%mod<<endl; //输出
    return 0;
}
​