一些记号:
- d ( x ) d(x) d(x) 表示 x x x 的因数个数。
- 如无特殊说明,以下记为 p p p 的变量的取值集合为质数集合。
- 为了方便,有时用 a / b a/b a/b 表示 ⌊ a b ⌋ \lfloor\dfrac{a}{b}\rfloor ⌊ba⌋。
- 记模数为 P P P。
有个加号不太好处理,我们分开两部分来求: ∏ i = 1 n d ( i ) i \prod\limits_{i=1}^nd(i)^i i=1∏nd(i)i 和 ∏ i = 1 n d ( i ) μ ( i ) \prod\limits_{i=1}^nd(i)^{\mu(i)} i=1∏nd(i)μ(i)。
首先看第一部分:
∏
i
=
1
n
d
(
i
)
i
=
∏
p
k
≤
n
(
k
+
1
)
c
n
t
\prod_{i=1}^nd(i)^i=\prod_{p^k\leq n}(k+1)^{cnt}
i=1∏nd(i)i=pk≤n∏(k+1)cnt
其中
c
n
t
=
p
k
S
u
m
(
⌊
n
p
k
⌋
)
−
p
k
+
1
S
u
m
(
⌊
n
p
k
+
1
⌋
)
cnt=p^kSum(\lfloor\dfrac{n}{p^k}\rfloor)-p^{k+1}Sum(\lfloor\dfrac{n}{p^{k+1}}\rfloor)
cnt=pkSum(⌊pkn⌋)−pk+1Sum(⌊pk+1n⌋),
S
u
m
(
n
)
=
∑
i
=
1
n
i
=
n
(
n
+
1
)
2
Sum(n)=\sum\limits_{i=1}^ni=\dfrac{n(n+1)}{2}
Sum(n)=i=1∑ni=2n(n+1)。
我们按 p p p 分类讨论:
-
对于 p ≤ n p\leq \sqrt n p≤n 的,我们直接暴力计算即可,时间复杂度 O ( n log n log n log P ) = O ( n log P ) O(\dfrac{\sqrt n}{\log \sqrt n}\log n\log P)=O(\sqrt n\log P) O(lognnlognlogP)=O(nlogP)。
-
对于 p > n p>\sqrt n p>n 的,有 k = 1 k=1 k=1,于是变形为
∏ n < p ≤ n 2 p × S u m ( n / p ) = pow ( 2 , ∑ n < p ≤ n p × S u m ( n / p ) ) \prod_{\sqrt n<p\leq n}2^{p\times Sum(n/p)}=\operatorname{pow}\left(2,\sum_{\sqrt n< p\leq n}p\times Sum(n/p)\right) n<p≤n∏2p×Sum(n/p)=pow⎝⎛2,n<p≤n∑p×Sum(n/p)⎠⎞
于是只需要求出 ( ∑ n < p ≤ n p × S u m ( n / p ) ) m o d ( P − 1 ) \left(\sum\limits_{\sqrt n< p\leq n}p\times Sum(n/p)\right)\bmod (P-1) (n<p≤n∑p×Sum(n/p))mod(P−1),注意是模 φ ( P ) = P − 1 \varphi(P)=P-1 φ(P)=P−1。直接整除分块,需要解决的问题是区间质数和,用 min_25 筛解决即可。
然后看第二部分:
∏
i
=
1
n
d
(
i
)
μ
(
i
)
\prod_{i=1}^n d(i)^{\mu(i)}
i=1∏nd(i)μ(i)
注意到只有当
i
i
i 没有平方因子时才有贡献,设
i
=
p
1
p
2
⋯
p
m
i=p_1p_2\cdots p_m
i=p1p2⋯pm,则
d
(
i
)
=
∏
j
=
1
m
(
1
+
1
)
=
2
m
d(i)=\prod\limits_{j=1}^m(1+1)=2^m
d(i)=j=1∏m(1+1)=2m,于是设
d
′
(
i
)
d'(i)
d′(i) 表示
i
i
i 的不同的质因数个数,有:
∏
i
=
1
n
d
(
i
)
μ
(
i
)
=
∏
i
=
1
n
(
2
d
′
(
i
)
)
μ
(
i
)
=
pow
(
2
,
∑
i
=
1
n
d
′
(
i
)
μ
(
i
)
)
\prod_{i=1}^nd(i)^{\mu(i)}=\prod_{i=1}^n\left(2^{d'(i)}\right)^{\mu(i)}=\operatorname{pow}\left(2,\sum_{i=1}^{n}d'(i)\mu(i)\right)
i=1∏nd(i)μ(i)=i=1∏n(2d′(i))μ(i)=pow(2,i=1∑nd′(i)μ(i))
于是只需求出
(
∑
i
=
1
n
d
′
(
i
)
μ
(
i
)
)
m
o
d
(
P
−
1
)
\left(\sum\limits_{i=1}^{n}d'(i)\mu(i)\right)\bmod (P-1)
(i=1∑nd′(i)μ(i))mod(P−1)。
设
f
(
n
,
i
)
=
∑
x
=
1
n
[
(
x
∈
P
)
∨
(
lpf
(
x
)
>
p
i
)
]
d
′
(
x
)
μ
(
x
)
f(n,i)=\sum\limits_{x=1}^n[(x\in \mathbb{P})\lor(\operatorname{lpf}(x)>p_i)]d'(x)\mu(x)
f(n,i)=x=1∑n[(x∈P)∨(lpf(x)>pi)]d′(x)μ(x),
g
(
n
,
i
)
=
∑
x
=
1
n
[
(
x
∈
P
)
∨
(
lpf
(
x
)
>
p
i
)
]
μ
(
x
)
g(n,i)=\sum\limits_{x=1}^n[(x\in \mathbb{P})\lor(\operatorname{lpf}(x)>p_i)]\mu(x)
g(n,i)=x=1∑n[(x∈P)∨(lpf(x)>pi)]μ(x),那么:
f
(
n
,
i
)
=
f
(
n
,
i
+
1
)
−
(
(
f
(
n
/
p
i
+
1
,
i
+
1
)
−
f
(
p
i
+
1
,
i
+
1
)
)
+
(
g
(
n
/
p
i
+
1
,
i
+
1
)
−
g
(
p
i
+
1
,
i
+
1
)
)
)
g
(
n
,
i
)
=
g
(
n
,
i
+
1
)
−
(
g
(
n
/
p
i
+
1
,
i
+
1
)
−
g
(
p
i
+
1
,
i
+
1
)
)
\begin{aligned} f(n,i)&=f(n,i+1)-\bigg(\big(f(n/p_{i+1},i+1)-f(p_{i+1},i+1)\big)+\big(g(n/p_{i+1},i+1)-g(p_{i+1},i+1)\big)\bigg)\\ g(n,i)&=g(n,i+1)-\bigg(g(n/p_{i+1},i+1)-g(p_{i+1},i+1)\bigg) \end{aligned}
f(n,i)g(n,i)=f(n,i+1)−((f(n/pi+1,i+1)−f(pi+1,i+1))+(g(n/pi+1,i+1)−g(pi+1,i+1)))=g(n,i+1)−(g(n/pi+1,i+1)−g(pi+1,i+1))
初始化
f
(
n
,
∣
P
∣
)
=
∑
x
=
1
n
[
x
∈
P
]
d
′
(
x
)
μ
(
x
)
=
−
∣
P
∣
f(n,|\mathbb{P}|)=\sum\limits_{x=1}^n[x\in \mathbb{P}]d'(x)\mu(x)=-|\mathbb{P}|
f(n,∣P∣)=x=1∑n[x∈P]d′(x)μ(x)=−∣P∣ 时(其中
P
\mathbb{P}
P 是小于等于
n
n
n 的质数集合),需要求区间质数个数,在第一部分 min_25 筛的时候顺便筛出来即可。
取模时要用 __int128
。
代码如下:
#include<bits/stdc++.h>
#define N 320000
#define ll long long
using namespace std;
namespace modular
{
const ll P=1e12+39;
inline ll add(ll x,ll y,ll mod){return x+y>=mod?x+y-mod:x+y;}
inline ll dec(ll x,ll y,ll mod){return x-y<0?x-y+mod:x-y;}
inline ll mul(ll x,ll y,ll mod){return (__int128)x*y%mod;}
}using namespace modular;
inline ll poww(ll a,ll b,ll mod)
{
ll ans=1;
while(b)
{
if(b&1) ans=mul(ans,a,mod);
a=mul(a,a,mod);
b>>=1;
}
return ans;
}
ll n;
ll b[N<<1];
ll g0[N<<1],gp1[N],g1[N<<1];
ll f[N<<1],g[N<<1];
int sn,tot,id1[N],id2[N];
int cnt,p[N];
bool notprime[N];
void init()
{
for(int i=2;i<=sn;i++)
{
if(!notprime[i])
{
p[++cnt]=i;
gp1[cnt]=add(gp1[cnt-1],i,P-1);
}
for(int j=1;j<=cnt&&i*p[j]<=sn;j++)
{
notprime[i*p[j]]=1;
if(!(i%p[j])) break;
}
}
for(ll l=1,r;l<=n;l=r+1)
{
r=n/(n/l);
ll x=n/l;
b[++tot]=x;
if(x<=sn) id1[x]=tot;
else id2[n/x]=tot;
}
}
ll Sum(ll n,ll mod)
{
return (n&1)?mul((n+1)/2,n,mod):mul(n+1,n/2,mod);
}
int get(ll x)
{
return x<=sn?id1[x]:id2[n/x];
}
ll work1()
{
ll ans=1;
for(int i=1;i<=cnt;i++)
{
ll x=p[i];
for(int k=1;x<=n;k++,x*=p[i])
{
ll y=x*p[i];
ll tmp=dec(mul(x,Sum(n/x,P-1),P-1),mul(y,Sum(n/y,P-1),P-1),P-1);
ans=mul(ans,poww(k+1,tmp,P),P);
}
}
for(int i=1;i<=tot;i++)
{
g0[i]=b[i]-1;
g1[i]=dec(Sum(b[i],P-1),1,P-1);
}
for(int i=1;i<=cnt;i++)
{
ll pi2=1ll*p[i]*p[i];
for(int j=1;j<=tot&&pi2<=b[j];j++)
{
int v=get(b[j]/p[i]);
g0[j]=g0[j]-(g0[v]-(i-1));
g1[j]=dec(g1[j],mul(p[i],dec(g1[v],gp1[i-1],P-1),P-1),P-1);
}
}
ll sum=0;
for(ll l=sn+1,r;l<=n;l=r+1)
{
r=n/(n/l);
sum=add(sum,mul(dec(g1[get(r)],g1[get(l-1)],P-1),Sum(n/l,P-1),P-1),P-1);
}
return mul(ans,poww(2,sum,P),P);
}
ll work2()
{
for(int j=1;j<=tot;j++)
f[j]=g[j]=dec(0,g0[j],P-1);
for(int i=cnt-1;i>=0;i--)
{
ll pi2=1ll*p[i+1]*p[i+1];
for(int j=1;j<=tot&&pi2<=b[j];j++)
{
int v=get(b[j]/p[i+1]);
f[j]=dec(f[j],add(add(f[v],i+1,P-1),add(g[v],i+1,P-1),P-1),P-1);
g[j]=dec(g[j],add(g[v],i+1,P-1),P-1);
}
}
return poww(2,f[get(n)],P);
}
int main()
{
scanf("%lld",&n);
sn=sqrt(n);
init();
ll ans1=work1();
ll ans2=work2();
printf("%lld\n",mul(ans1,ans2,P));
return 0;
}