不老不死的竹林引路人
特别感谢 icyM3tra 神提供的标程 OrzOrzOrz。
速度上总时间直接从 25s（我原来的程序） 飙到了 5s，单个测试点最慢也只有 800ms！
虽然有这么快的代码，但我还是不卡常了，测试数据能卡掉非线性做法对我出题的本意就足够了。
不过我的代码习惯确实太坏/qd，谢罪谢罪。
当然继续欢迎大家来虐爆这题！
贴一下 题面。
非常抱歉这个题的模数方面设置的不是太好，需要特殊处理一下负数特别小的情况。
题目背景是瞎写的，不用过多在意。
首先我们可以将 \(f\) 的贝尔级数写出来：
\[f_p(x)=1-(p^2+1)x+p^2x^2 \]将其因式分解得到：
\[f_p(x)=(1-x)(1-p^2x) \]我们考虑对 \(f\) 杜教筛，寻找函数 \(g\)。
考虑令 \(g=\mathbf{id^2}\)，则 \(g\) 的贝尔级数即为：
\[g_p(x)=\dfrac{1}{1-p^2x} \]显然 \(f_p(x)g_p(x)=1-x=\mu_p(x)\)，说明 \(f\ast g=\mu\)。
然后 \(g\) 的前缀和可以 \(O(1)\) 求，\(f\ast g=\mu\) 的前缀和可以再用一个杜教筛来求。
至于 \(f\) 的前 \(n^{2/3}\) 个值，实际上可以线性筛求出来，不过就是分类讨论一下最小质因子的次数。
然后就做完了。
时间复杂度 \(O(Tn^{2/3})\)。
奉劝大家选用其他 \(g\) 的时候不要用需要套杜教筛求的，不然理论复杂度就不是 \(O(n^{2/3})\) 了。。。
\(f\ast g\) 可以再套一层杜教筛是因为它的求取和 \(f\) 同级，也就是说调用一次 \(S_f\) 仅仅对应调用一次 \(S_{\mu}\)。因此理论复杂度还是 \(O(n^{2/3})\)。
然后笔者自己的程序跑得奇慢，觉得原因大概是用了 map 和自己代码习惯不好，所以把最后几个点的时限开到了 \(4s\)。
当然如果你愿意可以再根号平衡一下，把 \(n^c\) 适当调大来加快，不过我懒得去这么做了。
Wky 大师说这个可以 Min_25 筛，理论复杂度会更优秀，但我不会。。。
似乎因为构造的方式过于简单导致大师似乎可以一眼看出这个 \(g=\mathbf{id^2}\)。。。
那还要个鬼的贝尔级数。。。
但像我这样的凡人还是老老实实想一些比较通用的方法吧。。。
代码（来自 icyM3tra 的 std）：
#include<stdio.h>
#include<bitset>
#include<math.h>
#define rep(i,l,r) for (int i=l; i<=r; ++i)
const int N=1e7+5,M=205,K=666667,sq=45000;
const double eps=1e-14;
int T,P,mx,n,k,i6,t,sG,sF,mG[M],mF[M];
double I[sq+5]; std::bitset<N>v;
int p[K],mu[N],f[N],F[N],G[N];
inline int add(int x,int y) { return (x+=y)<P?x:x-P; }
inline int mul(int x,int y) { return 1ll*x*y%P; }
inline int power(int x,int y) {
int s=1;
while (y) {
if (y&1) s=mul(s,x);
x=mul(x,x),y>>=1;
}
return s;
}
int main() {
scanf("%d%d",&T,&P),i6=power(6,P-2);
const int MX=(T<=10?10000000:200001);
F[1]=f[1]=1,G[1]=mu[1]=1;
rep(i,1,sq) I[i]=1./i+eps;
rep(i,2,MX) {
if (!v[i]) p[++k]=i,mu[i]=-1,f[i]=P-mul(i,i)-1;
for (int j=1; j<=k&&(t=p[j]*i)<=MX; ++j) {
v[t]=1;
if (i%p[j]==0) {
if (i/p[j]%p[j]) f[t]=mul(f[i/p[j]],p[j]*p[j]);
break;
}
mu[t]=-mu[i],f[t]=mul(f[i],f[p[j]]);
}
G[i]=add(G[i-1],mu[i]+P);
F[i]=add(F[i-1],f[i]);
}
const auto dv=[](int x,int y) { return int(x*I[y]); };
const auto SS=[](int x) { return mul(mul(x,x+1),(2ull*x+1)*i6%P); };
const auto qryG=[](int x) { return x>mx?G[n/x]:mG[x]; };
const auto qryF=[](int x) { return x>mx?F[n/x]:mF[x]; };
while (T--) {
scanf("%d",&n),mx=n/MX;
for (int i=mx; i; --i) {
const int m=dv(n,i);
const int mid=sqrt(m);
sG=sF=0;
rep(j,2,mid) {
sG=add(sG,qryG(i*j));
sF=(1ll*j*j*qryF(i*j)+sF)%P;
}
int t=dv(m,mid),tx=SS(mid),ty;
for (int l=mid+1,r; l<=m; tx=ty,l=r+1) {
ty=SS(r=dv(m,--t));
sG=(1ll*(r-l+1)*qryG(dv(n,t))+sG)%P;
sF=(1ll*(ty-tx+P)*qryF(dv(n,t))+sF)%P;
}
mG[i]=add(1,P-sG);
mF[i]=add(mG[i],P-sF);
}
printf("%d\n",qryF(1));
}
return 0;
}
代码（笔者原本的代码）：
#include<iostream>
#include<cstdio>
#include<map>
#define ll __int128
using namespace std;
namespace Ehnaev{
inline ll read() {
ll ret=0,f=1;char ch=getchar();
while(ch<48||ch>57) {if(ch==45) f=-f;ch=getchar();}
while(ch>=48&&ch<=57) {ret=(ret<<3)+(ret<<1)+ch-48;ch=getchar();}
return ret*f;
}
inline void write(ll x) {
static char buf[22];static ll len=-1;
if(x>=0) {do{buf[++len]=x%10+48;x/=10;}while(x);}
else {putchar(45);do{buf[++len]=-(x%10)+48;x/=10;}while(x);}
while(len>=0) putchar(buf[len--]);
}
}using Ehnaev::read;using Ehnaev::write;
inline void writeln(ll x) {write(x);putchar(10);}
inline ll Fm(ll x,ll y) {return x>=0?x%y:(x+(-x/y+1)*y)%y;}
const ll M=5e6;
ll T,mo,n,cnt;
ll prime[M+5],f[M+5],mu[M+5],sf1[M+5],smu1[M+5];
bool ff[M+5];
map<ll,ll> sf2,smu2;
inline ll S2(ll x) {return x*(x+1)*(2*x+1)/6;}
inline ll Smu(ll x) {
if(x<=M) return smu1[x];
if(smu2.find(x)!=smu2.end()) return smu2[x];
ll r=1;
for(ll i=2,j;i<=x;i=j+1) {
j=x/(x/i);r=Fm(r-(j-i+1)*Smu(x/i),mo);
}
return smu2[x]=r;
}
inline ll Sf(ll x) {
if(x<=M) return sf1[x];
if(sf2.find(x)!=sf2.end()) return sf2[x];
ll r=Smu(x);
for(ll i=2,j;i<=x;i=j+1) {
j=x/(x/i);r=Fm(r-(S2(j)-S2(i-1))*Sf(x/i),mo);
}
return sf2[x]=r;
}
inline void Init() {
ff[1]=1;f[1]=1;mu[1]=1;
for(ll i=2;i<=M;i++) {
if(!ff[i]) {prime[++cnt]=i;f[i]=Fm(-i*i-1,mo);mu[i]=Fm(-1,mo);}
for(ll j=1;j<=cnt&&i*prime[j]<=M;j++) {
ff[i*prime[j]]=1;
if(i%prime[j]==0) {
mu[i*prime[j]]=0;
ll tmp=i/prime[j];
if(tmp%prime[j]==0) f[i*prime[j]]=0;
else f[i*prime[j]]=Fm(prime[j]*prime[j]*f[tmp],mo);
break;
}
f[i*prime[j]]=Fm(f[prime[j]]*f[i],mo);mu[i*prime[j]]=Fm(-mu[i],mo);
}
}
for(ll i=1;i<=M;i++) sf1[i]=Fm(sf1[i-1]+f[i],mo);
for(ll i=1;i<=M;i++) smu1[i]=Fm(smu1[i-1]+mu[i],mo);
}
int main() {
T=read();mo=read();Init();
while(T--) {
n=read();
writeln(Sf(n));
}
return 0;
}