\(\\\)

\(Description\)


定义二元函数\(F(x,y)\)表示,用 \(x\times y\) 的矩形不可旋转的铺成一个任意边长的正方形,所需要的最少的矩形个数。

现在\(T\)组询问,每次给出一个\(N\),求\(\prod_{i=1}^N\prod_{j=1}^N F(i,j)\) 模\(19260817\)的值。

  • \(N\in[1,10^6],T\in[1,10^3]\)

\(\\\)

下面的表述中用\([x,y]\)表示\(Lcm(x,y)\),用\((x,y)\)表示\(Gcd(x,y)\)。

\(\\\)

\(Solution1:\text O(N+T\sqrt NlogN)\)


出题人的做法。对于\(N\)更大一些,\(T\)较小的情况表现比较优秀。

考虑一个边长为\(x\times y\)的矩形铺成的最小正方形的边长,因为不可旋转,所以答案为\([x,y]\)。

那么用这种矩形铺出这个正方形需要的块数就是\(\frac{[x,y]^2}{xy}\),总面积除以单块面积嘛。

然后所求可以形式化的写出:

\[\prod_{i=1}^N\prod_{j=1}^N\frac{[i,j]^2}{ij}=\prod_{i=1}^N\prod_{j=1}^N\frac{(\frac{ij}{(i,j)})^2}{ij}=\prod_{i=1}^N\prod_{j=1}^N\frac{ij}{(i,j)^2}=\frac{\prod_{i=1}^N\prod_{j=1}^Nij}{\prod_{i=1}^N\prod_{j=1}^N(i,j)^2}
\]

上面的转化是基于以下定理:

\[i\times j=(i,j)\times[i,j]
\]

然后考虑对于每一次询问的\(N\),如何快速求出答案。

\(\\\)

首先将分母分子分开考虑,对于分子,有:

\[\prod_{i=1}^N\prod_{j=1}^Nij=\prod_{i=1}^Ni^N\times N\ !=(N\ !)^{2N}
\]

具体证明不太好说,大致是每一次循环到的\(i\)都要被乘上\(N\)次,手玩一下差不多就懂了。

这部分显然可以\(\text O(N)\)求\(N!\),然后再每一个数都算一个快速幂即可,预处理复杂度\(\text O(Nlog(2N))\)。

\(\\\)

然后考虑分子部分如何求,首先先将平方提出来:

\[\prod_{i=1}^N\prod_{j=1}^N(i,j)^2=\bigg(\prod_{i=1}^N\prod_{j=1}^N(i,j)\bigg)^2
\]

因为是连乘显然合理。

然后括号里的东西显然可以根据\(Gcd\)的值分情况讨论:

\[\prod_{i=1}^N\prod_{j=1}^N(i,j)=\prod_{d=1}^N d^{\sum_{i=1}^N\sum_{j=1}^N[(i,j)=d]}=\prod_{d=1}^N d^{\sum_{i=1}^{\lfloor\frac Nd\rfloor}\sum_{j=1}^{\lfloor\frac Nd\rfloor}[(i,j)=1]}
\]

然后根据仪仗队的解法,我们知道有:

\[\sum_{i=1}^{\lfloor\frac Nd\rfloor}\sum_{j=1}^{\lfloor\frac Nd\rfloor}[(i,j)=1]=\bigg(\sum_{i=1}^{\lfloor \frac Nd\rfloor}\varphi(i)\bigg)\times 2-1
\]

然后可以线性筛求\(\varphi\),预处理\(\varphi\)的前缀和,查询这个指数就是\(\text O(1)\)的了,预处理复杂度\(\text O(N)\)。

我们设\(g[x]\)表示:

\[g[x]=\sum_{i=1}^{x}\sum_{j=1}^{x}[(i,j)=1]=\bigg(\sum_{i=1}^{x}\varphi(i)\bigg)\times 2-1
\]

那么现在需要解决的式子是

