筛法

Vision271 · 2022-07-25 08:34:31 · 个人记录

概述

筛法原本是筛取质数的一种算法。在 OI 中，它被推广到了筛积性函数值。

埃氏筛

枚举每个数，筛去它们的整数倍。
复杂度为 O(n\log\log n)，证明非常困难，我们从心一下。
埃氏筛有很多 trick，例如只取 \leqslant \sqrt{n} 的质数来筛，分块筛（这个是时间换空间），etc.，但大都不实用。

欧拉筛（线性筛）

枚举所有数 i，筛去他们的质数 p_i 倍。
证明：
- 显然所有合数都可以被写成 i\times p（p 是质数）。
- 当 i\bmod p\neq 0,i\times p 在枚举至 i 时筛掉；
- 当遇到第一个 i\bmod p=0 ，显然这个 p 是 i 的最小质因子（也就是最小因子），记为 p_0 ,把 i\times p_0 筛掉然后 break；
- 对于 p>p_0 的 i\times p，i\times p=p_0\times k\times p=p_0\times (k\times p)。
- 其中 p_0<p,k\times p>i ，从而会在枚举到更大的 i 时删掉。
- 综上，每个数只会被筛掉一次，所以复杂度为线性。
这一证明有如下等价形式：将每个数写成质数-指数向量，则 n 由 \dfrac{n}{lowbit(n)} 筛出，这里 lowbit(n) 表示最低非零位的位权。
小技巧：可以通过记录每个数是被谁筛的来递归求解每个数的质因数分解。单次复杂度 O(\sum\limits_{p_i\mid x} k_i)。
欧拉筛也可以使用根号加速的优化，但通常很少用：
- 如果不根号加速会 T，那么所筛的范围根本存不下。
- 如果是在卡常...好吧，这还是有一点意义的。
- 不过事实上一般欧拉筛不是简单的筛质数，求某些积性函数的前缀和就已经 O(n) 了，别老惦记卡这点不在瓶颈上的常数了。
给出示范代码：

bool np[maxv];
int p[maxv],ptot;
il void init(){
    np[1]=1; int res;
    For(i,2,usev){
        if(!np[i]) p[++ptot]=i;
        for(int j=1;j<=ptot && i*(ll)p[j]<=usev;++j){
            res=i*p[j],np[res]=1;
            if(!(i%p[j])) break;
        }
    }
    return;
}

杜教筛

杜教筛可以在亚线性时间内求出某些特征积性函数的前缀和。
所谓某些特征积性函数，指的是存在一个对应的积性函数 g，使得 f\ast g=h，且 g,h 的前缀和可以 O(1) 得知的函数 f。
对于这种函数，我们考虑化式子，记 F,G,H 分别为 f,g,h 的前缀和：

\begin{aligned} H(n) & =\sum\limits_{i=1}^n \sum\limits_{d\mid i} f(\frac{i}{d})g(d) \\ & =\sum\limits_{d=1}^n g(d) \sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor} f(j) \\ & =\sum\limits_{d=1}^n g(d)F(\lfloor\frac{n}{d}\rfloor) \\ & =g(1)F(n)+\sum\limits_{d=2}^n g(d)F(\lfloor\frac{n}{d}\rfloor) \end{aligned}

我们这里用到的数论函数肯定都是 f(1)=1 的啦（积性嘛），所以事实上上式代表着

F(n)=H(n)-\sum\limits_{d=2}^n g(d)F(\lfloor\frac{n}{d}\rfloor)

注意到式右的 F 值我们是已知的（如果我们从小到大来筛的话），G,H 可以 O(1) 得知，则对 \lfloor\frac{n}{d}\rfloor 整除分块，单轮复杂度 O(\sqrt{n})。
当然事实上我们只关心一个点的 F(n)，故该式一般递归求解。对复杂度积分，发现在线筛预处理 n^\frac{2}{3} 内的 f 值时取得最优复杂度为 \Theta(n^\frac{2}{3})。

对 \mu

注意到 \mu\ast I=\epsilon，且 I,\epsilon 的前缀和都可以 O(1) 得知，于是是板子。

对 \varphi

