H y p e r l i n k Hyperlink Hyperlink
https://www.lydsy.com/JudgeOnline/problem.php?id=4816
https://www.luogu.org/problem/P3704
D e s c r i p t i o n Description Description
求
∏ i = 1 n ∏ j = 1 m f g c d ( i , j ) \huge\prod_{i=1}^n\prod_{j=1}^mf_{gcd(i,j)} i=1∏nj=1∏mfgcd(i,j)
f f f表示斐波那契数列
数据范围: T ≤ 1000 , 1 ≤ n , m ≤ 1 0 6 T\leq1000,1\leq n,m\leq 10^6 T≤1000,1≤n,m≤106
S o l u t i o n Solution Solution
问题转换成
∏
d
=
1
n
∏
i
=
1
n
∏
j
=
1
m
[
g
c
d
(
i
,
j
)
=
=
d
]
f
d
\huge \prod_{d=1}^n\prod_{i=1}^n\prod_{j=1}^m[gcd(i,j)==d]f_d
d=1∏ni=1∏nj=1∏m[gcd(i,j)==d]fd
把
f
d
f_d
fd拉出来,每个
f
d
f_d
fd乘的次数就是
∑
i
=
1
n
/
d
∑
j
=
1
m
/
d
[
g
c
d
(
i
,
j
)
=
=
1
]
\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[gcd(i,j)==1]
∑i=1n/d∑j=1m/d[gcd(i,j)==1],所以
∏ d = 1 n f d ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ m d ⌋ [ g c d ( i , j ) = = 1 ] \huge \huge \prod_{d=1}^nf_d^{\sum_{i=1}^{\lfloor \frac n d\rfloor}\sum_{j=1}^{\lfloor\frac m d\rfloor}[gcd(i,j)==1]} d=1∏nfd∑i=1⌊dn⌋∑j=1⌊dm⌋[gcd(i,j)==1]
到这里数论分块套数论分块就可以做到 O ( n ) O(n) O(n)了,但是会T飞
再根据之前莫比乌斯反演题的经验啊。。。
∑
i
=
1
n
/
d
∑
j
=
1
m
/
d
[
g
c
d
(
i
,
j
)
=
=
1
]
−
>
∑
i
=
1
n
/
d
μ
(
i
)
⌊
n
i
d
⌋
⌊
m
i
d
⌋
\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[gcd(i,j)==1]->\sum_{i=1}^{n/d}\mu(i)\lfloor\frac n{id}\rfloor\lfloor\frac m{id}\rfloor
∑i=1n/d∑j=1m/d[gcd(i,j)==1]−>∑i=1n/dμ(i)⌊idn⌋⌊idm⌋
继续用YY的GCD的方法,设
T
=
i
d
T=id
T=id,放入整个式子中
∏
d
=
1
n
f
d
∑
i
=
1
⌊
n
d
⌋
μ
(
i
)
⌊
n
i
d
⌋
⌊
m
i
d
⌋
\huge \huge \prod_{d=1}^nf_d^{\sum_{i=1}^{\lfloor\frac nd \rfloor}\mu(i)\lfloor\frac n{id}\rfloor\lfloor\frac m{id}\rfloor}
d=1∏nfd∑i=1⌊dn⌋μ(i)⌊idn⌋⌊idm⌋
−
>
∏
T
=
1
n
(
∏
d
∣
T
f
d
μ
(
T
/
d
)
)
⌊
n
T
⌋
⌊
m
T
⌋
->\huge \huge \prod_{T=1}^n(\prod_{d|T}f_d^{\mu(T/d)})^{\lfloor\frac n{T}\rfloor\lfloor\frac m{T}\rfloor}
−>T=1∏n(d∣T∏fdμ(T/d))⌊Tn⌋⌊Tm⌋
然后暴力求中间那一坨,然后整除分块即可
时间复杂度:
O
(
A
+
B
+
C
)
O(A+B+C)
O(A+B+C)
A
:
A:
A:线性筛+快速幂,
B
:
B:
B:暴力计算中间那一坨(跑不满),
C
:
C:
C:回答询问
O
(
n
l
o
g
n
+
n
l
o
g
2
n
+
T
n
l
o
g
n
)
O(nlogn+nlog^2n+T\sqrt{n}logn)
O(nlogn+nlog2n+Tnlogn)
C o d e Code Code
#include<cstdio>
#include<cctype>
#include<algorithm>
#define N 1000000
#define LL long long
#define mod 1000000007
using namespace std;
LL mu[N+10],f[N+10],F[N+10],g[N+10];int prime[N>>3],n,m,t;
bool vis[N+10];
inline LL read()
{
LL f=0,d=1;char c;
while(c=getchar(),!isdigit(c)) if(c=='-') d=-1;f=(f<<3)+(f<<1)+c-48;
while(c=getchar(),isdigit(c)) f=(f<<3)+(f<<1)+c-48;
return d*f;
}
inline LL ksm(LL x,LL y)
{
LL ans=1;
for(;y;y>>=1,(x*=x)%=mod) if(y&1) (ans*=x)%=mod;
return ans;
}
inline LL solve(int a,int b)//C
{
if(a>b) a^=b,b=a^b,a^=b;
LL res=1,inv=0;
for(register int l=1,r=0;l<=a;l=r+1)
{
r=min(a/(a/l),b/(b/l));
inv=1ll*F[r]*ksm(F[l-1],mod-2)%mod;
(res*=ksm(inv,1ll*(a/l)*(b/l)%(mod-1)))%=mod;//扩展欧拉定理a^b=a^{b mod phi(p)}
}
return res;
}
signed main()
{
mu[1]=1;f[1]=1;g[1]=1;F[0]=F[1]=1;
for(register int i=2;i<=N;i++)//A
{
f[i]=(f[i-1]+f[i-2])%mod;g[i]=ksm(f[i],mod-2);F[i]=1;
if(vis[i]==0) {prime[++m]=i;mu[i]=-1;}
for(register int j=1;j<=m&&prime[j]*i<=N;j++)
{
vis[prime[j]*i]=true;
if(i%prime[j]==0) break;
mu[prime[j]*i]=-mu[i];
}
}
for(register int i=1;i<=N;i++)//B
{
if(mu[i]==0) continue;
for(register int j=i;j<=N;j+=i)
F[j]=1ll*F[j]*(mu[i]==1?f[j/i]:g[j/i])%mod;
}
for(register int i=2;i<=N;i++) F[i]=1ll*F[i]*F[i-1]%mod;
for(t=read();t--;)
{
n=read();m=read();
printf("%lld\n",solve(n,m));
}
}