【51NOD1965】奇怪的式子（min_25筛）

一些记号：

• \(d(x)\) 表示 \(x\) 的因数个数。
• 如无特殊说明，以下记为 \(p\) 的变量的取值集合为质数集合。
• 为了方便，有时用 \(a/b\) 表示 \(\lfloor\dfrac{a}{b}\rfloor\)。
• 记模数为 \(P\)。

有个加号不太好处理，我们分开两部分来求：\(\prod\limits_{i=1}^nd(i)^i\) 和 \(\prod\limits_{i=1}^nd(i)^{\mu(i)}\)。

首先看第一部分：

\[\prod_{i=1}^nd(i)^i=\prod_{p^k\leq n}(k+1)^{cnt} \]

其中 \(cnt=p^kSum(\lfloor\dfrac{n}{p^k}\rfloor)-p^{k+1}Sum(\lfloor\dfrac{n}{p^{k+1}}\rfloor)\)，\(Sum(n)=\sum\limits_{i=1}^ni=\dfrac{n(n+1)}{2}\)。

我们按 \(p\) 分类讨论：

• 对于 \(p\leq \sqrt n\) 的，我们直接暴力计算即可，时间复杂度 \(O(\dfrac{\sqrt n}{\log \sqrt n}\log n\log P)=O(\sqrt n\log P)\)。

• 对于 \(p>\sqrt n\) 的，有 \(k=1\)，于是变形为

  \[\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) \]

  于是只需要求出 \(\left(\sum\limits_{\sqrt n< p\leq n}p\times Sum(n/p)\right)\bmod (P-1)\)，注意是模 \(\varphi(P)=P-1\)。

  直接整除分块，需要解决的问题是区间质数和，用 min_25 筛解决即可。

然后看第二部分：

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

注意到只有当 \(i\) 没有平方因子时才有贡献，设 \(i=p_1p_2\cdots p_m\)，则 \(d(i)=\prod\limits_{j=1}^m(1+1)=2^m\)，于是设 \(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) \]

于是只需求出 \(\left(\sum\limits_{i=1}^{n}d'(i)\mu(i)\right)\bmod (P-1)\)。

设 \(f(n,i)=\sum\limits_{x=1}^n[(x\in \mathbb{P})\lor(\operatorname{lpf}(x)>p_i)]d'(x)\mu(x)\)，\(g(n,i)=\sum\limits_{x=1}^n[(x\in \mathbb{P})\lor(\operatorname{lpf}(x)>p_i)]\mu(x)\)，那么：

\[\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,|\mathbb{P}|)=\sum\limits_{x=1}^n[x\in \mathbb{P}]d'(x)\mu(x)=-|\mathbb{P}|\) 时（其中 \(\mathbb{P}\) 是小于等于 \(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;
}