「SDOI2017」数字表格 莫比乌斯反演

题目链接:https://loj.ac/problem/2000
QQ截图20190402205913.jpg
到这里为止时间复杂度显然是$O(TN^{2})$可以获得30分。
那要AC我们就得继续对上式进行变形。设$i=dk$并带入下式并变形得
QQ截图20190402205923.jpg
显然这个表达式可以O(NlogN)的求出前缀积。然后我们根据这个表达式的性质再使用整数分块,然后就可以在$O(T\sqrt{N}logN)$内求出答案。
代码:

#include <bits/stdc++.h>

const int MAXN = 1000000;
const long long P = 1e9+7;

using namespace std;

int T,n,m,cnt;
int prime[MAXN+5],mu[MAXN+5];
long long s[MAXN+5],g[MAXN+5],invg[MAXN+5],f[MAXN+5],invf[MAXN+5];
bool vis[MAXN+5];

long long pow_mod(long long a,long long n){
    long long ans = 1;
    while (n) {
        if (n&1) ans = ans * a % P;
          a = a * a % P;
          n >>= 1;
     }
     return ans % P;
}

inline void init()
{
    mu[1] = f[1] = g[1] = invf[1] = g[0] = invg[0] = 1;
    for(int i = 2;i <= MAXN;i++){
        g[i] = 1;
        f[i] = (f[i-1] + f[i-2]) % P;
        invf[i] = pow_mod(f[i],P-2);
     }
    
    for(int i = 2;i <= MAXN;i++){
        if(!vis[i]){
            prime[++cnt] = i;
            mu[i] = -1;
        }
        for(int j = 1;j <= cnt && prime[j] * i <= MAXN;j++){
            vis[prime[j]*i] = true;
            if(i%prime[j] == 0){
                mu[i*prime[j]] = 0;
                break;    
            }
            else mu[i*prime[j]] = -mu[i];
        }
    }
        
    for(int i = 1;i <= MAXN;i++){
        for(int j = i;j <= MAXN;j += i){            
            if(mu[j/i] ==  1) g[j] = g[j]*f[i]%P;
            if(mu[j/i] == -1) g[j] = g[j]*invf[i]%P;
        }
    }

    for(int i = 1;i <= MAXN;i++){
        g[i] *= g[i-1];g[i] %= P;
        invg[i] = pow_mod(g[i],P-2);        
    }
}

inline long long cal()
{
    long long ret = 1;
    int last;
    if(n > m) swap(n,m);
    for(int i = 1;i <= n;i = last + 1){
        last = min(n/(n/i),m/(m/i));
        ret = ret * pow_mod(g[last]*invg[i-1]%P,1LL*(n/i)*(m/i)%(P-1)) % P;
    }
    return ret;
}

int main()
{    
    init();
        
    scanf("%d",&T);    
    while(T--){
        scanf("%d %d",&n,&m);
        printf("%lld\n",cal());
    } 
    return 0;
}
Last modification:April 2nd, 2019 at 09:00 pm
If you think my article is useful to you, please feel free to appreciate

Leave a Comment