zoukankan      html  css  js  c++  java
  • min25筛学习笔记

    4.min25筛

    听说这玩意能干杜教筛干不了的事?

    同杜教筛一样,这也是用来求积性函数前缀和的东西。其复杂度为 \(O(\dfrac{n^{0.75}}{\log n})\),大部分时候要略优于杜教筛。

    min25筛作用的积性函数,应保证对于一切质数 \(p\)\(f(p)\) 均是有关 \(p\) 的低阶多项式。同时,质数的正整数次幂的位置的 \(f\) 值要能够简单求出。

    既然是多项式,就可以计算令 \(f(p)=p^k\) 时函数的前缀和,然后分别加在一起即可。

    我们考虑将所有 \(i\) 按质数与合数分类:

    \[\sum\limits_{i=1}^nf(i)=\sum\limits_{p\in\mathtt{prime}}f(p)+\sum\limits_{i\notin\mathtt{prime}}f(i) \]

    接着,考虑枚举合数的最小质因子(该数一定 \(\leq\sqrt p\))及其次数。

    \[\sum\limits_{p\in\mathtt{prime}}f(p)+\sum\limits_{p\in\mathtt{prime}}\sum\limits_{p^a\leq n}\sum\limits_{i=1}^{\left\lfloor n/p^a\right\rfloor}[i\text{的最小质因子大于}p]f(ip^a) \]

    因为保证 \(f\) 是积性函数,所以可以变为

    \[\sum\limits_{p\in\mathtt{prime}}f(p)+\sum\limits_{p\in\mathtt{prime}}\sum\limits_{p^a\leq n}f(p^a)\sum\limits_{i=1}^{\left\lfloor n/p^a\right\rfloor}[i\text{的最小质因子大于}p]f(i) \]


    现在,我们考虑令 \(g(n,j)=\sum\limits_{i=1}^n[i\in\mathtt{prime}\lor i\text{的最小质因子大于第}j\text{个质数}]i^k\),其中 \(k\) 是一常数。明显,\(i^k\) 这个函数是完全积性的,但它仅在质数处与我们强制令 \(f(p)=p^k\) 得到的函数值相等。

    考虑由 \(g(n,j-1)\) 转移至 \(g(n,j)\)。显然,此时有些原本计入的 \(i\) 不能再被计入了。这样的 \(i\) 是所有最小质因子恰为第 \(j\) 个质数(设为 \(p_j\))的合数。则我们显然可以提出一个 \(p_j\) 来。于是

    \[g(n,j)=g(n,j-1)-(p_j)^k\Big(g(n/p_j,j-1)-g(p_{j-1},j-1)\Big) \]

    其中,最后又减掉的那个 \(g\) 是为了避免重复计算前 \(j-1\) 个质数的影响(它们会使得最小质因子并非 \(p_j\)

    因为 \(i^k\) 是完全积性函数,所以里面套着的 \(g(n/p_j,j-1)\) 可能还有的 \(p_j\) 因子也可继续相乘。

    因为 \(\left\lfloor\dfrac{\left\lfloor n/a\right\rfloor}{b}\right\rfloor=\left\lfloor\dfrac{n}{ab}\right\rfloor\),所以对于某个 \(n\),其所会访问到的状态一共只有所有的 \(\left\lfloor\dfrac{n}{x}\right\rfloor\),共 \(\sqrt n\) 个。可以对这 \(\sqrt n\) 个数建立值与下标的双射,这样就压缩了 \(g\) 所需要的下标范围。


    我们此时来设 \(h(n,j)=\sum\limits_{i=1}^n[i\text{的最小质因子大于第}j\text{个质数}]f(i)\)。则答案即为 \(h(n,0)\)

    考虑将 \(h(n,j)\) 划作两部分,即 \(>p_j\) 的质数的贡献以及最小质因子 \(>p_j\) 的合数的贡献。

    质数的贡献是 \(g(n,|\mathtt{prime}|)-\sum\limits_{i=1}^{j-1}f(p_i)\)

    合数的贡献是 \(\sum\limits_{k=j+1}^{|\mathtt{prime}|}\sum\limits_{(p_k)^a\leq n}f\Big((p_k)^a\Big)\sum\limits_{i=1}^{\left\lfloor n/(p_k)^a\right\rfloor}[i\text{的最小质因子大于}p_k]f(i)\)。后面一坨的定义刚好就是 \(h\Big(n/(p_k)^a,k\Big)\)。于是可以化为 \(\sum\limits_{k=j+1}^{|\mathtt{prime}|}\sum\limits_{(p_k)^a\leq n}f\Big((p_k)^a\Big)\Bigg(h\Big(n/(p_k)^a,k\Big)+[a\neq1]\Bigg)\)

    注意到后面还有一个 \([a\neq1]\),因为我们在min25筛的时候只考虑了质数与合数,没有考虑 \(1\),所以这里要补上;但是,当 \(a=1\) 时,如果考虑了 \(1\) 就会成为质数,已经在前面处理过了。这也意味着在min25筛执行完毕后,还要加上 \(1\) 的答案。

    总递归式为

    \[h(n,j)=g(n,|\mathtt{prime}|)-\sum\limits_{i=1}^{j-1}f(p_i)+\sum\limits_{k=j+1}^{|\mathtt{prime}|}\sum\limits_{(p_k)^a\leq n}f\Big((p_k)^a\Big)\Bigg(h\Big(n/(p_k)^a,k\Big)+[a\neq1]\Bigg) \]

    质数筛到 \(\sqrt n\) 就行了。

    I.【模板】Min_25筛

    代码:

    #include<bits/stdc++.h>
    using namespace std;
    typedef long long ll;
    const int mod=1e9+7;
    const int inv6=166666668;
    ll n;
    int m,pri[1001000],sp[2][1001000];//sp is the prefix sum of prime numbers' f, when 0 and 1 are the linear and quardratic terms respectively
    void sieve(){
    	for(int i=2;i<=m;i++){
    		if(!pri[i])pri[++pri[0]]=i,sp[0][pri[0]]=(sp[0][pri[0]-1]+i)%mod,sp[1][pri[0]]=(sp[1][pri[0]-1]+1ll*i*i)%mod;
    		for(int j=1;j<=pri[0]&&i*pri[j]<=m;j++){pri[i*pri[j]]=true;if(!(i%pri[j]))break;} 
    	}
    }
    ll kth[1001000];//the kth of n/i
    int tot,sml[1001000],lar[1001000];//sml is the bucket of elements no greater than m, while lar is for those greater than m.
    int g[2][1001000];//the g array can DP directly.
    int H(ll x,int y){
    	if(pri[y]>=x)return 0;
    	int X=(x<=m?sml[x]:lar[n/x]);
    	int ret=(0ll+g[1][X]-g[0][X]+mod-(sp[1][y]-sp[0][y])+mod)%mod;
    	for(int i=y+1;i<=pri[0]&&1ll*pri[i]*pri[i]<=x;i++){
    		ll pa=pri[i];
    		for(int a=1,tmp;pa<=x;a++,pa*=pri[i])tmp=pa%mod,(ret+=1ll*tmp*(tmp-1)%mod*(H(x/pa,i)+(a!=1))%mod)%=mod;
    	}
    //	printf("%d,%d:%d\n",x,y,ret);
    	return ret;
    }
    int main(){
    	scanf("%lld",&n),m=sqrt(n),sieve();
    //	for(int i=1;i<=pri[0];i++)printf("%d ",pri[i]);puts("");
    //	for(int i=1;i<=pri[0];i++)printf("%d ",sp[0][i]);puts("");
    //	for(int i=1;i<=pri[0];i++)printf("%d ",sp[1][i]);puts("");
    	for(ll i=1;i<=n;i=n/(n/i)+1){
    		kth[++tot]=n/i;
    		int I=kth[tot]%mod;//use the prefix sum formula of x and x^2 as the DP initial values of g(don't forget to ignore 1)
    		g[0][tot]=(1ll*I*(I+1)/2+mod-1)%mod,g[1][tot]=(1ll*I*(I+1)%mod*(2*I+1)%mod*inv6+mod-1)%mod;
    		if(kth[tot]<=m)sml[kth[tot]]=tot;else lar[n/kth[tot]]=tot;//because here we need to store values in buckets, so we can't use the moduloed values here.
    	}
    	for(int i=1;i<=pri[0];i++)for(int j=1;j<=tot&&1ll*pri[i]*pri[i]<=kth[j];j++){//cyclicize array g.
    		ll nexi=kth[j]/pri[i];int I=(nexi<=m?sml[nexi]:lar[n/nexi]);//get the id of n/p[i]
    		(g[0][j]+=mod-1ll*pri[i]*(g[0][I]+mod-sp[0][i-1])%mod)%=mod;
    		(g[1][j]+=mod-1ll*pri[i]*pri[i]%mod*(g[1][I]+mod-sp[1][i-1])%mod)%=mod;
    	}
    	printf("%d\n",(H(n,0)+1)%mod);
    	return 0;
    }
    

  • 相关阅读:
    PHP 实现 一致性哈希 算法(转的)
    一致性 hash 算法
    分布式设计与开发---宏观概述
    Lvs+keepalived+nginx+php的session 保持的算法
    从零搭建Web网站
    Java动态代理机制详解(JDK 和CGLIB,Javassist,ASM)
    哪个线程执行 CompletableFuture’s tasks 和 callbacks?
    HTTP 用户认证
    java 发送 HTTP 请求
    Http basic Auth 认证方式帮助类
  • 原文地址:https://www.cnblogs.com/Troverld/p/14619809.html
Copyright © 2011-2022 走看看