本文版权归ljh2000和博客园共有,欢迎转载,但须保留此声明,并给出原文链接,谢谢合作。
本文作者:ljh2000
作者博客:http://www.cnblogs.com/ljh2000-jump/
转载请注明出处,侵权必究,保留最终解释权!
题目链接:BZOJ4816
正解:莫比乌斯反演
解题报告:
这类题都是套路…
不过这题有所不同的是由于是连乘,所以莫比乌斯函数应该丢到指数去,保证不为1的时候指数为0就好了。
推到了一个看上去貌似是$O(Tnlogn)$的式子…
其实复杂度有保证,应该是$O(Tn^{\frac{3}{4}})$。
记忆化一下,卡卡常,能跑过去。
正解还能往下推推,很明显的就是可以设新变量代替原来的乘积,老套路了…
预处理逆元之后可以做到$O(T\sqrt{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\n",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;
}
//有志者,事竟成,破釜沉舟,百二秦关终属楚;苦心人,天不负,卧薪尝胆,三千越甲可吞吴。