注意到 \varphi\ast I=id，且 I,id 的前缀和都可以 O(1) 得知，于是还是板子。
不过，事实上我们有更好的方法。将 \varphi(i) 展开，得到如下式子：

\begin{aligned} \sum\limits_{i=1}^n \varphi(i) & =\sum\limits_{i=1}^n \sum\limits_{j=1}^i [gcd(i,j)=1] \\ & =\sum\limits_{i=1}^n \sum\limits_{j=1}^i \sum\limits_{d\mid i \And d\mid j} \mu(d) \\ & =\sum\limits_{d=1}^n \mu(d) \sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor} i \\ & =\sum\limits_{d=1}^n \mu(d) \dfrac{\lfloor\frac{n}{d}\rfloor(\lfloor\frac{n}{d}\rfloor+1)}{2} \end{aligned}

那么我们对这个东西整除分块，直接杜来查 \mu 的前缀和，就可以了。当然，常数未必更小，无非是杜 \mu 和杜 \varphi 的区别。另外也可以把这个式子拆一拆，把 2,1 拆出来什么的。
给出对这两个函数的示范代码（这里为了卡常使用了自定义哈希函数）：

namespace du{
    bool np[maxn23];
    int p[maxn23],ptot;
    int mu[maxn23],smu[maxn23];
    il void init(){
        np[1]=1,mu[1]=smu[1]=1; int res;
        For(i,2,usen23){
            if(!np[i]) p[++ptot]=i,mu[i]=-1;
            smu[i]=smu[i-1]+mu[i];
            for(int j=1;j<=ptot && i*(ll)p[j]<=usen23ll;++j){
                res=i*p[j],np[res]=1;
                if(i%p[j]) mu[res]=-mu[i];
                else break;
            }
        }
        return;
    }

    struct my_hash{
        static ull splitmix64(ull x){
            x+=0x9e3779b97f4a7c15;
            x=(x^(x>>30))*0xbf58476d1ce4e5b9;
            x=(x^(x>>27))*0x94d049bb133111eb;
            return x^(x>>31);
        }

        size_t operator()(ull x) const {
            static const ull FIXED_RANDOM = chrono::steady_clock::now().time_since_epoch().count();
            return splitmix64(x + FIXED_RANDOM);
        }
    };
    u_map<int,int,my_hash> smu2;
    il int getmu(int n){
        if(n<=usen23) return smu[n];
        if(smu2.find(n)!=smu2.end()) return smu2[n];
        int ret=1,res;
        for(int l=2,r;;l=r+1){
            res=n/l,r=n/res;
            ret=ret-(r-l+1)*getmu(res);
            if(r==n) break;
        }
        return smu2[n]=ret;
    }

    il ll getphi(int n){
        static ll ret; ret=0;
        static int res;
        for(int l=1,r;;l=r+1){
            res=n/l,r=n/res;
            ret=ret+(getmu(r)-getmu(l-1))*(res*(res+1ll)/2);
            if(r==n) break;
        }
        return ret;
    }
}

对 id\mu

直接瞪不太容易瞪出来合适的 g，我们先假设 g 已经找到，把卷积展开：

(f\ast g)(n)=\sum\limits_{d\mid n}d\mu(d)g(\frac{n}{d})

意识到 f=id\mu，而

(id\ast id)(n)=\sum\limits_{d\mid n} d\frac{n}{d}=nd(n)

考虑类似的手法，令 g=id，于是

\begin{aligned} (f\ast g)(n) & =\sum\limits_{d\mid n} d\mu(d)\frac{n}{d} \\ & =n\sum\limits_{d\mid n} \mu(d) \\ & =n(\mu \ast I)(n) \\ & =n\epsilon(n) \\ & =\epsilon(n) \end{aligned}

额...有点出乎意料。不过确实 g 和 h=\epsilon 都可以 O(1) 求前缀和，直接杜就好了。

对 id\varphi

故技重施，快进到

\begin{aligned} (f\ast g)(n) & =n\sum\limits_{d\mid n} \varphi(d) \\ & =n(\varphi\ast I)(n) \\ & =nid(n) \\ & =id_2(n) \end{aligned}

那么可以直接杜了。

P3768 简单的数学题