\[\prod_{i=1}^N\prod_{j=1}^N(i,j)=\prod_{d=1}^N d^{\ g[\lfloor\frac Nd \rfloor]}
\]

然后注意到这个东西可以除法分块,考虑对于一个\(\lfloor\frac Nd\rfloor\)相同的区间\([L,R]\),答案为

\[\prod_{i=L}^R i^{\ g[\lfloor\frac NL\rfloor]}=\bigg(\prod_{i=1}^N i\bigg)^{g[\lfloor\frac NL\rfloor]}=\bigg(\frac{R\ !}{(L-1)\ !}\bigg)^{g[\lfloor\frac NL\rfloor]}
\]

然后连乘积部分可以用阶乘相除的方法求解,因为模意义下我们还需要处理每一个阶乘的逆元。

然后快速幂求一下,在套回去平方一下,再求个逆元乘上\((N\ !)^{2N}\)就是答案了。

单次处理除法分块+快速幂复杂度\(\text O(\sqrt NlogN)\),总复杂度\(\text O(N+T\sqrt NlogN)\)

\(\\\)

\(Code1:\text O(N+T\sqrt NlogN)\)


#include<cmath>
#include<cctype>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
#define R register
#define gc getchar
#define N 1000010
#define mod 19260817ll
using namespace std;
typedef long long ll; ll g[N]={1,1},phi[N]={0,1},prm[N];
ll fac[N]={0,1},fac2[N]={0,1},inv[N]={1}; inline int rd(){
int x=0; bool f=0; char c=gc();
while(!isdigit(c)){if(c=='-')f=1;c=gc();}
while(isdigit(c)){x=(x<<1)+(x<<3)+(c^48);c=gc();}
return f?-x:x;
} inline ll qpow(ll x,ll t){
ll res=1ll;
while(t){
if(t&1) (res*=x)%=mod;
(x*=x)%=mod; t>>=1;
}
return res;
} inline void init(){
for(R int i=2;i<N;++i){
g[i]=1;
if(!phi[i]){phi[i]=i-1;prm[++prm[0]]=i;}
for(R int j=1,k;j<=prm[0]&&(k=prm[j]*i)<N;++j)
if(i%prm[j]) phi[k]=phi[i]*phi[prm[j]];
else{phi[k]=phi[i]*prm[j];break;}
}
for(R ll i=2;i<N;++i) fac[i]=fac[i-1]*i%mod,phi[i]+=phi[i-1];
for(R int i=2;i<N;++i) fac2[i]=qpow(fac[i],(i<<1));
for(R int i=1;i<N;++i) inv[i]=qpow(fac[i],mod-2);
} inline void work(){
ll n=rd(),ans=1ll;
for(R ll l=1,t,r;l<=n;l=r+1){
t=n/l; r=n/t;
(ans*=qpow(fac[r]*inv[l-1]%mod,phi[t]*2-1))%=mod;
}
ans=(qpow(ans,(mod-2)<<1)*fac2[n])%mod;
printf("%lld\n",ans);
} int main(){
init();
int t=rd();
while(t--) work();
return 0;
}

\(\\\)

\(Solution2:\text O(N\sqrt NlogN+T)\)


我的做法。卡着时限也就过了,对\(T\)大的时候表现会很优秀。

继续考虑这个式子:

\[\prod_{i=1}^N\prod_{j=1}^N\frac{[i,j]^2}{ij}=\prod_{i=1}^N\prod_{j=1}^N\frac{(\frac{ij}{(i,j)})^2}{ij}=\prod_{i=1}^N\prod_{j=1}^N\frac{ij}{(i,j)^2}=\frac{\prod_{i=1}^N\prod_{j=1}^Nij}{\prod_{i=1}^N\prod_{j=1}^N(i,j)^2}
\]

分子部分不变,还是那么预处理,分母我有一个不同的做法。

首先还是把平方套在外面:

\[\prod_{i=1}^N\prod_{j=1}^N(i,j)^2=\bigg(\prod_{i=1}^N\prod_{j=1}^N(i,j)\bigg)^2
\]

