X000011

P1890 gcd区间

\(\gcd\) 是满足结合律的，所以考虑用 ST 表解决

时间复杂度 \(O((n\log n+m)\log a_i)\)

考虑到 \(n\) 很小，你也可以直接算出所有的区间 \(\gcd\)

时间复杂度 \(O(n^2\log a_i+m)\)

实现的时候我写了 ST 表

加了 fread&fwrite IO优化，跑得还是非常快的

#include<stdio.h>

#include<ctype.h>

#define gc (l==r&&(r=(l=c)+fread(c,1,1<<21,stdin),l==r)?EOF:*l++)

int n,m,t,x,y,p=-1,p1,a[10],lg[1002],f[1002][10];

char c[1<<21],*l=c,*r=c,cw[(1<<21)+10];

inline int read() { //快读

	int x=0; char ch=gc;

	while (!isdigit(ch)) ch=gc;

	while (isdigit(ch)) x=x*10+(ch^48),ch=gc;

	return x;

}

inline void flush() { fwrite(cw,1,p+1,stdout),p=-1; }

inline void write(int x) { //快写

	(p>>21)&&(flush(),0);

	do { a[++p1]=x%10; }while (x/=10);

	do { cw[++p]=a[p1]^48; }while (--p1);

}

int gcd(int x,int y) { return y?gcd(y,x%y):x; }

int main() {

	n=read(),m=read(),lg[0]=-1;

	for (int i=1; i<=n; ++i) lg[i]=lg[i>>1]+1,f[i][0]=read();

	for (int j=0,t=1; j<lg[n]; ++j,t<<=1)

		for (int i=1,mx=n-(t<<1)+1; i<=mx; ++i)

			f[i][j+1]=gcd(f[i][j],f[i+t][j]);

	while (m--) {

		x=read(),y=read(),t=lg[y-x+1];

		write(gcd(f[x][t],f[y-(1<<t)+1][t])),cw[++p]='\n';

	}

	return flush(),0;

}

P3327 [SDOI2015]约数个数和

先证明这么一个公式：

\[d(ij)=\sum_{x\mid i}\sum_{y\mid j}[(x,y)==1]
\]

考虑一个质因子 \(p\) ，它在 \(i\) 中出现次数为 \(a\) ，在 \(j\) 中出现次数为 \(b\)

那么它在 \(d(ij)\) 中的贡献为 \(a+b+1\)（ \(p^{0\sim a+b}\) ）

在右边的式子中，产生贡献需要保证 \(x,y\) 互质

那么 \(x,y\) 中必定有一个数不包含质因子 \(p\) ，而另一个随便包含多少个

简单计算得贡献为 \(a+b+1\)

发现不同质因子产生的贡献互相独立，那么式子得证

（这段看起来可能会比较不知所云，如果你不想了解直接记公式就行

然后开始推式子吧

\(x/y\) 表示 \(\left\lfloor \dfrac xy\right\rfloor\)

\[\begin{aligned}
& {\color{white}=}\sum_{i=1}^n\sum_{j=1}^m\sum_{x\mid i}\sum_{y\mid j}[(x,y)==1] \\
& = \sum_{x=1}^n\sum_{y=1}^m[(x,y)==1](n/x)(m/y) \\
& = \sum_{i=1}^n\sum_{j=1}^m[(i,j)==1](n/i)(m/j) \\
& = \sum_{i=1}^n\sum_{j=1}^m(n/i)(m/j)
\sum_{d\mid (i,j)}\mu(d) \\
& = \sum_{d=1}^{\min(n,m)}\mu(d)\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}\left\lfloor \dfrac {n/d}i\right\rfloor\left\lfloor \dfrac {m/d}j\right\rfloor \\
\end{aligned}\]

考虑预处理 \(\sum\limits_{i=1}^n(n/i)\)

不会的话请右转P1403

由这道题我们知道这个式子其实就是 \(1\sim n\) 所有数的约数个数之和

约数个数可以线性筛然后处理个前缀和（ \(S\) ）即可

现在只需要求 \(\sum\limits_{d=1}^{\min(n,m)}\mu(d)S(n/d)S(m/d)\)

这个就比较ez了直接筛 \(\mu\) 算前缀和，然后数论分块

时间复杂度 \(O(\sqrt n+\sqrt m)\)

