【Bzoj4816】【Sdoi2017】数字表格
一眼推公式+莫比乌斯反演
下面来推公式：
(借用dalao的图 http://www.cnblogs.com/SiriusRen/p/6703062.html)
然后就可以分块啦！！
P.S. 被卡TLE了好多次。。。。
后来又新预处理了一些东西才AC的。。
跑得好慢。。
1 #include<cstdio> 2 #include<algorithm> 3 #define ll long long 4 using namespace std; 5 const ll maxn=1000010; 6 const ll mod=1000000007; 7 ll mu[maxn],prime[maxn],p,n,m,T,x[maxn],f[maxn],num[maxn]; 8 bool is_prime[maxn]; 9 ll qpow(ll x,ll y){ 10 if (y<0) y+=(mod-1); 11 ll now=x; 12 ll res=1; 13 while (y) { 14 if (y & 1) { 15 res=res*now; 16 if (res>=mod) res%=mod; 17 } 18 now=now*now; 19 if (now>=mod) now%=mod; 20 y >>= 1; 21 } 22 return res; 23 } 24 void get_mu_f(){ 25 f[0]=0; 26 f[1]=1; 27 num[0]=0; 28 num[1]=1; 29 for (ll i=0; i<=1000000/3; i++) is_prime[i*3]=is_prime[i*3+1]=is_prime[i*3+2]=true; 30 is_prime[0]=is_prime[1]=false; 31 mu[1]=1; 32 x[0]=x[1]=1; 33 for (ll i=2; i<maxn; i++) { 34 f[i]=(f[i-1]+f[i-2]) % mod; 35 num[i]=qpow(f[i],mod-2); 36 x[i]=1; 37 if (is_prime[i]) { 38 prime[++p]=i; 39 mu[i]=-1; 40 } 41 for (ll j=1; j<=p; j++) { 42 ll now=prime[j]; 43 if (now*i>=maxn) break; 44 is_prime[now*i]=false; 45 if (i % now==0){ 46 mu[i*now]=0; 47 break; 48 } 49 else mu[i*now]=-mu[i]; 50 } 51 } 52 } 53 void get_x(){ 54 for (ll i=1; i<maxn; i++) 55 for (ll j=i; j<maxn; j+=i) { 56 if (mu[j/i]==1) x[j]=(x[j]*f[i]) % mod; 57 else if (mu[j/i]==-1) x[j]=(x[j]*num[i]) % mod; 58 //x[j]=(x[j]*qpow(f[i],mu[j / i])) % mod; 59 } 60 for (ll i=1; i<maxn; i++) x[i]=(x[i]*x[i-1]) % mod; 61 } 62 int main(){ 63 scanf("%lld",&T); 64 get_mu_f(); 65 get_x(); 66 while (T--) { 67 scanf("%lld%lld",&n,&m); 68 ll j=0; 69 ll ans=1; 70 if (n>m) swap(n,m); 71 for (ll i=1; i<=n; i=j+1){ 72 j=min(n/(n/i),m/(m/i)); 73 ans=(ans*qpow((x[j]*qpow(x[i-1],mod-2))%mod,((n/i)*(m/i))%(mod-1)))%mod; 74 } 75 printf("%lld\n",ans); 76 } 77 return 0; 78 }
↑不要介意这丑陋的代码。。