题意：求 \sum\limits_{i=1}^n \sum\limits_{j=1}^n ij(i,j)\pmod{p}。
数据范围：n\leqslant 10^{10}，4s。
考虑莫反：

\begin{aligned} \sum\limits_{i=1}^n \sum\limits_{j=1}^n ij(i,j) & =\sum\limits_{d=1}^n d \sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor} i\times d \sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor} j\times d\times [(i,j)=1] \\ & =\sum\limits_{d=1}^n d^3 \sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor} i \sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor} j \times \epsilon((i,j)) \\ & =\sum\limits_{d=1}^n d^3 \sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor} i \sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor} j \sum\limits_{t\mid i\And t\mid j} \mu(t) \\ & =\sum\limits_{d=1}^n d^3 \sum\limits_{t=1}^{\lfloor\frac{n}{d}\rfloor} \mu(t) \sum\limits_{i=1}^{\lfloor\frac{n}{td}\rfloor} i\times t \sum\limits_{j=1}^{\lfloor\frac{n}{td}\rfloor} j\times t \\ & =\sum\limits_{d=1}^n d^3 \sum\limits_{t=1}^{\lfloor\frac{n}{d}\rfloor} t^2\mu(t) \sum\limits_{i=1}^{\lfloor\frac{n}{td}\rfloor} i \sum\limits_{j=1}^{\lfloor\frac{n}{td}\rfloor} j \\ & =\sum\limits_{d=1}^n d^3 \sum\limits_{t=1}^{\lfloor\frac{n}{d}\rfloor} t^2\mu(t) \frac{\lfloor\frac{n}{td}\rfloor^2(\lfloor\frac{n}{td}\rfloor+1)^2}{4} \\ & =\sum\limits_{td=1}^n td^2\frac{\lfloor\frac{n}{td}\rfloor^2(\lfloor\frac{n}{td}\rfloor+1)^2}{4} \sum\limits_{d\mid td} d\mu(\frac{td}{d}) \\ & =\sum\limits_{td=1}^n td^2 \varphi(td)\frac{\lfloor\frac{n}{td}\rfloor^2(\lfloor\frac{n}{td}\rfloor+1)^2}{4} \end{aligned}

...真（炎国粗口）长。
事实上，我们可以用欧拉反演：

\begin{aligned} \sum\limits_{i=1}^n \sum\limits_{j=1}^n ij(i,j) & =\sum\limits_{i=1}^n \sum\limits_{j=1}^n ij\times id((i,j)) \\ & =\sum\limits_{i=1}^n \sum\limits_{j=1}^n ij \sum\limits_{d\mid i\And d\mid j} \varphi(d) \\ & =\sum\limits_{d=1}^n d^2\varphi(d) \sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor} i \sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor} j \\ & =\sum\limits_{d=1}^n d^2\varphi(d)\frac{\lfloor\frac{n}{d}\rfloor^2(\lfloor\frac{n}{d}\rfloor+1)^2}{4} \end{aligned}

就是这么简洁。
问题变成求 f(n)=n^2\varphi(n) 的前缀和 F，考虑怎么杜：取 g=id_2，有

(f\ast g)(n)=\sum\limits_{d\mid n} d^2\varphi(d)(\frac{n}{d})^2=n^2\sum\limits_{d\mid n}\varphi(d)=n^3=id_3(n)

那么 g=id_2,h=id_3，问题解决。应当指出 G(n)=\frac{n(n+1)(2n+1)}{6},H(n)=ID^2(n)，建议背过。推一下杜的式子：

\begin{aligned} ID^2(n) & =\sum\limits_{i=1}^n \sum\limits_{d\mid i}f(\frac{n}{i})id_2(i) \\ & =\sum\limits_{d=1}^n id_2(d)F(\lfloor\frac{n}{d}\rfloor) \\ & =F(n)+\sum\limits_{d=2}^n id_2(d)F(\lfloor\frac{n}{d}\rfloor) \end{aligned}

所以有

F(n)=(\frac{n(n+1)}{2})^2-\sum\limits_{d=2}^n id_2(d)F(\lfloor\frac{n}{d}\rfloor)

终于结束了。呼~

Powerful Number 筛

