1. 程式人生 > >bzoj 4816 數字表格 —— 反演

bzoj 4816 數字表格 —— 反演

題目:https://www.lydsy.com/JudgeOnline/problem.php?id=4816

推導過程同:http://www.cnblogs.com/zhouzhendong/p/8666106.html

可以利用階乘處理逆元。

程式碼如下:

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
int const xn=1e6+5,mod=1e9+7;
int f[xn],sf[xn],mu[xn],cnt,pri[xn],g[xn],sg[xn],invf[xn],invg[xn];
bool vis[xn]; int rd() { int ret=0,f=1; char ch=getchar(); while(ch<'0'||ch>'9'){if(ch=='-')f=0; ch=getchar();} while(ch>='0'&&ch<='9')ret=ret*10+ch-'0',ch=getchar(); return f?ret:-ret; } ll pw(ll a,ll b) { ll ret=1; a=a%mod; b=b%(mod-1); if(b<0)b+=(mod-1); for(;b;b>>=1ll,a=(a*a)%mod)if
(b&1)ret=(ret*a)%mod; return ret; } int upt(int x){while(x>=mod)x-=mod; while(x<0)x+=mod; return x;} void init() { int mx=1e6; mu[1]=1; for(int i=2;i<=mx;i++) { if(!vis[i])pri[++cnt]=i,mu[i]=-1; for(int j=1;j<=cnt&&(ll)i*pri[j]<=mx;j++) { vis[i
*pri[j]]=1; if(i%pri[j])mu[i*pri[j]]=-mu[i]; else break; } } f[1]=1; sf[0]=1; for(int i=2;i<=mx;i++)f[i]=upt(f[i-1]+f[i-2]); for(int i=1;i<=mx;i++)sf[i]=(ll)f[i]*sf[i-1]%mod; int jcn=pw(sf[mx],mod-2); for(int i=mx;i;i--) invf[i]=(ll)jcn*sf[i-1]%mod,jcn=(ll)jcn*f[i]%mod;//inv(f) for(int i=0;i<=mx;i++)g[i]=1; for(int i=1;i<=mx;i++) for(int j=i;j<=mx;j+=i) { if(mu[j/i]==1)g[j]=(ll)g[j]*f[i]%mod; else if(mu[j/i]==-1)g[j]=(ll)g[j]*invf[i]%mod; } sg[0]=1; for(int i=1;i<=mx;i++)sg[i]=(ll)g[i]*sg[i-1]%mod; invg[mx]=pw(sg[mx],mod-2); for(int i=mx-1;i>=0;i--)invg[i]=(ll)invg[i+1]*g[i+1]%mod;//inv(sg) //>=0 } int main() { int T=rd(); init(); while(T--) { int n=rd(),m=rd(); int mn=min(n,m),ans=1; for(int i=1,j;i<=mn;i=j+1) { j=min(n/(n/i),m/(m/i)); int tmp=(ll)sg[j]*invg[i-1]%mod; tmp=pw(tmp,(ll)(n/i)*(m/i)); ans=(ll)ans*tmp%mod; } printf("%d\n",ans); } return 0; }