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

一些记号：

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

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

首先看第一部分：
$$\prod _{i=1}^{n}d(i{)}^i=\prod _{p^k\le n}(k+1{)}^{cnt}$$
其中 $cnt=p^{k}Sum(\lfloor \frac{n}{p^k}\rfloor )-p^{k+1}Sum(\lfloor \frac{n}{p^{k+1}}\rfloor )$，$Sum(n)=\sum _{i=1}^{n}i=\frac{n(n+1)}{2}$。

我们按 $p$ 分类讨论：

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

• 对于 $p>\sqrt{n}$ 的，有 $k=1$，于是变形为
  $$\prod _{\sqrt{n}<p\le n}{2}^{p\times Sum(n/p)}=\mathrm{pow}\left(2,\sum _{\sqrt{n}<p\le n}p\times Sum(n/p)\right)$$
  于是只需要求出 $\left(\sum _{\sqrt{n}<p\le n}p\times Sum(n/p)\right)\mathrm{m}\mathrm{o}\mathrm{d}(P-1)$，注意是模 $\phi (P)=P-1$。

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

然后看第二部分：
$$\prod _{i=1}^{n}d(i{)}^{\mu (i)}$$
注意到只有当 $i$ 没有平方因子时才有贡献，设 $i=p_1{p}_2\cdots p_m$，则 $d(i)=\prod _{j=1}^m(1+1)=2^m$，于是设 $d^{\prime }(i)$ 表示 $i$ 的不同的质因数个数，有：
$$\prod _{i=1}^{n}d(i{)}^{\mu (i)}=\prod _{i=1}^n{\left(2^{d^{\prime }(i)}\right)}^{\mu (i)}=\mathrm{pow}\left(2,\sum _{i=1}^n{d}^{\prime }(i)\mu (i)\right)$$
于是只需求出 $\left(\sum _{i=1}^n{d}^{\prime }(i)\mu (i)\right)\mathrm{m}\mathrm{o}\mathrm{d}(P-1)$。

设 $f(n,i)=\sum _{x=1}^n[(x\in \mathbb{P})\vee (\mathrm{lpf}(x)>p_i)]d^{\prime }(x)\mu (x)$，$g(n,i)=\sum _{x=1}^n[(x\in \mathbb{P})\vee (\mathrm{lpf}(x)>p_i)]\mu (x)$，那么：
$$\begin{array}{rl}f(n,i)& =f(n,i+1)-\left(\right(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)\left)\right)\\ g(n,i)& =g(n,i+1)-(g(n/p_{i+1},i+1)-g(p_{i+1},i+1))\end{array}$$
初始化 $f(n,|\mathbb{P}|)=\sum _{x=1}^n[x\in \mathbb{P}]d^{\prime }(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;
}