[模板]Min_25筛
用途
快速($O(\frac{n^{3/4}}{logn})$)地计算一些函数f的前缀和，以及（作为中间结果的）只计算质数的前缀和
一般要求$f(p)$是积性函数，$f(p)$是多项式的形式，且$f(p^k)$可以快速计算
做法
首先考虑求出范围内的质数的取值的和
如果有$f(p)=\sum{a_ip^i}$
那么我们构造$h_i(x)=x^i$，不难发现$h_i$是完全积性的
就是说，把f在质数的时候的式子拆开，然后让它在不是质数的时候也成立
考虑求其中的一个h，接下来设$pri_j$是第j个质数
设$g(n,j)=\sum\limits_{i=1}^{n}[i \in Prime OR min\_divisor(i)>pri_j]h(i)$，即n范围内质数或者最小质因数$>pri_j$的h值之和
可以当成是埃氏筛法的过程，已经筛掉了前j个质数
那么$g(n,inf)$就是要求的质数的和，g(n,0)就是所有数（1以外的）的h的和
（关于计算g(n,0)，次数要是小的话可以直接代公式，大的话就要斯特林数或者别的什么来求自然数幂和了..不过别忘了减掉1）
考虑转移，有:
$g(n,j)=g(n,j-1) , pri_j*pri_j>n$
$g(n,j)=g(n,j-1)-h(pri_j)(g(\lfloor\frac{n}{pri_j}\rfloor,j-1)-\sum\limits_{i=1}^{j-1}h(pri_i)) , pri_j*pri_j<=n$
考虑$pri_j*pri_j<=n$的情况，它在从$g(n,j-1)$转移过来的时候，需要减掉那些最小质因子是$pri_j$的，而h又是完全积性的，所以可以直接乘。但这样减会多减掉那些比$pri_j$小的质数乘$pri_j$，需要再加回来
于是我们就求出来质数的h的和啦！然而这有什么卵用呢..
设$S(n,j)=\sum\limits_{i=1}^{n}[min\_divisor(i)>=pri_j]f(i)$，那么$S(N,1)+f(1)$就是要求的答案
考虑转移，有$S(n,j)=\sum a_ig_i(n,inf) - \sum\limits_{i=1}^{j-1}f(pri_i) + \sum\limits_{k=j}^{pri_k*pri_k<=n}\sum\limits_{e=1}^{pri_k^{e+1}<=n}f(pri_k^e)S(\lfloor \frac{n}{pri_k^e} \rfloor,k+1)+f(pri_k^{e+1})$
就是说，既然我们已经求出了质数的情况，那就先把符合要求的质数加上，然后再加上合数，枚举它的质因子然后都除下去就好了（因为这里已经没有完全积性了，所以一次要把某个质因子都除掉），因为S里是没有1的，所以还要额外加上$p^{e+1}$的情况
例题
loj6053 简单的函数
注意到除了f(2)=3以外，f(p)=p-1
于是拆成p和1去做即可。最后算S的时候，如果j<=1，那再加个2
1 #include<bits/stdc++.h> 2 #include<tr1/unordered_map> 3 #define CLR(a,x) memset(a,x,sizeof(a)) 4 #define MP make_pair 5 using namespace std; 6 typedef long long ll; 7 typedef unsigned long long ull; 8 typedef pair<int,int> pa; 9 const int maxn=1e6+10,P=1e9+7; 10 11 inline ll rd(){ 12 ll x=0;char c=getchar();int neg=1; 13 while(c<'0'||c>'9'){if(c=='-') neg=-1;c=getchar();} 14 while(c>='0'&&c<='9') x=x*10+c-'0',c=getchar(); 15 return x*neg; 16 } 17 18 int pri[maxn],sump[maxn],cnt,tot,sqrtn; 19 bool np[maxn]; 20 ll N,w[maxn]; 21 tr1::unordered_map<ll,int> id; 22 int g[maxn],h[maxn]; 23 24 inline int S(ll n,int j){ 25 if(n<=1||pri[j]>n) return 0; 26 int re=(0ll+g[id[n]]-h[id[n]]-sump[j-1]+j-1)%P; 27 if(j==1) re=(re+2)%P; 28 for(int i=j;i<=cnt&&1ll*pri[i]*pri[i]<=n;i++){ 29 ll tmp=pri[i]; 30 for(int k=1;tmp*pri[i]<=n;k++){ 31 re=(re+1ll*S(n/tmp,i+1)*(pri[i]^k)+(pri[i]^(k+1)))%P; 32 tmp=tmp*pri[i]; 33 } 34 }return re; 35 } 36 37 int main(){ 38 //freopen("","r",stdin); 39 ll i,j; 40 N=rd();sqrtn=sqrt(N); 41 for(i=2;i<=sqrtn;i++){ 42 if(!np[i]){ 43 pri[++cnt]=i; 44 sump[cnt]=(sump[cnt-1]+i)%P; 45 } 46 for(j=1;j<=cnt&&pri[j]*i<=sqrtn;j++){ 47 np[pri[j]*i]=1; 48 if(i%pri[j]==0) break; 49 } 50 } 51 for(i=1;i<=N;i=j+1){ 52 j=N/(N/i); 53 w[++tot]=N/i;id[w[tot]]=tot; 54 g[tot]=1ll*(w[tot]+2)%P*((w[tot]-1)%P)%P; 55 if(g[tot]&1) g[tot]+=P; 56 g[tot]=g[tot]/2%P; 57 h[tot]=(w[tot]-1)%P; 58 } 59 for(j=1;j<=cnt;j++){ 60 for(i=1;i<=tot&&1ll*pri[j]*pri[j]<=w[i];i++){ 61 g[i]=(g[i]-1ll*pri[j]*(g[id[w[i]/pri[j]]]-sump[j-1]))%P; 62 h[i]=(0ll+h[i]-(h[id[w[i]/pri[j]]]-j+1))%P; 63 } 64 } 65 printf("%d\n",((S(N,1)+1)%P+P)%P); 66 67 return 0; 68 }
