HDU 6125 Free from square

从1到n中取出最多k个数，它们的乘积不被某个除了1的完全平方数整除，求方案数

很显然每个质数最多只能出现一次。

大于sqrt（n）的质数在每个选取的数中只能取一次，小于sqrt(n)的质数最多有8个，就可以状压了。

ps:在外面除以排列数做法，我的好像只能用c++过

#include<stdio.h>
#include<string.h>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;
const int MOD=(int)1e9+7;
int p[555],vis[555],pnum;
ll dp[96][96][300];
int bg[96][300];
int cj[300],lg[300];
void init(){
    for(int i=2;i<=500;i++){
        if(!vis[i]){
            p[pnum++]=i;
            for(int j=i+i;j<=500;j+=i)vis[j]=1;
        }
    }
    for(int i=0;i<8;i++)lg[1<<i]=i;
    cj[0]=1;
    for(int i=1;i<(1<<8);i++){
        int lowbit=i&(-i);
        cj[i]=min(501,cj[i-lowbit]*p[lg[lowbit]]);
    }
}
int a[10],b[10];
int main()
{
    init();
    int t;
    scanf("%d",&t);
    while(t--){
        memset(dp,0,sizeof(dp));
        memset(bg,0,sizeof(bg));
		int n,k;
        scanf("%d%d",&n,&k);
        k=min(k,95);
        int sa,sb;
        sa=sb=0;
        for(int i=0;i<pnum;i++){
            if(p[i]*p[i]<=n)a[sa++]=p[i];
            else if(p[i]<=n)b[sb++]=p[i];
        	else break;
		}
        for(int i=0;i<sb;i++){
            for(int j=0;j<(1<<sa);j++){   
                if(b[i]*cj[j]<=n)bg[i][j]=1;
            }
        }
        dp[0][0][0]=1;
        for(int i=0;i<sa;i++){
            for(int j=1;j<(1<<sa);j++){   
            	if(cj[j]>n)continue;
                for(int l=0;l<j;l++){  //这里不是l<(1<<sa)是防止重复，也可以写成l<(1<<sa)，在外面对选取x个物品的方案除以x!也就是排列数
                    if((j&l)==0&&cj[j]<=n){
                        dp[i+1][0][l|j]=(dp[i+1][0][l|j]+dp[i][0][l])%MOD;;
                    }
                }
            }
        }
        for(int i=1;i<=sb;i++){
            for(int j=0;j<=k;j++)for(int l=0;l<(1<<sa);l++)dp[j][i][l]=dp[j][i-1][l];    
            for(int j=0;j<k;j++){
                for(int l=0;l<(1<<sa);l++){
                    if(bg[i-1][l]==0)continue;
                    for(int x=0;x<=((1<<sa)-1);x++){
                        if(x&l)continue;
                        dp[j+1][i][x|l]=(dp[j+1][i][x|l]+dp[j][i-1][x])%MOD;
                    }
                }
            }
        }
        ll ans=0;
        for(int x=1;x<=k;x++)for(int i=0;i<(1<<sa);i++)ans=(ans+dp[x][sb][i]+dp[x-1][sb][i])%MOD;
        printf("%lld\n",ans);
    }
    return 0;
}