【BZOJ4816】【SDOI2017】数字表格 [莫比乌斯反演]

数字表格

Time Limit: 50 Sec  Memory Limit: 128 MB
[Submit][Status][Discuss]

Description

　　Doris刚刚学习了fibonacci数列。用f[i]表示数列的第i项，那么

　　f[0]=0

　　f[1]=1

　　f[n]=f[n-1]+f[n-2],n>=2

　　Doris用老师的超级计算机生成了一个n×m的表格，第i行第j列的格子中的数是f[gcd(i,j)]，其中gcd(i,j)表示i,j的最大公约数。Doris的表格中共有n×m个数，她想知道这些数的乘积是多少。答案对10^9+7取模。

Input

　　第一个一个数T，表示数据组数。

　　接下来T行，每行两个数n,m

Output

　　输出T行，第i行的数是第i组数据的结果

Sample Input

　　3
　　2 3
　　4 5
　　6 7

Sample Output

　　1
　　6
　　960

HINT

　　T<=1000,1<=n,m<=10^6

Solution

　　运用莫比乌斯反演，得到式子：

　　这样我们对于内外分块即可，复杂度为O(n^(0.75)*T)。

　　然后我们会发现在BZOJ上过不去，怎么办呢？卡常！BearChild运用了如下的卡常技巧：

　　　　1.  读入优化； 2. O(n)预处理逆元； 3. 内嵌汇编实现乘和取模； 4. 记录n/i，避免多次除法。

Code

#include<iostream>  
#include<string>  
#include<algorithm>  
#include<cstdio>  
#include<cstring>  
#include<cstdlib>  
#include<cmath>
using namespace std; 
typedef long long s64;
 
const int ONE = 1e6+5;
const int MOD = 1e9+7;
const int PHI = 1e9+6;
 
int T;
int n,m;
int prime[ONE],p_num,miu[ONE];
int F[ONE];
bool isp[ONE];
s64 Ans;
 
int get() 
{
        int res=1,Q=1;  char c;
        while( (c=getchar())<48 || c>57)
        if(c=='-')Q=-1;
        if(Q) res=c-48; 
        while((c=getchar())>=48 && c<=57) 
        res=res*10+c-48; 
        return res*Q; 
}
 
int Quickpow(int a,int b)
{
        int res = 1;
        while(b)
        {
            if(b&1) res = (s64)res*a%MOD;
            a = (s64)a*a%MOD;
            b>>=1;
        }
        return res;
}
 
void Deal_first(int MaxN)
{
        F[1]=1;
        F[0]=0; for(int i=2; i<=MaxN; i++) F[i] = ((s64)F[i-1]+F[i-2]) % MOD;
        F[0]=1; for(int i=2; i<=MaxN; i++) F[i] = (s64)F[i]*F[i-1] % MOD;
         
        miu[1] = 1;
        for(int i=2; i<=MaxN; i++)
        {
            if(!isp[i])
                prime[++p_num] = i, miu[i] = -1;
            for(int j=1; j<=p_num, i*prime[j]<=MaxN; j++)
            {
                isp[i * prime[j]] = 1;
                if(i % prime[j] == 0)
                {
                    miu[i * prime[j]] = 0;
                    break;
                }
                miu[i * prime[j]] = -miu[i];
            }
            miu[i] += miu[i-1];
        }
}
 
int f(int n,int m)
{
        if(n > m) swap(n,m);
        s64 Ans = 0;
        for(int i=1,j=0; i<=n; i=j+1)
        {
            j = min(n/(n/i), m/(m/i));
            Ans += (s64)(n/i) * (m/i)%PHI * ((s64)(miu[j] - miu[i-1] + PHI)%PHI) % PHI;
            Ans %= PHI;
        }
        return Ans;
}
 