总体时间复杂度 \(O(N+T(\sqrt n+\sqrt m))\) ，这里 \(N=50000\)

#include<stdio.h>

typedef long long ll;

const int N=50001; int T,t,t1,t2,n,m,mx,cnt;

int p[6000],v[N],mu[N],s[N],f[N]; double a[N]; ll ans,S[N];

inline int min(int x,int y) { return x<y?x:y; }

int main() {

	scanf("%d",&T),mu[1]=S[1]=1;

	for (int i=1; i<N; ++i) a[i]=1./i*(1+1e-14);

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

		v[i]||(p[++cnt]=i,mu[i]=-1,s[i]=2,f[i]=1);

		for (int j=1; (t=p[j]*i)<N; ++j) {

			v[t]=1;

			if (i%p[j]) mu[t]=-mu[i],s[t]=(s[i]<<1),f[t]=1;

			else { f[t]=f[i]+1,s[t]=s[i]*a[f[t]]*(f[t]+1); break; }

		}

		mu[i]+=mu[i-1],S[i]=S[i-1]+s[i];

	}

	while (T--) {

		scanf("%d%d",&n,&m),ans=0;

		for (int l=1,r,mx=min(n,m); l<=mx; l=r+1) {

			r=min(n*a[t1=n*a[l]],m*a[t2=m*a[l]]);

			ans+=S[t1]*S[t2]*(mu[r]-mu[l-1]);

		}

		printf("%lld\n",ans);

	}

	return 0;

}

代码加了一些卡常优化（参考SP26045），效果显著，然而还是只排在了最优解第二页/kk

P4213 【模板】杜教筛（Sum）

杜教筛板子题。我已经写了一篇博客讲这个，所以直接丢链接了：here

这里提供一份无注释的代码有注释的可以在上面那篇文章里找到

#include<stdio.h>

#include<unordered_map>

typedef long long ll;

const int N=7800000; ll ph[N+1];

int T,n,t,cnt,p[600000],v[N+1],mu[N+1],phi[N+1];

std::unordered_map<int,ll>smu,sph;

ll Sphi(int n) {

	if (n<=N) return ph[n]; if (sph[n]) return sph[n];

	int t; ll s=0;

	for (unsigned l=2,r; l<=n; l=r+1)

		r=n/(t=n/l),s+=Sphi(t)*(r-l+1);

	return sph[n]=((1ll*n+1)*n>>1)-s;

}

int Smu(int n) {

	if (n<=N) return mu[n]; if (smu[n]) return smu[n];

	int t; ll s=0;

	for (unsigned l=2,r; l<=n; l=r+1)

		r=n/(t=n/l),s+=Smu(t)*(r-l+1);

	return smu[n]=1ll-s;

}

int main() {

	scanf("%d",&T),ph[1]=mu[1]=1;

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

		v[i]||(p[++cnt]=i,mu[i]=-1,phi[i]=i-1);

		for (int j=1; j<=cnt&&(t=p[j]*i)<=N; ++j) {

			v[t]=1;

			if (i%p[j]==0) { phi[t]=phi[i]*p[j]; break; }

			mu[t]=-mu[i],phi[t]=phi[i]*(p[j]-1);

		}

		mu[i]+=mu[i-1],ph[i]=ph[i-1]+phi[i];

	}

	while (T--) {

		scanf("%d",&n);

		printf("%lu %d\n",Sphi(n),Smu(n));

	}

	return 0;

}

P1447 [NOI2010]能量采集

一道傻逼紫题

看图口胡易知题目要求的是 \(\sum\limits_{i=1}^n\sum\limits_{j=1}^m2\gcd(i,j)-1\)

考虑求出 \(S=\sum\limits_{i=1}^n\sum\limits_{j=1}^m\gcd(i,j)\) ，那么答案就是 \(2S-nm\)

\[S=\sum_{i=1}^n\sum_{j=1}^m\sum_{d\mid (i,j)}\varphi(d)=\sum_{d=1}^{\min(n,m)}\!\!\!\varphi(d)\sum_{i=1}^n\sum_{j=1}^m[d\mid i,d\mid j]=\sum_{d=1}^{\min(n,m)}\!\!\!\varphi(d)(n/d)(m/d)
\]

然后就随便做了我写的是线性筛预处理+整除分块

#include<stdio.h>

