Min

Min

从.html习得的一种算法

这个算法可以认为是洲阁筛的一种简易替代品(另一种写法),时空常数上和代码复杂度上都比洲阁筛优秀,时间复杂度基本和洲阁筛是一样的.
使用条件和洲阁筛一样,对积性函数求和.下面是具体流程

假设F(x)为积性函数,且如果x为质数, F(x)=∑xk F ( x ) = ∑ x k $F(x)=sum x^k$
按照洲阁筛的方法,对于每个 k k $k$分开统计贡献
对于每个不超过n" role="presentation">n−−√n$sqrt n$的质数 p p $p$,设g(x)=∑[x是质数或x的每个质因子都不小于p]xk" role="presentation">g(x)=∑[x是质数或x的每个质因子都不小于p]xkg(x)=∑[x是质数或x的每个质因子都不小于p]xk$g(x)=sum [x是质数或x的每个质因子都不小于p]x^k$
从小到大枚举 p p $p$,再从大到小枚举所有要求的x>=p2" role="presentation">x>=p2x>=p2$x>=p^2$,令 g(x)−=(g(x/p)−g(p−1))∗pk g ( x ) − = ( g ( x / p ) − g ( p − 1 ) ) ∗ p k $g(x)-=(g(x/p)-g(p-1))*p^k$
上面这个式子的意义大概是枚举合数中的最小质因子,再把它的贡献减去.这个部分复杂度和洲阁筛是一样的.

再定义S(n,x)代表 ∑ni=2[i的最小质因子不小于x]F(i) ∑ i = 2 n [ i 的 最 小 质 因 子 不 小 于 x ] F ( i ) $sum_{i=2} ^n[i的最小质因子不小于x]F(i)$
求S的时候现将质数的贡献加进来,这个可以通过g算出来
接下来从p开始往上枚举质数p
如果 p2 p 2 $p^2$大于 n n $n$就break,否则枚举正整数j" role="presentation">jj$j$,满足 pj+1<=n p j + 1 <= n $p^{j+1}<=n$,把 S(n/pj,p的下一个质数)∗F(pj)+F(pj+1) S ( n / p j , p 的 下 一 个 质 数 ) ∗ F ( p j ) + F ( p j + 1 ) $S(n/p^j,p的下一个质数)*F(p^j)+F(p^{j+1})$算入贡献
这部分的复杂度感觉不好证,不过据说是和洲阁筛差不多的

以SPOJDIVCNT3为例,具体实现见代码.这个做法还是比较容易写错的,而且姿势不好的话常数会大很多.
运行时间上最朴素写法12s,加上优化后4.7s,网上随便找到一个洲阁筛21.9s

#include<bits/stdc++.h>
using namespace std;
#define maxn 2000100
long long Min(long long x,long long y){return x<y?x:y;
}
int mysqrt(long long n){int x=sqrt(n);while(x*(long long)x<n) x++;while(x*(long long)x>n) x--;//必须满足S是最大的满足S*S<=n的数return x;
}
int S;
long long g[maxn];
int prim[maxn],ilem;
long long n;
void getg(){S=mysqrt(n),ilem=0;long long *l=g+S,*s=g;for(int i=1;i<=S;i++)s[i]=i,l[i]=n/i;//s[i]->g(i),l[i]->g(n/i),此时g(i)代表不超过i的质数的数量(包括1)for(int p=2;p<=S;p++){if(s[p]==s[p-1]) continue;//p不是质数long long r=p*(long long)p,v=s[p-1];int ed1=S/p,ed2=Min(n/r,(long long)S);for(int i=1;i<=ed1;i++)l[i]-=l[i*p]-v;if(n/p<INT_MAX){int m=n/p;for(int i=ed1+1;i<=ed2;i++)l[i]-=s[m/i]-v;}else{long long m=n/p;for(int i=ed1+1;i<=ed2;i++)l[i]-=s[m/i]-v;}
/*** 如果不这么写的话,常数大概会增加一倍左右.把int和long long分开在n=1e11的数据下有0.4s的优化(开O2) ***/
/*** 论减少除法数量的重要性... ***/for(int i=S;i>=r;i--)s[i]-=s[i/p]-v;prim[++ilem]=p;}prim[++ilem]=S+1;//这句很容易被遗漏for(int i=1;i<=S;i++)s[i]=s[i]*4,l[i]=l[i]*4;
}
long long F(long long x,int y){if(prim[y]>x) return 0;long long res=g[x<=S?x:S+n/x]-g[prim[y]-1],st;for(int i=y;i<=ilem;i++){st=prim[i];if(st*(long long)prim[i]>x) break;for(int j=1;;j++){if(st*prim[i]>x) break;res+=F(x/st,i+1)*(j*3+1)+(j*3+4);st=st*prim[i];}}return res;
}
int main(){
//  freopen("DIVCNT3.in","r",stdin);int T;scanf("%d",&T);while(T--){scanf("%lld",&n);getg();printf("%lldn",F(n,1)+1);//注意把1的贡献算进来.}return 0;
}