然后平方里面的部分解法就不太一样了,设一个\(g[n]\)表示:

\[g[n]=\prod_{i=1}^n(i,n)=\prod_{d=1}^nd^{\ \varphi(\lfloor\frac{N}{d}\rfloor)}
\]

第一个等号是定义,第二个等号理解为,设\((i,n)=d\),那么这样的数对有\(\varphi(\lfloor\frac{N}{d}\rfloor)\)个,因为除掉\(d\)两数应该互质。

然后类似仪仗队的做法,把求和扩展到求积,那么考虑将下三角矩阵的积求出来,答案是

\[\prod_{i=1}^Ng[i]
\]

然后整个矩阵的积是下三角的平方除以对角线,对角线的积是\(\prod(i,i)=N\ !\),所以整个矩阵的积等于

\[\frac{(\prod_{i=1}^Ng[i])^2}{N\ !}
\]

回到所求,有

\[\frac{\prod_{i=1}^N\prod_{j=1}^Nij}{\prod_{i=1}^N\prod_{j=1}^N(i,j)^2}=\frac{(N\ !)^{2N}}{\bigg(\frac{(\prod_{i=1}^Ng[i])^2}{N\ !}\bigg)^2}=\frac{(N\ !)^{2N+2}}{(\prod_{i=1}^Ng[i])^4}
\]

复杂度预处理 \(g\) 数组\(\text O(N\sqrt NlogN)\),分别代表质因数分解和快速幂的复杂度。

其他线性筛、求阶乘、求逆元等的复杂度都不会太高,之后显然每一个数的答案都可以\(\text O(1)\)得到了。

\(\\\)

\(Code2:\text O(N\sqrt NlogN+T)\)


#include<cmath>
#include<cctype>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
#define R register
#define gc getchar
#define N 1000010
#define mod 19260817ll
using namespace std;
typedef long long ll; ll g[N]={1},fac[N]={0,1}; int phi[N]={0,1},prm[N]; inline int rd(){
int x=0; bool f=0; char c=gc();
while(!isdigit(c)){if(c=='-')f=1;c=gc();}
while(isdigit(c)){x=(x<<1)+(x<<3)+(c^48);c=gc();}
return f?-x:x;
} inline ll qpow(ll x,ll t){
ll res=1ll;
while(t){
if(t&1) (res*=x)%=mod;
(x*=x)%=mod; t>>=1;
}
return res;
} inline void init(){
for(R int i=2;i<N;++i){
if(!phi[i]){phi[i]=i-1;prm[++prm[0]]=i;}
for(R int j=1,k;j<=prm[0]&&(k=prm[j]*i)<N;++j)
if(i%prm[j]) phi[k]=phi[i]*phi[prm[j]];
else{phi[k]=phi[i]*prm[j];break;}
}
for(R ll i=2;i<N;++i) fac[i]=fac[i-1]*i%mod;
for(R int i=2;i<N;++i) fac[i]=qpow(fac[i],(i<<1)+2);
for(R int i=1,r;i<N;++i){
r=sqrt(i); g[i]=1ll;
for(R int j=1;j<=r;++j)
if(i%j==0){
(g[i]*=qpow(j,phi[i/j]))%=mod;
if(i/j!=j) (g[i]*=qpow(i/j,phi[j]))%=mod;
}
(g[i]*=g[i-1])%=mod;
}
for(R int i=1;i<N;++i) g[i]=qpow(g[i],(mod-2)*4);
for(R int i=1;i<N;++i) (g[i]*=fac[i])%=mod;
} int main(){
init();
int t=rd();
while(t--) printf("%lld\n",g[rd()]);
return 0;
}

\(\\\)

\(Solution3:\text O(N\times lnN\times logN+T)\)


优化一下第二个做法。

