U214268 不老不死的竹林引路人

不老不死的竹林引路人

特别感谢 icyM3tra 神提供的标程 OrzOrzOrz。

速度上总时间直接从 25s（我原来的程序） 飙到了 5s，单个测试点最慢也只有 800ms！（UPD: AD20230402 看到了 947ms 的标程，只能说鲍队太猛了Orzzzzzzzzzzz）

虽然有这么快的代码，但我还是不卡常了，测试数据能卡掉非线性做法对我出题的本意就足够了。

不过我的代码习惯确实太坏/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;
}