【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 } 

↑不要介意这丑陋的代码。。

posted @ 2017-04-27 23:14  MoerBlack  阅读(226)  评论(0编辑  收藏  举报