给定积性函数 $f(x)$,求 $F(n)=\sum_{i=1}^nf(i)$。

如果没有特殊说明,忽略掉所有函数在 $1$ 处的取值。

第二部分

设 $G(n)=\sum_{i=1}^n f(i)[i\in Prime]$,假设我们现在对于所有 $i$,求出了 $G(\left\lfloor\frac{n}{i}\right\rfloor)$ 的值,考虑如何计算 $F(n)$。

考虑设 $F(i,j)$ 表示 $\sum_{k=1}^j f(k)[\mathrm{min}_k\ge p_i]$,其中 $\mathrm{min}_k$ 表示 $k$ 的最小质因子,而 $p_i$ 表示第 $i$ 个质数。

我们考虑如何计算 $F(i,j)$。

我们考虑枚举这个最小的质因子,设为 $k$,显然 $k\ge i$,接下来我们考虑计数满足 $\mathrm{min}j=p_k$ 的 $j$,显然有 $p_k^2\le j$,否则的话 $j$ 一定有更小的质因子(我们先不考虑质数)。接下来要考虑的是 $j$ 里面 $p_k$ 出现了几次,不难发现答案应该是 $F(k+1,\left\lfloor\frac{j}{p_k}\right\rfloor)+F(k+1,\left\lfloor\frac{j}{p_k^2}\right\rfloor)+\cdots$。注意我们的 $F(i,j)$ 里面是不考虑 $1$ 的,所以 $p_k^e$ 都没有被考虑,我们考虑把它们加上,这样我们有:$\sum{p_{k}^{e+1}\le j,e\ge 1} f(p_k^{e+1})+F(k+1,\frac{j}{p_k^e})$。我们这里先不考虑 $f(p_k)$。不难发现我们还有一些包括 $f_{p_k}$ 在内的一些质数的取值没有计算,所以我们需要加上它们,于是我们有:

$$
F(i,j)=G(j)-\sum\limits_{k=1}^{i-1}f(p_k)+\sum\limits_{p_k^2\le j,k\ge i}\sum\limits_{p_k^{e+1}\le j,e\ge 1} f(p_k^{e+1})+F(k+1,\left\lfloor\frac{j}{p_k^e}\right\rfloor)*f(p_k^e)
$$

综上,我们可以通过爆搜来计算 $F(i,j)$,最后答案就是 $F(1,n)$。

这段的复杂度并没有被证明,但是实测跑的并不慢。复杂度只被证明是 $O(n^{1-\epsilon})$

注意,上面的部分不需要记忆化,而且只能求出 $F(n)$,而不能求出所有的 $F(\left\lfloor\frac{n}{i}\right\rfloor)$。

由 $\sum$ 的限制可以知道,我们预处理质数需要预处理 $\sqrt{n}$ 以内的质数和第一个比 $\sqrt{n}$ 大的质数,这是因为我们递归的停止条件是 $p_k^2\le j$,所以是需要最后一个质数满足 $p^2>n$ 的。

第一部分

关键在于如何求出我们的 $G$。我们找一个完全积性函数 $h$ 使得对于所有的 $p\in Prime$,有 $h(p)=f(p)$,且 $h$ 便于计算前缀和。设 $G(i,n)=\sum_{k=1}^n h(k)[k\in Prime \lor min_k> p_i]$,这里我们递推考虑 $p_i$ 筛掉了哪些数,就是那些最小质因子为 $p_i$ 的合数。则有:

$$
G(i,n)=G(i-1,n)-\sum\limits_{k=1}^n h(k)[k\notin Prime \land min_k=p_i]
$$

假设 $p_i^2\le n$

考虑 $G(i-1,\lfloor\frac{n}{p_i}\rfloor)$ 里面包含了哪些数,首先它包含 $p_1,p_2,\cdots,p_{i-1}$,其次,它包含所有最小质因子 $\ge p_i$ 的数给它乘上一个 $h(p_i)$,那么我们所有最小质因子大于等于 $p_i$ 的就变成了所有最小质因子为 $p_i$ 的数,但是前面那部分 $p_1p_i,\cdots,p_{i-1}p_{i}$ 我们是不要的,所以再把它们加回来,于是我们有:

