【51nod1965】奇怪的式子

Portal --> 51nod1965

Solution

　　怎么说呢。。这题。。做的有点痛苦。。

　　首先看这个式子长得。。比较奇怪，指数里面那个加号有点烦人，而且这个函数不是一个积性函数也有点烦人，那就考虑把这个函数拆成两部分来算

[prodlimits_{i=1}^{n}sigma_0 (i)^{mu (i)+i}=prodlimits_{i=1}^{n}sigma_0(i)^icdotprodlimits_{i=1}^{n}sigma_0(i)^{mu (i)} ]

　　首先看第一部分的处理，我们记(Sum(i)=frac{n(n+1)}{2})，那么有：

[egin{aligned} &prodlimits_{i=1}^{n}sigma_0(i)^i=prodlimits_{p}prodlimits_{k}(k+1)^W\ &W=p^kcdot Sum(lfloor frac{n}{p^k} floor)-p^{k+1}cdot Sum(lfloor frac{n}{p^k+1} floor)\ end{aligned} ]

　　具体一点的话就是，因为这里是连乘，然后(sigma_0(x))是一个积性函数，所以我们考虑将每一个(i)分解质因数，然后按照分解的质因数把每个(sigma_0)写开

　　那么，枚举素数(p)，枚举这个素数的指数(k)，然后统计分解质因数后(p)的指数恰好为(k)的数的个数，然后算贡献就好了

　　然而这个式子还是不能直接用来求解第一部分的答案，所以我们还要进一步处理

　　考虑分成(p<=sqrt n)和(p>sqrt n)两类来看：

　　首先对于(p<=sqrt n)的部分，是可以直接计算而不会T掉的

　　然后对于(p>sqrt n)的部分，(p)的指数(k)只能是(1)（因为(p>sqrt n)），所以：

