Min_25 篩用來求一類特殊數論函式的字首和。

要求的函式 \(f\) 必須滿足以下條件：

1. \(f\) 是積性函式。
2. 對於質數 \(p\)，\(f(p)\) 是一個關於 \(p\) 的簡單低次多項式。
3. 對於質數 \(p\) 和正整數 \(e\)，\(f(p^e)\) 可以快速計算。

例題 Luogu P5325

轉化題目

題目中的函式 \(f(p^k)=p^k(p^k-1)\) 顯然滿足要求。

即對於 \(x=p^k\)，\(f(x)=x^2-x\)，其中 \(p\) 是質數，\(k\) 是正整數。

顯然我們可以將這樣的函式拆成兩個，即對於 \(x=p^k\)，把 \(f(x)\) 拆成 \(f_1(x)=x\)

和 \(f_2(x)=x^2\)，分別計算它們，求和後就是原來的 \(f\)。

在接下來的討論中，我們設 \(f'(x)=x^k\)。分別將 \(k=1,2\) 帶入並求和，就得到了答案。

觀察到 \(f'\) 是完全積性函式，且在質數處和 \(f_k\) 取值相同。

來做一些奇怪的工作

設 \(\mathbb P\) 表示質數集， \(P_i\) 表示第 \(i\) 個質數。特別地，設 \(P_0=0\)。設 \(L(i)\) 表示 \(i\) 的最小質因子。考慮有這樣一個數組 \(g\)，它滿足以下定義：

\[g(n,j)=\sum \limits_{1 \leq i \leq n \and (i \in \mathbb P \or L(i)>P_j)} f'(i) \]

其中 \(j \leq Q\)

，\(Q\) 是最大的滿足 \(P_Q \leq \sqrt n\) 的正整數。

即，在 \(1\) 到 \(n\) 的所有 \(i\) 中，對所有滿足 \(i\) 是質數或 \(i\) 的最小質因子大於 \(P_j\) 的 \(i\) 的 \(f'\) 值求和。

考慮 \(g(n,j)\) 的轉移。如果 \(n<P_j^2\)，那麼 \(n\) 的最小質因子必然不會為 \(P_j\)，所以這個限制不會對貢獻產生影響，於是有 \(g(n,j)=g(n,j-1)\)。

如果 \(n \geq P_j^2\) 怎麼辦呢？這種情況下，\(g(n,j)\) 會在 \(g(n,j-1)\) 的基礎上減少。觀察到 \(P_{j-1},P_j\)

是連續的兩個質數，於是從貢獻中消失的 \(i\) 的最小質因子一定是 \(P_j\)。這一部分 \(i\) 的貢獻即是 \(f'(P_j) \times( g(\left \lfloor \dfrac{n}{P_j} \right \rfloor,j-1)-g(P_{j-1},j-1))\)。為什麼呢？觀察式子開頭的 \(f'(P_j)\)，這是因為 \(f'\) 是完全積性函式，且從貢獻中消失的 \(i\) 的最小質因子一定是 \(P_j\)，我們可以將這一部分的貢獻提出。\(g(\left \lfloor \dfrac{n}{P_j} \right \rfloor,j-1)\) 中的 \(i\)（姑且記作 \(i'\)）滿足 \(L(i')>P_{j-1}\) 或 \(i' \in \mathbb P\)。

對於 \(L(i')>P_{j-1}\) 這一部分，\(i'\) 的最小質因子大於 \(P_{j-1}\)，那麼 \(i' \times P_j\)（為什麼要乘上 \(P_j\)？因為我們將它提到外面去了）的最小質因子只可能為 \(P_j\)，\(i= i' \times P_j\) 正是貢獻中應該消失的 \(i\)。

為什麼要加上 \(g(P_{j-1},j-1)\) 呢？這對應著第二部分 \(i' \in \mathbb P\)。對於質數 \(i'\)，僅當 \(i' \times P_j\) 的最小質因子為 \(P_j\) 時才被刪掉。觀察到對於小於 \(P_j\) 的質數 \(i'\)，我們不應該刪去 \(i =i'\times P_j\) 的貢獻。同時，觀察到 \(i'\) 要小於等於 \(\left \lfloor \dfrac{n}{P_j} \right \rfloor\) 才會錯誤地刪掉貢獻。因為 \(P_j\) 一定不大於 \(\sqrt n\)，所以有 \(P_j \leq \left \lfloor \dfrac{n}{P_j} \right \rfloor\)，這個限制可以忽略。我們要加回來的就是小於 \(P_j\) 的質數的貢獻，即前 \(j-1\) 個質數的 \(f'\) 值之和。由 \(g\) 的定義可知，這和 \(g(P_{j-1},j-1)\) 相等。

於是我們得到了 \(g\) 的遞推式：

\[g(n,j) = \begin{cases} g(n,j-1) & P_j^2 > n \\ g(n,j-1)-f'(P_j) \times( g(\left \lfloor \dfrac{n}{P_j} \right \rfloor,j-1)-g(P_{j-1},j-1)) & P_j^2 \leq n\end{cases} \]

求解答案

設 \(S(n,x)\) 表示 \(\sum \limits_{1 \leq i \leq n} [L(i)> P_x] f_k(i)\)，即 \(1 \sim n\) 中所有最小質因子大於 \(P_x\) 的數的 \(f_k\) 值之和，注意是 \(f_k\) 而非 \(f'\)。

分類討論：

• 如果符合條件的 \(i\) 是質數

  那麼這一部分 \(i\) 之和為 \(g(n,Q)- \sum \limits_{i=1}^{x} f_k(P_i)\)。後面的和式可以預處理的時候做一個字首和算出來。

• 如果符合條件的 \(i\) 是合數

  列舉最小質因子 \(p_j^e\)。這一部分的貢獻即為

  \[\sum \limits_{p_j^e \leq n \and j > i }f_k(p^e) \times (S(\left \lfloor \dfrac{n}{p_j^e}\right \rfloor,j)+[e>1]) \]

  \(+[e>1]\) 是因為 \(p_j^e\) 在 \(e>1\) 的時候本身就是合數，應該被算進去。而 \(e=1\) 的時候就是質數，不在這個分類裡面。

綜合以下，我們可以得到

\[S(n,x) = g(n,Q)- \sum \limits_{i=1}^{x} f_k(P_i) \sum \limits_{p_j^e \leq n \and j > i }f_k(p^e) \times (S(\left \lfloor \dfrac{n}{p_j^e}\right \rfloor,j)+[e>1]) \]

這樣問題就解決啦，答案就是 \(S(n,0)\)。值得一提的是，我們沒有計算 \(f(1)\) 的值，因此還要加上 \(1\)。

時間複雜度據說是 \(O(\dfrac{n^{\frac34}}{\log n})\)，大概是 \(1000 \text{ms}\) 跑 \(n=10^{10}\)。

# include <bits/stdc++.h>

const int N=1000010,INF=0x3f3f3f3f,MOD=1e9+7,INV6=166666668,INV2=500000004;
typedef long long ll;

ll prime[N],psum;
bool vis[N];
ll g1[N],g2[N],prsum1[N],prsum2[N],w[N],wtot,idx1[N],idx2[N];

ll MAXN;

ll n;

inline int read(void){
	int res,f=1;
	char c;
	while((c=getchar())<'0'||c>'9')
		if(c=='-')f=-1;
	res=c-48;
	while((c=getchar())>='0'&&c<='9')
		res=res*10+c-48;
	return res*f;
}
inline void init(void){ // 預處理
	for(int i=2;i<=MAXN;++i){
		if(!vis[i]){
			prime[++psum]=i;
			prsum1[psum]=(prsum1[psum-1]+i)%MOD;
			prsum2[psum]=(prsum2[psum-1]+1ll*i*i)%MOD;
		}
		for(int j=1;i*prime[j]<=MAXN&&j<=psum;++j){
			vis[i*prime[j]]=true;
			if(i%prime[j]==0){
				break;
			}
		}
	}
	return;
}
inline ll getsum1(ll x){ // 記得取模
	x%=MOD;
	return x*(x+1)%MOD*INV2%MOD;
}
inline ll getsum2(ll x){
	x%=MOD;
	return x*(x+1)%MOD*(2ll*x+1)%MOD*INV6%MOD;
}
inline int getid(ll x){ // trick: n/x 只有 O(\sqrt n) 種取值，於是我們開大小為 2\sqrt n 的陣列即可，idx1[x] 儲存 x <= \sqrt n 時 n/x 的取值，idx2[n/x] 儲存 x > \sqrt n 時 n/x 的取值
	return (x<=MAXN)?idx1[x]:idx2[n/x];
}
ll S(ll x,ll i){ // 遞迴計算 S
	if(prime[i]>=x){
		return 0;
	}
	int nowid=getid(x);
	ll res=(((g2[nowid]-g1[nowid])-(prsum2[i]-prsum1[i]))%MOD+MOD)%MOD;
	for(int j=i+1;j<=psum&&prime[j]*prime[j]<=x;++j){
		for(ll k=1,nowp=prime[j];nowp<=x;++k,nowp=nowp*prime[j]){
			ll truenowp=nowp%MOD;
			res=(res+(truenowp-1)*truenowp%MOD*(S(x/nowp,j)+(k>1)))%MOD;
		}
	}
	return res;
}
int main(void){
	scanf("%lld",&n);
	MAXN=1ll*sqrt(n);
	init();
	for(ll l=1,r;l<=n;l=r+1){
		ll val=n/l;
		r=n/val,w[++wtot]=val;
		g1[wtot]=getsum1(val)-1,g2[wtot]=getsum2(val)-1; // 預處理出 g(...,0)
		if(val<=MAXN){
			idx1[val]=wtot;
		}else{
			idx2[n/val]=wtot;
		}
	}
    // g1,g2 分別代表 k=1,2 時 g 陣列的值
	for(int i=1;i<=psum;++i){ // 滾動計算 g
		for(int j=1;prime[i]*prime[i]<=w[j]&&j<=wtot;++j){
			int nowid=getid(w[j]/prime[i]);
			g1[j]=(g1[j]-prime[i]*(g1[nowid]-prsum1[i-1])%MOD+MOD)%MOD;
			g2[j]=(g2[j]-prime[i]*prime[i]%MOD*(g2[nowid]-prsum2[i-1])%MOD+MOD)%MOD;
		}
	}
	printf("%lld",(S(n,0)+1ll)%MOD);
	return 0;
}