LGP5702题解
感觉 Min_25 相当于是给出了一个 \(O(f(M)\sqrt{n}\log n)\) 求 \(\prod_{i=1}^{n}f(i)\) 的方法（\(f(M)\) 是 \(f\) 的返回类型做一次“乘法”运算需要的复杂度）
上面基本上是在瞎扯，只不过是尝试在总结这一类的问题应该怎么做（
比如在快速阶乘算法中这个函数就是：
\[f(i)=i
\]
在此题中，我们知道：
\[H_n=H_{n-1}+\frac{1}{n}
\]
那么就是 \(f(i)=\frac{1}{i}=\frac{F(i)}{G(i)}\)。
和快速阶乘算法一样，我们设 \(g_d(x)=\prod_{i=1}^{d}f(x+i)\)，\(B=\lfloor\sqrt{n}\rfloor\)
假设目前我们知道了 \(g_d(iB)\)（\(i\in[0,d]\)），想要知道 \(g_{d+1}(iD)\)（\(i\in[0,d+1]\)）或 \(g_{2d}(iB)\)（\(i\in[0,2d]\)）
\[g_d(iB)(i\in[0,d])\to g_{d+1}(iB)(i\in[0,d+1])
\]
很显然有 \(g_{d+1}(iB)=g_d(iB)f(iB+d+1)\)，对于 \(g_{d+1}((d+1)B)=\prod_{i=1}^{n}f((d+1)B+i)\) 直接算就行。该部分复杂度 \(O(f(M)d)\)。
\[g_d(iB)(i\in[0,d])\to g_{2d}(iB)(i\in[0,2d])
\]
先设 \(h_i(i)=g_d(iB)\)，然后可以通过点值平移求出 \(h(i)=g(iB)\)（\(i\in[d,2d]\)），再一个点值平移求出 \(g(iB+d)\)（\(i\in[0,2d]\)）。
对于分式 \(\frac{F(x)}{G(x)}\)，显然只能对于分子和分母分别点值平移，然后通分乘起来了。
#include<cstdio>
#include<cmath>
#define IMP(lim,act) for(int qwq=(lim),i=0;i^qwq;++i)act
const int M=1<<19|5;
struct Barrett{
typedef unsigned long long ull;
typedef __uint128_t LL;
ull m,B;
Barrett(const ull&m=2):m(m),B((LL(1)<<64)/m){}
friend inline ull operator%(const ull&a,const Barrett&mod){
ull r=a-mod.m*(LL(mod.B)*a>>64);return r>=mod.m?r-mod.m:r;
}
}mod;
int P,ifac[M];
int buf[M<<2],*w[25];
inline int Getlen(const int&n){
int len(0);while((1<<len)<n)++len;return len;
}
inline int Add(const int&a,const int&b){
return a+b>=P?a+b-P:a+b;
}
inline int Del(const int&a,const int&b){
return b>a?a-b+P:a-b;
}
inline void swap(int&a,int&b){
int c=a;a=b;b=c;
}
inline int qpow(int a,int b=P-2){
int ans(1);for(;b;b>>=1,a=1ll*a*a%mod)if(b&1)ans=1ll*ans*a%mod;return ans;
}
inline void init(const int&n,const int&g){
const int&m=Getlen(n);int*now=buf;w[m]=now;now+=1<<m;
w[m][0]=1;w[m][1]=qpow(g,P-1>>m+1);for(int i=2;i^1<<m;++i)w[m][i]=1ll*w[m][i-1]*w[m][1]%mod;
for(int k=m-1;k>=0&&(w[k]=now,now+=1<<k);--k)IMP(1<<k,w[k][i]=w[k+1][i<<1]);
}
inline void DFT(int*f,const int&M){
const int&n=1<<M;
for(int len=n>>1,d=M-1;d>=0;--d,len>>=1)for(int k=0;k^n;k+=len<<1){
int*W=w[d],*L=f+(k),*R=f+(k|len),x,y;IMP(len,(x=*L,y=*R)),*L++=Add(x,y),*R++=1ll**W++*Del(x,y)%mod;
}
}
inline void IDFT(int*f,const int&M){
const int&n=1<<M;
for(int len=1,d=0;d^M;++d,len<<=1)for(int k=0;k^n;k+=len<<1){
int*W=w[d],*L=f+(k),*R=f+(k|len),x,y;IMP(len,(x=*L,y=1ll**W++**R%mod)),*L++=Add(x,y),*R++=Del(x,y);
}
const int&k=qpow(n);IMP(n,f[i]=1ll*f[i]*k%mod);for(int i=1;(i<<1)<n;++i)swap(f[i],f[n-i]);
}
inline void Getinv(int*f,const int&n){
static int g[M];g[0]=f[0];for(int i=1;i<n;++i)g[i]=1ll*g[i-1]*f[i]%mod;
int t,c=qpow(g[n-1]);for(int i=n-1;i>=1;--i)t=f[i],f[i]=1ll*g[i-1]*c%mod,c=1ll*c*t%mod;f[0]=c;
}
inline void PT(int*f,int*g,const int&n,const int&m){
static int F[M],G[M],H[M];H[0]=1;IMP(n,H[0]=1ll*H[0]*(m-i)%mod);IMP(n+n,G[i]=m-n+i);G[0]=1;Getinv(G,n+n);
IMP(n,F[i]=1ll*ifac[i]*(n-i-1&1?P-ifac[n-i-1]:ifac[n-i-1])%mod*f[i]%mod);
for(int i=1;i<n;++i)H[i]=1ll*(m+i)*G[i]%mod*H[i-1]%mod;
const int&len=Getlen(n+n);DFT(F,len);DFT(G,len);IMP(1<<len,F[i]=1ll*F[i]*G[i]%mod);IDFT(F,len);
IMP(n,g[i]=1ll*F[n+i]*H[i]%mod);IMP(1<<len,F[i]=G[i]=H[i]=0);
}
inline int GetH(const int&n,const int&p,const int&g){
static int F1[M],G1[M],F2[M],G2[M];const int&B=sqrt(n),&m=Getlen(B+1)-2;
mod=Barrett(P=p);init(B<<1,g);ifac[0]=ifac[1]=1;for(int i=2;i<=B;++i)ifac[i]=1ll*(P-P/i)*ifac[P%i]%mod;
for(int i=1;i<=B;++i)ifac[i]=1ll*ifac[i-1]*ifac[i]%mod;F1[0]=1;F1[1]=B+1;F2[0]=1;F2[1]=1;
for(int i=m;i>=0;--i){
const int&q=B>>i+1,&p=B>>i;
PT(F1,G1,q+1,q+1);IMP(q,F1[q+i+1]=G1[i]);PT(F1,G1,q*2+1,1ll*q*qpow(B)%mod);
PT(F2,G2,q+1,q+1);IMP(q,F2[q+i+1]=G2[i]);PT(F2,G2,q*2+1,1ll*q*qpow(B)%mod);
IMP(q*2+1,(F2[i]=(1ll*F1[i]*G2[i]+1ll*F2[i]*G1[i])%mod,F1[i]=1ll*F1[i]*G1[i]%mod)),G1[i]=G2[i]=0;
if(q*2+1==p){
IMP(q*2+1,(F2[i]=(F1[i]+1ll*F2[i]*(i*B+p))%mod,F1[i]=1ll*F1[i]*(i*B+p)%mod));
F1[p]=1;F2[p]=0;IMP(p,(F2[p]=(F1[p]+1ll*F2[p]*(p*B+i+1))%mod,F1[p]=1ll*F1[p]*(p*B+i+1)%mod));
}
}
int ans(0);IMP(B,ans=(ans+1ll*F2[i]*qpow(F1[i]))%mod),F1[i]=F2[i]=0;F1[B]=F2[B]=0;
for(int i=B*B+1;i<=n;++i)ans=(ans+qpow(i))%mod;return ans;
}
signed main(){
int T,N,P,G;scanf("%d",&T);while(T--)scanf("%d%d%d",&N,&P,&G),printf("%d\n",GetH(N,P,G));
}
