题目链接:
题目大意:
快速求:
[prod_{i=1}^{n}prod_{j=1}^{m}f_{gcd(i,j)}
]
其中 (f_i) 表示斐波那契数列第 (i) 个数。
正文:
在写本题之前,建议先拿 【Luogu P2257】 YY的GCD 练练手。
将原式化为我们能够接受的时限:
[egin{aligned}prod_{i=1}^{n}prod_{j=1}^{m}f_{gcd(i,j)}&=prod_{d=1}^{min(n,m)}prod_{i=1}^{n}prod_{j=1}^{m}f_dleft[d==gcd(i,j)
ight]\&=prod_{d=1}^{min(n,m)}{f_d}^{left(sum_{i=1}^{lfloorfrac{n}{d}
floor}sum_{j=1}^{lfloorfrac{m}{d}
floor}left[gcd(i,j)==1
ight]
ight)}\&=prod_{d=1}^{min(n,m)}{f_d}^{left(sum_{k|d}mu(k)leftlfloorfrac{n}{kd}
ight
floorleftlfloorfrac{m}{kd}
ight
floor
ight)}\&=prod_{T=1}^{min(n,m)}left(prod_{d|T}{f_d}^{mu(frac{T}{d})}
ight)^{leftlfloorfrac{n}{T}
ight
floorleftlfloorfrac{m}{T}
ight
floor}end{aligned}
]
设 (F_i) 括号内的数,即 (F_i=prod_{d|i}{f_d}^{mu(frac{i}{d})}),这个函数用 (O(nlog n)) 预处理出来,在外面再套一个整除分块就是答案。
代码:
inline ll qpow(ll a, ll b)
{
if(b < 0) b = p - 2;
ll ans = 1;
for (; b; b >>= 1)
{
if (b & 1) ans = (ans * a % p + p) % p;
a = (a * a % p + p) % p;
}
return ans;
}
inline void prework()
{
miu[1] = 1;
for (ll i = 2; i <= N - 10; i++)
{
if(!vis[i]) {pri[++cnt] = i, miu[i] = -1;}
for (int j = 1; j <= cnt && pri[j] * i <= N - 10; j++)
{
vis[pri[j] * i] = 1;
if (i % pri[j] == 0)
{
miu[i * pri[j]] = 0;
break;
}
else
miu[i * pri[j]] = -miu[i];
}
}
sum[1] = fib[1] = fib_inv[1] = 1LL; //求斐波那契数列及其逆元(这里用fib代替)
for (int i = 2; i <= N - 10; i++)
fib[i] = (fib[i - 1] + fib[i - 2]) % p,
fib_inv[i] = qpow(fib[i], p - 2),
sum[i] = 1LL;
for (int i = 1; i <= N - 10; i++) //求F(n)及其前缀和(积)
for (int j = i; j <= N - 10; j += i)
{
int t = j / i;
if(miu[t] == 1) sum[j] = (sum[j] * fib[i] % p + p) % p;
else if(miu[t] == -1) sum[j] = (sum[j] * fib_inv[i] % p + p) % p;
}
inv[0] = sum[0] = 1;
for (int i = 1; i <= N - 10; i++)
sum[i] = (sum[i] * sum[i - 1] % p + p) % p,
inv[i] = qpow(sum[i], p - 2);
}
int main()
{
prework();
for (read(t); t--; )
{
ans = 1LL;
read(n);read(m);
if(n > m)
{
ll c = n; n = m; m = c;
}
for (register int l = 1, r; l <= n; l = r + 1)
{
r = min (n / (n / l), m / (m / l));
ans = ans * qpow((sum[r] * inv[l - 1] % p + p) % p, (1ll * (n / r) * (m / r))) % p;
}
print(ans);
}
return 0;
}