int n,m,t,tt,cnt; long long ans;

int p[40000],v[100010],phi[100010];

inline int min(int x,int y) { return x<y?x:y; }

int main() {

	scanf("%d%d",&n,&m); n>m&&(t=n,n=m,m=t),phi[1]=1;

	for (int i=2; i<=n; ++i) {

		v[i]||(p[++cnt]=i,phi[i]=i-1);

		for (int j=1; j<=cnt&&(t=p[j]*i)<=n; ++j) {

			v[t]=1;

			if (i%p[j]==0) { phi[t]=phi[i]*p[j]; break; }

			phi[t]=phi[i]*(p[j]-1);

		}

		phi[i]+=phi[i-1];

	}

	for (int l=1,r; l<=n; l=r+1) {

		r=min(n/(t=n/l),m/(tt=m/l));

		ans+=1ll*t*tt*(phi[r]-phi[l-1]);

	}

	printf("%lld\n",(ans<<1)-1ll*n*m);

	return 0;

}

P4450 双亲数

老套路了

要求的就是 \(\sum\limits_{i=1}^A\sum\limits_{j=1}^B[(i,j)==d]\)

假设 \(A\le B\)

\[\begin{aligned}
\sum_{i=1}^A\sum_{j=1}^B[(i,j)==d]
&=\sum_{i=1}^{A/d}\sum_{j=1}^{B/d}[(i,j)==1] \\
&=\sum_{i=1}^{A/d}\sum_{j=1}^{B/d}\sum_{d'\mid (i,j)}\mu(d') \\
&=\sum_{d'=1}^{A/d}\sum_{d'\mid i}\sum_{d'\mid j}\mu(d') \\
&=\sum_{d'=1}^{A/d}\mu(d')(A/d/d')(B/d/d')
\end{aligned}\]

这题没了

#include<stdio.h>

const int N=1000002; long long ans;

int n,m,d,cnt,t,tt,p[500000],v[N],mu[N];

inline int min(int x,int y) { return x<y?x:y; }

int main() {

	scanf("%d%d%d",&n,&m,&d),mu[1]=1;

	(n/=d)>(m/=d)&&(t=n,n=m,m=t);

	for (int i=2; i<=n; ++i) { //线性筛

		v[i]||(p[++cnt]=i,mu[i]=-1);

		for (int j=1; j<=cnt&&1ll*p[j]*i<=n; ++j) {

			v[t=p[j]*i]=1;

			if (i%p[j]==0) break; mu[t]=-mu[i];

		}

		mu[i]+=mu[i-1];

	}

	for (int l=1,r; l<=n; l=r+1) { //整除分块

		r=min(n/(t=n/l),m/(tt=m/l));

		ans+=1ll*(mu[r]-mu[l-1])*t*tt;

	}

	printf("%lld\n",ans);

	return 0;

}

P3768 简单的数学题

\[\begin{aligned}
\sum_{i=1}^n\sum_{j=1}^nij(i,j)
&= \sum_{i=1}^n\sum_{j=1}^nij\sum_{d\mid(i,j)}\varphi(d) \\
&= \sum_{d=1}^n\varphi(d)\sum_{d\mid i}\sum_{d\mid j}ij \\
&= \sum_{d=1}^n\varphi(d)d^2\left(\sum_{i=1}^{n/d}i\right)^2
\end{aligned}\]

然后用杜教筛可以 \(O(n^{2/3})\) 计算

考虑用杜教筛处理 \(f(n)=\varphi(n)n^2\) 的前缀和

发现这个是积性函数我们再整一个完全积性函数 \(g(n)=n^2\)

\[\begin{aligned}
\sum_{i=1}^n(f\ast g)(i)
&=\sum_{i=1}^n\sum_{d\mid i}f(d)g(\frac id) \\
&=\sum_{i=1}^n\sum_{d\mid i}\varphi(d)d^2\times\left(\dfrac id\right)^2 \\
&=\sum_{i=1}^ni^2\sum_{d\mid i}\varphi(d) \\
&=\sum_{i=1}^ni^3
\end{aligned}\]

用小学奥数的知识可得 \(\sum\limits_{i=1}^ni^3=\left(\dfrac{n(n+1)}2\right)^2\)

这样就好了

然后注意一下取模问题这题细节有点多

#include<stdio.h>

