「51nod 1965」奇怪的式子

为什么模数是偶数我都敢先取模后除以2了。不知道是不是傻了

51nod 1965

题意

\[\prod_{i=1}^n \sigma_0(i)^{\mu(i)+i} \mod(10^{12}+39)\]

其中 \(\sigma_0(i)\) 表示 \(i\) 的正约数个数,\(10^{12}+39\) 是质数

\(n\le 10^{11}\)


做法

\[ \prod_{i=1}^n \sigma_0(i)^{\mu(i)+i} = \prod_{i=1}^n \sigma_0(i)^i \prod_{i=1}^n \sigma_0(i)^{\mu(i)} \]

两部分可以分开来做

考虑每个质数的贡献

\[ \prod_{i=1}^n \sigma_0(i)^i = \prod_{p\text{是质数}} \prod_{p^i\le n} (i+1)^{f(\lfloor \frac{n}{p^i} \rfloor) * p^i - f(\lfloor \frac{n}{p^{i+1}} \rfloor) * p^{i+1}} \]

其中 \(f(x)=\sum\limits_{i=1}^x i=\frac{x(x+1)}{2}\)

  • 对于 \(p\le \sqrt{n}\) 枚举 \(p\) 计算贡献

  • 对于 \(p>\sqrt{n}\),指数只能为 \(1\)

    于是贡献为

    \[ \prod_{p>\sqrt{n},p\text{是质数}} 2^{f(\lfloor \frac{n}{p} \rfloor) * p} \]

    用Min_25筛出对于每一个 \(\lfloor \frac{n}{x} \rfloor\)\(\sum\limits_{i=1}^{\lfloor \frac{n}{x} \rfloor} [i\text{是质数}] * i\)

    分块计算即可

要求 \(\prod\limits_{i=1}^n \sigma_0(i)^{\mu(i)}\)

可以发现对于有平方因子的 \(i\)\(\mu(i)=0\),所以没有贡献

\(g(i)\) 表示 \(i\) 的质因子数

\[ \begin{align} & \prod_{i=1}^n \sigma_0(i)^{\mu(i)} \\ = & \prod_{i=1}^n 2^{\mu(i)g(i)} \\ = & 2^{\left(\sum_{i=1}^n \mu(i)g(i)\right)} \end{align} \]

考虑怎么求 \(\sum\limits_{i=1}^n \mu(i)g(i)\)

\[ \begin{align} F(a,b)&=\sum_{i=2}^a [i\text{是质数 或 } pmin_p\ge prime_b] * \mu(i) \\ S(a,b)&=\sum_{i=2}^a [i\text{是质数 或 } pmin_p\ge prime_b] * \mu(i)g(i) \end{align} \]

其中 \(pmin_i\) 表示 \(i\) 的最小质因子,\(prime_i\) 表示第 \(i\) 个质数

于是可以用Min_25筛转移

\[ \begin{align} F(a,b)&=F(a,b+1)-\left(F(\lfloor \frac{a}{prime_b} \rfloor,b+1)+b\right) \\ S(a,b)&=S(a,b+1)-\left(S(\lfloor \frac{a}{prime_b} \rfloor,b+1)+F(\lfloor \frac{a}{prime_b},b+1)+2b\right) \end{align} \]

注意这里的模数有点大,需要用快速乘

总复杂度 \(\mathcal O(\frac{n^{\frac{3}{4}}}{\log n})\)


代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
#include<cstdio>
#include<algorithm>
#include<ctype.h>
#include<string.h>
#include<math.h>

using namespace std;
#define ll long long

const int N = 316230;
const ll P = 1000000000039;
int T, sn, cnt, id, prime[N];
ll n, f[45], a[N<<1], g[N<<1], h[N<<1], S[N<<1];
inline ll Mul(ll x, ll y, ll mod){
ll ans=(x*y-(ll)((long double)x*y/mod)*mod)%mod;
return ans<0?ans+mod:ans;
}
inline ll calc(ll x){ return (x&1)?Mul(x, (x+1)/2, P-1):Mul(x/2, x+1, P-1);}
inline int Id(ll x){ return x<=sn?x:id-n/x+1;}
inline ll Pow(ll x, ll y){
ll ans=1;
for(; y; y>>=1, x=Mul(x, x, P)) if(y&1) ans=Mul(ans, x, P);
return ans;
}
int main() {
scanf("%d", &T);
while(T--){
scanf("%lld", &n), sn=sqrt(n);
cnt=id=0;
for(ll i=1; i<=n; i=a[id]+1)
a[++id]=n/(n/i), g[id]=calc(a[id])-1, h[id]=a[id]-1;
for(int i=2; i<=sn; ++i) if(h[i]!=h[i-1]){
prime[++cnt]=i;
for(int j=id, tmp; (ll)i*i<=a[j]; --j)
(g[j]-=Mul(i, g[tmp=Id(a[j]/i)]-g[i-1], P-1)-P+1)%=P-1, h[j]-=h[tmp]-(cnt-1);
}
memset(f, 0, sizeof f);
for(int i=1; i<=cnt; ++i)
for(ll x=prime[i], j=2; x<=n; x*=prime[i], ++j)
f[j]+=Mul(calc(n/x), x, P-1)-Mul(calc(n/(x*prime[i])), x*prime[i], P-1);
ll ans=1;
for(int i=2; i<=40; ++i) ans=Mul(ans, Pow(i, (f[i]%(P-1)+(P-1))%(P-1)), P);
for(int i=1; i<=sn; ++i) ans=Mul(ans, Pow(2, Mul(calc(a[i]), g[id-i+1]-g[max(sn, id-i)]+P-1, P-1)), P);

for(int i=1; i<=id; ++i) S[i]=h[i]=-h[i];
for(int i=cnt; i; --i) for(int j=id; (ll)prime[i]*prime[i]<=a[j]; --j){
int t=Id(a[j]/prime[i]);
h[j]-=h[t]+i;
S[j]-=S[t]+h[t]+2*i;
}
printf("%lld\n", Mul(ans, Pow(2, (S[id]%(P-1)+P-1)%(P-1)), P));
}
return 0;
}