给定积性函数 $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; 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(){ 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; 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(){ 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; int Ans=GetF(); printf("%lld\n",((Ans+1)%mod+mod)%mod); return 0; }
|