#include<unordered_map>

#include<math.h>

typedef long long ll;

const int M=4700000;

int p,m,inv,t,ans,cnt,pr[330000],phi[M],v[M];

ll n,tt; std::unordered_map<ll,int>mps;

int _power(int x,int y) {

	int ans=1;

	while (y) (y&1)&&(ans=1ll*ans*x%p),x=1ll*x*x%p,y>>=1;

	return ans;

}

int s2(int x) { return 1ll*x*(x+1)%p*(x<<1|1)%p*inv%p; }

int s3(int x) { x=(1ll*x*(x+1)>>1)%p; return 1ll*x*x%p; }

int S(ll x) {

	if (x<=m) return phi[x]; if (mps[x]) return mps[x];

	int s=0; ll t;

	for (ll l=2,r; l<=x; l=r+1)

		r=x/(t=x/l),s=(1ll*S(t)*(s2(r%p)-s2((l-1)%p)+p)+s)%p;

	return (s=s3(x%p)-s)<0&&(s+=p),mps[x]=s;

}

int main() {

	scanf("%d%lld",&p,&n),m=pow(n,2./3);

	phi[1]=1,inv=_power(6,p-2);

	for (int i=2; i<=m; ++i) {

		v[i]||(pr[++cnt]=i,phi[i]=i-1);

		for (int j=1; j<=cnt&&1ll*pr[j]*i<=m; ++j) {

			v[t=pr[j]*i]=1;

			if (i%pr[j]==0) { phi[t]=phi[i]*pr[j]; break; }

			phi[t]=phi[i]*(pr[j]-1);

		}

		phi[i]=(1ll*i*i%p*phi[i]+phi[i-1])%p;

	}

	for (ll l=1,r; l<=n; l=r+1)

		r=n/(tt=n/l),ans=((1ll*S(r)-S(l-1)+p)*s3(tt%p)+ans)%p;

	printf("%d\n",ans);

	return 0;

}

P1891 疯狂LCM

\[\begin{aligned}
\sum_{i=1}^n[i,n]
&=\sum_{i=1}^n\dfrac{ni}{(n,i)} \\
&=n\sum_{i=1}^n\dfrac i{(n,i)} \\
&=n\sum_{d\mid n}\sum_{i=1}^n\dfrac id[(i,n)==d] \\
&=n\sum_{d\mid n}\sum_{i=1}^{n/d}i[(i,n/d)==1] \\
&=n\sum_{d\mid n}\sum_{i=1}^di[(i,d)==1] \\
&=n\sum_{d\mid n}\dfrac{\varphi(d)d}2
\end{aligned}\]

讲一下最后一步

首先满足 \((i,d)==1\) 的 \(i\) 有 \(\varphi(d)\) 个

然后对于每一个满足 \((i,d)==1\) 的 \(i\) ，有 \((d-i,d)==1\)

进而可以证明所有满足条件的 \(i\) 的平均数是 \(\dfrac d2\)

则它们的和为 \(\dfrac{\varphi(d)d}2\)

然后这里注意一个特殊情况， \(d=1\) 时这个是不成立的，答案就是 1

筛完欧拉函数直接算时间复杂度是 \(O(n+T\sqrt n)\)

我们设 \(f_n=\sum_{d\mid n}\dfrac{\varphi(d)d}2\)

可以 \(O(n\log n)\) 完成预处理，然后 \(O(1)\) 处理询问

时间复杂度 \(O(n\log n+T)\)

#include<stdio.h>

const int N=1000001;

int T,t,cnt,n,p[100000],v[N],phi[N]; long long tt,s[N];

int main() {

	scanf("%d",&T),s[1]=phi[1]=1;

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

		++s[i],v[i]||(p[++cnt]=i,phi[i]=i-1);

		for (int j=1; j<=cnt&&p[j]*i<N; ++j) {

			v[t=p[j]*i]=1;

			if (i%p[j]==0) { phi[t]=phi[i]*p[j]; break; }

			phi[t]=phi[i]*(p[j]-1);

		}

		tt=(1ll*phi[i]*i>>1);

		for (int j=i; j<N; j+=i) s[j]+=tt;

	}

	while (T--) scanf("%d",&n),printf("%lld\n",s[n]*n);

	return 0;

}

P5221 Product