下简称 PN 筛。PN 筛可以在亚线性时间内求出某些特征积性函数的前缀和。
所谓某些特征积性函数，指的是存在一个对应的积性函数 g，使得 g(p)=f(p)，且求 g 的前缀和不比杜更复杂。称 g 为 f 的质数拟合函数。
对于此种函数，我们构造一个函数 h 使得 f=g\ast h。h 的存在性显然，因为 h 实际上是 g^{-1}\ast f，而逆元存在性我们在狄利克雷卷积处已经证过了。从这里也能得出 h 的积性。
将上面的卷积式展开，得到：

\begin{aligned} F(n) & =\sum\limits_{i=1}^n \sum\limits_{d\mid i} g(d)h(\frac{n}{d}) \\ & =\sum\limits_{d=1}^n h(d)G(\lfloor\frac{n}{d}\rfloor) \end{aligned}

这里 g 的前缀和我们可以 O(1) 得知，于是只要解决 h 的问题，就可以使用和杜教筛一样的整除分块办法来求解。关键是 h 要满足一个性质：其只在 Powerful Nubmer 处有值。

Powerful Number

定义：记正整数 n 的质数-指数向量为 n=\prod\limits_{i=1}^{m} p_i^{k_i}。n 是 Powerful Number（下简称 PN），当且仅当 \forall i\in [1,m],k_i>1。
性质 1：所有 PN 都可以写成 a^2b^3 的形式。
- 对于 k_i\bmod 2=0，将 p_i^{k_i} 合并入 a^2；
- 对于 k_i\bmod 2=1，先将 p_i^3 合并入 b^3，再将剩余的 p_i^{k_i-3} 合并入 a^2。
性质 2：[1,n] 中的 PN 只有 O(\sqrt{n}) 个。
- 考虑枚举 a，于是 PN 的数量为
\sum\limits_{a=1}^{\sqrt{n}} \sqrt[3]{\frac{n}{a^2}}=O(\sqrt{n})
- 这个等号是积分求得的，我暂时不会。
- 怎么求出？线筛出 \sqrt{n} 内的所有质数，然后搜所有质数的指数即可。
性质 3：对任何的 f=g\ast h\And f(p)=g(p)，h 函数仅在 PN 处可能有值。
- 考虑将卷积式在 p 处展开：
\begin{aligned} f(p) & =g(1)h(p)+g(p)h(1) \\ & =h(p)+g(p) \end{aligned}
- 故显然 h(p)=0。因为 h 是积性函数，则 h 在非 PN 处一定无值。
- 注意 1 是 PN，因为 1 的质因数分解形式中 m=0。
于是我们可以在筛 PN 的过程中直接计算 F，唯一的问题是 h 在 PN 处的值到底是多少。两种求法：
- 推出 h(p^k) 的公式，或者说，h(p^k) 仅关于 p 和 k 的表达式。比较难，但有一个很优的复杂度为 O(\dfrac{\sqrt{n}}{\ln \sqrt{n}})。
- 由 f=g\ast h，展开得
\begin{aligned} f(p^k) & =\sum\limits_{i=0}^kg (p^i)h(p^{k-i}) \\ & =h(p^k)+\sum\limits_{i=1}^k g(p^i)h(p^{k-i}) \end{aligned}
- 于是有
h(p^k)=f(p^k)-\sum\limits_{i=1}^k g(p^i)h(p^{k-i})
- 枚举 p,k 计算之，若已知 f(p^k)，则复杂度为
O(\frac{\sqrt{n}}{\ln \sqrt{n}})\times \log^2 n\approx O(\sqrt{n}\log)
- 分别是质数个数、枚举 k、枚举 i。事实上，这是一个很宽的上界，远远不满（较大质数的次数很小，远不及一般意义上的 \log_2 n）。
- 后来在 Min_25 筛的板题中的实验也表明，纯暴力求（理论复杂度 O(\sqrt{n}\log)）和公式法（理论复杂度 O(\dfrac{\sqrt{n}}{\log}) 的用时一样，可见瓶颈似乎不在这部分（该题需要杜）。
- 是的，按道理讲这里我们还没有筛出 f(p^k)，但积性函数就是通过给出 f(p^k) 的取值定义的，故我们在题面中就知道了。注意这里的 g 也应当用这种方式求，如果杜的话...感觉上显然复杂度不对，试了一下也没跑出来。
然后搜索所有 PN 直接计算 F，式子为

