51nod1986 Jason曾不想做的数论题

http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1986

方便起见，公式中的区间内只考虑整数，X的gcd,lcm定义为每个元素的gcd,lcm，d|X表示X中元素均为d的倍数

$prodlimits_{{Xin[1,m]^n}}(lcmX)^{gcdX} \=prodlimits_{g=1}^mprodlimits_{Xin[1,m/g]^n}(lcmX)^{gcdot[gcdX=1]} \= prodlimits_{g=1}^m prodlimits_{Xin[1,m/g]^n}(gcdot lcmX)^{gsumlimits_{d|gcdX}mu(d)} \= prodlimits_{g=1}^m prodlimits_{d=1}^{m/g} ig( prodlimits_{Xin[1,m/(gd)]^n}(gdcdot lcmX) ig)^{gcdotmu(d)} \= prodlimits_{t=1}^m ig( prodlimits_{Xin[1,m/t]^n}(tcdot lcmX) ig)^{phi(t)} \= prodlimits_k ig( prodlimits_{Xin[1,k]^n}lcmX ig)^{sumlimits_{lfloor m/t floor=k}phi(t)} cdot (prodlimits_{lfloor m/t floor=k}t^{phi(t)})^{k^n} \= prodlimits_k F^{G} cdot H^{k^n}$

线性筛预处理1~m，可以 $ O(1) $ 回答欧拉函数区间和G。
F可以枚举质数算贡献，较小的质数暴力容斥计算lcm中这个质数为每个幂次时的方案数，超过$sqrt{k}$ 的质数幂次不超过1，可以按k/p的取值分类计算，于是F可以在 $ O(sqrt kcdot log(10^9+7)) $ 内计算。
由于H包含指数部分，不方便高效地预处理。如果能求出1~m模 $ (10^9+7) $ 的离散对数，即可$O(m) $预处理 $ O(logn) $ 询问H。 $ 10^9+7 $ 的最小正原根为5，把 $ [0,10^9+7) $ 分为一些小区间，将区间作为整体参与运算，估算出一个下界t，使得区间内的数至少乘上原根的t次方才可能<=m。于是我们可以从小到大枚举原根的幂$ 5^0,5^1,5^2...5^{10^9+4},5^{10^9+5} $，利用之前估计的下界，快速跳过原根的幂>m的情况，再加上一些底层优化，可以在几秒内完成这部分预处理（这也是整个算法最耗时的部分），比朴素算法快了几倍。由于估算的下界t并不是准确值，这部分的渐进时间复杂度不太好分析，估计在 $ O(frac{10^9+7}{log_5(10^9+7)}+m) $ 附近。
回答询问时，枚举m/t的取值k，查询F,G,H等然后乘起来，询问的时间复杂度为 $ O(m^{3/4}log(10^9+7)) $。

#include<bits/stdc++.h>
typedef long long i64;
typedef unsigned long long u64;
typedef unsigned u32;
namespace Z_P1d2{
const u32 MOD=500000003u,iMOD=3707490901u,RR=29087838u,ONE=294967272u;
inline u32 mul(u32 a,u32 b){
    u64 c=(u64)a*b;
    u32 v=c+u32(c)*iMOD*u64(MOD)>>32;
    return v>=MOD?v-MOD:v;
}
inline u32 $(u32 x){return mul(x,RR);}
inline u32 i$(u32 x){return mul(x,1);}
inline u32 power(u32 a,i64 n){
    u32 v=ONE;
    for(;n;n>>=1,a=mul(a,a))if(n&1)v=mul(v,a);
    return v;
}
}

namespace Z_P{
const u32 MOD=1000000007u,iMOD=2226617417u,RR=582344008u,ONE=294967268u;
inline u32 mul(u32 a,u32 b){
    u64 c=(u64)a*b;
    u32 v=c+u32(c)*iMOD*u64(MOD)>>32;
    return v>=MOD?v-MOD:v;
}
inline u32 $(u32 x){return mul(x,RR);}
inline u32 i$(u32 x){return mul(x,1);}
inline u32 power(u32 a,i64 n){
    u32 v=ONE;
    for(;n;n>>=1,a=mul(a,a))if(n&1)v=mul(v,a);
    return v;
}
}
const int N=1e8,cP=6e6+7,P=1e9+7,P1=P-1;
int ps[cP],log5p[cP],pp=0;
std::bitset<N+64>np;
int t[10007],phi[N+7],log5[N+7];

template<int P>
inline void adds(int&a,int b){if((a+=b)>=P)a-=P;}
template<int P>
inline void subs(int&a,int b){if((a-=b)<0)a+=P;}
template<int P>
inline int add(int a,int b){return a+=b,a>=P?a-P:a;}
template<int P>
inline int sub(int a,int b){return a-=b,a<0?a+P:a;}
template<int P>
inline int mul(int a,int b){return i64(a)*b%P;}
template<int P>
int pw(int a,i64 n){
    if(P==::P){
        using namespace Z_P;
        return i$(power($(a),n));
    }
    if(!n)return 1;
    using namespace Z_P1d2;
    int v=i$(power($(a%(P1/2)),n));
    if((v^a)&1)v+=P1/2;
    return v;
}

