题解 P5071 【[Ynoi2015] 此时此刻的光辉】

$\text{Link}$

我写的第 12 道 $\text{Ynoi}$ ，为了防止以后忘了写篇题解。

题意

区间乘积因数个数。

$n,m\le10^5,a_i\le10^9$

有关时间复杂度的约定

本片题解中， $n,m$ 默认同阶， $w$ 代表值域， $f(x)$ 表示 $x$ 内质数个数， $g(x)$ 表示 $x$ 以内数质因数个数最多为多少。

思路

显然可以分解质因数+因数个数定理搞一次操作。区间查询可以用莫队。

这样直接做时间复杂度是 $O(n(\sqrt{w}+g(w)\sqrt{n}))$ 。

这样复杂度过高，分解质因数可以用 $\text{Pollard-Rho}$ 优化，做到 $O(n(w^{\frac 1 4}+g(w)\sqrt{n}))$ 。

但这样还是不够快，考虑设置阈值 $S$ ，小于 $S$ 的质因数暴力分解出来做前缀和，大的 $\text{Pollard-Rho}$ 莫队直接转移，于是做到了时间复杂度 $O(n(w^{\frac 1 4}+\sqrt{n}(\log_Sw+f(S))))$ ，空间复杂度 $O(S+n(\log_Sw+f(S)))$ 。

$S$ 取 $1000$ 很不错，时间复杂度 $O(n(w^{\frac 1 4}+\sqrt{n}))$ ，空间复杂度 $O(n)$ ，都带个 $170$ 左右的常数。

$\text{Pollard-Rho}$ 的倍增上界为 $31$ 时表现比较优秀。

如果复杂度分析有误，请在讨论区中指出。

#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int N=1e5+10,size=170,mod=19260817;
int pri[]={1,2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,997};
int num[1010];
namespace PR{
    vector<int>ans;
    inline int qpow(int x,int y,const int mod){
        int res=1;
        while(y){
            if(y&1) res=1ll*res*x%mod;
            x=1ll*x*x%mod;
            y>>=1;
        }
        return res;
    }
    inline bool Miller_Rabin(int x,int a){
        int y=x-1;
        while(y){
            int res=qpow(a,y,x);
            if(res!=1&&res!=x-1) return 0;
            if(y&1||res==x-1) return 1;
            y>>=1;
        }
        return 1;
    }
    inline bool Isprime(int x){
        if(x<2||x==46856248255981) return 0;
        if(x==2||x==3||x==7||x==61||x==24251) return 1;
        if(x%6!=1&&x%6!=5) return 0;
        return Miller_Rabin(x,2)&&Miller_Rabin(x,61);
    }
    inline int f(int x,int c,int n){
        return (1ll*x*x+c)%n;
    }
    inline int Pollard_Rho(int x){
        int s=0,t=0,c=1ll*rand()%(x-1)+1,val=1;
        int stp,g=1;
        for(;;g<<=1,s=t,val=1){
            for(stp=1;stp<=g;stp++){
                t=f(t,c,x);
                val=1ll*val*abs(t-s)%x;
                if(!(stp&31)){
                    int q=__gcd(val,x);
                    if(q>1) return q;
                }
            }
            int q=__gcd(val,x);
            if(q>1) return q;
        }   
    }
    inline void fj(int x,int cnt){
        if(x<2) return ;
        if(Isprime(x)) {for(int i=1;i<=cnt;i++)ans.push_back(x);/*printf("-%d\n",x);*/return ;}
        int pri=x;
        for(;pri>=x;pri=Pollard_Rho(x));
        int s=0;
        for(;x%pri==0;x/=pri)s++;
        fj(x,cnt),fj(pri,cnt*s);
    }
}
int inv[N*2],sum[N][size];
inline void Prefix(int n){
    inv[1]=1;
    for(int i=2;i<=n;i++){
        inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
    }
    for(int i=1;i<=168;i++){
        num[pri[i]]=i;
    }
}
int n,m,curL,curR,sz,cnt[N*2],p[N][3],a[N],now=1,divv[N],f[N],ss[N];
int g[N];
vector<int>appear;
struct que{
    int id,l,r; 
    friend bool operator<(const que &a,const que &b){
        return divv[a.l]^divv[b.l]?a.l<b.l:a.r<b.r;
    }
}ask[N];
inline void add(int x){
    for(int i=1;i<=ss[x];i++){
        now=1ll*now*inv[cnt[p[x][i]]]%mod*(cnt[p[x][i]]+1)%mod;
        cnt[p[x][i]]++;
    }
}
inline void del(int x){
    for(int i=1;i<=ss[x];i++){
        now=1ll*now*inv[cnt[p[x][i]]]%mod*(cnt[p[x][i]]-1)%mod;
        cnt[p[x][i]]--;
    }
}
inline int getans(int l,int r){
    int ans=now;
    for(int i=1;i<=168;i++){
        ans=1ll*ans*(sum[r][i]-sum[l-1][i]+1)%mod;
    }
    return ans;
}
int main(){
//  freopen("test.in","r",stdin);
    n=read(),m=read();
    Prefix(n*2+10);
    sz=n/sqrt(m*2/3)+1;
    for(int i=1;i<=n;i++){
        a[i]=read();
        divv[i]=i/sz;
        memcpy(sum[i],sum[i-1],sizeof(sum[0]));
        for(int j=1;j<=168;j++){
            while(a[i]%pri[j]==0){
                sum[i][j]++;
                a[i]/=pri[j];
            }
        }
        if(a[i]>1){
            PR::ans.clear();
            PR::fj(a[i],1);
            for(auto j:PR::ans){
                p[i][++ss[i]]=j;
                appear.push_back(j);
            }
        }
    }
    sort(appear.begin(),appear.end());
    appear.erase(unique(appear.begin(),appear.end()),appear.end());
    for(int i=1;i<=n;i++){
        for(int j=1;j<=ss[i];j++){
            p[i][j]=lower_bound(appear.begin(),appear.end(),p[i][j])-appear.begin();
        }
    }
    for(int i=1;i<=m;i++){
        int l=read(),r=read();
        ask[i]={i,l,r};
    }
    sort(ask+1,ask+m+1);
    curL=1,curR=0;
    for(int i=0;i<n*2;i++) cnt[i]=1;
    for(int i=1;i<=m;i++){
        int L=ask[i].l,R=ask[i].r;
        while(curL>L) add(--curL);
        while(curR<R) add(++curR);
        while(curL<L) del(curL++);
        while(curR>R) del(curR--);
        g[ask[i].id]=getans(L,R);
    }
    for(int i=1;i<=m;i++){
        write(g[i]);
        putc('\n');
    }
    flush();
    return 0;
}