「算法笔记」Min_25 筛
修改于 2023.1.26。
一、Min_25 筛
对于积性函数 \(f(x)\)，\(f(p^e)\) 能快速求，\(f(p)\) 可以表示成完全积性函数的线性组合，求 \(\sum_{i=1}^n f(i)\)。
\(n\leq 10^{10}\)。
-
找到 完全 积性函数 \(f'\) 使得在质数处的取值 \(f'(p)=f(p)\)。
一般 \(f'(p)\) 形如 \(p^c\)。
区别：比如要提出 \(p\)，对 \(f\) 需要提出整个 \(p^e\) 再 \(f(p^e)f(\frac i{p^e})\)，对 \(f'\) 只需提出一个 \(p\) 再 \(f'(p)f'(\frac i p)\)。利用 \(f'\) 能方便求质数处取值和，合数处取值和再交给 \(f\)，而合数的最小质因子 \(\leq \sqrt n\) 所以枚举起来很快。
-
将 \(1\)、质数、合数分别处理。
\[\sum_{i=1}^n f(i)=f(1)+\sum_{p\leq n}f(p)+\sum_{p,e}f(p^e)(\sum_{i=1}^{n/p^e}[minp_i>p]f(i)) \]最后一部分是枚举合数的最小质因子及其次数，这样就利用积性函数的性质拆出最小质因子来递归子问题。
-
处理 \(\sum_{p\leq n}f(p)=\sum_{p\leq n}f'(p)\)：
设 \(g(n,j)\) 表示运行 \(j\) 次埃氏筛后，\(\leq n\) 仍未被筛掉的数（注意不考虑 \(1\)）的 \(f'\) 之和，即 \(\sum_{i=2}^n[i\in \text{prime} \lor minp_i>p_j]f'(i)\)。用 \(\leq\sqrt n\) 以内的质数们 \(P\) 筛后合数就被筛完了，那么 \(\sum_{p\leq n}f'(p)=g(n,|P|)\)。
边界 \(g(n,0)=\sum_{i=2}^n f'(i)\)。转移考虑 \(p_j\) 筛掉的，应该是 \(p_j\times\) 一个 \(minp\geq p_j\) 的数。注意 \(g(\frac{n}{p_j},j-1)\) 中包含的 \(<p_j\) 的质数未被筛掉。
\[g(n,j)= \begin{cases} g(n,j-1)&(p_j^2>n)\\ g(n,j-1)-f'(p_j)(g(\frac{n}{p_j},j-1)-\sum_{i=1}^{j-1}f'(p_i)))&(p_j^2\leq n) \end{cases} \]第一维很大：只需存 \(\frac n x\) 这些有用的，那么要么 \(x\leq \sqrt n\) 要么 \(\frac n x<\sqrt n\)，开两个数组分给编号。
-
递归处理整个：
设 \(S(n,j)=\sum_{i=2}^n [minp_i>p_j]f(i)\)。
\[S(n,j)=(\sum_{p\leq n}f'(p)-\sum_{i=1}^j f'(p_i))+(\sum_{i>j,p_i^e\leq n} f(p_i^e)\times (S(\frac{n}{p_i^e},i)+[e>1])) \]因为 \(S\) 里不算 \(1\)，所以若 \(e>1\) 那么 \(f(p_i^e)\) 没算进去，要额外加上。
\(ans=S(n,0)+f(1)=S(n,0)+1\)。
-
时间复杂度 \(\mathcal O(\frac{n^{0.75}}{\log n})\)。
二、例题
栗子：（注意常数除了 \(1,0\) 都不是完全积性函数）
- 筛 \(\mu\)：\(\mu(p^e)=(-1)^e\)。\(\mu(p)=-1\)，筛质数处取值时取 \(f'=1\) 最后再 \(\times (-1)\)。
- 筛 \(\varphi\)：\(\varphi(p^e)=(p-1)p^{e-1}\)。\(\varphi(p)=p-1\)，取 \(f'=\text{Id}\) 与 \(f'=1\) 相减。
- 筛 \(d^k\)：\(d(p^e)^k=(e+1)^k\)。\(d(p)^k=2^k\)（常数），取 \(f'=1\) 最后再 \(\times 2^k\)。
有些时候不会简单告诉你要求积性函数的前缀和，需要自己去推式子/找性质来发现。
技巧：
-
非严格次大质因子产生贡献：\(f(i)=val(smxp_i)\)，并规定 \(smxp_1=smxp_p=0\)。求 \(\sum_{i=1}^n f(i)\)。
对于 \(p_i^e x\)（其中 \(p_i\) 为 \(p_i^ex\) 的最小质因子），分别考虑 \(x\) 是合数、\(x\) 是质数、\(x\) 是 \(1\) 时的次大质因子。
\[S(n,j)=\sum_{i>j,p_i^e\leq n}(S(\frac{n}{p_i^e},i)+val(p_i)\times (\sum_{p_i^ex\leq n,x>p_i}[x\in\text{prime}]+[e>1])) \]（由于 \(S\) 里质数处无贡献，所以不用管 \(S(\frac{n}{p_i^e},i)\) 中只算合数的限制）
区间质数个数用 Min_25 筛。\(\sum_{p_i^ex\leq n,x>p_i}[x\in\text{prime}]=\max(g(\frac{n}{p_i^e},|{P|})-i,0)\)。
-
同上，只不过 \(smxp_1=0,smxp_p=1\)：\(S\) 式子不变，最后额外算是质数贡献。
-
分别求模 \(4\) 余 \(1\) 和模 \(4\) 余 \(3\) 的质数个数。
设 \(g_1(n,j),g_3(n,j)\) 分别表示运行 \(j\) 次埃氏筛后，\(\leq n\) 仍未被筛掉的模 \(4\) 余 \(1\) / 模 \(4\) 余 \(3\) 的数的个数。\(g_1(n,0)=\lfloor\frac{n-1}{4}\rfloor,g_3(n,0)=\lfloor\frac{n+1}{4}\rfloor\)。对于 \(g_r\)，每次在剩余模 \(4\) 余 \(r\) 中继续筛。
转移考虑 \(p_j\) 新筛掉的：对于 \(g_1(\frac{n}{p_j},j-1)\) 或 \(g_3(\frac{n}{p_j},j-1)\) 中留下的数 \(a\)，当 \(a\times p_j\equiv 1\pmod 4\) 时会对 \(g_1(n,j)\) 产生影响，\(\equiv 3\) 时对 \(g_3(n,j)\) 产生影响。
于是分类讨论：若 \(p_j\equiv 1\pmod 4\)，\(g_1(n,j)\gets g_1(n,j-1)-(g_1(\frac{n}{p_j},j-1)-sp_1(j-1))\)，\(g_3\gets g_3\) 同理；若 \(p_j\equiv 3\pmod 4\)，\(g_1(n,j)\gets g_1(n,j-1)-(g_3(\frac{n}{p_j},j-1)-sp_3(j-1))\)，\(g_3\gets g_1\) 同理。
1. LOJ#6235. 区间素数个数
2021.2.28
求 \(1\sim n\) 之间的质数个数。
\(2\leq n\leq 10^{11}\)。
\(f'=1\)。
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=7e5+5;
ll n,x,val[N],g[N];
int s,vis[N],cnt,p[N],tot,id1[N],id2[N];
int id(ll x){return x<=s?id1[x]:id2[n/x];}
signed main(){
scanf("%lld",&n),s=sqrt(n);
for(int i=2;i<=s;i++){
if(!vis[i]) p[++cnt]=i;
for(int j=1;j<=cnt&&i*p[j]<=s;j++){
vis[i*p[j]]=1;
if(i%p[j]==0) break;
}
}
for(ll l=1,r;l<=n;l=r+1){
r=n/(n/l),val[++tot]=x=n/l;
(x<=s?id1[x]:id2[n/x])=tot,g[tot]=x-1;
}
for(int j=1;j<=cnt;j++)
for(int i=1;i<=tot&&1ll*p[j]*p[j]<=val[i];i++) //某 sb 把 val[i] 写成 i 了
g[i]-=g[id(val[i]/p[j])]-(j-1);
printf("%lld\n",g[id(n)]);
return 0;
}
2. P5325 【模板】Min_25筛
2021.3.1
积性函数 \(f(x)\) 满足 \(f(p^k)=p^k(p^k-1)\)，求 \(\sum_{i=1}^n f(i)\pmod{10^9+7}\)。
\(1\leq n\leq 10^{10}\)。
\(f(p)=p(p-1)=p^2-p\)，分别用 \(f'(p)=p^2\) 和 \(f'(p)=p\) 算质数处的取值相减。
\(1^2+2^2+3^2+\cdots+n^2=\frac{n(n+1)(2n+1)}{6}\)。
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=2e5+5,mod=1e9+7;
ll n,x,val[N];
int s,vis[N],cnt,p[N],tot,id1[N],id2[N],sp1[N],sp2[N],g1[N],g2[N],iv2=(mod+1)/2,iv6=(mod+1)/6;
int id(ll x){return x<=s?id1[x]:id2[n/x];}
int f(ll x){return x%=mod,x*(x-1)%mod;}
int S(ll x,int y){
if(x<=p[y]) return 0;
int ans=((g1[id(x)]-sp1[y])%mod-(g2[id(x)]-sp2[y])%mod)%mod;
for(int i=y+1;i<=cnt&&1ll*p[i]*p[i]<=x;i++){
ll pw=1;
for(int e=1;pw*p[i]<=x;e++)
pw*=p[i],ans=(ans+1ll*f(pw)*(S(x/pw,i)+(e>1))%mod)%mod;
}
return ans;
}
signed main(){
scanf("%lld",&n),s=sqrt(n);
for(int i=2;i<=s;i++){
if(!vis[i])
p[++cnt]=i,sp1[cnt]=(sp1[cnt-1]+1ll*i*i%mod)%mod,
sp2[cnt]=(sp2[cnt-1]+i)%mod;
for(int j=1;j<=cnt&&i*p[j]<=s;j++){
vis[i*p[j]]=1;
if(i%p[j]==0) break;
}
}
for(ll l=1,r;l<=n;l=r+1){
r=n/(n/l),val[++tot]=x=n/l;
(x<=s?id1[x]:id2[n/x])=tot;
x%=mod,g1[tot]=x*(x+1)%mod*(2*x+1)%mod*iv6%mod-1,g2[tot]=(1+x)*x%mod*iv2%mod-1; //注意要 x%=mod
}
for(int j=1;j<=cnt;j++)
for(int i=1;i<=tot&&1ll*p[j]*p[j]<=val[i];i++)
g1[i]=(g1[i]-1ll*p[j]*p[j]%mod*(g1[id(val[i]/p[j])]-sp1[j-1])%mod)%mod,
g2[i]=(g2[i]-1ll*p[j]*(g2[id(val[i]/p[j])]-sp2[j-1])%mod)%mod;
printf("%d\n",(S(n,0)+1+mod)%mod);
return 0;
}
3. LOJ#6053. 简单的函数
2021.3.1
积性函数 \(f(x)\) 满足 \(f(p^c)=p\oplus c\)，求 \(\sum_{i=1}^n f(i)\pmod{10^9+7}\)。
\(1\leq n\leq 10^{10}\)。
对于 \(2\) 之外的质数都满足 \(f(p)=p-1\)，拆成 \(f'(p)=p\) 与 \(f'(p)=1\) 相减。
先都按 \(f(p)=p-1\) 算，最后若 \(n\geq 2\) 则将答案加上 \((2\oplus 1)-(2-1)=2\) 即可。
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=2e5+5,mod=1e9+7;
ll n,x,val[N];
int s,vis[N],cnt,p[N],tot,id1[N],id2[N],sp1[N],sp2[N],g1[N],g2[N],iv2=(mod+1)/2;
int id(ll x){return x<=s?id1[x]:id2[n/x];}
int S(ll x,int y){
if(x<=p[y]) return 0;
int ans=((g1[id(x)]-sp1[y])%mod-(g2[id(x)]-sp2[y])%mod)%mod;
for(int i=y+1;i<=cnt&&1ll*p[i]*p[i]<=x;i++){
ll pw=1;
for(int e=1;pw*p[i]<=x;e++)
pw*=p[i],ans=(ans+1ll*(p[i]^e)*(S(x/pw,i)+(e>1))%mod)%mod;
}
return ans;
}
signed main(){
scanf("%lld",&n),s=sqrt(n);
for(int i=2;i<=s;i++){
if(!vis[i])
p[++cnt]=i,sp1[cnt]=(sp1[cnt-1]+i)%mod,sp2[cnt]=cnt;
for(int j=1;j<=cnt&&i*p[j]<=s;j++){
vis[i*p[j]]=1;
if(i%p[j]==0) break;
}
}
for(ll l=1,r;l<=n;l=r+1){
r=n/(n/l),val[++tot]=x=n/l;
(x<=s?id1[x]:id2[n/x])=tot;
x%=mod,g1[tot]=(1+x)*x%mod*iv2%mod-1,g2[tot]=x-1;
}
for(int j=1;j<=cnt;j++)
for(int i=1;i<=tot&&1ll*p[j]*p[j]<=val[i];i++)
g1[i]=(g1[i]-1ll*p[j]*(g1[id(val[i]/p[j])]-sp1[j-1])%mod)%mod,
g2[i]=(g2[i]-1ll*(g2[id(val[i]/p[j])]-sp2[j-1])%mod)%mod;
printf("%d\n",(S(n,0)+1+(n>=2?2:0)+mod)%mod);
return 0;
}
4. 求和
2021.4.1
给出 \(n,k\)，求：
\[\sum_{1\leq a_1,a_2,\cdots,a_k\leq n}\left\lfloor\frac{n}{\text{lcm}(a_1,a_2,\cdots,a_k)}\right\rfloor\pmod{10^9+7} \]\(1\leq n\leq 10^{10}\)，\(1\leq k\leq 10^9\)。
首先要转化一下：\(\lfloor\frac n x\rfloor=\sum_{i=1}^n [x\mid i]\)。因为整除与 \(\text{lcm}\) 比较相配。
\(d(p^e)^k=(e+1)^k\)，\(d(p)^k=2^k\)（常数），注意 \(f'=2^k\) 不是完全积性函数，需要取 \(f'=1\) 最后再 \(\times 2^k\)。
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=2e5+5,mod=1e9+7;
ll n,x,val[N];
int k,s,vis[N],cnt,p[N],tot,id1[N],id2[N],g[N],mp[N];
int id(ll x){return x<=s?id1[x]:id2[n/x];}
int qpow(int x,int n){
if(mp[x]) return mp[x]; //记忆化
int tx=x,ans=1;
for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
return mp[tx]=ans;
}
int S(ll x,int y){
if(x<=p[y]) return 0;
int ans=1ll*(g[id(x)]-y)%mod*qpow(2,k)%mod;
for(int i=y+1;i<=cnt&&1ll*p[i]*p[i]<=x;i++){
ll pw=1;
for(int e=1;pw*p[i]<=x;e++)
pw*=p[i],ans=(ans+1ll*qpow(e+1,k)*(S(x/pw,i)+(e>1))%mod)%mod;
}
return ans;
}
signed main(){
scanf("%lld%d",&n,&k),s=sqrt(n);
for(int i=2;i<=s;i++){
if(!vis[i]) p[++cnt]=i;
for(int j=1;j<=cnt&&i*p[j]<=s;j++){
vis[i*p[j]]=1;
if(i%p[j]==0) break;
}
}
for(ll l=1,r;l<=n;l=r+1){
r=n/(n/l),val[++tot]=x=n/l;
(x<=s?id1[x]:id2[n/x])=tot,g[tot]=(x-1)%mod;
}
for(int j=1;j<=cnt;j++)
for(int i=1;i<=tot&&1ll*p[j]*p[j]<=val[i];i++)
g[i]=(g[i]-(g[id(val[i]/p[j])]-(j-1))%mod)%mod;
printf("%d\n",(S(n,0)+1+mod)%mod);
return 0;
}
5. UOJ#188. 【UR #13】Sanrd
2023.1.26
设 \(f(i)\) 表示 \(i\) 的非严格次大质因子（重复的质因子算多次，比如 \(f(4)=2\)），若 \(i\) 为质数或 \(1\) 则为 \(0\)。求 \(\sum_{i=l}^r f(i)\)。
\(r\leq 10^{11}\)。
因为质数时没贡献，不用管。还是枚举合数的最小质因子 \(p_i\) 及其出现次数 \(e\)，注意到对于 \(p_i^ex\)，若 \(x\) 是合数则其次大质因子为 \(x\) 的次大质因子可以递归下去，若 \(x\) 是质数则为 \(p_i\)，若 \(x=1\) 且 \(e>1\) 则为 \(p_i\)。故有：
（由于 \(S\) 里质数处无贡献，所以不用管 \(S(\frac{n}{p_i^e},i)\) 中只算合数的限制）
区间质数个数用 Min_25 筛。\(\sum_{p_i^ex\leq n,x>p_i}[x\in\text{prime}]=\max(g(\frac{n}{p_i^e},|{P|})-i,0)\)。
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=7e5+5;
ll n,l,r,x,val[N],g[N];
int s,vis[N],cnt,p[N],tot,id1[N],id2[N];
int id(ll x){return x<=s?id1[x]:id2[n/x];}
ll S(ll x,int y){
if(x<=p[y]) return 0;
ll ans=0;
for(int i=y+1;i<=cnt&&1ll*p[i]*p[i]<=x;i++){
ll pw=1;
for(int e=1;pw*p[i]<=x;e++)
pw*=p[i],ans+=S(x/pw,i)+p[i]*(max(g[id(x/pw)]-i,0ll)+(e>1));
}
return ans;
}
ll calc(ll tn){
s=sqrt(n=tn),cnt=tot=0;
for(int i=2;i<=s;i++){
if(!vis[i]) p[++cnt]=i;
for(int j=1;j<=cnt&&i*p[j]<=s;j++){
vis[i*p[j]]=1;
if(i%p[j]==0) break;
}
}
for(ll l=1,r;l<=n;l=r+1){
r=n/(n/l),val[++tot]=x=n/l;
(x<=s?id1[x]:id2[n/x])=tot,g[tot]=x-1;
}
for(int j=1;j<=cnt;j++)
for(int i=1;i<=tot&&1ll*p[j]*p[j]<=val[i];i++)
g[i]-=g[id(val[i]/p[j])]-(j-1);
return S(n,0);
}
signed main(){
scanf("%lld%lld",&l,&r),
printf("%lld\n",calc(r)-calc(l-1));
return 0;
}
6. LOJ#572.「LibreOJ Round #11」Misaka Network 与求和
2023.1.26 莫比乌斯反演 + 杜教筛 + Min_25 筛
给出 \(n,k\)，设 \(f(i)\) 表示 \(i\) 的非严格次大质因子的 \(k\) 次方，并规定 \(f(1)=0,f(p)=1\)，求：
\[\sum_{i=1}^n\sum_{j=1}^n f(\gcd(i,j))\pmod{2^{32}} \]\(1\leq n,k\leq 2\times 10^9\)。
先套路莫反：
外面可以套整除分块，里面需要算 \(sum(T)=\sum_{d\mid T}f(d)\mu(\frac T d)=f*\mu\) 的前缀和。
杜教筛，\(1\) 和 \(sum*1\) 的前缀和都能求，其中 \(sum*1=f*\varepsilon=f\)，\(f\) 的前缀和类似 UOJ#188 Min_25 筛。
#include<bits/stdc++.h>
#define ui unsigned int
using namespace std;
const int N=1e5+5;
int n,k,s,x,vis[N],cnt,p[N],tot,val[N],id1[N],id2[N],g[N];
ui pk[N],vS[N],vD[N],ans;
int id(int x){return x<=s?id1[x]:id2[n/x];}
ui qpow(ui x,int n){
ui ans=1;
for(;n;n>>=1,x*=x) if(n&1) ans*=x;
return ans;
}
ui S(int x,int y){
if(!y&&vS[id(x)]) return vS[id(x)]; //记忆化
if(x<=p[y]) return 0;
ui ans=0;
for(int i=y+1;i<=cnt&&1ll*p[i]*p[i]<=x;i++){
int pw=1;
for(int e=1;1ll*pw*p[i]<=x;e++)
pw*=p[i],ans+=S(x/pw,i)+pk[i]*(max(g[id(x/pw)]-i,0)+(e>1));
}
if(!y) vS[id(x)]=ans;
return ans;
}
ui D(int n){ //杜教筛
if(vD[id(n)]) return vD[id(n)];
ui ans=0;
for(int l=2,r;l<=n;l=r+1)
r=n/(n/l),ans+=(r-l+1)*D(n/l);
return vD[id(n)]=(S(n,0)+g[id(n)])-ans; //+g 处理 f(p)=1：额外算上质数贡献
}
signed main(){
scanf("%d%d",&n,&k),s=sqrt(n);
for(int i=2;i<=s;i++){
if(!vis[i]) p[++cnt]=i,pk[cnt]=qpow(i,k);
for(int j=1;j<=cnt&&i*p[j]<=s;j++){
vis[i*p[j]]=1;
if(i%p[j]==0) break;
}
}
for(int l=1,r;l<=n;l=r+1){
r=n/(n/l),val[++tot]=x=n/l;
(x<=s?id1[x]:id2[n/x])=tot,g[tot]=x-1;
}
for(int j=1;j<=cnt;j++)
for(int i=1;i<=tot&&1ll*p[j]*p[j]<=val[i];i++)
g[i]-=g[id(val[i]/p[j])]-(j-1);
for(int l=1,r;l<=n;l=r+1)
r=n/(n/l),ans+=1u*(n/l)*(n/l)*(D(r)-D(l-1));
printf("%u\n",ans);
return 0;
}
7. 数
2023.1.26 找性质
给出 \(n\)，求有多少 \(k\in[1,n]\) 满足：
- 使得 \(a^2+b^2+k=c^2\) 成立且 \(a,b,c\) 为等差数列的正整数 \(a,b,c\) 有唯一解。
\(1\leq n\leq 10^{11}\)。
由于“等差数列”，\(c=2b-a\)，代入得 \(k=(2b-a)^2-a^2-b^2=b(3b-4a)\)。
考虑 \(k\) 的一个分解 \(k=xy\)，对应一组解 \(b=x,a=\frac{3x-y}{4}\)。由于“正整数”，故该解合法等价于 \(3x>y\) 且 \(3x\equiv y\pmod 4\)，而 \(3x\equiv -x\pmod 4\)，后面一个条件可以改写为 \(x+y\equiv 0\pmod 4\)。我们只关心是否只有一个 \(k=xy\) 满足 \(3x>y\) 且 \(x+y\equiv 0\pmod 4\)。
设 \(k=u\cdot 2^v\,(2\not\mid u )\)，发现当 \(v=0\) 时，只有 \(k\equiv 3\pmod 4\) 才存在解，并且 \(k\) 是质数才存在唯一解（否则有更多分配方式。显然 \(k\equiv 3\pmod 4\) 分解起来模 \(4\) 肯定是一个 \(1\) 一个 \(3\)）；当 \(v>0\) 时，只有 \(v=2/4\) 且 \(u\) 为质数时才有唯一解。
综上，合法的 \(k\) 满足：\(k\) 为模 \(4\) 余 \(3\) 的质数；\(k\) 为奇质数 \(\times 4\) 或 \(k=4\)；\(k\) 为奇质数 \(\times 16\) 或 \(k=16\)。
#include<bits/stdc++.h>
#define ll long long
#define g(x) g[x][0]+g[x][1]
using namespace std;
const int N=7e5+5;
ll n,x,val[N],g[N][2];
int s,vis[N],cnt,p[N],tot,id1[N],id2[N],sp[N][2];
int id(ll x){return x<=s?id1[x]:id2[n/x];}
signed main(){
scanf("%lld",&n),s=sqrt(n);
for(int i=2;i<=s;i++){
if(!vis[i])
p[++cnt]=i,sp[cnt][0]=sp[cnt-1][0]+(i%4==1),
sp[cnt][1]=sp[cnt-1][1]+(i%4==3);
for(int j=1;j<=cnt&&i*p[j]<=s;j++){
vis[i*p[j]]=1;
if(i%p[j]==0) break;
}
}
for(ll l=1,r;l<=n;l=r+1){
r=n/(n/l),val[++tot]=x=n/l;
(x<=s?id1[x]:id2[n/x])=tot;
g[tot][0]=(x-1)/4,g[tot][1]=(x+1)/4;
}
for(int j=2;j<=cnt;j++) //2%4 既不是 1 也不是 3，要从 p=3 开始筛
for(int i=1;i<=tot&&1ll*p[j]*p[j]<=val[i];i++){
int k=id(val[i]/p[j]),o=p[j]%4==3; //考虑用来筛的质数模 4 余几
g[i][0]-=g[k][o]-sp[j-1][o];
g[i][1]-=g[k][o^1]-sp[j-1][o^1];
}
printf("%lld\n",g[id(n)][1]+(g(id(n/4))+(n>=4))+(g(id(n/16))+(n>=16)));
return 0;
}
浙公网安备 33010602011771号