F(n)=\sum\limits_{i\in PN\And i\leqslant n} h(i)G(\lfloor\frac{n}{i}\rfloor)

若 G 可以 O(1) 得知，则复杂度为搜索的 O(\sqrt{n})。
若 G 需要杜，则复杂度为 O(n^\frac{2}{3})，因若先杜一次 G(n)，则访问到的所有 G 都是已经筛出的值。
理由：对于 i=1，G 是 G(n)；否则 G 在杜中是求 G(n) 时被减去的项，整除分块过了。
综上，两部分复杂度加合，可以认为 PN 筛的复杂度上界是 O(n^\frac{2}{3})。
最后，再度给出便于实现的一般过程：
1. 构造 f 的质数拟合函数 g。
2. 构造快速计算 G 的方法。
3. 计算 h(p^k)，通常是预处理把表打出来，这样第 4 步直接查表就好。
4. 搜索所有 PN，过程中累加答案。
给出示范代码，这里 f(p^k)=p^k(p^k-1)，即 P5325 【模板】Min_25 筛（h 部分使用纯暴力计算）：

namespace du{
    //f=iphi(i),g=id,h=id_2
    const int maxn23=4641593,maxpn23=325167;
    const ll usen23=4641590;
    ll n23;
    bool np[maxn23];
    ll p[maxpn23]; int ptot;
    ll f[maxn23],sf[maxn23];
    il void init(){
        n23=ceil(pow(n,2/(long double)3));
        np[1]=1,f[1]=sf[1]=1; ll res;
        for(ll i=2;i<=n23;++i){
            if(!np[i]) p[++ptot]=i,f[i]=i*(i-1)%mod;
            sf[i]=addr(sf[i-1],f[i]);//printf("phi[%lld]:%lld f[%lld]:%lld sf[%lld]:%lld\n",i,phi[i],i,f[i],i,sf[i]);
            for(int j=1;i*p[j]<=n23 && j<=ptot;++j){
                res=i*p[j],np[res]=1;
                if(i%p[j]) f[res]=f[i]*p[j]%mod*(p[j]-1)%mod;
                else{
                    f[res]=f[i]*p[j]%mod*p[j]%mod;
                    break;
                }
            }
        }
        return;
    }

    il ll getid(ll nnow){nnow%=mod; return nnow*(nnow+1)%mod*inv[2]%mod;}
    il ll getid2(ll nnow){nnow%=mod; return nnow*(nnow+1)%mod*(2*nnow+1)%mod*inv[6]%mod;}

    struct my_hash{
        static ull splitmix64(ull x){
            x+=0x9e3779b97f4a7c15;
            x=(x^(x>>30))*0xbf58476d1ce4e5b9;
            x=(x^(x>>27))*0x94d049bb133111eb;
            return x^(x>>31);
        }

        size_t operator()(ull x) const {
            static const ull FIXED_RANDOM = chrono::steady_clock::now().time_since_epoch().count();
            return splitmix64(x + FIXED_RANDOM);
        }
    };
    u_map<ll,ll,my_hash> sf2;
    il ll getf(ll nnow){
        if(nnow<=n23) return sf[nnow];
        if(sf2.find(nnow)!=sf2.end()) return sf2[nnow];
        ll ret=getid2(nnow); static ll res;
        for(ll l=2,r;;l=r+1){
            res=nnow/l,r=nnow/res;
            minu(ret,minur(getid(r),getid(l-1))*getf(res)%mod);
            if(r==nnow) break;
        }
        return sf2[nnow]=ret;
    }
}

