51nod 1221矩阵中不重复的元素V3

普通代码:

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

效率高代码:

#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=(a);i<=(b);i++)
#define per(i,a,b) for(int i=(a);i>=(b);i--)
#define REP(i,n) for(int i=0;i<(n);i++)
#define fi first
#define se second
#define pb push_back
#define mp make_pair
#define SZ(x) ((int)(x).size())
#define ALL(x) (x).begin(),(x).end()
using namespace std;
typedef pair<int,int> pii;
typedef vector<int> vi;
typedef long double ld;
typedef long long ll;
typedef ll LL;
typedef pair<LL,int> pli;

namespace{
	const ll inf=3e18;
	const int mod=1e9+7;
	inline int mul(const int x,const int y){
		return (ll)x*y%mod;
	}
	inline void add(int &x,const int y){
		x=(x+y>=mod?x+y-mod:x+y);
	}
	inline void sub(int &x,const int y){
		x=(x>=y?x-y:x+mod-y);
	}
	inline ll Pow(ll x,ll y){
		ll res=1;
		for(;y&&res<=inf;y--) res*=x;
		return res;
	}
	int GCD[64][64];
	void init(){
		REP(i,64) rep(j,1,63){
			GCD[i][j]=__gcd(i,j);
		}
	}
	inline LL lcm(const LL x,const int y){
		return x/GCD[x%y][y]*y;
	}
}

const int SZ=1<<21;
namespace Hash{
	int val[SZ],head[SZ],link[SZ],stk[SZ],top,sz;
	LL key[SZ];
	inline void add(const LL a,const int b){
		int t=a&(SZ-1),k=head[t];
		for(;k;k=link[k]){
			if(key[k]==a){
				::add(val[k],b);
				return;
			}
		}
		val[++sz]=b,key[sz]=a;
		link[sz]=head[t],head[t]=sz;
		stk[++top]=t;
	}
	void solve(pli res[],int &cnt){
		cnt=sz=0;
		for(;top;top--){
			for(int k=head[stk[top]];k;k=link[k]){
				res[++cnt]=mp(key[k],val[k]);
			}
			head[stk[top]]=0;
		}
	}
}

const int N=64,M=5e5;
int wl[N],wr[N],f[N][N],g[N][N],isp[N],K,ans;
ll vl[N],vr[N],a,b,n,m;

void part1(){
	vl[0]=vr[0]=inf;
	for(K=1;vr[K-1]>=2;K++){
		vl[K]=pow(a,(ld)1/K);
		for(;Pow(vl[K]+1,K)<=a;vl[K]++);
		for(;Pow(vl[K],K)>a;vl[K]--);
		vr[K]=pow(a+n,(ld)1/K);
		for(;Pow(vr[K]+1,K)<=a+n;vr[K]++);
		for(;Pow(vr[K],K)>a+n;vr[K]--);
	}
	K--;
	per(i,K,1){
		wl[i]=(vl[i]-1)%mod;
		wr[i]=(vr[i]-1)%mod;
		for(int j=i+i;j<=K;j+=i){
			sub(wl[i],wl[j]);
			sub(wr[i],wr[j]);
		}
	}
	per(i,K,1) per(j,K,1){
		if(vl[i-1]<vr[j]) f[i][j]=wl[i-1];
		else f[i][j]=wr[j];
		rep(x,i,K) rep(y,j+(x==i),K){
			sub(f[i][j],f[x][y]);
		}
	}
	rep(i,1,K) rep(j,i,K){
		rep(x,1,i) rep(y,j,K){
			add(g[i][j],f[x][y]);
		}
	}
	rep(i,K/2+1,K){
		isp[i]=1;
		rep(j,2,i-1){
			if(i%j==0) isp[i]=0;
		}
	}
	init();
}

void rebuild(vi &p,int prod,int tax[]){
	fill(tax,tax+prod,0);
	for(auto x:p){
		for(int i=0;i<prod;i+=x){
			tax[i]=1;
		}
	}
	REP(i,prod){
		tax[i]=1-tax[i];
		if(i) tax[i]+=tax[i-1];
	}
}

void part2(){
	static int tax[M];
	static pli dp[M];
	rep(u,1,K){
		int prod=1,sz=0;
		vi p; tax[0]=1;
		dp[++sz]=mp(u,1);
		LL lim=(LL)u*(b+m);
		rep(v,u,K){
			// [v*(b+1),u*(b+m)]
			if((LL)v*(b+1)>lim||!g[u][v]) break;
			// cerr<<"#"<<u<<" "<<v<<"   "<<sz<<" "<<prod<<endl;
			rep(i,1,sz){
				const pli &x=dp[i];
				LL w=lcm(x.fi,v);
				int val=mul(g[u][v],(u==v?x.se:mod-x.se));
				LL R=lim/w,L=((LL)v*(b+1)-1)/w;
				int tmp=mul(R/prod%mod,tax[prod-1]);
				add(tmp,tax[R%prod]);
				add(ans,mul(tmp,val));
				tmp=mul(L/prod%mod,tax[prod-1]);
				add(tmp,tax[L%prod]);
				sub(ans,mul(tmp,val));
			}
			int skip=0;
			rep(i,u+1,v-1){
				if(v%i==0) skip=1;
			}
			if(v==u||skip) continue;
			if(isp[v]&&(u<=K/2?prod*v<M:prod*v<M/5)){
				p.pb(v),prod*=v;
				rebuild(p,prod,tax);
			}
			else{
				rep(i,1,sz){
					const pli &x=dp[i];
					if(x.fi%v==0) continue;
					Hash::add(x.fi,x.se);
					LL nw=lcm(x.fi,v);
					if(nw>lim) continue;
					Hash::add(nw,mod-x.se);
				}
				Hash::solve(dp,sz);
			}
		}
	}
}

int main(){
	cin>>m>>n>>a>>b;
	a--,b--,part1(),part2();
	printf("%d\n",ans);
	return 0;
}

点赞——收藏——评论——关注

一键四连