$$
G(i,n)=G(i-1,n)-G(i-1,\lfloor\frac{n}{p_i}\rfloor)h(p_i)+(\sum_{j=1}^{i-1}h(p_j))h(p_i)
$$

那么最后的答案为 $G(n)=G(m,n)$,其中 $m$ 可以设为 $n$ 内的质数个数。

考虑这部分的复杂度,这部分实际上是一个可以进行 $O(1)$ 转移的 DP,因此只需要考虑其状态数,第一维显然是 $O(\sqrt n)$ 的,第二维需要满足 $n\ge p_i^2$,所以状态数等同于下面的式子:

$$
\sum\limits_{i=1}^{\sqrt{n}}[i\in Prime]\sum\limits_{j,\not\exist j’<j,s.t.\lfloor\frac{n}{j’}\rfloor=\lfloor\frac{n}{j}\rfloor}[\lfloor\frac{n}{j}\rfloor\ge i^2]
$$

注意第二个 $\sum$,$\lfloor\frac{n}{j}\rfloor$ 相同的 $j$ 我们只会计算一次。

其中,我们认为当 $i\le n^\frac{1}{4}$ 的时候代价为 $\sqrt{n}$,其余情况认为是 $\frac{n}{i^2}$,前面直接加起来,后面用积分积起来,而质数的限制我们直接除以 $\ln n$ 即可,这是因为 $\pi(n)=\frac{n}{n\ln n}$,等价于给质数每 $\frac{1}{\ln n}$ 出现一次,而相邻 $\frac{1}{\ln n}$ 个数以上函数取值差异有限,所以可以认为一段 $\ln n$ 中的数函数值相等,于是我们就可以这样计算:

$$
\frac{1}{\ln n}(\sum\limits_{i=1}^{n^{\frac{1}{4}}}\sqrt{n}+\int_{n^{\frac{1}{4}}}^{\sqrt n}\frac{n}{x^2}dx)
$$

前面 $\sum$ 显然是 $n^{\frac{3}{4}}$,后面是 $n(n^{-\frac{1}{4}}-n^{-\frac{1}{2}})$,可以认为是 $n^{\frac{3}{4}}$,所以总复杂度为 $\frac{n^{\frac{3}{4}}}{\ln n}$。

注意我们暂且认为 $f(1)=h(1)=0$

同时,第一部分启发我们,枚举所有数除了最大质因子以外的部分是很快的。

整个算法在数据范围是 $10^{10}$ 时消耗的时间约为 $2s$ 左右。

例题

洛谷模板题

链接

如果给定的函数 $f$ 找不到对应的完全积性函数 $h$ 满足要求的话,可以尝试把 $f$ 拆成 $a_1f_1+a_2f_2+\cdots a_nf_n$ 的形式,然后如果后者每个都可以找到对应的 $h$,然后对后面每个 $f_i$ 都做一遍第一部分,然后就可以得到 $G$。

这提示我们,$G$ 是可以拆开来计算的。

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
#include<bits/stdc++.h>
#define mset(a,b) memset((a),(b),sizeof((a)))
#define rep(i,l,r) for(int i=(l);i<=(r);i++)
#define dec(i,l,r) for(int i=(r);i>=(l);i--)
#define cmax(a,b) (((a)<(b))?(a=b):(a))
#define cmin(a,b) (((a)>(b))?(a=b):(a))
#define Next(k) for(int x=head[k];x;x=li[x].next)
#define vc vector
#define ar array
#define pi pair
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define N 1000010
#define M number
using namespace std;

typedef double dd;
typedef long double ld;
typedef long long ll;
typedef unsigned int uint;
typedef unsigned long long ull;
#define int long long
typedef pair<int,int> P;
typedef vector<int> vi;

const int INF=0x3f3f3f3f;
const dd eps=1e-9;
const int Len=100010;
const int mod=1e9+7;