namespace pn{
    //f(p^k)=p^k(p^k-1),g=idphi,h=? 这里 g 是 du::f
    const int maxpsq=9597;
    ll nsq; int uptot;
    vector<ll> h[maxpsq];//h[i][j] 表示 h(p[i]^j) 的值
    il void init(){//计算 h(p^k)。
        //注意到至少 h(p^2) 处才有值，故只需要考虑 sqrt(n) 以内的质数 
        nsq=ceil(sqrt(n)); if(n==1) return;//这是因为 n=1 时 p[1]=0，但我不想在循环里加判断 
        ll powp[36]={1}; int top=0;//2^35>10^10，足够了 
        for(int i=1;du::p[i]<=nsq;++i){
            while(powp[top]*du::p[i]<=n) ++top,powp[top]=powp[top-1]*du::p[i];
            h[i].resize(top+1),h[i][0]=1,h[i][1]=0;
            For(k,2,top){
                powp[k]%=mod;
                h[i][k]=powp[k]*(powp[k]-1)%mod;
                For(j,1,k) minu(h[i][k],(powp[j]*minur(powp[j],powp[j-1])%mod)*h[i][k-j]%mod);
            }
            top=0,uptot=i;
        }
        return;
    }

    ll nnow,lim[maxpsq],sum;
    il void dfs(int now,ll vn,ll hn){//考虑到第几个质数，当前的值是多少，当前的 h 是多少 
        if(now>uptot || vn*du::p[now]>=lim[now]){//计算贡献。这个剪枝非常重要。 
            sum=(sum+hn*du::getf(nnow/vn))%mod;
            return;
        }
        dfs(now+1,vn,hn);//k 为 0
        int k=2;
        for(ll res=vn*du::p[now];res<lim[now];res=res*du::p[now],++k) dfs(now+1,res*du::p[now],hn*h[now][k]%mod);
        return; 
    }

    il ll getf(ll nnowin){
        nnow=nnowin;
        For(i,1,uptot) lim[i]=nnow/du::p[i]+1;//这是一个防止 10^10 相乘爆 ll 的 trick，记 v*p=x，则 x*p<=nnow->x<=nnow/p->x<floor(nnow/p)+1 
        sum=0;
        dfs(1,1,1);
        return sum;
    }
}

P5325 【模板】Min_25 筛

题意：求 F(n)，这里积性函数 f(p^k)=p^k(p^k-1)。
数据范围：n\leqslant 10^{10}。
首先我们令 k=1 看看 f 的表现，发现 f(p)=p(p-1)，注意到 \varphi(p)=p-1，故考虑构造 f 的质数拟合函数 g=id\varphi。
嗯？这个函数好像似曾相识。是的，这个函数可以杜，参看杜教筛的最后一个例子。
于是是板子，O(\sqrt{n}\log+n^\frac{2}{3})。
但本题有特殊的一点在于 h 是可以求公式的，思路类似杜，过程如下：

\begin{aligned} f(p^k) & =(g\ast h)(p^k) \\ & =\sum\limits_{i=0}^k g(p^{k-i})h(p^i) \\ & =h(p^k)+\sum\limits_{i=1}^k p^{k-i}\varphi(p^{k-i})h(p^i) \\ & =h(p^k)+\sum\limits_{i=1}^k p^{k-i}(p-1)p^{k-i-1}h(p^i) \\ & =h(p^k)+\sum\limits_{i=1}^k p^{2k-2i-1}(p-1)h(p^i) \end{aligned}

（未完成）
又 h(p)=0，于是通过累加法可以得到 h(p^k)=(k-1)(p-1)p^k。

HDU7186 Jo loves counting

题意：
- 对一个数 n，我们定义 good_n 为所有质数-指数向量的质数部分和 n 的质数-指数向量的质数部分相同的数构成的集合。
- 形式化地，记 n=\prod\limits_{i=1}^m p_i^{k_i}，则 i=\prod\limits_{i=1}^m p_i^{k_i}\And k_i\geqslant 0 \leftrightarrow i\in good_n。
- 定义函数 f(n)=\frac{n}{|good_n|}，求 \frac{1}{m}F(m)，对 4179340454199820289（一个大质数，=29\times 2^{57}+1）取模。
数据范围：T\leqslant 12,m\leqslant 10^{12}，但 m>10^6 的 T\leqslant 6。
我们先考察一下 |good_n|，发现其其实是 \prod\limits_{i=1}^m k_i，因为 good 中元素的质数-指数向量中，p_i 的次数可以取 1\sim k_i。
于是发现 f 在 p 处的表现为 f(p)=\frac{p}{1}=p，考虑取 id 为其质数拟合函数。
后面就是板子咯，而且还不用杜，O(T\sqrt{m}\log)。