题解 P3704 【[SDOI2017]数字表格】
qwaszx
·
·
题解
哭了，卡了两个小时常了，Kelin为什么这么nb啊明明我复杂度比他低还没他跑得快
现在我比他快了
首先化式子.设n\leq m，除法均为整除.
\begin{array}{lcc}\prod\limits_{i=1}^n\prod\limits_{j=1}^mf[gcd(i,j)]=\prod\limits_{d=1}^n\prod\limits_{i=1}^n\prod\limits_{j=1}^m(f[d])^{[gcd(i,j)]}\\=\prod\limits_{d=1}^n(f[d])^{\sum\limits_{i=1}^n\sum\limits_{j=1}^m[gcd(i,j)=d]}\\=\prod\limits_{d=1}^n(f[d])^{\sum\limits_{k=1}^\frac{n}{d}\mu(k)\frac{n}{kd}\frac{m}{kd}}\\=\prod\limits_{T=1}^n\left(\prod\limits_{d|T}(f[d])^{\mu(\frac{T}{d})}\right)^{\frac{n}{T}\frac{m}{T}}\end{array}
于是我们预处理\prod\limits_{d|T}(f[d])^{\mu(\frac{T}{d})}即可O(\sqrt{n}\log n)回答询问了
如何预处理呢?
绝大多数人都是枚举d，然后更新T
然而事实上我们有更好的做法.参考rqy的blog.
设g_{i,n}=\prod\limits_{d|T,d\text{只含前i种质因子}}(f[\frac{T}{d}])^{\mu(d)}
和inv_{i,n}=\dfrac{1}{g_{i,n}}=\prod\limits_{d|T,d\text{只含前i种质因子}}(f[\frac{T}{d}])^{-\mu(d)}
于是枚举p_i.
如果p_i\not|n，那么g_{i,n}=g_{i-1,n},inv_{i,n}=inv_{i-1,n};
否则，g_{i,n}=g_{i-1,n}\times\prod\limits_{dp_i|T}(f[\frac{\frac{T}{p_i}}{d}])^{\mu(dp_i)}=g_{i-1,n}\times\prod\limits_{dp_i|T}(f[\frac{\frac{T}{p_i}}{d}])^{-\mu(d)}=g_{i-1,n}\times inv_{i-1,\frac{n}{p_i}}
同理有inv_{i,n}=inv_{i-1,n}\times g_{i-1,\frac{n}{p_i}}
于是我们就可以枚举素数来更新了.滚动数组消掉第一维即可.
于是我们把跑g和跑inv得时间从n\log mod降到了n\log \log n
可是还有一个问题!g和inv的初值分别是f和inv(f).inv(f)如果硬算还有n\log mod的复杂度，这远远大于刚才的复杂度.
于是我观摩了Kelin的代码冥思苦想一番之后想起了一种神仙的做法.
我们来考虑f的前缀积sf，以及invf的前缀积si.
那么我们有invf(n)=si(n)\times sf(n-1)以及si(n-1)=si(n)\times f(n)
于是只要知道了si(N)我们就可以O(n)递推出所有invf，而si(N)=sf(N)^{-1}，跑一次逆元求出来即可.
这个技巧还有很多用处emm它可以线性求数论函数逆元.
于是预处理的复杂度只剩n\log \log n了我就不信哪个神仙能线筛
```cpp
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#define mod 1000000007
using namespace std;
const int N=1e6+5;
int s[N],invs[N];
int m,n,T,p[N],prime[100000],mu[N],cnt,nn[1005],mm[1005],sf[N];
inline int qpower(int a,int b)
{
int ans=1;for(;b;b>>=1,a=1ll*a*a%mod)if(b&1)ans=1ll*ans*a%mod;return ans;
}
void make(int n)
{
mu[1]=p[1]=invs[0]=s[1]=invs[1]=sf[1]=1;
for(int i=2;i<=n;i++)
{
s[i]=(s[i-1]+s[i-2])%mod,sf[i]=1ll*sf[i-1]*s[i]%mod;if(!p[i])prime[++cnt]=i,mu[i]=-1;
for(int j=1,x;j<=cnt&&(x=i*prime[j])<=n;j++)
{
p[x]=1;if(i%prime[j])mu[x]=-mu[i];
else break;
}
}
int sis=qpower(sf[n],mod-2);
for(int i=n;i>=2;i--)invs[i]=1ll*sis*sf[i-1]%mod,sis=1ll*sis*s[i]%mod;
for(int i=1;i<=cnt;i++)
{
for(int jj=n/prime[i],j=jj*prime[i];jj>=1;jj--,j=jj*prime[i])
{
s[j]=1LL*s[j]*invs[jj]%mod;
invs[j]=1ll*invs[j]*s[jj]%mod;
}
}
for(int i=2;i<=n;i++)s[i]=1ll*s[i]*s[i-1]%mod,invs[i]=1ll*invs[i]*invs[i-1]%mod;;
}
int main()
{
scanf("%d",&T);int maxn=0;
for(int i=1;i<=T;i++)scanf("%d%d",nn+i,mm+i),maxn=max(maxn,min(nn[i],mm[i]));
make(maxn);
for(int TTT=1;TTT<=T;TTT++)
{
int n=nn[TTT],m=mm[TTT];
if(n>m)swap(n,m);int ans=1;
int i=1,lt=sqrt(n);
for(;i<=lt;i++)
{
ans=(1ll*ans*qpower(1ll*s[i]*invs[i-1]%mod,1ll*(n/i)*(m/i)%(mod-1)))%mod;
}
for(;i<=n;i=lt+1)
{
lt=min(n/(n/i),m/(m/i));
ans=(1ll*ans*qpower(1ll*s[lt]*invs[i-1]%mod,(n/i)*(m/i)))%mod;
}
printf("%d\n",ans);
}
}
```
~~卡常是无止境的~~
```