【BZOJ 4816】【SDOI 2017】数字表格

考虑到 gcd(i,j) gcd(i,j)的这个形式是一个常见的莫比乌斯反演，尝试构造函数。

1、枚举gcd，转化为 ans=∏nd=1f(d)h(d) ans = \prod_{d=1}^{n}f(d)^{h(d)}，其中 h(d)=∑ni=1∑mj=1[gcd(i,j)==d] h(d)=\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)==d]；
2、 h(d) h(d)是莫比乌斯反演的模板题，直接上结论： h(d)=∑⌊nd⌋k=1μ(k)⌊ndk⌋⌊mdk⌋ h(d)=\sum _{k=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\mu(k)\left \lfloor \frac{n}{dk} \right \rfloor\left \lfloor \frac{m}{dk} \right \rfloor，所以转化为 ans=∏nd=1f(d)h(d)=∑⌊nd⌋k=1μ(k)⌊ndk⌋⌊mdk⌋ ans = \prod_{d=1}^{n}f(d)^{h(d)=\sum _{k=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\mu(k)\left \lfloor \frac{n}{dk} \right \rfloor\left \lfloor \frac{m}{dk} \right \rfloor}；
做到这里通过枚举约数的方法就能够得到答案，这个方法的时间复杂度是 O((n+m)34) O((n+m)^{\frac{3}{4}})，但是乘上一个 T T就有可能超时。这样至少就是60分。
3、设T=dkT=dk，并且改变一下运算顺序（因为全部都是乘法），转换为 ans=∏nT=1(∏d|Tf(d)μ(Td))⌊nT⌋⌊mT⌋ ans=\prod _{T=1}^{n}(\prod _{d|T}f(d)^{\mu(\frac{T}{d})})^{\left \lfloor \frac{n}{T} \right \rfloor\left \lfloor \frac{m}{T} \right \rfloor}；
4、记 g(x)=∏d|Tf(d)μ(Td) g(x)=\prod _{d|T}f(d)^{\mu(\frac{T}{d})}，现在如果我们已经知道了 g(x) g(x)的前缀积，我们就可以在 O(N−−√logN) O(\sqrt {N}logN)的复杂度内求出每次询问，其中根号是分块的近似复杂度， log log是快速幂的复杂度；
5、那么如何求 g(x) g(x)呢？可以通过枚举d，对于每个满足 d|T d|T的 g(T) g(T)都乘上 f(d)μ(Td) f(d)^{\mu(\frac{T}{d})}，只需提前预处理 μ(x) \mu(x)、 f(x) f(x)和 f(x) f(x)的逆元即可，根据调和级数的相关证明可以得到这部分的时间复杂度是 O(lnNlogN) O(lnNlogN)；
6、那么如何求 f(x) f(x)的逆元呢？显然不能 O(P) O(P)求出所有数的逆元。可以选择扩欧来做，时间复杂度为 O(NlnN) O(NlnN)，不过网上有 O(N) O(N)的做法：设 Si=∏i1ai S_{i}=\prod_{1}^{i}a_i，然后求出 Sn S_n的逆元 invSn invS_n，然后通过递推式 invSi−1=ai∗invSi invS_{i-1}=a_{i}*invS_{i}求出 invS invS数组，再通过 invai−1=invSi∗Si−1 inva_{i-1}=invS_{i}*S_{i-1}就可以求出所有的逆元；
7、至此所有的问题都已经解决了，注意最后快速幂求答案的时候指数模一下mod-2不然会有神奇的错误。。。

#include<cmath>
#include<cstdio>
#include<vector>
#include<queue>
#include<cstring>
#include<iomanip>
#include<stdlib.h>
#include<iostream>
#include<algorithm>
#define ll long long
#define inf 1000000000
#define mo 1000000007
#define N 1000000
#define fo(i,a,b) for(i=a;i<=b;i++)
#define fd(i,a,b) for(i=a;i>=b;i--)
using namespace std;
int mbs[N+5],p[N],vis[N+5];
int f[N+5],Sf[N+5],invf[N+5],invS[N+5];
int g[N+5],Sg[N+5],invg[N+5];
int T,n,m,res,i,last;
int quick_multi(int x,int a)
{int res = 1,t = x;while (a) {if (a&1) res=(1ll*res*t)%mo; t=(1ll*t*t)%mo; a>>=1;}return res;
}
void mobius()
{int i,j; mbs[1] = 1;fo(i,2,N){if (!vis[i]) {p[++p[0]] = i; mbs[i] = -1;}for (j = 1;j <= p[0] && i * p[j] <= N; j++){vis[i*p[j]] = 1;if (i % p[j] == 0) {mbs[i*p[j]] = 0; break;}mbs[i*p[j]] = - mbs[i];}}
}
void fibonacci()
{int i,j; f[0] = 0; f[1] = 1; Sf[1] = 1;fo(i,2,N){f[i] = f[i-1] + f[i-2];f[i] -= f[i]>mo?mo:0;Sf[i] = 1ll*Sf[i-1]*f[i]%mo; }invS[N] = quick_multi(Sf[N],mo-2); invf[1] = 1;fd(i,N,2) invS[i-1] = 1ll*invS[i]*f[i] % mo;fd(i,N,2) invf[i] = 1ll*Sf[i-1]*invS[i] % mo;invf[0] = 1;
}
void pre_g()
{int i,d;fo(i,1,N) g[i] = 1;fo(d,1,N){for (i = 1;i * d <= N; i++)if (mbs[i] == 1) g[i*d] = 1ll*g[i*d]*f[d]%mo; elseif (mbs[i] == -1) g[i*d] = 1ll*g[i*d]*invf[d]%mo;}Sg[1] = 1;fo(i,2,N) Sg[i] = 1ll*Sg[i-1]*g[i]%mo;invg[N] = quick_multi(Sg[N],mo-2);fd(i,N,2) invg[i-1] = 1ll*invg[i]*g[i]%mo;invg[0] = 1;
}
int main()
{mobius();fibonacci();pre_g();/*n = 100;fo(i,1,n) cout<<invg[i]<<" "; cout<<endl;*/scanf("%d",&T);while (T--){scanf("%d%d",&n,&m);if (n > m) swap(n,m);res = 1;for (i = 1;i <= n; i=last+1){last = min(n/(n/i),m/(m/i));res = 1ll*res*quick_multi(1ll*Sg[last]*invg[i-1]%mo,1ll*(n/i)*(m/i)%(mo-1))%mo;}printf("%d\n",res);}return 0;
}

【BZOJ 4816】【SDOI 2017】数字表格相关推荐

BZOJ 2154 Crash的数字表格 (莫比乌斯反演)
Crash的数字表格今天的数学课上,Crash小朋友学习了最小公倍数(Least Common Multiple).对于两个正整数a和b,LCM(a, b)表示能同时被a和b整除的最小正整数.例如, ...
[BZOJ 4819] [SDOI 2017] 新生舞会
Description 学校组织了一次新生舞会,Cathy作为经验丰富的老学姐,负责为同学们安排舞伴. 有 \(n\) 个男生和 \(n\) 个女生参加舞会买一个男生和一个女生一起跳舞,互为舞伴. C ...
【BZOJ】【2154】Crash的数字表格
莫比乌斯反演 PoPoQQQ讲义第4题题解:http://www.cnblogs.com/jianglangcaijin/archive/2013/11/27/3446169.html 感觉两次sq ...
BZOJ 2154 [国家集训队]Crash的数字表格 / JZPTAB（莫比乌斯反演，经典好题）（Luogu P1829）
整理的算法模板合集: ACM模板点我看算法全家桶系列!!! 实际上是一个全新的精炼模板整合计划 P1829 [国家集训队]Crash的数字表格 / JZPTAB(反演,经典好题) Problem S ...
bzoj2154 Crash的数字表格
2154: Crash的数字表格 Time Limit: 20 Sec Memory Limit: 259 MB Submit: 4549 Solved: 1643 [Submit][Status ...
[BZOJ 1879][SDOI 2009]Bill的挑战题解(状压DP)
[BZOJ 1879][SDOI 2009]Bill的挑战 Description Solution 1.考虑状压的方式. 方案1:如果我们把每一个字符串压起来,用一个布尔数组表示与每一个字母的匹配关 ...
bzoj2154: Crash的数字表格
好神的莫比乌斯函数然后O(sqrt(n)*sqrt(n))好神的优化啊. #include<cstdio> #include<cstring> #include<ccty ...
P1829 [国家集训队]Crash的数字表格（推了好久的mobius反演）
P1829 [国家集训队]Crash的数字表格 / JZPTAB 推导过程 ∑i=1n∑j=1mlcm(i,j)\sum_{i = 1} ^{n} \sum_{j = 1} ^{m} lcm(i, j ...
P1829 [国家集训队]Crash的数字表格 / JZPTAB
P1829 [国家集训队]Crash的数字表格 / JZPTAB 题意: 求∑i=1n∑j=1mlcm(i,j)\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j)∑i=1n∑j ...
获得N^2个往返接力数字表格的算法
在描述算法之前,先看看下面的5*5的表格: 1 3 4 10 11 2 5 9 12 19 6 8 13 18 20 7 14 17 21 24 15 ...

【BZOJ 4816】【SDOI 2017】数字表格

【BZOJ 4816】【SDOI 2017】数字表格相关推荐

最新文章

热门文章