UOJ#188. 【UR #13】Sanrd（min25篩）

題目連結：https://uoj.ac/problem/188
題目大意：定義數論函式 f ( x ) f(x) f(x) 的值為 x x x 的第二大質因子，其中質因子可以重複貢獻，例如 f ( 36 ) = 3 f(36)=3 f(36)=3 ，且當且僅當 x ∈ p r i m e o r x = 1 , f ( x ) = 0 x\in prime\space or\space x=1,f(x)=0 x∈primeorx=1,f(x)=0 。給定區間 [ l , r ] [l,r] [l,r] ，求 ∑ i = l r f ( i ) \sum\limits_{i=l}^rf(i)

i=l∑r​f(i) 的值。
題解：對於這樣的一個數論函式，顯然我們幾乎是無法構造合適的數論函式 g ( x ) g(x) g(x) ，使得能夠利用杜教篩去求 f ( x ) f(x) f(x) 的字首和。那我們要往min25篩的解決方向去想，但這個函式並不是積性函式，要如何魔改我們的min25篩，從而可以求出這樣一個怪異的函式的字首和呢？回想一下min25篩中 S ( i , j ) S(i,j) S(i,j) 的轉移式，我們有 S ( i , j ) = g [ i ] − s u m [ j ] + ∑ k = j + 1 n u m ∑ r = 1 p k r ≤ i f ( p k r ) ( S ( ⌊ i p k r ⌋ , k ) + [ r > 1 ] ) S(i,j)=g[i]-sum[j]+\sum\limits_{k=j+1}^{num}\sum\limits_{r=1}^{p_{k}^{r}\leq i}f(p_k^r)(S(\lfloor\frac{i}{p_{k}^{r}}\rfloor,k)+[r>1])

S(i,j)=g[i]−sum[j]+k=j+1∑num​r=1∑pkr​≤i​f(pkr​)(S(⌊pkr​i​⌋,k)+[r>1]) ，其中前半部分的 g [ i ] − s u m [ j ] g[i]-sum[j] g[i]−sum[j] 為 f ( i ) f(i) f(i) 在質數點時的總和，後半部分為 f ( i ) f(i) f(i) 在合數點時的總和。不難發現min25篩的原理就是從每個數的最小質因數開始篩起，一直整除到其為質數或者質數的乘方。那麼這就好辦了，對與被篩到為質數的情況，只要知道上一次所被篩掉的質因數的值，再將 g ( x ) g(x) g(x) 定義為 [ 1 , x ] [1,x]

[1,x] 內質數的個數，把這個值乘上 g [ i ] − s u m [ j ] g[i]-sum[j] g[i]−sum[j] ，就是質數情況的答案了，而對於被篩到為質數乘方的情況，顯然答案一定為 p p p ，因為根據 f ( x ) f(x) f(x) 的定義， ∀ p ∈ p r i m e , r > 2 , f ( p r ) = p \forall p\in prime,r>2,f(p^r)=p ∀p∈prime,r>2,f(pr)=p 。此時我們篩法的轉移式便可以寫為 S ( i , j ) = p j ∗ ( g [ i ] − s u m [ j ] ) + ∑ k = j + 1 n u m ∑ r = 1 p k r ≤ i ( S ( ⌊ i p k r ⌋ , k ) + p k ∗ [ r > 1 ] ) S(i,j)=p_j*(g[i]-sum[j])+\sum\limits_{k=j+1}^{num}\sum\limits_{r=1}^{p_{k}^{r}\leq i}(S(\lfloor\frac{i}{p_{k}^{r}}\rfloor,k)+p_k*[r>1]) S(i,j)=pj​∗(g[i]−sum[j])+k=j+1∑num​r=1∑pkr​≤i​(S(⌊pkr​i​⌋,k)+pk​∗[r>1]) 。特別地，令 p 0 = 0 p_0=0 p0​=0 。然後套用min25篩的程式碼即可。
下面是AC程式碼。

#include<bits/stdc++.h>
using namespace std;
#define debug(x) cerr<<#x<<":"<<x<<endl
#define debug2(x,y) cerr<<#x<<":"<<x<<"  "<<#y<<":"<<y<<endl
#define debug3(x,y,z) cerr<<#x<<":"<<x<<"  "<<#y<<":"<<y<<"  "<<#z<<":"<<z<<endl
#define int unsigned long long
struct fastio{
    fastio(){
        ios::sync_with_stdio(0);
        cin.tie(0);
        cout.tie(0);
    }
}fio;
int N,sqN,num,cnt;
const int MAX=2e6+5;
int prime[MAX],part[MAX],sum[MAX],g[MAX],ind1[MAX],ind2[MAX],w[MAX];
inline int getind(int x){
    if(x>sqN) return ind2[N/x];
    return ind1[x];
}
void init(){
    sqN=sqrt(N);
    num=cnt=0;
    memset(prime,0,MAX);
    memset(part,0,MAX);
    memset(sum,0,MAX);
    memset(g,0,MAX);
    memset(ind1,0,MAX);
    memset(ind2,0,MAX);
    memset(w,0,MAX);
    for(int i=2;i<=sqN;++i){
        if(!part[i]){
            prime[++num]=i;
            part[i]=i;
            sum[num]=sum[num-1]+1;
        }
        for(auto j:prime){
            if(!j) continue;
            if(i*j>sqN) break;
            part[i*j]=j;
            if(i%j==0) break;
        }
    }
    for(int l=1,r;l<=N;l=r+1){//notice! g[cnt] must decrease 1!(because f(1) is out of range)
        w[++cnt]=N/l;
        r=N/w[cnt];
        if(w[cnt]>sqN) ind2[r]=cnt;
        else ind1[w[cnt]]=cnt;
        g[cnt]=w[cnt]-1;
    }
    for(int j=1;j<=num;++j){
        for(int i=1;i<=cnt&&prime[j]*prime[j]<=w[i];++i){
            int ind=getind(w[i]/prime[j]);
            g[i]-=g[ind]-sum[j-1];
        }
    }
}
int S(int i,int j){
    if(prime[j]>=i) return 0;
    int res=prime[j]*(g[getind(i)]-sum[j]);
    for(int k=j+1;k<=num&&prime[k]*prime[k]<=i;++k){
        for(int r=1,t=prime[k];t<=i;t*=prime[k],++r){//notice! t cannot be moduled!
            res+=S(i/t,k)+prime[k]*(r>1);
        }
    }
    return res;
}
signed main(){
    int L,R,ans;
    cin>>L>>R;
    N=L-1;
    init();
    ans=S(N,0);
    N=R;
    init();
    cout<<S(N,0)-ans<<endl;
}

第一次寫min25篩題，用了約3個小時吧（蒟蒻線上爬爬爬）。估計因為兩次開了預處理的原理，時間有點慢了。