依次考虑每一个数对每一个 \(g\) 的贡献,所以一个数 \(i\) 能做贡献的数有\(\lfloor\frac Ni\rfloor\)个,直接暴力乘更新 \(g\) 就好。

此时复杂度神奇的降到了调和级数,也就是\(ln N\)。

\(\\\)

\(Code3:\text O(N\times lnN\times logN+T)\)


#include<cmath>
#include<cctype>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
#define R register
#define gc getchar
#define N 1000010
#define mod 19260817ll
using namespace std;
typedef long long ll; ll g[N]={1,1},fac[N]={0,1}; int phi[N]={0,1},prm[N]; inline int rd(){
int x=0; bool f=0; char c=gc();
while(!isdigit(c)){if(c=='-')f=1;c=gc();}
while(isdigit(c)){x=(x<<1)+(x<<3)+(c^48);c=gc();}
return f?-x:x;
} inline ll qpow(ll x,ll t){
ll res=1ll;
while(t){
if(t&1) (res*=x)%=mod;
(x*=x)%=mod; t>>=1;
}
return res;
} inline void init(){
for(R int i=2;i<N;++i){
g[i]=1;
if(!phi[i]){phi[i]=i-1;prm[++prm[0]]=i;}
for(R int j=1,k;j<=prm[0]&&(k=prm[j]*i)<N;++j)
if(i%prm[j]) phi[k]=phi[i]*phi[prm[j]];
else{phi[k]=phi[i]*prm[j];break;}
}
for(R ll i=2;i<N;++i) fac[i]=fac[i-1]*i%mod;
for(R int i=2;i<N;++i) fac[i]=qpow(fac[i],(i<<1)+2);
for(R int i=1;i<N;++i)
for(R ll j=1;j*(ll)i<(ll)N;++j) (g[j*i]*=qpow(i,phi[j]))%=mod;
for(R int i=1;i<N;++i) (g[i]*=g[i-1])%=mod;
for(R int i=1;i<N;++i) g[i]=qpow(g[i],(mod-2)*4);
for(R int i=1;i<N;++i) (g[i]*=fac[i])%=mod;
} int main(){
init();
int t=rd();
while(t--) printf("%lld\n",g[rd()]);
return 0;
}

最新文章

  1. Redsi和Memcached区别总结
  2. Objective-C 桥接模式 -- 简单实用和说明
  3. ZOJ Problem Set - 1049 I Think I Need a Houseboat
  4. [No000077]打造自己的Eclipse
  5. ZeroMQ接口函数之 :zmq_socket – 创建ZMQ套接字
  6. hdu 5676 ztr loves lucky numbers
  7. Disconf
  8. 使用dynamic类型改进反射
  9. MyEclipse安装插件的三种方法和使用心得
  10. jquery滚动条
  11. NWERC 2012 Problem A Admiral
  12. 解决SQL Server的TEXT、IMAGE类型字段的长度限制
  13. 由Lucnene 对于预治疗的文字,全角半角转换器(可执行)
  14. linux c编程获得当前进程的进程名和执行路径
  15. jvisualvm
  16. 01-01_环境准备_pyenv
  17. JSON基础(JavaScript)
  18. ASP.NET平台下从浏览器地址栏输入之后发生的事
  19. docker (2) 通用/镜像命令
  20. Cannot locate BeanDefinitionParser for element [scoped-proxy]

热门文章

  1. 介绍css 的3D 变换(3D transform)
  2. CCNP路由实验之八 路由重公布
  3. mybatis中批量插入数据
  4. Android MaoZhuaWeiBo开发Service抓取个人信息-2
  5. SQLyog软件里无法插入中文(即由默认的latin1改成UTF8编码格式)
  6. 2015-2016 ACM-ICPC Pacific Northwest Regional Contest (Div. 2) S Surf
  7. RK3288的gpio设置【转】
  8. YTU 2625: B 构造函数和析构函数
  9. textView设置按下和焦点改变时让字体颜色发生变化
  10. django 数据库连接模块解析及简单长连接改造