zoukankan      html  css  js  c++  java
  • [莫比乌斯反演][三元环计数]「SDOI2018」旧试题

    • 题目传送门

    • 一个结论:(d(ijk)=sum_{a|i,b|j,c|k}[(a,b)=(a,c)=(b,c)=1])

    • 证明:分每个质因子独立考虑,因为 (a,b,c) 两两互质,故对于任意一个质因子 (p) 这三个数中最多存在一个数出现过 (p) 这个质因子,这样一个质因子 (p) 造成的乘积贡献为 (i,j,k)(p) 的个数之和再加 (1),也就是 (ijk)(p) 的个数之和再加 (1)

    • 于是开始反演

    • [sum_{i=1}^Asum_{j=1}^Bsum_{k=1}^Csum_{a|i,b|j,c|k}sum_{x|a,x|b}sum_{y|a,y|c}sum_{z|b,z|c}mu(x)mu(y)mu(z) ]

    • [=sum_{x,y,z}mu(x)mu(y)mu(z)sum_{lcm(x,y)|a}sum_{lcm(x,z)|b}sum_{lcm(y,z)|c}lfloorfrac Aa floorlfloorfrac Bb floorlfloorfrac Cc floor ]

    • (f(n,x)=sum_{x|i}lfloorfrac ni floor)(对于一个 (n) 可以 (O(nlog n)) 预处理出所有 (f(n,x))):

    • [sum_{x,y,z}mu(x)mu(y)mu(z)f(A,lcm(x,y))f(B,lcm(x,z))f(C,lcm(y,z)) ]

    • 这看上去还是三方的,但我们注意到 (f(n,>n)=0),而这个式子中 (f) 的第二个参数又是 (lcm),我们可以直观感觉到贡献不为 (0)((x,y,z)) 是很少的

    • 由于 ((x,y)(x,z)(y,z)) 同时出现,考虑建一个图,点集为 ([1,max(A,B,C)]) 中所有 (mu) 不为 (0) 的数,如果 (lcm(x,y)lemax(A,B,C)) 就连边 ((x,y)),我们要考虑的就是这个图所有的三元环,注意自环的特殊处理

    • 先考虑如何建图:设 (n=max(A,B,C)),先枚举 (d=(x,y)) 再枚举 (i=frac xd,j=frac yd),满足 (frac{xy}d=ijdle n) 也就是 (ijlefrac nd),如果 ((i,j)=1) 就连边 ((id,jd))

    • 这样的复杂度为 (sum_{d=1}^nsum_{i=1}^{lfloorfrac nd floor}lfloorfrac n{id} floor=sum_{d=1}^nO(frac ndlogfrac nd)<O(nlog^2n))

    • 建完这张图之后发现 (n=10^5) 时边数只有 (760741),于是可以枚举这张图的三元环

    • 三元环计数黑科技:把所有的点按照度数从大到小排序,度数相同情况下任意。然后把所有边定向,排序后在前面(度数大)的点(连向在后面(度数小)的点,形成一个 DAG

    • 枚举三元环上排在最前面的点 (u),把 (u) 的所有出点标记,然后枚举 (u) 的出边 (u ightarrow v)(v) 的出边 (v ightarrow w),如果 (w) 被标记则找到了一个三元环

    • 复杂度为 (O(msqrt m))(m) 为边数。分析((d_u) 为点 (u) 的度数,(cnt_u) 为点 (u) 的出度数):

    • 复杂度相当于对每条边的出点 (u)(cnt_u) 的和

    • 若一条边的出点 (u) 满足 (d_ulesqrt m),则 (cnt_u) 必然也不超过 (sqrt m),这一部分的复杂度 (O(msqrt m))

    • 否则设入点为 (v),则必然 (d_vge d_u>sqrt_m),由于图中 (d>sqrt m) 的点不超过 (sqrt m) 个,故每个 (cnt_u) 贡献不超过 (sqrt m) 次,又因为所有点的 (cnt) 之和为 (m),故这一部分复杂度也是 (O(msqrt m))

    Code

    #include <bits/stdc++.h>
    
    template <class T>
    inline void read(T &res)
    {
    	res = 0; bool bo = 0; char c;
    	while (((c = getchar()) < '0' || c > '9') && c != '-');
    	if (c == '-') bo = 1; else res = c - 48;
    	while ((c = getchar()) >= '0' && c <= '9')
    		res = (res << 3) + (res << 1) + (c - 48);
    	if (bo) res = ~res + 1;
    }
    
    typedef long long ll;
    
    const int N = 1e5 + 5, M = 2e6 + 5, djq = 1e9 + 7;
    
    int n, A, B, C, tot, pri[N], miu[N], ecnt, nxt[M], adj[N], go[M], val[M], d[N],
    fA[N], fB[N], fC[N], X[M], Y[M], Z[M], ans, m, a[N], id[N], vis[N];
    bool mark[N];
    
    void sieve()
    {
    	mark[0] = mark[miu[1] = 1] = 1;
    	for (int i = 2; i <= 100000; i++)
    	{
    		if (!mark[i]) miu[pri[++tot] = i] = -1;
    		for (int j = 1; j <= tot; j++)
    		{
    			if (1ll * i * pri[j] > 100000) break;
    			mark[i * pri[j]] = 1;
    			if (i % pri[j] == 0) break;
    			else miu[i * pri[j]] = -miu[i];
    		}
    	}
    }
    
    void add_edge(int u, int v, int w)
    {
    	nxt[++ecnt] = adj[u]; adj[u] = ecnt; go[ecnt] = v; val[ecnt] = w;
    }
    
    inline bool comp(int a, int b) {return d[a] > d[b];}
    
    void work()
    {
    	read(A); read(B); read(C);
    	n = std::max(std::max(A, B), C);
    	ecnt = m = ans = 0; int ec = 0;
    	for (int i = 1; i <= n; i++) adj[i] = 0;
    	for (int d = 1; d <= n; d++)
    		for (int i = 1; i <= n / d; i++) if (miu[i * d])
    			for (int j = 1; j <= n / d / i; j++)
    				if (miu[j * d] && std::__gcd(i, j) == 1)
    					X[++ec] = i * d, Y[ec] = j * d, Z[ec] = i * j * d;
    	for (int i = 1; i <= n; i++)
    	{
    		fA[i] = fB[i] = fC[i] = d[i] = 0;
    		for (int j = i; j <= n; j += i)
    			fA[i] += A / j, fB[i] += B / j, fC[i] += C / j;
    		if (miu[i]) a[++m] = i;
    	}
    	std::sort(a + 1, a + m + 1);
    	for (int i = 1; i <= m; i++) id[a[i]] = i;
    	for (int i = 1; i <= ec; i++) if (id[X[i]] < id[Y[i]])
    		add_edge(X[i], Y[i], Z[i]);
    	for (int i = 1; i <= m; i++)
    	{
    		int u = a[i];
    		for (int e = adj[u], v = go[e]; e; e = nxt[e], v = go[e])
    			vis[v] = val[e];
    		for (int e = adj[u], v = go[e]; e; e = nxt[e], v = go[e])
    		{
    			ans = (ans + 1ll * fA[u] * fB[val[e]] % djq * fC[val[e]] % djq
    				* miu[v] + djq) % djq;
    			ans = (ans + 1ll * fB[u] * fA[val[e]] % djq * fC[val[e]] % djq
    				* miu[v] + djq) % djq;
    			ans = (ans + 1ll * fC[u] * fA[val[e]] % djq * fB[val[e]] % djq
    				* miu[v] + djq) % djq;
    			ans = (ans + 1ll * fA[v] * fB[val[e]] % djq * fC[val[e]] % djq
    				* miu[u] + djq) % djq;
    			ans = (ans + 1ll * fB[v] * fA[val[e]] % djq * fC[val[e]] % djq
    				* miu[u] + djq) % djq;
    			ans = (ans + 1ll * fC[v] * fA[val[e]] % djq * fB[val[e]] % djq
    				* miu[u] + djq) % djq;
    			for (int x = adj[v], w = go[x]; x; x = nxt[x], w = go[x])
    				if (vis[w])
    				{
    					ll sum = 0;
    					sum += 1ll * fA[val[e]] * fB[val[x]] % djq * fC[vis[w]];
    					sum += 1ll * fA[val[e]] * fC[val[x]] % djq * fB[vis[w]];
    					sum += 1ll * fB[val[e]] * fA[val[x]] % djq * fC[vis[w]];
    					sum += 1ll * fB[val[e]] * fC[val[x]] % djq * fA[vis[w]];
    					sum += 1ll * fC[val[e]] * fA[val[x]] % djq * fB[vis[w]];
    					sum += 1ll * fC[val[e]] * fB[val[x]] % djq * fA[vis[w]];
    					if (sum %= djq, miu[u] * miu[v] * miu[w] == 1)
    						ans = (ans + sum) % djq;
    					else ans = (ans - sum + djq) % djq;
    				}
    		}
    		for (int e = adj[u], v = go[e]; e; e = nxt[e], v = go[e])
    			vis[v] = 0;
    		if (miu[u] == 1) ans = (1ll * fA[u] * fB[u] % djq * fC[u] + ans) % djq;
    		else ans = (ans - 1ll * fA[u] * fB[u] % djq * fC[u] % djq + djq) % djq;
    	}
    	printf("%d
    ", ans);
    }
    
    int main()
    {
    	sieve();
    	int T; read(T);
    	while (T--) work();
    	return 0;
    }
    
  • 相关阅读:
    PAIRING WORKFLOW MANAGER 1.0 WITH SHAREPOINT 2013
    Education resources from Microsoft
    upgrade to sql server 2012
    ULSViewer sharepoint 2013 log viewer
    Top 10 Most Valuable Microsoft SharePoint 2010 Books
    讨论 Setsockopt选项
    使用 Alchemy 技术编译 C 语言程序为 Flex 可调用的 SWC
    Nagle's algorithm
    Nagle算法 TCP_NODELAY和TCP_CORK
    Design issues Sending small data segments over TCP with Winsock
  • 原文地址:https://www.cnblogs.com/xyz32768/p/12422811.html
Copyright © 2011-2022 走看看