BZOJ4816 [Sdoi2017]数字表格

本文版权归ljh2000和博客园共有，欢迎转载，但须保留此声明，并给出原文链接，谢谢合作。

本文作者：ljh2000
作者博客：http://www.cnblogs.com/ljh2000-jump/
转载请注明出处，侵权必究，保留最终解释权！

题目链接：BZOJ4816

正解：莫比乌斯反演

解题报告：

　　这类题都是套路…

　　不过这题有所不同的是由于是连乘，所以莫比乌斯函数应该丢到指数去，保证不为1的时候指数为0就好了。

　　推到了一个看上去貌似是$O(Tnlogn)$的式子…
　　其实复杂度有保证，应该是$O(Tn^{frac{3}{4}})$。

　　记忆化一下，卡卡常，能跑过去。

　　正解还能往下推推，很明显的就是可以设新变量代替原来的乘积，老套路了…

　　预处理逆元之后可以做到$O(Tsqrt{n})$，预处理复杂度是$O(nlogn)$的，每次暴力$for$倍数就好了。

//It is made by ljh2000
//有志者，事竟成，破釜沉舟，百二秦关终属楚；苦心人，天不负，卧薪尝胆，三千越甲可吞吴。
#include <algorithm>
#include <iostream>
#include <cstring>
#include <vector>
#include <cstdio>
#include <string>
#include <queue>
#include <cmath>
#include <ctime>
#define lc root<<1
#define rc root<<1|1
#define rep(i,j,k) for(int i=j;i<=k;i++)
#define reg(i,x) for(int i=first[x];i;i=next[i])
using namespace std;
typedef long long LL;
const int mod = 1000000007;
const int MAXN = 1000011;
int n,m,savn[1011],savm[1011],maxn,prime[MAXN],cnt;
LL mobius[MAXN],fib[MAXN],g[MAXN],nf[MAXN],ng[MAXN],ans;
bool vis[MAXN];
inline LL fast_pow(LL x,LL y){ LL r=1; while(y>0) { if(y&1) r=r*x%mod; x=x*x%mod; y>>=1; } return r; }
inline int getint(){
    int w=0,q=0; char c=getchar(); while((c<'0'||c>'9') && c!='-') c=getchar();
    if(c=='-') q=1,c=getchar(); while (c>='0'&&c<='9') w=w*10+c-'0',c=getchar(); return q?-w:w;
}

inline void Init(){
	mobius[1]=1;
	for(int i=2;i<=maxn;i++) {
		if(!vis[i]) { prime[++cnt]=i; mobius[i]=-1; }
		for(int j=1;j<=cnt && prime[j]*i<=maxn;j++) {
			vis[i*prime[j]]=1;
			if(i%prime[j]==0) { mobius[i*prime[j]]=0; break; }
			mobius[i*prime[j]]=-mobius[i];
		}
	}
	fib[0]=0; fib[1]=1;
	for(int i=2;i<=maxn;i++) fib[i]=fib[i-1]+fib[i-2],fib[i]%=mod;
	for(int i=1;i<=maxn;i++) nf[i]=fast_pow(fib[i],mod-2),g[i]=1;
	for(int i=1;i<=maxn;i++)
		for(int j=i;j<=maxn;j+=i) {
			//乘号...
			if(mobius[j/i]==1) g[j]*=fib[i];
			else if(mobius[j/i]==-1) g[j]*=nf[i];
			g[j]%=mod;
		}
	g[0]=1; for(int i=1;i<=maxn;i++) g[i]=g[i-1]*g[i]%mod;
	ng[0]=1; for(int i=1;i<=maxn;i++) ng[i]=fast_pow(g[i],mod-2);
}

inline void solve(int n,int m){
	ans=1; int nexn,nexm,nex; LL now;
	for(int x=1;x<=n;x=nex+1) {
		nexn=n/(n/x); nexm=m/(m/x);
		nex=min(nexn,nexm); nex=min(nex,n);
		now=1LL*g[nex]%mod*ng[x-1]%mod;
		ans*=fast_pow(now,1LL*(n/x)*(m/x));//是乘呀...
		ans%=mod;
	}
	printf("%lld
",ans);
}

inline void work(){
	int T=getint(); for(int o=1;o<=T;o++) { savn[o]=getint(),savm[o]=getint(); if(savn[o]>savm[o]) swap(savn[o],savm[o]); }
	for(int o=1;o<=T;o++) maxn=max(maxn,savn[o]);
	Init(); 
	for(int o=1;o<=T;o++)
		solve(savn[o],savm[o]);
}

int main()
{
#ifndef ONLINE_JUDGE
	freopen("product.in","r",stdin);
	freopen("product.out","w",stdout);
#endif
    work();
    return 0;
}
//有志者，事竟成，破釜沉舟，百二秦关终属楚；苦心人，天不负，卧薪尝胆，三千越甲可吞吴。