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)>>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\n",solve(qs[i][1],qs[i][0])); return 0; }