\[\begin{aligned}
\prod_{i=1}^n\prod_{j=1}^n\dfrac{[i,j]}{(i,j)}
&= \prod_{i=1}^n\prod_{j=1}^n\dfrac{ij}{(i,j)^2}\\
&= \left(\prod_{i=1}^n\prod_{j=1}^nij\right)\left(\prod_{i=1}^n\prod_{j=1}^n(i,j)\right)^{-2}
\end{aligned}\]

分别计算

\[\prod_{i=1}^n\prod_{j=1}^nij=\prod_{i=1}^ni^nn!=(n!)^{2n}
\]

\[\prod_{i=1}^n\prod_{j=1}^n(i,j)=\prod_{d=1}^nd^{\sum_{i=1}^{n/d}\sum_{j=1}^{n/d}[(i,j)=1]}=\prod_{d=1}^nd^{\sum_{d'=1}^{n/d}\mu(d)(n/dd')^2}
\]

两重整除分块直接上，时间复杂度应该是 \(O(n^{3/4})\) 的

但这里涉及到连续自然数积的幂的计算，所以还是要把 1~n 扫一遍，时间复杂度是 \(O(n)\)

加上线性筛预处理还是 \(O(n)\)

#include<stdio.h>

const int P=104857601,P1=P-1; long long s; bool v[1000010];

int n,cnt,F=1,ans=1,t,t1,f,pr[100000],mu[1000010];

inline int _power(int x,int y,int p) {

	int s=1;

	while (y) (y&1)&&(s=1ll*s*x%p),x=1ll*x*x%p,y>>=1;

	return s;

}

int main() {

	scanf("%d",&n),mu[1]=1;

	for (int i=2; i<=n; ++i) {

		v[i]||(pr[++cnt]=i,mu[i]=-1);

		for (int j=1; j<=cnt&&(t=pr[j]*i)<=n; ++j) {

			v[t]=1;

			if (i%pr[j]==0) break; mu[t]=-mu[i];

		}

		mu[i]+=mu[i-1];

	}

	for (int l=1,r; l<=n; l=r+1) {

		r=n/(t=n/l),s=0,f=1;

		for (int i=l; i<=r; ++i) f=1ll*f*i%P;

		for (int l1=1,r1; l1<=t; l1=r1+1)

			r1=t/(t1=t/l1),s+=1ll*t1*t1*(mu[r1]-mu[l1-1]);

		ans=1ll*ans*_power(f,s%P1,P)%P,F=1ll*F*f%P;

	}

	printf("%d\n",1ll*_power(F,n<<1,P)*_power(ans,P-3,P)%P);

	return 0;

}

P3312 [SDOI2014]数表

数论+数据结构，一道比较好的题

设 \(\sigma(i)\) 为 \(i\) 的约数和，我们可以用线性筛预处理 \(\sigma\)

而我们要求的是 \(\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\sigma((i,j))\le a]\sigma((i,j))\)

假设 \(n\le m\)

设 \(f(i)=\sum\limits_{d\mid i}[\sigma(d)\le a]\sigma(d)\mu(i/d)\) ，推式子可得答案为 \(\sum\limits_{i=1}^n(n/i)(m/i)f(i)\)

发现不同的询问中 \(a\) 是不相等的，在线做就很困难

考虑把所有询问离线，按 \(a\) 从小到大排序，然后我们通过暴力枚举 \(d\) 和 \(i/d\) 更新所有 \(f(d*(i/d))\) 即 \(f(i)\) 的值

通过调和级数可以证明这个过程是 \(O(n\log n)\) 的

但是我们发现回答每个询问的复杂度瓶颈不再是整除分块而是求 \(f\) 的前缀和，也就是线性的

所以我们不维护 \(f(i)\) ，用 BIT 维护 \(f(i)\) 的前缀和

时间复杂度由 \(O(n\log n+T(\sqrt n+n))\) 变成了 \(O(n\log^2n+T\sqrt n\log n)\)

于是就可以通过了

然后这里有一些奇妙的卡常技巧

实际上是我在做完这题后查了别人的代码，然后学会的

但第一次接触这种技巧是在SP26045 虽然当时没看懂

所以这个坑留到SP26045的题解再填

所以代码没写注释就别看了吧

#include<stdio.h>

#include<algorithm>

const int N=100001,inf=2147483647; int T,n,m,A,id,t,cnt,S;

