[BZOJ4816][Sdoi2017]数字表格数学

最新推荐文章于 2019-04-24 23:00:30 发布

HbFS-

最新推荐文章于 2019-04-24 23:00:30 发布

阅读量417

点赞数

分类专栏：数学

本文链接：https://blog.csdn.net/di4CoveRy/article/details/72821661

版权

数学专栏收录该内容

23 篇文章 0 订阅

订阅专栏

考虑每个 $f_i$ 对答案的贡献，就能得到式子

\sum k = 1 n f \sum d | k μ (d) ⌊ n k d ⌋ ⌊ m k d ⌋ k

$\sum_{k=1}^nf_k^{\sum_{d|k}\mu (d)\lfloor \frac{n}{kd}\rfloor \lfloor \frac{m}{kd}\rfloor}$
转化成枚举

kd $kd$ ，令

T=kd $T=kd$

\prod T = 1 n (\prod d | T f μ (T d) d) ⌊ n T ⌋ ⌊ m T ⌋

$\prod_{T=1}^n(\prod_{d|T}f_d^{\mu(\frac{T}{d})})^{\lfloor \frac{n}{T}\rfloor \lfloor \frac{m}{T}\rfloor}$
预处理小括号内的部分，令

gT=∏d|Tfμ(Td)d $g_T=\prod_{d|T}f_d^{\mu(\frac{T}{d})}$
单次询问只需要回答

\prod T = 1 n g ⌊ n T ⌋ ⌊ m T ⌋ T

$\prod_{T=1}^ng_T^{\lfloor \frac{n}{T}\rfloor \lfloor \frac{m}{T}\rfloor}$
易证

⌊nT⌋ $\lfloor \frac{n}{T}\rfloor$ 和

⌊mT⌋ $\lfloor \frac{m}{T}\rfloor$ 只有

O(N−−√) $O(\sqrt{N})$ 种取值，分段计数即可
预处理

O(NN−−√) $O(N\sqrt{N})$ ，询问

O(TN−−√log2N) $O(T\sqrt{N}log_2N)$ ，log来自快速幂

#include <bits/stdc++.h>
#define mod 1000000007
#define N 1000050
#define tp 1000000
using namespace std;
typedef long long LL;
inline int rd() {int r;scanf("%d",&r);return r;}

bool np[N];
int eu[N],pr[N],g[N],fib[N],sum[N],siv[N],cnt,n,m,ans;
inline int qp(int a,long long b) {
    int ret = 1;
    while (b) {
        if (b&1) ret = 1LL * ret * a % mod;
        b >>= 1, a = 1LL * a * a % mod;
    }
    return ret;
}

void preset() {
    eu[1] = 1;
    for (int i=2;i<=tp;i++) {
        if (!np[i]) pr[++cnt] = i, eu[i] = -1;
        for (int j=1;j<=cnt && pr[j]*i<=tp;j++) {
            np[i*pr[j]] = 1;
            if (i%pr[j] == 0) {
                eu[i*pr[j]] = 0;
                break;
            }
            eu[i*pr[j]] = -eu[i];
        }
    }   


    fib[0] = 0, fib[1] = 1;
    for (int i=2;i<=tp;i++) fib[i] = (fib[i-1] + fib[i-2]) % mod;
    for (int i=1;i<=tp;i++) g[i] = 1;   

    sum[0] = siv[0] = 1;
    for (int i=1;i<=tp;++i) {
        int cur = fib[i];
        int inv = qp(cur, mod-2);
        for (int j=1;i*j<=tp;++j) {
            if (eu[j] == 0) continue;
            g[i*j] = eu[j] == 1 ?
                1LL * g[i*j] * cur % mod:
                1LL * g[i*j] * inv % mod;
        }
        sum[i] = 1LL * sum[i-1] * g[i] % mod;
        siv[i] = 1LL * siv[i-1] * qp(g[i], mod-2) % mod;
    }

}

void solve() {
    n = rd(), m = rd(), ans = 1;
    if (n>m) swap(n,m);
    for (int i=1,lst,cur;i<=n;i=lst+1) {
        lst = min(n/(n/i), m/(m/i));
        cur = 1LL * sum[lst] * siv[i-1] % mod;
        ans = 1LL * ans * qp(cur, 1LL*(n/i)*(m/i)) % mod;
    }
    printf("%d\n",ans);
}

int main() {
    freopen("product.in","r",stdin);
    freopen("product.out","w",stdout);
    preset();
    for (int T=rd();T;T--) solve();
    return 0;
}