zoukankan      html  css  js  c++  java
  • P4619 [SDOI2018]旧试题

    P4619 [SDOI2018]旧试题

    hoho!!!终于做掉了!!!Orz (color{black}{ exttt{c}}color{red}{ exttt{yn2006}})

    这个 (d(ijk)) 一看就是要拆的。不然那根本没法下手。

    (d(S)) 的本质入手。(S) 的一部分来自 (i) ,一部分来自 (j) ,一部分来自 (k)

    考虑直接直接 (d(i)d(j)d(k)) 算重的原因,发现是来自 (i,j,k) 的部分有些不互质,互质就能保证 (ijk) 的每一个因数只算一次。

    所以我们有 (d(ijk)=sum_{x|i}sum_{y|j}sum_{z|k}[gcd(x,y)=1][gcd(y,z)=1][gcd(z,x)=1])

    带回去推式子

    [sum_{i=1}^{A}sum_{j=1}^{B}sum_{k=1}^{C}d(ijk)\ =sum_{i=1}^{A}sum_{j=1}^{B}sum_{k=1}^{C}sum_{x|i}sum_{y|j}sum_{z|k}[gcd(x,y)=1][gcd(y,z)=1][gcd(z,x)=1]\ =sum_{x=1}^{A}sum_{y=1}^{B}sum_{z=1}^{C}[gcd(x,y)=1][gcd(y,z)=1][gcd(z,x)=1]dfrac{A}{x}dfrac{B}{y}dfrac{C}{z}\ =sum_{x=1}^{A}sum_{y=1}^{B}sum_{z=1}^{C}sum_{a|gcd(x,y)}mu(a)sum_{b|gcd(y,z)}mu(b)sum_{c|gcd(z,x)}mu(c)dfrac{A}{x}dfrac{B}{y}dfrac{C}{z}\ =sum_{a=1}^{A}sum_{b=1}^{B}sum_{c=1}^{C}mu(a)mu(b)mu(c)sum_{x=1}^{A}sum_{y=1}^{B}sum_{z=1}^{C}sum_{a|x,a|y}sum_{b|y,b|z}sum_{c|x,c|z}dfrac{A}{x}dfrac{B}{y}dfrac{C}{z}\ =sum_{a=1}^{A}sum_{b=1}^{B}sum_{c=1}^{C}mu(a)mu(b)mu(c)sum_{x=1}^{frac{A}{operatorname{lcm}(a,c)}}sum_{y=1}^{frac{B}{operatorname{lcm}(a,b)}}sum_{z=1}^{frac{C}{operatorname{lcm}(b,c)}}dfrac{A}{xoperatorname{lcm}(a,c)}dfrac{B}{yoperatorname{lcm}(a,b)}dfrac{C}{zoperatorname{lcm}(b,c)}\ ]

    最后面那个 (sum_{x=1}^{frac{A}{operatorname{lcm}(a,c)}}sum_{y=1}^{frac{B}{operatorname{lcm}(a,b)}}sum_{z=1}^{frac{C}{operatorname{lcm}(b,c)}}dfrac{A}{xoperatorname{lcm}(a,c)}dfrac{B}{yoperatorname{lcm}(a,b)}dfrac{C}{zoperatorname{lcm}(b,c)}) 其实是可以预处理之后 (O(log)) 的。

    [sum_{i=1}^{frac{A}{operatorname{lcm}(a,c)}}dfrac{A}{xoperatorname{lcm}(a,c)}=sum_{i=1}^{frac{A}{operatorname{lcm}(a,c)}}dfrac{dfrac{A}{operatorname{lcm}(a,c)}}{x} ]

    我们预处理 (f_x=sum_{i=1}^{x}dfrac{x}{i}) 就可以 (O(log)) 搞上面那个东西。这个东西水爆了,不会建议回炉重造。

    但是千万注意这里有个 (log) ,不要像我一样傻乎乎的忽略了 (operatorname{lcm}) 的复杂度当成 (O(1)) ,不然后面会死的很惨。

    发现完全没法往下做了。。。

    这时候需要一点数学(数论)的直觉。

    首先 (mu(a)) 不能为 (0) ,其次 (operatorname{lcm}(a,b) le max(A,B,C)) ,感觉这个东西特别少的样子。

    写个暴力输出一下 (10^5) 以内 ((a,b)) 的个数,你发现纯暴力直接T没。

    于是我们开始优化这个暴力求个数的过程。

    首先套路的枚举一下 (gcd) 。我们在对于每一个数 (x) 开vector存它倍数,如果 (mu(x)) 不为 (0) 那么我们平方遍历它所有的倍数, (gcd(d_j,d_k)=x) 的时候计数器++。

    发现这个东西跑的飞快,而且,计数器只有 821535 !!!

    有一个特别牛逼的结论,那就是三元环的数量上限是 (O(msqrt m)) 级别的。这个可以参考网上一大堆三元环计数的博客,我们是通过枚举来统计答案的。

    我们惊奇的发现可以枚举 ((a,b,c)) 这种三元环了,只要把上面统计个数直接改成拉边即可,而且数量不会很大。上面的暴力没有白写,随便改改就是正解。。。

    对了有个显然的事情,邻接表拉边快遍历慢,vector 拉边慢遍历快,而且有O2加持。这题的瓶颈显然是遍历吧,拉边不过8e5次,应该用 vector 存图。不过邻接表会慢成什么样子没试过,小概率比 vector 快。

    可是 (msqrt m) 上界不是 (10^9) 吗?

    别忘了自然数是很神奇的,这又不是出题人手造恶意卡的图。。。

    事实上上限大概 (10^7)

    吼啊,直接枚举!就完结了???

    有个细节:这个三元环可以自环。于是我们要分:三个点都不同、三个点都相同、三个点有两个相同三种情况来讨论(常数。。。)

    极限数据应该是

    2
    100000 100000 100000
    100000 100000 100000
    

    因为这个时候那个 (sqrt m) 最大。

    之前提过了 (operatorname{lcm}) 是有 (log) 的。。。于是直接T飞,但是只用8s,没有特别离谱

    首先一波乱搞发现 (operatorname{lcm}) 不用开 long long ,就是预处理的时候把那些 (operatorname{lcm}) 大于 (10^5) 的直接舍掉,然后 先除以 (gcd) 再乘,即 lcm(x,y)=x/gcd(x,y)*y ,这样就不用 long long 了。

    还有上面说可以预处理 (O(log)) 的那个式子,我们发现他可以不用取模,只不过极限情况有可能炸 long long (甚至可能根本炸不掉)。于是我们可以先拿一个大的数存下来,最后再取模,减少取模次数

    发现两个优化很管用,本地直接跑到了 (5.2s)

    勇了一发,最大点5.00s,AC了!!!喜提次劣解。woc怎么还有更慢的???

    upd:更慢的那个不知道哪里去了,最劣解坐稳了。。。

    注释代码来测瓶颈。发现瓶颈还是 (operatorname{lcm}) 。发现这些边的 (operatorname{lcm}) 是可以直接预处理的。于是大大减少了 (operatorname{lcm}) 的次数。

    成功跑进了 (3s) ,稳了。

    //Orz cyn2006
    #include<bits/stdc++.h>
    using namespace std;
    #define fi first
    #define se second
    #define mkp(x,y) make_pair(x,y)
    #define pb(x) push_back(x)
    #define sz(v) (int)v.size()
    typedef long long LL;
    typedef double db;
    template<class T>bool ckmax(T&x,T y){return x<y?x=y,1:0;}
    template<class T>bool ckmin(T&x,T y){return x>y?x=y,1:0;}
    #define rep(i,x,y) for(int i=x,i##end=y;i<=i##end;++i)
    #define per(i,x,y) for(int i=x,i##end=y;i>=i##end;--i)
    inline int read(){
        int x=0,f=1;char ch=getchar();
        while(!isdigit(ch)){if(ch=='-')f=0;ch=getchar();}
        while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
        return f?x:-x;
    }
    const int N=100005;
    const int M=821535+5;
    #define mod 1000000007
    void fmod(int&x){x+=x>>31&mod,x-=mod,x+=x>>31&mod;}
    int A,B,C,mx,ans;
    int mu[N],pri[N],pct,sum[N],tot,tag[N];
    bool vis[N];
    vector<int>d[N];
    int U[M],V[M],ord[N],deg[N],tu[M],tv[M],cnt,LCM[M],id[M];
    vector<pair<int,int> >e[N];
    pair<int,int>p[N];
    void init(const int&N){
    	mu[1]=1;
    	for(int i=2;i<=N;++i){
    		if(!vis[i])pri[++pct]=i,mu[i]=-1;
    		for(int j=1;j<=pct&&i*pri[j]<=N;++j){
    			vis[i*pri[j]]=1;
    			if(i%pri[j]==0){mu[i*pri[j]]=0;break;}
    			mu[i*pri[j]]=-mu[i];
    		}
    	}
    	for(int i=1;i<=N;++i){
    		if(!mu[i])continue;
    		for(int j=1;i*j<=N;++j)if(mu[i*j])d[i].pb(i*j);
    	}
    	for(int i=1;i<=N;++i)
    		for(int j=0,up=sz(d[i]);j<up;++j)
    			for(int k=j;k<up&&1ll*d[i][j]*d[i][k]/i<=N;++k)
    				if(__gcd(d[i][j],d[i][k])==i)U[++tot]=d[i][j],V[tot]=d[i][k],LCM[tot]=d[i][j]/i*d[i][k];
    	for(int i=1;i<=N;++i){
    		int res=0;
    		for(int l=1,r;l<=i;l=r+1)
    			r=i/(i/l),res+=(i/l)*(r-l+1);
    		sum[i]=res;
    	}
    }
    int lcm(int x,int y){return x/__gcd(x,y)*y;}
    LL f(int x,int y,int z){return 1ll*sum[x]*sum[y]*sum[z];}
    void Main(){
    	A=read(),B=read(),C=read(),mx=max(A,max(B,C)),fill(deg+1,deg+mx+1,0),fill(tag+1,tag+mx+1,0),ans=0,cnt=0;
    	rep(i,1,tot)if(LCM[i]<=mx)++deg[U[i]],++deg[V[i]],++cnt,tu[cnt]=U[i],tv[cnt]=V[i],id[cnt]=i;
    	rep(i,1,mx)p[i]=mkp(deg[i],i);
    	sort(p+1,p+mx+1);
    	rep(i,1,mx)ord[p[i].se]=i,e[i].clear();
    	rep(i,1,cnt)
    		if(ord[tu[i]]>ord[tv[i]])e[tu[i]].pb(mkp(tv[i],LCM[id[i]]));
    		else e[tv[i]].pb(mkp(tu[i],LCM[id[i]]));
    	rep(u,1,mx){
    		for(int i=0,up=sz(e[u]);i<up;++i)tag[e[u][i].fi]=u;
    		for(int i=0,up1=sz(e[u]);i<up1;++i){
    			int v=e[u][i].fi;
    			for(int j=0,up2=sz(e[v]);j<up2;++j){
    				int w=e[v][j].fi;
    				if(tag[w]==u){
    					int x=e[u][i].se,y=e[v][j].se,z=lcm(w,u);unsigned long long res=0;
    					if(u!=v&&u!=w&&v!=w)
    						res=f(A/x,B/y,C/z)+f(A/x,B/z,C/y)+f(A/y,B/x,C/z)+f(A/y,B/z,C/x)+f(A/z,B/x,C/y)+f(A/z,B/y,C/x);
    					else if(u==v&&u==w&&v==w)res=f(A/x,B/y,C/z);
    					else res=f(A/x,B/y,C/z)+f(A/y,B/z,C/x)+f(A/z,B/x,C/y);
    					int t=res%mod;
    					fmod(ans+=1ll*mu[u]*mu[v]*mu[w]*t);
    				}
    			}
    		}
    	}
    	printf("%d
    ",ans);
    }
    signed main(){
    	init(100000);
    	for(int T=read();T;--T)Main();
    }
    /*
    5
    10 10 10
    100 100 100
    1000 1000 1000
    10000 10000 10000
    100000 100000 100000
    
    2
    100000 100000 100000
    100000 100000 100000
    
    */
    
  • 相关阅读:
    时间计算
    DateTime
    C# trim split dataGrid
    something
    生活
    如何导入外部的源码到eclipse中
    java类中获取ServletContext的方法
    获取spring容器上下文(webApplicationContext)的几种方法
    java反射(转)
    mysql常见命令
  • 原文地址:https://www.cnblogs.com/zzctommy/p/14050558.html
Copyright © 2011-2022 走看看