1221 矩阵中不重复的元素V3

一个m*n的矩阵。

该矩阵的第一列是a^b,(a+1)^b,.....(a + n - 1)^b

第二列是a^(b+1),(a+1)^(b+1),.....(a + n - 1)^(b+1)

.......

第m列是a^(b + m - 1),(a+1)^(b + m - 1),.....(a + n - 1)^(b + m - 1)

(a^b表示a的b次方）

下面是一个4*4的矩阵：

2^2=4, 2^3=8, 2^4=16, 2^5=32

3^2=9, 3^3=27, 3^4=81, 3^5=243

4^2=16, 4^3=64, 4^4=256, 4^5=1024

5^2=25, 5^3=125, 5^4=625, 5^5=3125

问这个矩阵里有多少不重复的数（比如4^3 = 8^2，这样的话就有重复了)

2^2=4, 2^3=8, 2^4=16, 2^5=32

3^2=9, 3^3=27, 3^4=81, 3^5=243

4^2=16, 4^3=64, 4^4=256, 4^5=1024

m = 4, n = 3, a = 2, b = 2。其中2^4与4^2是重复的元素。

由于数据会很大，输出Mod 1000000007(10^9 + 7)的结果。

 收起

输入

输入数据包括4个数：m,n,a,b。m，n为矩阵的长和宽（2 <= m,n <= 5 * 10^15)。a，b为矩阵的第1个元素，a^b（2 <= a , b <= 10^15）。

输出

输出不重复元素的数量Mod 1000000007(10^9 + 7)的结果。

输入样例

4 3 2 2

输出样例

11

题解：

先考虑 xa=ybxa=yb 意味着什么

记 p(x)p(x) 为最小的 pp, 使得 ∃k,pk=x∃k,pk=x。 那么 xa=p(x)k(x)×axa=p(x)k(x)×a

则 xa=ybxa=yb 等价于 p(x)k(x)×a=p(y)k(y)×bp(x)k(x)×a=p(y)k(y)×b 等价于 p(x)=p(y)∧k(x)×a=k(y)×bp(x)=p(y)∧k(x)×a=k(y)×b

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

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

于是问题转化为，求 card{i×j|i∈[l1,r1],j∈[l2,r2]}card{i×j|i∈[l1,r1],j∈[l2,r2]}其中l1l1,r1r1为O(log2(n))O(log2(n))级别，l2l2,r2r2为 O(n)O(n) 级别

容斥。

先加上在 [l1,r1][l1,r1] 中任选一个数,并在 [l2,r2][l2,r2] 任选一个数的方案数

再减去在 [l1,r1][l1,r1] 任选两个数 x,yx,y，并在 [l2,r2][l2,r2] 任选两个数 a,ba,b ,使得x×a=y×bx×a=y×b的方案数

再加上在 [l1,r1][l1,r1] 任选三个数 x,y,zx,y,z，并在 [l2,r2][l2,r2] 任选三个数 a,b,ca,b,c ,使得x×a=y×b=z×cx×a=y×b=z×c的方案数 ......

可以发现 当 [l1,r1][l1,r1] 中的数确定之后, [l2,r2][l2,r2] 的方案数是可以算出来的 因为乘积一定为 lcmlcm 的倍数

设为 xx 倍 则l2≤x×lcm/max&&x×lcm/min≤r2l2≤x×lcm/max&&x×lcm/min≤r2 于是 l2×max≤x×lcm≤r2×minl2×max≤x×lcm≤r2×min

所以 xx 的取值数=r2×min/lcm−l2×max/lcm+1r2×min/lcm−l2×max/lcm+1

于是暴力容斥复杂度 =O(2r1−l1+1)=O(n)=O(2r1−l1+1)=O(n)

l1l1,r1r1一共有 log22(n)log22(n)种

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

当然还过不了

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

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

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

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

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

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

我是筛出n−−√n以内所有p(x)=xp(x)=x的xx 然后枚举 求更优的方法。

#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;
}