BZOJ2988 : DIVISORS

令$NX=PQ$，其中$N<Pleq B$，则$1leq Q<X$，对于每个满足条件的$N$，在最大的$Q$处统计它的贡献。

枚举最大的$Q$，那么$NX$要是$Q$和$X$的倍数，且不是$Q+1,Q+2,...,X-1$的倍数。

考虑容斥原理，指数级枚举$NX$一定是$Q+1$到$X-1$这些数中哪些数的倍数，剩下的不管，则对答案的贡献为$(-1)^{打破的限制数量}lfloorfrac{BQ}{lcm} floor$。

注意到$1..60$在$6 imes 10^{13}$范围内的$lcm$数量不多(只有大约一百多万个)，将容斥中的指数级枚举换成DP：

设$f[i][j]$表示考虑了$i$到$X-1$，被打破的限制的$lcm=j$的方案数乘以容斥系数之和，然后直接转移即可。

注意精细地实现做到$O(1)$转移，比如通过质因数分解来去掉求$lcm$的$O(log X)$复杂度。

#include<cstdio>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N=65,M=219,K=8197;
int Case,X,S,i,j,k,x,y,z;ll B,lim,ans;
int id[6][4][3][3],cnt,g[M],va[M],p2[9],p3[9],p5[9],p7[9];
int p[N],tot,all;
ll f[M][K],vb[K],_vb[K];
int q[K],at[K],to[K],en[N][M],r;
inline bool isprime(int x){
  for(int i=2;i<x;i++)if(x%i==0)return 0;
  return 1;
}
inline int get(int x,int y){
  int t=0;
  while(x%y==0)x/=y,t++;
  return t;
}
inline int mask(int x){
  int t=0;
  for(int i=0;i<tot;i++)if(x%p[i]==0)t|=1<<i;
  return t;
}
inline void gen(int o){
  int t=0;
  for(int i=0;i<6;i++)for(int j=0;j<4;j++)for(int x=0;x<3;x++)for(int y=0;y<3;y++)g[t++]=id[max(i,get(o,2))][max(j,get(o,3))][max(x,get(o,5))][max(y,get(o,7))];
}
inline bool cmp(int x,int y){return vb[x]<vb[y];}
int main(){
  for(p2[0]=p3[0]=p5[0]=p7[0]=i=1;i<6;i++){
    p2[i]=p2[i-1]*2;
    p3[i]=p3[i-1]*3;
    p5[i]=p5[i-1]*5;
    p7[i]=p7[i-1]*7;
  }
  for(i=0;i<6;i++)for(j=0;j<4;j++)for(x=0;x<3;x++)for(y=0;y<3;y++){
    va[cnt]=p2[i]*p3[j]*p5[x]*p7[y];
    id[i][j][x][y]=cnt++;
  }
  scanf("%d",&Case);
  while(Case--){
    scanf("%lld%d",&B,&X);
    tot=0;
    lim=B*X;
    for(i=X;i>=11;i--)if(isprime(i))p[tot++]=i;
    for(S=0;S<1<<tot;S++){
      ll t=1;
      for(i=0;i<tot;i++)if(S>>i&1){
        t*=p[i];
        if(t>lim)break;
      }
      vb[S]=min(t,lim+1);
    }
    all=(1<<tot)-1;
    for(i=0;i<=all;i++)q[i]=i;
    sort(q,q+all+1,cmp);
    for(i=0;i<=all;i++)at[q[i]]=i,_vb[i]=vb[q[i]];
    for(j=0;j<cnt;j++)en[X+1][j]=all;
    for(i=X;~i;i--){
      lim=B*i;
      for(j=0;j<cnt;j++){
        k=en[i+1][j];
        while(~k){
          if(va[j]>lim/_vb[k])k--;
          else break;
        }
        en[i][j]=k;
      }
    }
    ans=0;
    S=mask(X);
    for(i=0;i<cnt;i++){
      r=en[X][i];
      for(j=0;j<=r;j++)f[i][j]=0;
    }
    gen(X);
    f[g[0]][at[S]]=1;
    for(i=X-1;i;i--){
      lim=B*i;
      S=mask(i);
      gen(i);
      for(k=0;k<=all;k++)to[k]=at[q[k]|S];
      for(j=cnt-1;~j;j--){
        int v=va[j],nj=g[j],r=en[i][j];
        for(k=r;~k;k--)if(f[j][k]){
          ll A=v*_vb[k],W=f[j][k];
          int nk=to[k];
          A=va[nj]*_vb[nk];
          ans+=lim/A*W;
          if(A<=lim-B)f[nj][nk]-=W;
        }
      }
    }
    printf("%lld
",ans);
  }
  return 0;
}