int p[10000],mu[N],d[N],s[N],ans[N],v[N],pw[N]; double inv[N];

inline int min(int x,int y) { return x<y?x:y; }

inline int lb(int x) { return -x&x; }

inline void upd(int i,int x) { for ( ; i<=N; i+=lb(i)) s[i]+=x; }

inline int qry(int i) { for (S=0; i; i-=lb(i)) S+=s[i]; return S; }

struct query { int n,m,A,id; }a[20010];

struct node { int x,d; }b[N];

int operator < (query p,query q) { return p.A<q.A; }

int operator < (node p,node q) { return p.d<q.d; }

int main() {

	scanf("%d",&T),mu[1]=d[1]=pw[1]=1;

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

		v[i]||(p[++cnt]=i,mu[i]=-1,d[i]=i+1,pw[i]=i*i);

		for (int j=1,t1; j<=cnt&&(t=(t1=p[j])*i)<N; ++j) {

			v[t]=1;

			if (i%t1) mu[t]=-mu[i],pw[t]=1ll*t1*t1,d[t]=d[i]*(t1+1);

			else { d[t]=1ll*d[i]*((pw[t]=pw[i]*t1)-1)/(pw[i]-1); break; }

		}

	}

	for (int i=1; i<N; ++i) b[i]={i,d[i]},inv[i]=1./i*(1+1e-9);

	for (int i=1; i<=T; ++i) scanf("%d%d%d",&n,&m,&A),a[i]={n,m,A,i};

	std::sort(a+1,a+T+1),std::sort(b+1,b+N);

	for (int i=1,j=1; i<=T; ++i) {

		for (n=a[i].n,m=a[i].m,A=a[i].A,id=a[i].id; b[j].d<=A&&j<=N; ++j)

			for (int t1=b[j].x,k=1; (t=t1*k)<N; ++k) upd(t,b[j].d*mu[k]);

		for (int l=1,r,t1,mx=min(n,m); l<=mx; l=r+1) {

			r=min(n*inv[t=n*inv[l]],m*inv[t1=m*inv[l]]);

			ans[id]+=t*t1*(qry(r)-qry(l-1));

		}

	}

	for (int i=1; i<=T; ++i) printf("%d\n",ans[i]&inf);

	return 0;

}

SP5971 LCMSUM - LCM Sum

这波啊，这波是双倍经验

上面有讲P1891，两题完全一样

不说了

P4917 天守阁的地板

题意：求 \(\sum\limits_{i=1}^n\sum\limits_{j=1}^n\dfrac{[i,j]^2}{ij}\)

\[\sum_{i=1}^n\sum_{j=1}^n\dfrac{[i,j]^2}{ij}=\sum_{i=1}^n\sum_{j=1}^n\dfrac{ij}{(i,j)^2}
\]

然后这题就变成了 P5221

但是照着那题那么写会跑得很慢，所以我们换个方法推式子

分别计算

\[\prod_{i=1}^n\prod_{j=1}^nij=\prod_{i=1}^ni^nn!=(n!)^{2n}
\]

\[\prod_{i=1}^n\prod_{j=1}^n(i,j)=\prod_{d=1}^nd^{\sum_{i=1}^{n/d}\sum_{j=1}^{n/d}[(i,j)=1]}=\prod_{d=1}^nd^{2\left(\sum_{d'=1}^{n/d}\varphi(d')\right)-1}
\]

预处理 \(\varphi\) 的前缀和，我们就有了简单的 \(O(T\sqrt n\log n)\) 做法

#include<stdio.h>

const int P=19260817,P1=P-1,N=1000001; bool v[N];

int T,n,cnt,ans,t,t1,pr[100000],phi[N],fac[N],inv[N];

inline int _power(int x,int y,int p) {

	int s=1;

	while (y) (y&1)&&(s=1ll*s*x%p),x=1ll*x*x%p,y>>=1;

	return s;

}

