题解 P3704 【[SDOI2017]数字表格】

原题地址

一句话题意

求\(\begin{aligned}\prod_{i=1}^n\prod_{j=1}^mf_{\gcd(i,j)}\end{aligned}\)，多组询问

解题思路

考虑更换枚举\(gcd\)的方式，由于是连乘\(\prod\)处理起来可能有点抽象，不像\(\sum\)一样非常好理解

\(\begin{aligned}\prod_{i=1}^n\prod_{j=1}^mf_{\gcd(i,j)}\end{aligned}\)

\(\begin{aligned}=\prod_{p=1}^nf_p^{\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=p]}\end{aligned}\)，这里就是在更换方式枚举\(gcd\)，对于每一个\(n\)，由于是\(\prod\)，对于之前的\(n\)也要再次进行计算

\(\begin{aligned}=\prod_{p=1}^nf_p^{\sum_{i=1}^\frac np\sum_{j=1}^\frac mp[\gcd(i,j)=1]}\end{aligned}\)，反演老套路，对幂上的柿子进行反演

\(\begin{aligned}=\prod_{p=1}^nf_p^{\sum_{i=1}^{\frac np}\sum_{j=1}^{\frac mp}\sum_{d|i,d|j}\mu(d)}\end{aligned}\)

\(\begin{aligned}=\prod_{p=1}^nf_p^{\sum_{d=1}^n\mu(d)\lfloor\frac n{pd}\rfloor\lfloor\frac m{pd}\rfloor}\end{aligned}\)

然后推到到这里按照反演老套路该\(dp\)化\(Q\)然后整除分块直接直播切题了，但是怎么化\(Q\)？？？？？？？歪？

\(\begin{aligned}\prod_{Q=1}^n\prod_{d|Q}\end{aligned}\)按照老套路，一定会有这个东西，在思考如何去处理\(f\)以及\(\mu\)我们更换枚举之后对于枚举\(d\)的时候不会对\(Q\)产生影响，把分块的指数提出来，然后改一下下标就可以了

\(\begin{aligned}=\prod_{Q=1}^n\left(\prod_{d|Q}f_{d/Q}^{\mu(d)}\right)^{\lfloor\frac n{Q}\rfloor\lfloor\frac m{Q}\rfloor}\end{aligned}\)

对于里面那一堆东西\(\begin{aligned}\prod_{d|Q}f_{d/Q}^{\mu(d)}\end{aligned}\)我们直接预处理然后分块计算就行了，复杂度为\(O(n+Tlog n\sqrt n)\)

#include<bits/stdc++.h>
    
#define LL long long
    
#define _ 0

#define int long long

using namespace std;
    
/*Grievous Lady*/
    
template <typename T> void read(T & t){
    t = 0;int f = 1;char ch = getchar();
    while(ch < '0' || ch > '9'){if(ch == '-')f =- 1;ch = getchar();}
    do{t = t * 10 + ch - '0';ch = getchar();}while(ch >= '0' && ch <= '9');t *= f;
}
    
#define mod 1000000007

const int kato = 1e7 + 10;

const int luna = 1e6;

inline LL quick_pow(LL a , LL b){
    LL res = 1;
    for(; b ; b >>= 1 , a = a * a % mod){
        if(b & 1){
            res = res * a % mod;
        }
    }
    return res % mod;
}

int prime[kato] , cnt , f[kato] , invf[kato] , yuni[kato] , mu[kato];

bool ispri[kato];

inline void get_yuni(){
    for(int i = 2;i <= luna;i ++){
        if(!ispri[i]){
            prime[++ cnt] = i;
            mu[i] = -1;
        }
        for(int j = 1;j <= cnt && i * prime[j] <= luna;j ++){
            ispri[i * prime[j]] = 1;
            if(i % prime[j] == 0){
                break;
            }
            mu[i * prime[j]] = -mu[i];
        }
    }
    for(int i = 2;i <= luna;i ++){
        f[i] = (f[i - 1] + f[i - 2]) % mod;
        invf[i] = quick_pow(f[i] , mod - 2);
        yuni[i] = 1;
    }
    for(int d = 1;d <= luna;d ++){
        if(mu[d] == 1){
            for(int Q = d , cnt = 1;Q <= luna;Q += d , cnt ++){
                yuni[Q] = yuni[Q] * f[cnt] % mod;
            }
        }
        if(mu[d] == -1){
            for(int Q = d , cnt = 1;Q <= luna;Q += d , cnt ++){
                yuni[Q] = yuni[Q] * invf[cnt] % mod;
            }
        }
    }
    for(int i = 1;i <= luna;i ++){
        yuni[i] = yuni[i] * yuni[i - 1] % mod;
    }
}

LL n , m , t;

inline int Ame_(){
    read(t);
    mu[1] = f[1] = invf[1] = yuni[1] = yuni[0] = 1;
    get_yuni();
    // for(int i = 1;i <= 10;i ++){
    //     cerr << mu[i] << ' '; 
    // }
    for(; t --> 0 ;){
        LL ans = 1;
        read(n) , read(m);
        if(n > m) swap(n , m);
        for(LL l = 1 , r;l <= n;l = r + 1){
            r = min(n / (n / l) , m / (m / l));
            ans = ans * quick_pow(yuni[r] * quick_pow(yuni[l - 1] , mod - 2) % mod , (n / l) * (m / l)) % mod;
        }
        printf("%lld\n" , ans);
    }
    // fclose(stdin);
    // fclose(stdout);
    return ~~(0^_^0);
}
    
int Ame__ = Ame_();
    
signed main(){;}