void Solve()
{
        n=get();    m=get();
        if(n > m) swap(n,m);
        Ans = 1;
        for(int i=1,j=0; i<=n; i=j+1)
        {
            j = min(n/(n/i), m/(m/i));
            Ans = Ans * Quickpow( (s64)F[j] * Quickpow(F[i-1],MOD-2) % MOD , f(n/i,m/i) % PHI) % MOD;
        }
        printf("%lld
",Ans);
}
 
int main()
{
        Deal_first(ONE-1);
        T = get();
        while(T--)
            Solve();
        return 0;
}

#include<iostream>  
#include<string>  
#include<algorithm>  
#include<cstdio>  
#include<cstring>  
#include<cstdlib>  
#include<cmath>
using namespace std; 
typedef long long s64;
 
const int ONE = 1e6+2;
const int MOD = 1e9+7;
const int PHI = 1e9+6;
 
int T;
int n,m;
int prime[ONE],p_num,miu[ONE];
int Niyu[ONE];
int F[ONE];
bool isp[ONE];
int Ans;
 
inline int get() 
{
        int res=1,Q=1;  char c;
        while( (c=getchar())<48 || c>57)
        if(c=='-')Q=-1;
        if(Q) res=c-48; 
        while((c=getchar())>=48 && c<=57) 
        res=res*10+c-48; 
        return res*Q; 
}

inline int modmul(const int &a, const int &b,const int &M)
{
        int ret;
        __asm__ __volatile__("	mull %%ebx
	divl %%ecx
" : "=d"(ret) : "a"(a), "b"(b), "c"(M));
        return ret;
}

inline int Quickpow(int a,int b)
{
        int res = 1;
        while(b)
        {
            if(b&1) res = modmul(res,a,MOD);
            a = modmul(a,a,MOD);
            b>>=1;
        }
        return res;
}
 
inline void Deal_first(int MaxN)
{
        F[0]=0; F[1]=1;
        int val=1;
        for(int i=2; i<=MaxN; i++)
        {
            F[i] = F[i-1]+F[i-2];
            if(F[i] >= MOD) F[i] -= MOD;
            val = modmul(val,F[i],MOD);
        }
        Niyu[MaxN] = Quickpow(val, MOD-2);
        for(int i=MaxN-1;i>=1;i--) Niyu[i] = modmul(Niyu[i+1],F[i+1],MOD);
        Niyu[0] = Niyu[1]; F[0]=1;
        for(int i=1; i<=MaxN; i++) F[i] = modmul(F[i],F[i-1],MOD);
         
        miu[1] = 1;
        for(int i=2; i<=MaxN; i++)
        {
            if(!isp[i])
                prime[++p_num] = i, miu[i] = -1;
            for(int j=1; j<=p_num, i*prime[j]<=MaxN; j++)
            {
                isp[i * prime[j]] = 1;
                if(i % prime[j] == 0)
                {
                    miu[i * prime[j]] = 0;
                    break;
                }
                miu[i * prime[j]] = -miu[i];
            }
            miu[i] += miu[i-1];
        }
}
 
inline int f(int n,int m)
{
        if(n > m) swap(n,m);
        int Ans = 0;
        for(int i=1,j=0; i<=n; i=j+1)
        {
            int x=n/i, y=m/i;
            j = min(n/x, m/y);
            Ans = ((s64)Ans + modmul(modmul(x,y,PHI) , ((s64)miu[j] - miu[i-1] + PHI), PHI) )%PHI;
        }
        return Ans;
}
 
inline void Solve()
{
        n=get();    m=get();
        if(n > m) swap(n,m);
        Ans = 1;
        for(int i=1,j=0; i<=n; i=j+1)
        {
            int x=n/i, y=m/i;
            j = min(n/x, m/y);
            Ans = (s64)modmul(Ans , Quickpow( modmul(F[j],Niyu[i-1],MOD) , f(x,y)), MOD);
        }
        printf("%d
",Ans);
}
 
int main()
{
        Deal_first(ONE-1);
        T = get();
        while(T--)
            Solve();
        return 0;
}