这题还是挺值得一做的。
简要题意：
求 (1000) 次
[prod_{i=1}^n prod_{j=1}^n frac{lcm(i,j)^2}{i*j}
]
(i、j) 规模均为 (1e6)。
淦。
[prod_{i=1}^n prod_{j=1}^n frac{lcm(i,j)^2}{i*j}
]
[= prod_{i=1}^n prod_{j=1}^n frac{i*j}{gcd(i,j)^2}
]
由于是累乘， 所以拆开来算
[= frac{ prod_{i=1}^n prod_{j=1}^n i*j } { prod_{i=1}^n prod_{j=1}^n gcd(i,j)^2}
]
上面的直接预处理， 重点是下面的。
[prod_{i=1}^n prod_{j=1}^n gcd(i,j)^2
]
[= prod_{g=1}^n g^{2 * sum_{i=1}^{lfloor frac{n}g
floor} sum_{j=1}^{lfloor frac{n}g
floor} [gcd(i,j)==1]}
]
设(f(n) = sum_{i=1}^{n} sum_{j=1}^{n} [gcd(i,j)==1])
于是
[prod_{g=1}^n g^{2 * sum_{i=1}^{lfloor frac{n}g
floor} sum_{j=1}^{lfloor frac{n}g
floor} [gcd(i,j)==1]}
]
[= prod_{g=1}^n g^{2 * f(lfloor frac{n}g
floor)}
]
可以分块套快速幂算。
(f(n))怎么算？
[f(n) = sum_{i=1}^{n} sum_{j=1}^{n} [gcd(i,j)==1]
]
[= 2*sum_{i=1}^n varphi(i)-1
]
也可以预处理。
总复杂度就是 (O(n)) 预处理 加上 (O( T sqrt n * 快速幂的log))
Luogu 数据 AC 代码：
#include<bits/stdc++.h>
using namespace std;
#define li long long
const int mod = 19260817; // prime
const int maxn = 1e6 + 15;
li ksm(li a, li b, li p) {
li res = 1ll;
for(;b;b>>=1, a=(a*a)%p)
if(b&1) res = (res*a)%p;
return res%p;
}
li jc[maxn], phi[maxn];
int pr_cnt, prime[maxn], v[maxn];
void euler(int n) {
jc[0] = 1ll;
phi[1] = 1ll;
for(int i=2; i<=n; ++i) {
if(!v[i]) {
v[prime[++pr_cnt] = i] = i;
phi[i] = i-1ll;
}
for(int j=1; j<=pr_cnt; ++j) {
if(prime[j] > n/i || prime[j] > v[i]) break;
v[prime[j] * i] = prime[j];
phi[prime[j] * i] = phi[i] * (i%prime[j] ? prime[j]-1 : prime[j]);
}
}
for(int i=1; i<=n; ++i) {
jc[i] = 1ll * jc[i-1] * i % mod;
phi[i] += phi[i-1];
}
}
li f(int n) {
return 2ll*phi[n] - 1ll;
}
signed main()
{
euler(1000000);
int T,n; cin >> T; while(T--) {
scanf("%d", &n);
li ans = ksm(jc[n],2*n,mod);
li downans = 1ll;
for(int i=1,j;i<=n;i=j+1) {
j = min(n, n/(n/i));
li b = 2 * f(n/i);
li a = jc[j] * ksm(jc[i-1],mod-2,mod) % mod;
downans *= ksm(a,b,mod);
downans = (downans%mod+mod)%mod;
}
ans *= ksm(downans,mod-2,mod);
ans = (ans%mod+mod)%mod;
cout << ans << '
';
}
return 0;
}