[egin{aligned} &prodlimits_{p>sqrt n}prodlimits_{k}(k+1)^W=prodlimits_{p>sqrt n}2^{W'}\ &W'=Sum(lfloor frac{n}{p} floor)\ end{aligned} ]

　　这个部分的话，因为(lfloor frac{n}{p} floor)的取值只有(2sqrt n)个，所以我们可以直接用莫比乌斯反演中的那个常用技巧，分块一下求就可以了

　　然后整个过程中我们需要预先处理出来的东西有两个：约数个数和约数和

　　这两个东西用min_25筛一下就好了

　　然而！

　　这样我们求出来的是(2)的指数部分，这个数在计算的过程中如果不加取模操作的话会奇大无比，这就很尴尬

　　注意到题目中的模数是一个素数，那么我们可以用一下由欧拉定理得到的一个结论：

[x^yequiv x^{y mod phi(p)}(mod p) ]

　　也就是说我们可以对指数部分的计算(\%phi(p))（也就是(p-1)），并不会影响结果，问题解决ovo

　　
　　
　　

　　然后看第二部分，这一个部分看上去会比较麻烦（实际上也确实有点麻烦）

　　把式子单独列出来：

[prodlimits_{i=1}^{n}sigma_0(i)^{mu(i)} ]

　　首先来分析一下什么样的(i)才会被算进去

　　由于(mu(i))有一个很好的性质：当(i)有平方因子的时候(mu(i)=0)，所以我们接下来对(i)的讨论可以多一个前提就是(i)分解质因数之后每一个素数的指数都是(1)

　　那么沿用处理第一部分的时候的想法，我们将每个(i)分解质因数拆开来写

　　用(d(i))表示(i)的质因数个数，那么可以得到：

[prodlimits_{i=1}^{n}sigma_0(i)^{mu(i)}=prodlimits_{i=1}^{n}2^{d(i)cdotmu(i)} ]

　　有了上面的欧拉定理得到的结论支撑，我们就可以放心大胆地先求指数再求结果啦

　　用类似min_25筛中设状态的方式：

　　记(minp(i))表示(i)的最小质因子，(d(i))表示(i)的质因数个数，(P)为素数列表

　　首先大力定义

[f(i,j)=sumlimits_{k=1}^{i}[iin Prime或minp(i)>P_j]mu(i)cdot d(i) ]

　　然后为了方便转移我们还需要定义一个辅助数组

[f1(i,j)=sumlimits_{k=1}^{i}[iin Prime或minp(i)>P_j]mu(i) ]

　　那么我们可以得到转移：

[egin{aligned} &f1(i,j)=f1(i,j+1)+(-1)*(f1(lfloor frac{n}{P_j} floor,j+1)-f1(P_j,j+1))\ &f(i,j)=f(i,j+1)+(-1)*(f(lfloor frac{n}{P_j} floor,j+1)-f(P_j,j+1))+(-1)*(f1(lfloor frac{n}{P_j} floor,j+1)-f1(P_j,j+1)) end{aligned} ]

　　初始化(f1(i,max)=f(i,max)=(-1)*[1,i])范围内素数个数，(max)表示的是(>=sqrt n)的第一个素数在列表中的下标

　　具体一点的话就是，首先所有的((-1))其实都是(mu (p))的值，因为是素数的(mu)值所以是(-1)

　　然后转移的话其实就是要加上([1,i])范围内的(minp(x)=P_j)的数的贡献，为了保证有(P_j)的存在并且是最小的，所以调用的是(f(lfloorfrac{n}{P_j} floor,j+1))和(f1(lfloorfrac{n}{P_j} floor,j+1))，然后还要保证素数不会被多算一次所以要去掉素数部分也就是(f(P_j,j+1))和(f1(P_j,j+1))

　　然后就可以快乐dp得出答案啦

​　　

　　一些实现的细节问题：

　　　　1.这题的模数比较巨大，两个数的乘积可能爆long long，所以要。。用一个奇妙黑科技，也就是用long double 模拟取模这样来取模qwq

　　　　2.注意在算的过程中的模数到底是用(phi(MOD))还是(MOD)要看清楚。。调了好久。。

　　　　3.在求的过程中因为数字都比较巨大所以。。要积极取模

（调试调到后来直接放弃思考见到乘号就取模导致代码长得十分不可描述并且常数巨大qwq）

​ 　　

　　代码大概长这个样子

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define ll long long
#define ld long double 
using namespace std;
const int MAXN=1e6+10;
const ll MOD=1e12+39,PHIMOD=MOD-1;
ll pcnt[MAXN*2],psum[MAXN*2],g[MAXN*2],f[MAXN*2],f1[MAXN*2];
ll P[MAXN],rec[MAXN*2];
int loc1[MAXN],loc2[MAXN];
bool vis[MAXN];
ll n,m,sq,cnt,cntval,ans,T;
void prework(int n);
void init_loc();
void get_g();
int Pos(ll x){return x<=sq?loc1[x]:loc2[n/x];}
ll Sum(ll l,ll r);
ll ksm(ll x,ll y);
ll mul(ll x,ll y,ll mod){
	ld tmp=(ld)x*(ld)y;
	return (x*y-((ll)(tmp/(ld)mod)*mod)+mod)%mod;
}
void solve1();
void solve2();

int main(){
#ifndef ONLINE_JUDGE
	freopen("a.in","r",stdin);
#endif
	scanf("%d",&T);
	prework(MAXN);
	for (int o=1;o<=T;++o){
		scanf("%lld",&n);
		sq=sqrt(n)+0.5;
		m=upper_bound(P+1,P+1+cnt,sq)-P-1;
		init_loc();
		get_g();
		ans=1;
		solve1();
		solve2();
		printf("%lld
",ans);
	}
}

void prework(int n){
	cnt=0;
	for (int i=2;i<=n;++i){
		if (!vis[i])
			P[++cnt]=i;
		for (int j=1;j<=cnt&&P[j]*i<=n;++j){
			vis[i*P[j]]=true;
			if (i%P[j]==0) break;
		}
	}
}

void init_loc(){
	cntval=0;
	for (ll i=1,pos;i<=n;i=pos+1){
		rec[++cntval]=n/i;
		pos=n/(n/i);
	}
	reverse(rec+1,rec+1+cntval);
	for (int i=1;i<=cntval;++i)
		if (rec[i]<=sq) loc1[rec[i]]=i;
		else loc2[n/rec[i]]=i;
}

void get_g(){
	for (int i=1;i<=cntval;++i)
		pcnt[i]=rec[i]-1,psum[i]=Sum(1,rec[i])-1;
	for (int j=1;j<=m;++j)
		for (int i=cntval;i>=1&&rec[i]>=1LL*P[j]*P[j];--i){
			pcnt[i]-=pcnt[Pos(rec[i]/P[j])]-pcnt[Pos(P[j]-1)];
			psum[i]=(psum[i]+PHIMOD-P[j]*(psum[Pos(rec[i]/P[j])]+PHIMOD-psum[Pos(P[j]-1)])%PHIMOD)%PHIMOD;
		}
}

ll Sum(ll l,ll r){
	ll tmp=l+r,tmp1=r-l+1;
	if (tmp&1) tmp1/=2;
	else tmp/=2;
	return mul(tmp,tmp1,PHIMOD);
}

ll ksm(ll x,ll y){
	ll ret=1,base=x;
	for (;y;y>>=1,base=mul(base,base,MOD))
		if (y&1) ret=mul(ret,base,MOD);
	return ret;
}

void solve1(){//prod sigma i^i
	ll tmp;
	for (int i=1;i<=m;++i){
		tmp=P[i];
		for (int k=1;tmp<=n;tmp*=P[i],++k)
			ans=mul(ans,ksm(k+1,(tmp*Sum(1,n/tmp)%PHIMOD-tmp*P[i]%PHIMOD*Sum(1,n/(tmp*P[i]))%PHIMOD)+PHIMOD)%PHIMOD,MOD);
	}
	tmp=0;
	for (int i=Pos(sq)+1;i<=cntval;++i)
		tmp+=mul(Sum(1,n/rec[i]),(psum[i]-psum[i-1]+PHIMOD)%PHIMOD,PHIMOD);
	ans=mul(ans,ksm(2,tmp),MOD);
}

void solve2(){//prod sigma i^{mu i}
	for (int i=2;i<=cntval;++i) f1[i]=f[i]=PHIMOD-pcnt[i];
	for (int j=m;j>=1;--j){
		for (int i=cntval;i>=2&&rec[i]>=P[j]*P[j];--i){
			f1[i]=(f1[i]+(-1LL)*(f1[Pos(rec[i]/P[j])]-f1[Pos(P[j])]))%PHIMOD;
			f[i]=(f[i]+(-1LL)*(f[Pos(rec[i]/P[j])]-f[Pos(P[j])])+(-1LL)*(f1[Pos(rec[i]/P[j])]-f1[Pos(P[j])]))%PHIMOD;
			if (f1[i]<0) f1[i]+=PHIMOD;
			if (f[i]<0) f[i]+=PHIMOD;
		}
	}
	ans=mul(ans,ksm(2,f[Pos(n)]),MOD);
}