zoukankan      html  css  js  c++  java
  • 【luogu P4213】【模板】杜教筛(Sum)(数学)(整除分块)

    【模板】杜教筛(Sum)

    题目链接:luogu P4213

    题目大意

    要你求 φ 函数的前缀和和 μ 函数的前缀和。
    (分别是欧拉函数和莫比乌斯函数)

    思路

    前置知识(们)

    积性函数:对于两个互质的数 (x,y)(f(xy)=f(x)f(y)),那 (f) 就是积性函数。
    完全积性函数:对于任意两个整数 (x,y)(f(xy)=f(x)f(y)),那 (f) 就是完全积性函数。

    一些积性函数:(varphi,mu,d,sigma)
    (分别是欧拉函数,莫比乌斯函数,约数个数,约数个数和)
    一些完全积性函数:(epsilon,I,id)
    (epsilon(n)=[n=1],I(n)=1,id(n)=n)


    狄利克雷卷积:
    有两个函数 (f,g),它们的狄利克雷卷积:((f*g)(n)=sumlimits_{d|n}f(d)g(dfrac{n}{d}))

    不难看出,它满足交换律和结合律。
    然后单位元就是 (epsilon)

    然后又一些性质:(根据定义简单推一推都不难看出)
    (mu*I=epsilon)
    (varphi*I=id)
    (mu*id=varphi)


    莫反:

    如果 (g(n)=sumlimits_{d|n}f(d))
    那么 (f(n)=sumlimits_{d|n}mu(d)g(dfrac{n}{d}))

    这个地方的证明其实可以用 (mu*I=epsilon)

    给出条件相当于 (g=f*I)

    然后 (g*mu=f*I*mu=f*epsilon=f)

    杜教筛

    杜教筛就是用来求积性函数的前缀和 (sum(n)=sumlimits_{i=1}^nf(i))

    考虑再找一个积性函数 (g),求他们的狄利克雷卷积:
    (sumlimits_{i=1}^n(f*g)(i))
    (=sumlimits_{i=1}^nsumlimits_{d|n}f(d)g(dfrac{i}{d}))
    (=sumlimits_{d=1}^ng(d)sumlimits_{i=1}^{leftlfloorfrac{n}{d} ight floor}f(i))
    (=sumlimits_{d=1}^ng(d)sum(leftlfloordfrac{n}{d} ight floor))

    然后考虑 (g(1)sum(n)),根据前缀和:
    (=sumlimits_{i=1}^ng(i)sum(leftlfloordfrac{n}{i} ight floor)-sumlimits_{i=2}^ng(i)sum(leftlfloordfrac{n}{i} ight floor))
    (=sumlimits_{i=1}^n(f*g)(i)-sumlimits_{i=2}^ng(i)sum(leftlfloordfrac{n}{i} ight floor))

    那这个式子就是我们杜教筛要用的式子了。

    有什么用呢,大概就是对于每个你要求的 (f),如果你能找到一个合适的 (g),它和 (f) 乘起来很好算而且 (g) 也好搞的话就可以拿来用。

    就这题的例子,(mu*I=epsilon),所以左边部分就 (epsilon) 的前缀和是 (1),然后右边你可以用整除分块 (O(sqrt{n})) 搞。
    然后 (varphi*I=id),然后左边部分就 (id) 的前缀和是 (dfrac{(1+n)n}{2}),右边也是整除分块过去。
    (然后这两个右边 (g(i)) 因为是整除分块也可以同前缀和求出,因为是 (I) 所以就是 (r-l+1),即块的大小)

    然后你可以预处理出 (n^{frac{2}{3}}) 以内的结果,复杂度就可以压到 (O(n^{frac{2}{3}}))
    (然后你还可以用 (map) 弄个小小的记忆化)

    代码

    #include<map>
    #include<cstdio>
    #define ll long long
    
    using namespace std;
    
    //const int Maxn = 1664511;
    const int Maxn = 2000000;
    int T, np[Maxn + 1], prime[Maxn + 1];
    int miu[Maxn + 1], x;
    ll phi[Maxn + 1];
    map <int, int> ans_miu;
    map <int, ll> ans_phi;
    
    void init() {//预处理 n^{2/3} 的部分
    	miu[1] = phi[1] = 1;
    	for (int i = 2; i <= Maxn; i++) {
    		if (!np[i]) {
    			prime[++prime[0]] = i;
    			miu[i] = -1;
    			phi[i] = i - 1;
    		}
    		for (int j = 1; j <= prime[0] && i * prime[j] <= Maxn; j++) {
    			if (i % prime[j]) miu[i * prime[j]] = -miu[i], phi[i * prime[j]] = phi[i] * (prime[j] - 1), np[i * prime[j]] = prime[j];
    				else {
    					miu[i * prime[j]] = 0, phi[i * prime[j]] = phi[i] * prime[j], np[i * prime[j]] = prime[j];
    					break;
    				}
    		}
    	}
    	for (int i = 1; i <= Maxn; i++)//前缀和起来
    		miu[i] += miu[i - 1], phi[i] += phi[i - 1];
    }
    
    ll get_phi(int x) {
    	if (x <= Maxn) return phi[x];
    	if (ans_phi[x]) return ans_phi[x];
    	
    	ll re = 1ll * (1ll + x) * x / 2;//id 函数的前缀和
    	for (ll l = 2, r; l <= x; l = r + 1) {
    		r = x / (x / l);
    		re -= 1ll * (r - l + 1) * get_phi(x / l);
    	}
    	return ans_phi[x] = re;
    }
    
    int get_miu(int x) {
    	if (x <= Maxn) return miu[x];
    	if (ans_miu[x]) return ans_miu[x];
    	
    	ll re = 1;//ϵ 函数的前缀和
    	for (ll l = 2, r; l <= x; l = r + 1) {
    		r = x / (x / l);
    		re -= 1ll * (r - l + 1) * get_miu(x / l);
    	}
    	return ans_miu[x] = re;
    }
    
    int main() {
    	init();
    	
    	scanf("%d", &T);
    	while (T--) {
    		scanf("%d", &x);
    		printf("%llu %d
    ", get_phi(x), get_miu(x));
    	}
    	
    	return 0;
    } 
    
  • 相关阅读:
    CentOS挂载ntfs
    搬水果
    CentOS安装完显卡启动进不去系统
    二叉树遍历
    拓展的异常类
    Visual Studio 2008 IDE 出现 “帮助更新正在进行”错误的解决办法
    JSON工具SuperObject编译时"整数溢出"异常的解决办法
    引用命名区域
    ruby操作excel的几段代码
    Delphi2009&2007中,测试工程的初级使用(DUnit)
  • 原文地址:https://www.cnblogs.com/Sakura-TJH/p/luogu_P4213.html
Copyright © 2011-2022 走看看