51nod1986 Jason曾不想做的数论题

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

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

$\prod\limits_{{X\in[1,m]^n}}(lcmX)^{gcdX} \\=\prod\limits_{g=1}^m\prod\limits_{X\in[1,m/g]^n}(lcmX)^{g\cdot[gcdX=1]} \\= \prod\limits_{g=1}^m \prod\limits_{X\in[1,m/g]^n}(g\cdot lcmX)^{g\sum\limits_{d|gcdX}\mu(d)} \\= \prod\limits_{g=1}^m \prod\limits_{d=1}^{m/g} \big( \prod\limits_{X\in[1,m/(gd)]^n}(gd\cdot lcmX) \big)^{g\cdot\mu(d)} \\= \prod\limits_{t=1}^m \big( \prod\limits_{X\in[1,m/t]^n}(t\cdot lcmX) \big)^{\phi(t)} \\= \prod\limits_k \big( \prod\limits_{X\in[1,k]^n}lcmX \big)^{\sum\limits_{\lfloor m/t\rfloor=k}\phi(t)} \cdot (\prod\limits_{\lfloor m/t\rfloor=k}t^{\phi(t)})^{k^n} \\= \prod\limits_k F^{G} \cdot H^{k^n}$

线性筛预处理1~m，可以 $ O(1) $ 回答欧拉函数区间和G。
F可以枚举质数算贡献，较小的质数暴力容斥计算lcm中这个质数为每个幂次时的方案数，超过$\sqrt{k}$ 的质数幂次不超过1，可以按k/p的取值分类计算，于是F可以在 $ O(\sqrt k\cdot 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)>>;

    return v>=MOD?v-MOD:v;

}

inline u32 $(u32 x){return mul(x,RR);}

inline u32 i$(u32 x){return mul(x,);}

inline u32 power(u32 a,i64 n){

    u32 v=ONE;

    for(;n;n>>=,a=mul(a,a))if(n&)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)>>;

    return v>=MOD?v-MOD:v;

}

inline u32 $(u32 x){return mul(x,RR);}

inline u32 i$(u32 x){return mul(x,);}

inline u32 power(u32 a,i64 n){

    u32 v=ONE;

    for(;n;n>>=,a=mul(a,a))if(n&)v=mul(v,a);

    return v;

}

}

const int N=1e8,cP=6e6+,P=1e9+,P1=P-;

int ps[cP],log5p[cP],pp=;

std::bitset<N+>np;

int t[],phi[N+],log5[N+];

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)<)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<?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 ;

    using namespace Z_P1d2;

    int v=i$(power($(a%(P1/)),n));

    if((v^a)&)v+=P1/;

    return v;

}

int pri(int x){

    return std::upper_bound(ps+,ps+pp+,x)-ps-;

}

int ttt=,pw_x[],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=;

    for(int i=,p;p=ps[i],p*p<=m;++i){

        int u=;

        for(int j=,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+,r,c;l<=m;l=r+){

        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>(,v1));

    return v;

}

int solve(int m,int n){

    int sq=sqrt(m),v=;

    SQ=sq;

    for(int i=;i<=sq;++i)pw_x[i]=pw<P1>(i,n);

    for(int i=;i<=sq;++i)t[i]=;

    for(int l=,r,c,s,sx;l<=m;l=r+){

        r=m/(c=m/l);

        s=sub<P1>(phi[r],phi[l-]);

        sx=sub<P1>(log5[r],log5[l-]);

        v=mul<P>(v,cal(c,n,s,sx));

    }

    for(int i=;i<=sq;++i)v=mul<P>(v,pw<P>(ps[i],t[i]));

    return v;

}

const int D=;

int nx[P>>D|][];

void cal_pri(const int N){

    int tt=clock();

    phi[]=;

    for(int i=;i<=N/;++i){

        if(!np.test(i))ps[++pp]=i,phi[i]=i-;

        for(int j=,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]-);

        }

    }

    for(int i=N/+;i<=N/;++i){

        if(!np.test(i))ps[++pp]=i,phi[i]=i-;

        np.set(i*);

        if(i&){

            phi[i*]=phi[i];

            np.set(i*);

            phi[i*]=phi[i]*(i%?:);

        }else phi[i*]=phi[i]*;

    }

    for(int i=N/+;i<=N/;++i){

        if(!np.test(i))ps[++pp]=i,phi[i]=i-;

        np.set(i*);

        phi[i*]=phi[i]*(-(i&));

    }

    for(int i=N/+;i<=N;++i)if(!np.test(i))ps[++pp]=i,phi[i]=i-;

}

void cal_log5(const int N){

    for(int l=,r;l<P;l+=<<D){

        r=l+(<<D)-;

        int L=l,R=r,w=(1ll<<)%P,t=;

        do{

            L=mul<P>(L,);

            R=mul<P>(R,);

            w=mul<P>(w,);

            ++t;

        }while(L<=R&&L>N&&t<=);

        nx[l>>D][]=t;

        nx[l>>D][]=w;

    }

    int tt=clock();

    int gp=,gt=;

    do{

        do{

            gt+=nx[gp>>D][];

            gp=Z_P::mul(gp,nx[gp>>D][]);

        }while(gp>N);

        if(>>gp%&)log5[gp]=gt;

    }while(gp!=);

    log5[]=;

    log5[]=;

    log5[]=;

    log5[]=;

    log5[]=;

    for(int i=;i<=N;i+=){

        log5[i]=add<P1>(log5[i/],);

        log5[i+]=add<P1>(log5[i/+],);

        log5[i+]=add<P1>(log5[i/+],);

        log5[i+]=add<P1>(log5[i/+],);

    }

}

void pre(const int N){

    int tt=clock();

    cal_pri(N);

    cal_log5(N);

    for(int i=;i<=pp;++i)log5p[i]=add<P1>(log5p[i-],log5[ps[i]]);

    for(int i=;i<=N;++i){

        log5[i]=add<P1>(log5[i-],mul<P1>(phi[i],log5[i]));

        adds<P1>(phi[i],phi[i-]);

    }

}

int qs[][],qp;

int main(){

    int mx=1e7;

    scanf("%d",&qp);

    for(int i=;i<qp;++i){

        scanf("%d%d",qs[i],qs[i]+);

        mx=std::max(mx,qs[i][]);

    }

    pre(mx);

    for(int i=;i<qp;++i)printf("%d\n",solve(qs[i][],qs[i][]));

    return ;

}

巴特西

51nod1986 Jason曾不想做的数论题

最新文章

热门文章