• BZOJ3309 DZY Loves Maths 莫比乌斯反演、线性筛


    传送门


    推式子(默认(N leq M)):

    (egin{align*} sumlimits_{i=1}^N sumlimits_{j=1}^Mf(gcd(i,j)) & = sumlimits_{d=1}^N f(d) sumlimits_{i=1}^frac{N}{d} sumlimits_{j=1}^frac{M}{d} sumlimits_{p | gcd(i,j)} mu(p) \ &= sumlimits_{d=1}^N f(d) sumlimits_{p = 1}^frac{N}{d} mu(p) frac{N}{dp} frac{M}{dp} \ &= sumlimits_{T=1}^N frac{N}{T} frac{M}{T} sumlimits_{d | T} f(d) mu(frac{T}{d}) end{align*})

    看这个数据范围毫无疑问要将(g(T) = sumlimits_{d | T} f(d) mu(frac{T}{d}))线性筛出来,然后就可以数论分块求解。

    因为(f(x))比较难处理,所以直接下手似乎不太好做。

    注意到因为每一项有(mu(d)),所以如果(T = p_1^{e_1}p_2^{e_2}...p_k^{e_k} , d = p_1^{e_1'}p_2^{e_2'}...p_k^{e_k'}),那么(e_i - 1 leq e_i' leq e_i)(f(d) mu(frac{T}{d}))才会有贡献。所以所有有贡献的(f(d))一定为(max e_i)或者(max e_i - 1),考虑这两种贡献对(g(T))的影响就可以得到(g(T))的值。

    假设(e_i)中能够取到(max e_i)的值有(x)个,那么:

    (x = k)时,只有(d = frac{T}{prodlimits_{i=1}^k p_i})时贡献为(max e_i - 1),此时(mu(frac{T}{d}) = (-1)^k)。而与此对应的,取到(max e_i)的所有(f(d) mu(frac{T}{d}))的和就是((-1)^{k+1} max e_i)。所以总贡献就是((-1)^{k+1})

    而当(x eq k)时,不难得到取到(max e_i - 1)(max e_i)的所有值的贡献都是(0)

    所以我们只有所有质因子指数相等的数有贡献。在线性筛的时候记录一下最小质因子的个数和去掉最小质因子之后会变成哪个数,就可以计算出所有会有贡献产生的(T)。总复杂度(O(n+qsqrt{n}))

    #include<iostream>
    #include<cstdio>
    #include<cstring>
    #include<algorithm>
    //This code is written by Itst
    using namespace std;
    
    inline int read(){
    	int a = 0;
    	char c = getchar();
    	bool f = 0;
    	while(!isdigit(c) && c != EOF){
    		if(c == '-')
    			f = 1;
    		c = getchar();
    	}
    	if(c == EOF)
    		exit(0);
    	while(isdigit(c)){
    		a = a * 10 + c - 48;
    		c = getchar();
    	}
    	return f ? -a : a;
    }
    
    const int MAXN = 1e7 + 7;
    int prime[MAXN] , num[MAXN] , pre[MAXN] , cntP[MAXN] , val[MAXN] , cnt;
    
    void init(){
    	for(int i = 2 ; i <= 1e7 ; ++i){
    		if(!cntP[i]){
    			prime[++cnt] = i;
    			num[i] = cntP[i] = 1;
    		}
    		for(int j = 1 ; j <= cnt && i * prime[j] <= 1e7 ; ++j){
    			if(i % prime[j] == 0){
    				cntP[i * prime[j]] = cntP[i];
    				num[i * prime[j]] = num[i] + 1;
    				pre[i * prime[j]] = pre[i];
    				break;
    			}
    			num[i * prime[j]] = 1;
    			pre[i * prime[j]] = i;
    			cntP[i * prime[j]] = cntP[i] + 1;
    		}
    	}
    	for(int i = 2 ; i <= 1e7 ; ++i)
    		if(!pre[i] || num[pre[i]] == num[i])
    			val[i] = (cntP[i] & 1 ? 1 : -1);
    		else num[i] = 0;
    	for(int i = 2 ; i <= 1e7 ; ++i)
    		val[i] += val[i - 1];
    }
    
    int main(){
    #ifndef ONLINE_JUDGE
    	freopen("in","r",stdin);
    	//freopen("out","w",stdout);
    #endif
    	init();
    	for(int T = read() ; T ; --T){
    		long long N = read() , M = read() , sum = 0;
    		if(N > M) N ^= M ^= N ^= M;
    		for(int i = 1 , pi ; i <= N ; i = pi + 1){
    			pi = min(N / (N / i) , M / (M / i));
    			sum += (N / i) * (M / i) * (val[pi] - val[i - 1]);
    		}
    		cout << sum << '
    ';
    	}
    	return 0;
    }
    
  • 相关阅读:
    安装包安装服务,点修复出现的错误”Error 1001:指定的服务已存在“ 解决办法
    C# tostring 格式化输出 (转)
    WPF button 如何区分click和doubleclick
    二叉树算法
    关于OA流程相关数据表的设计
    没事干写写流程审批数据库的设计
    挂载system.img并提取文件
    使Checkbox禁用,不可选
    Android显示网速的另一种方法
    C语言 解压华为固件
  • 原文地址:https://www.cnblogs.com/Itst/p/10502948.html
Copyright © 2020-2023  润新知