题目链接

先考虑 $x^a=y^b$ 意味着什么

记 $p(x)$ 为最小的 $p$, 使得 $\exists k , p^k=x$。 那么 $x^a=p(x)^{k(x) \times a}$

则 $ x^a=y^b $ 等价于 $ p(x)^{k(x) \times a}=p(y)^{k(y) \times b}$ 等价于 $ p(x) = p(y) \land k(x) \times a = k(y) \times b $

考虑枚举 $p$ ，计算所有 $p(x)=p$ 的 $x$ 带来的重复

那么 使 $x$ 在给定区间 的 $k(x)$ 在一个区间中

于是问题转化为，求 $ card \{ i \times j \, | \, i \in [l_1,r_1] \, , j \in [l_2,r_2] \} $其中$l_1$,$r_1$为$O(log_2(n))$级别，$l_2$,$r_2$为 $O(n)$ 级别

容斥。

先加上在 $[l_1,r_1]$ 中任选一个数,并在 $[l_2,r_2]$ 任选一个数的方案数

再减去在 $[l_1,r_1]$ 任选两个数 $x,y$，并在 $[l_2,r_2]$ 任选两个数 $a,b$ ,使得$x \times a=y \times b$的方案数

再加上在 $[l_1,r_1]$ 任选三个数 $x,y,z$，并在 $[l_2,r_2]$ 任选三个数 $a,b,c$ ,使得$x \times a=y \times b=z \times c$的方案数 ......

可以发现 当 $[l_1,r_1]$ 中的数确定之后, $[l_2,r_2]$ 的方案数是可以算出来的 因为乘积一定为 $lcm$ 的倍数

设为 $x$ 倍 则$l_2 \le x \times lcm/max \, \& \& \, x \times lcm/min \le r_2$ 于是 $l_2 \times max \le x \times lcm \le r_2 \times min$

所以 $x$ 的取值数=$r_2 \times min/lcm - l_2 \times max/lcm +1$

于是暴力容斥复杂度 $=O(2^{r_1-l_1+1})=O(n)$

$l_1$,$r_1$一共有 $ log_2^2(n) $种

对于所有的一样的$l_1$,我们可以从小到大枚举$r_1$计算。 总复杂度还是 $O(n)$

当然还过不了

注意到 如果一个数选了之后对于$lcm,max,min$没有影响 那么他选和不选的贡献就会抵消

所以先枚举$max,min$ 之后$dfs$ 如果碰到对$lcm$无贡献的数，就直接退出

同时，如果选了它导致前面的数成为对$lcm$无贡献的数，那么它就不选

这样优化之后复杂度就很好了 李陶冶说"比较接近于光滑数的数量（其实比这个还要低一个数量级）"

但我不知道光滑数是什么东西 总之O(跑的过)

还有个问题，就是如何求出有多少个$p(x)$会使得$k(x)$在 $[l_1,r_1]$ 中

我是筛出$ \sqrt{n} $以内所有$p(x)=x$的$x$ 然后枚举 求更优的方法？

#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
const int D=1e9+7,Lo=52;
int cnt[Lo+2][Lo+2];
ll L1,R1,l2,r2,ans,now;
void init()
{
    ll m,n,a,b;
    cin>>m>>n>>a>>b;
    ans=(m%D)*(n%D)%D;
    L1=a;R1=a+n-1;l2=b;r2=b+m-1;
}

int Gcd[Lo+2][Lo+2];
int gcd(int x,int y)
{
    while(y)swap(x%=y,y);
    return x;
}

int l1,r1,mn;
ll L,R;//L(=l2*mx)<=x*Lcm<=R(=r2*mn)   
int mx_p[Lo+5],st[Lo+5],*top;
#define len (R/Lcm - (L+Lcm-1)/Lcm +1)
void dfs(int *h,const ll &Lcm,bool ji)
{
    if(Lcm>R)return ;
    if(h>top)
    {
        if(ji)now-=len;
        else now+=len;
        return ;
    }
    int x=*h;
    if(Lcm%x==0) return ;
    //if(len<=0) return ;
    dfs(h+1,Lcm,ji);
    dfs(h+1,Lcm*(x/Gcd[x][Lcm%x]),ji^1);
}

const int U=8e7;
bool mark[U];

int main()
{
    //freopen("1.in","r",stdin);freopen("1.out","w",stdout);
    init();
    int mp=sqrt(R1);
    for(int p=2;p<=mp;++p)
    if(!mark[p])
    {
        int k=0;
        ll x=1;
        while(x<L1&&x<=R1/p) { x*=p;if(x<=mp)mark[x]=1;++k; }
        if(x<L1||x>R1)continue;
        int l1=k; 
        while(x<=R1/p) { x*=p;if(x<=mp)mark[x]=1;++k; }
        if(k>l1)
        ++cnt[l1][k];
    }

    for(int i=1;i<=Lo;++i)
    for(int j=0;j<i;++j)
     Gcd[i][j]=gcd(i,j);
    for(int i=1;i<=Lo;++i)
    for(int j=1;j<i;++j)
    if(i%j==0) mx_p[i]=j; 

    for(l1=1;l1<=Lo;++l1)
    {
        int mr=0;
        for(r1=Lo;r1>l1;--r1)
        if(cnt[l1][r1]){mr=r1;break;}

        now=0;
        for(r1=l1+1;r1<=mr;++r1)
        {        
            L=l2*r1;
            for(mn=l1;mn<r1;++mn)
            {
               R=r2*mn;    
               if(L>R)continue;

               int Lcm=r1*(mn/gcd(r1,mn));
               int i=mn+1;
               for(;i<r1;++i)
               if(Lcm%i==0) break;
               if(i<r1)continue;

               top=st-1;
               for(int i=mn+1;i<r1;++i)
               if(mx_p[i]<=mn) 
                *++top=i;

               dfs(st,Lcm,0); 
            }
            now%=D;
            ans-=cnt[l1][r1]*now%D;
        }
    }
    cout<<(ans%D+D)%D;
}