template<typename T> inline void read(T &x) {
x=0; int f=1;
char c=getchar();
for(;!isdigit(c);c=getchar()) if(c == '-') f=-f;
for(;isdigit(c);c=getchar()) x=x*10+c-'0';
x*=f;
}

int pr[N],ta,h[N][2],sumh[N][2],id[N],idtot,n,g[N][2],inv6,inv2,G[N],pp[N][41],sumf[N];
bool no[N];
unordered_map<int,int> rk;

inline int ksm(int a,int b,int mod){
int res=1;while(b){if(b&1)res=1ll*res*a%mod;a=1ll*a*a%mod;b>>=1;}return res;
}
inline int inv(int x){return ksm(x,mod-2,mod);}
inline int f(int p,int k){
return pp[p][k]*(pp[p][k]-1)%mod;
}
inline void sieve(){
no[1]=1;
rep(i,2,Len){
if(!no[i]) pr[++ta]=i;
for(int j=1;j<=ta&&1ll*pr[j]*i<=Len;j++){
no[pr[j]*i]=1;if(i%pr[j]==0) break;
}
}
rep(i,1,ta){
h[i][0]=pr[i]*pr[i]%mod;sumh[i][0]=(sumh[i-1][0]+h[i][0])%mod;
h[i][1]=pr[i];sumh[i][1]=(sumh[i-1][1]+h[i][1])%mod;
}
rep(i,1,ta) pp[i][0]=1;
rep(i,1,ta){
rep(j,1,40) pp[i][j]=pp[i][j-1]*pr[i]%mod;
}
rep(i,1,ta){
sumf[i]=(sumf[i-1]+f(i,1))%mod;
}
}
inline int SumH(int n,int op){
n%=mod;
if(op==0) return n*(n+1)%mod*(2*n+1)%mod*inv6%mod;
else return n*(n+1)%mod*inv2%mod;
}
inline void Init(){
sieve();
for(int l=1,r;l<=n;l=r+1){
r=n/(n/l);id[++idtot]=n/l;
}
reverse(id+1,id+idtot+1);
rep(i,1,idtot) rk[id[i]]=i;
inv6=inv(6);inv2=inv(2);
}
inline void GetG(int op){
rep(i,1,idtot) g[i][op]=SumH(id[i],op)-1;
rep(i,1,ta){
dec(j,1,idtot){
if(pr[i]*pr[i]>id[j]) break;
// if(pr[i-1]*pr[i-1]>id[j]/pr[i]) break;
g[j][op]=((g[j][op]-g[rk[id[j]/pr[i]]][op]*h[i][op]%mod+sumh[i-1][op]*h[i][op]%mod)%mod+mod)%mod;
}
}
}
inline int dfs(int i,int j){
if(i>ta) return 0;
if(pr[i]>id[j]) return 0;
int ans=((G[j]-sumf[i-1])%mod+mod)%mod;
for(int k=i;k<=ta&&pr[k]*pr[k]<=id[j];k++){
int now=pr[k];
for(int e=1;now<=id[j]/pr[k];e++,now*=pr[k]){
ans=((ans+f(k,e+1))%mod+dfs(k+1,rk[id[j]/(now)])*f(k,e)%mod)%mod;
}
}
return ans;
}
inline int GetF(){return dfs(1,idtot);}

signed main(){
// freopen("my.in","r",stdin);
// freopen("std.out","w",stdout);
read(n);Init();GetG(0);GetG(1);
rep(i,1,idtot) G[i]=((g[i][0]-g[i][1])%mod+mod)%mod;
int Ans=GetF();
printf("%lld\n",((Ans+1)%mod+mod)%mod);
return 0;
}

犯的错误:

  • 在预处理 $G$ 的时候注意不要加上限制 $p_{i-1}^2\le \lfloor\frac{j}{p_i}\rfloor$,这是因为即使是 $p_i^2<j$ 的时候,$G(i,j)$ 也是有意义的,只不过其没必要再继续计算,直接从 $G(i-1,j)$ 传递过来即可。

求莫比乌斯函数前缀和