int main() {

	scanf("%d",&T),phi[1]=fac[1]=1;

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

		v[i]||(pr[++cnt]=i,phi[i]=i-1);

		for (int j=1; j<=cnt&&(t=pr[j]*i)<N; ++j) {

			v[t]=1;

			if (i%pr[j]) phi[t]=phi[i]*(pr[j]-1);

			else { phi[t]=phi[i]*pr[j]; break; }

		}

		(phi[i]+=phi[i-1])>=P1&&(phi[i]-=P1);

		fac[i]=1ll*fac[i-1]*i%P;

	}

	inv[N-1]=_power(fac[N-1],P-2,P);

	for (int i=N-1; i; --i)

		inv[i-1]=1ll*inv[i]*i%P,(--(phi[i]<<=1))>=P1&&(phi[i]-=P1);

	while (T--) {

		scanf("%d",&n),ans=1;

		for (int l=1,r; l<=n; l=r+1)

			r=n/(t=n/l),

			ans=1ll*ans*_power(1ll*fac[r]*inv[l-1]%P,phi[t],P)%P;

		ans=1ll*_power(fac[n],n<<1,P)*_power(ans,P-3,P)%P;

		printf("%d\n",ans);

	}

	return 0;

}

P3704 [SDOI2017]数字表格

假设 \(n\le m\) ，答案是：

\[\begin{aligned}
&~~~~~\prod_{d=1}^nf[d]^{\sum_{i=1}^n\sum_{j=1}^m[(i,j)=d]} \\
&=\prod_{d=1}^nf[d]^{\sum_{d'=1}^{n/d}\mu(d')(n/dd')(m/dd')} \\
&=\prod_{T=1}^n\left(\prod_{d\mid T}f[d]^{\mu(T/d)}\right)^{(n/T)(m/T)}
\end{aligned}\]

外层整除分块，中间 \(O(n(\log n+\log mod))\) 暴力预处理即可

经过测试，我们发现 \(f[1]\sim f[10^6]\) 都不能被 \(1e9+7\) 整除

这意味着我们还可以使用线性求多个数逆元的科技（这里不说了，右转P5431题解区）

那么预处理的时间复杂度是 \(O(n\log n)\)

时间复杂度 \(O(n\log n+T\sqrt n\log n)\)

#include<stdio.h>

#include<math.h>

const int P=1e9+7,P1=P-1,N=1000001;

int T,n,m,k,t,t1,cnt,S,p[100000]; double a[N];

int mu[N],s[N],f[N],v[N],pre[N],inv[N];

inline int power(int x,int y) {

	int s=1;

	while (y) (y&1)&&(s=1ll*s*x%P),x=1ll*x*x%P,y>>=1;

	return s;

}

inline int min(int x,int y) { return x<y?x:y; }

inline void mul(int &x,int y) { x=1ll*x*y%P; }

int main() {

	scanf("%d",&T),a[1]=mu[1]=f[1]=pre[0]=pre[1]=s[1]=1;

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

		v[i]||(p[++cnt]=i,mu[i]=-1);

		for (int j=1; j<=cnt&&(t=p[j]*i)<N; ++j) {

			v[t]=1;

			if (i%p[j]==0) break; mu[t]=-mu[i];

		}

		(f[i]=f[i-1]+f[i-2])>=P&&(f[i]-=P),s[i]=1;

		pre[i]=1ll*pre[i-1]*f[i]%P,a[i]=1./i*(1+1e-9);

	}

	S=power(pre[N-1],P-2),t=1;

	for (int i=N-1; i; --i) inv[i]=1ll*S*t%P*pre[i-1]%P,mul(t,f[i]);

	for (int i=1; i<N; ++i) for (int j=1; (t=i*j)<N; ++j)

		mul(s[t],mu[j]?((~mu[j])?f[i]:inv[i]):1);

	for (int i=1; i<N; ++i) pre[i]=1ll*pre[i-1]*s[i]%P;

	inv[N-1]=power(pre[N-1],P-2);

	for (int i=N-1; i; --i) inv[i-1]=1ll*inv[i]*s[i]%P;

	while (T--) {

		scanf("%d%d",&n,&m),n>m&&(t=n,n=m,m=t),k=min(n,sqrt(m)),S=1;

		for (int i=1; i<=k; ++i)

			mul(S,power(s[i],(long long)(n*a[i])*(int)(m*a[i])%P1));

		for (int l=k+1,r; l<=n; l=r+1)

			r=min(n*a[t=n*a[l]],m*a[t1=m*a[l]]),

			mul(S,power(1ll*pre[r]*inv[l-1]%P,1ll*t*t1%P1));

		printf("%d\n",S);

	}

	return 0;

}

巴特西