一些约定
- \(p_k\) 表示第 \(k\) 个质数,特别地,\(p_0=1\);
- 令 \(n/k\) 为 \(\left \lfloor \frac{n}{k}\right \rfloor\);
- 令 \(\text{lpf}(x)\) 为数 \(x\) 的最小质因数。
要干啥?
求 积性函数 \(F(i)\) 的前缀和。
推柿子
令 \(h(i,j)=\sum_{x=2}^i[\text{lpf}(x)> p_j]F(x)\),那么答案就是 \(h(n,0)+F(1)=h(n,0)+1\)(积性函数的第一项一定为 \(1\),\(1\) 和任何数互质嘛)。求解的普遍思路似乎是递推,不妨先按照定义展开(注意,本次展开只包含可被分解成多个质数的合数):
\[h(i,j)=\sum_{k=j+1}\sum_{e=1} F(p_k^e)\cdot h(i/p_{k}^e,k) \]其实就是枚举最小质因子及其幂次,然后进行递归(易得它们是互质的)。
另外就还剩下质数、单个质数的幂的情况,后者也是比较容易写的:
\[h(i,j)=\sum_{k=j+1}\sum_{e=1} F(p_k^e)\cdot h(i/p_{k}^e,k)+F(p_k^{e+1})\ \ \ \ (p_{k}^{e+1}\le i) \]需要注意后面的限制,可以发现,当 \(p_k^{e+1}>i\) 时,\(F(p_k^e)\cdot h(i/p_{k}^e,k)\) 一定也为零,所以这个限制是紧凑且合理的。
令 \(P(n)\) 为小于等于 \(n\) 的质数的 \(F\) 值的前缀和,那么最终就有(注意到质数只需要枚举到 \(\sqrt n\) 级别):
\[h(i,j)=P(i)-P(p_j)+\sum_{k=j+1}\sum_{e=1} F(p_k^e)\cdot h(i/p_{k}^e,k)+F(p_k^{e+1})\ \ \ \ (p_{k}^{e+1}\le i) \]于是我们只用计算 \(P(n)\).
首先,我们需要 \(F(n)\) 是一个 次数较小 的多项式,这样拆开的每一项就可以视作一个 完全积性函数 —— 更一般地,不妨设 \(f_k(n)=n^k\),只需要计算这个函数的 \(P_k(n)\),乘上系数再相加即可得到 \(P(n)\).
令 \(g(i,j)\) 为在 \([2,i]\) 之间且 \(\text{lpf}>p_j\) 的数字的 \(f\) 值之和,特别地,这一坨恒包含小于等于 \(i\) 的质数的 \(f\) 值。
考虑递推地求取 \(g(i,j)\) —— 我们可以将 \(g(i,j)\) 理解为在 \(g(i,j-1)\) 的基础上,去掉 最小质因子 为 \(p_j\) 的 合数 的 \(f\) 值:
\[g(i,j)=g(i,j-1)-f_k(p_j)\cdot\left (g(i/p_j,j-1)-\sum_{q=1}^{j-1}f_k(p_q)\right)\ \ \ \ \ (i\ge p_j^2) \]由于这是完全积性函数,那么有 \(f(r\cdot p_j)=f(r)\cdot f(p_j)\),后面的一大坨相当于计算合法的 \(r\). 可以发现完全积性函数的性质省去了枚举质数幂次的操作。
注意到后面 \(i\ge p_j^2\) 的条件。实际上,若 \(i<p_j^2\),那么就不可能找到合法的 \(r\),所以直接 \(g(i,j)\leftarrow g(i,j-1)\) 即可。
容易发现,\(P_k(n)=g(n,cnt)\)(\(cnt\) 是筛出质数的个数)。
关于复杂度:\(\text{It is said that}\) 计算 \(h(i,j)\) 的复杂度为 \(\mathcal O(n^{1-\epsilon})\),计算 \(g(i,j)\) 的复杂度为 \(\mathcal O\left(\frac{n^{\frac{3}{4}}}{\log n}\right)\). 证明不会,可能会一直咕下去。
代码实现
参照某谷 模板题。
在实现过程中,函数的 \(i,j\) 参数的 个数 都是小于 \(\sqrt n\) 的。\(j\) 的原因就不阐述,对于 \(i\),可以观察到它是一个类似 \(\left \lfloor \frac{\left \lfloor \frac{n}{a}\right \rfloor}{b}\right \rfloor\) 的向下嵌套的结构,可以证明它的取值属于 \(\left \lfloor \frac{n}{a}\right \rfloor\) 的取值集合,这样不仅可以证明它的个数在 \(\sqrt n\) 之内,而且还可以进行整除分块预处理。以下给出感性证明:不妨将 \(\left \lfloor \frac{n}{a}\right \rfloor\) 中的分母替换成 \(a'\) 使得 \(a'\mid n\) 且 \(\frac{n}{a'}=\left \lfloor \frac{n}{a}\right \rfloor\),这样嵌套就可以换成 \(\left \lfloor \frac{n}{a'b}\right \rfloor\),结论得证。
于是我们将合理的 \(i\) 存起来,需要注意的是,\(i\) 的取值可能很大,所以对于 \(i\le \sqrt n\) 的部分用 \(i\) 来查找,另一部分用 \(n/i\) 来查找(这个玩意实际上就是整除分块时有 \(i=n/l\),对应的 \(r\) 就是 \(n/i\))。
回到此题,可知 \(F(n)=n^2-n\),可以分成 \(f_1,f_2\) 来计算,然后合并即可。
#include <cstdio>
#define print(x,y) write(x),putchar(y)
template <class T>
inline T read(const T sample) {
T x=0; char s; bool f=0;
while((s=getchar())>'9' || s<'0')
f |= (s=='-');
while(s>='0' && s<='9')
x = (x<<1)+(x<<3)+(s^48),
s = getchar();
return f?-x:x;
}
template <class T>
inline void write(T x) {
static int writ[50],tp=0;
if(x<0) putchar('-'),x=-x;
do writ[++tp] = x-x/10*10, x/=10; while(x);
while(tp) putchar(writ[tp--]^48);
}
#include <cmath>
#include <iostream>
using namespace std;
typedef long long ll;
const int maxn = 1e6+5;
const int mod = 1e9+7;
const int inv = 166666668;
bool is[maxn];
ll n,id[maxn];
int S,f[2][maxn],pc,p[maxn];
int g[2][maxn],tot,ind[2][maxn];
inline int inc(int x,int y) {
return x+y>=mod?x+y-mod:x+y;
}
inline int dec(int x,int y) {
return x-y<0?x-y+mod:x-y;
}
void sieve(int n) {
for(int i=2;i<=n;++i) {
if(!is[i]) {
p[++pc] = i;
f[0][pc] = inc(f[0][pc-1],i);
f[1][pc] = inc(f[1][pc-1],1ll*i*i%mod);
}
for(int j=1; j<=pc && i*p[j]<=n; ++j) {
is[i*p[j]] = 1;
if(i%p[j]==0) break;
}
}
}
inline int getSum(ll n,bool d) {
int ret; n %= mod; // qwq
if(!d) ret = (n*(n+1)>>1)%mod;
else ret = n*(n+1)%mod*(n*2+1)%mod*inv%mod;
return (ret-1<0?ret-1+mod:ret-1);
}
void getf() {
ll l,r,x;
for(l=1;l<=n;l=r+1) {
r = min(n,n/(n/l));
id[++tot] = (x=n/l);
if(x<=S) ind[0][x]=tot;
else ind[1][n/x]=tot; // for better restoration
g[0][tot] = getSum(x,0);
g[1][tot] = getSum(x,1);
}
// 用滚动数组递推 g(i,j)
for(int i=1;i<=pc;++i) {
for(int j=1; j<=tot && 1ll*p[i]*p[i]<=id[j]; ++j) {
int pre = (id[j]/p[i]<=S)?
ind[0][id[j]/p[i]]:ind[1][n/(id[j]/p[i])];
g[0][j] = dec(g[0][j],1ll*p[i]*dec(g[0][pre],f[0][i-1])%mod);
g[1][j] = dec(g[1][j],1ll*p[i]*p[i]%mod*dec(g[1][pre],f[1][i-1])%mod);
}
}
}
int h(ll x,int y) {
if(x<=p[y]) return 0;
int ID = (x<=S)?ind[0][x]:ind[1][n/x];
int ans = dec(dec(g[1][ID],g[0][ID]),dec(f[1][y],f[0][y]));
for(int i=y+1; i<=pc && 1ll*p[i]*p[i]<=x; ++i) {
for(ll tmp=p[i];;tmp=tmp*p[i]) {
if(tmp*p[i]>x) break;
ans = inc(ans,inc(
tmp*(tmp-1)%mod*h(x/tmp,i)%mod,
tmp*p[i]%mod*((tmp*p[i]-1)%mod)%mod
));
}
}
return ans;
}
int main() {
n=read(9ll);
sieve(S=sqrt(n)); getf();
print(inc(h(n,0),1),'\n');
return 0;
}
题目小结
然鹅现在只做过一道……