$$
G(n)=\sum\limits_{i=1}^n-[i\in Prime]=-\sum\limits_{i=1}^n [i\in Prime]
$$

取 $h=I$ 即可。

欧拉函数前缀和

$$
G(n)=\sum\limits_{i=1}^n(i-1)[i\in Prime]=\sum\limits_{i=1}^ni[i\in Prime]+\sum\limits_{i=1}^n 1[i\in Prime]
$$

取 $h_1(x)=x,h_2=I$ 即可。

LOJ6053

有 $f(p)=p\oplus1$,实际上除了 $p=2$ 的情况其它情况都是 $f(p)=p-1$,我们暂且认为 $f(p)=p-1$,这可以直接用上面的做法做,然后特殊处理 $2$,即把所有 $G$ 的位置加上 $2$。然后做第二部分即可。

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
int pr[N],ta,h[N][2],sumh[N][2],id[N],idtot,n,g[N][2],inv6,inv2,G[N],sumf[N];
bool no[N];
unordered_map<int,int> rk;

inline int ksm(int a,int b,int mod){
int res=1;while(b){if(b&1)res=1ll*res*a%mod;a=1ll*a*a%mod;b>>=1;}return res;
}
inline int inv(int x){return ksm(x,mod-2,mod);}
inline int f(int p,int k){return pr[p]^k;}
inline void sieve(){
no[1]=1;
rep(i,2,Len){
if(!no[i]) pr[++ta]=i;
for(int j=1;j<=ta&&1ll*pr[j]*i<=Len;j++){
no[pr[j]*i]=1;if(i%pr[j]==0) break;
}
}
rep(i,1,ta){
h[i][0]=pr[i]%mod;sumh[i][0]=(sumh[i-1][0]+h[i][0])%mod;
h[i][1]=1;sumh[i][1]=(sumh[i-1][1]+h[i][1])%mod;
}
rep(i,1,ta){
sumf[i]=(sumf[i-1]+f(i,1))%mod;
}
}
inline int SumH(int n,int op){
n%=mod;
if(op==0) return n*(n+1)%mod*inv2%mod;
else return n;
}
inline void Init(){
sieve();
for(int l=1,r;l<=n;l=r+1){
r=n/(n/l);id[++idtot]=n/l;
}
reverse(id+1,id+idtot+1);
rep(i,1,idtot) rk[id[i]]=i;
inv6=inv(6);inv2=inv(2);
}
inline void GetG(int op){
rep(i,1,idtot) g[i][op]=SumH(id[i],op)-1;
rep(i,1,ta){
dec(j,1,idtot){
if(pr[i]*pr[i]>id[j]) break;
// if(pr[i-1]*pr[i-1]>id[j]/pr[i]) break;
g[j][op]=((g[j][op]-g[rk[id[j]/pr[i]]][op]*h[i][op]%mod+sumh[i-1][op]*h[i][op]%mod)%mod+mod)%mod;
}
}
}
inline int dfs(int i,int j){
if(i>ta) return 0;
if(pr[i]>id[j]) return 0;
int ans=((G[j]-sumf[i-1])%mod+mod)%mod;
for(int k=i;k<=ta&&pr[k]*pr[k]<=id[j];k++){
int now=pr[k];
for(int e=1;now<=id[j]/pr[k];e++,now*=pr[k]){
ans=((ans+f(k,e+1))%mod+dfs(k+1,rk[id[j]/(now)])*f(k,e)%mod)%mod;
}
}
return ans;
}
inline int GetF(){return dfs(1,idtot);}

signed main(){
// freopen("my.in","r",stdin);
// freopen("std.out","w",stdout);
read(n);Init();GetG(0);GetG(1);
rep(i,1,idtot) G[i]=((g[i][0]-g[i][1])%mod+mod)%mod;
rep(i,2,idtot) G[i]=(G[i]+2)%mod;
// rep(i,1,idtot) printf("G[%d]=%d\n",i,G[i]);
int Ans=GetF();
printf("%lld\n",((Ans+1)%mod+mod)%mod);
return 0;
}