int pri(int x){
    return std::upper_bound(ps+1,ps+pp+1,x)-ps-1;
}
int ttt=0,pw_x[10007],SQ;
inline int pwP1(int m,int n){
    return m<=SQ?pw_x[m]:pw<P1>(m,n);
}
int cal(int m,int n,int k,int lr){
    int tt0=clock();
    int pw_m=pwP1(m,n);
    int v=1;
    for(int i=1,p;p=ps[i],p*p<=m;++i){
        int u=0;
        for(int j=1,mp=m;;++j){
            mp/=p;
            int tmp=pw<P1>(m-mp,n);
            subs<P1>(u,tmp);
            if(mp<p){
                adds<P1>(u,mul<P1>(j,pw_m));
                break;
            }
        }
        adds<P1>(t[i],mul<P1>(u,k));
    }
    v=pw<P>(v,k);
    ttt+=clock()-tt0;
    int sq=sqrt(m),s00=pri(sq),s0=s00,s1;
    int v1=mul<P1>(sub<P1>(log5p[pri(m)],log5p[s00]),pw_m);
    for(int l=sq+1,r,c;l<=m;l=r+1){
        r=m/(c=m/l);
        s1=pri(r);
        subs<P1>(v1,mul<P1>(sub<P1>(log5p[s1],log5p[s0]),pwP1(m-c,n)));
        s0=s1;
    }
    v1=add<P1>(mul<P1>(v1,k),mul<P1>(lr,pw_m));
    v=mul<P>(v,pw<P>(5,v1));
    return v;
}

int solve(int m,int n){
    int sq=sqrt(m),v=1;
    SQ=sq;
    for(int i=0;i<=sq;++i)pw_x[i]=pw<P1>(i,n);
    for(int i=1;i<=sq;++i)t[i]=0;
    for(int l=1,r,c,s,sx;l<=m;l=r+1){
        r=m/(c=m/l);
        s=sub<P1>(phi[r],phi[l-1]);
        sx=sub<P1>(log5[r],log5[l-1]);
        v=mul<P>(v,cal(c,n,s,sx));
    }
    for(int i=1;i<=sq;++i)v=mul<P>(v,pw<P>(ps[i],t[i]));
    return v;
}

const int D=20;
int nx[P>>D|11][2];
void cal_pri(const int N){
    int tt=clock();
    phi[1]=1;
    for(int i=2;i<=N/5;++i){
        if(!np.test(i))ps[++pp]=i,phi[i]=i-1;
        for(int j=1,k;j<=pp&&(k=i*ps[j])<=N;++j){
            np.set(k);
            if(!(i%ps[j])){
                phi[k]=phi[i]*ps[j];
                break;
            }
            phi[k]=phi[i]*(ps[j]-1);
        }
    }
    for(int i=N/5+1;i<=N/3;++i){
        if(!np.test(i))ps[++pp]=i,phi[i]=i-1;
        np.set(i*2);
        if(i&1){
            phi[i*2]=phi[i];
            np.set(i*3);
            phi[i*3]=phi[i]*(i%3?2:3);
        }else phi[i*2]=phi[i]*2;
    }
    for(int i=N/3+1;i<=N/2;++i){
        if(!np.test(i))ps[++pp]=i,phi[i]=i-1;
        np.set(i*2);
        phi[i*2]=phi[i]*(2-(i&1));
    }
    for(int i=N/2+1;i<=N;++i)if(!np.test(i))ps[++pp]=i,phi[i]=i-1;
}
void cal_log5(const int N){
    for(int l=0,r;l<P;l+=1<<D){
        r=l+(1<<D)-1;
        int L=l,R=r,w=(1ll<<32)%P,t=0;
        do{
            L=mul<P>(L,5);
            R=mul<P>(R,5);
            w=mul<P>(w,5);
            ++t;
        }while(L<=R&&L>N&&t<=4);
        nx[l>>D][0]=t;
        nx[l>>D][1]=w;
    }
    int tt=clock();
    int gp=1,gt=0;
    do{
        do{
            gt+=nx[gp>>D][0];
            gp=Z_P::mul(gp,nx[gp>>D][1]);
        }while(gp>N);
        if(34>>gp%6&1)log5[gp]=gt;
    }while(gp!=1);
    log5[1]=0;
    log5[2]=381838282;
    log5[3]=884237698;
    log5[4]=763676564;
    log5[5]=1;
    for(int i=6;i<=N;i+=6){
        log5[i]=add<P1>(log5[i/2],381838282);
        log5[i+2]=add<P1>(log5[i/2+1],381838282);
        log5[i+4]=add<P1>(log5[i/2+2],381838282);
        log5[i+3]=add<P1>(log5[i/3+1],884237698);
    }
}
void pre(const int N){
    int tt=clock();
    cal_pri(N);
    cal_log5(N);
    for(int i=1;i<=pp;++i)log5p[i]=add<P1>(log5p[i-1],log5[ps[i]]);
    for(int i=1;i<=N;++i){
        log5[i]=add<P1>(log5[i-1],mul<P1>(phi[i],log5[i]));
        adds<P1>(phi[i],phi[i-1]);
    }
}
int qs[10007][2],qp;
int main(){
    int mx=1e7;
    scanf("%d",&qp);
    for(int i=0;i<qp;++i){
        scanf("%d%d",qs[i],qs[i]+1);
        mx=std::max(mx,qs[i][1]);
    }
    pre(mx);
    for(int i=0;i<qp;++i)printf("%d
",solve(qs[i][1],qs[i][0]));
    return 0;
}