【51nod1965】奇怪的式子
Portal --> 51nod1965
Solution
怎麽說呢。。這題。。做的有點痛苦。。
首先看這個式子長得。。比較奇怪,指數裏面那個加號有點煩人,而且這個函數不是一個積性函數也有點煩人,那就考慮把這個函數拆成兩部分來算
\[
\prod\limits_{i=1}^{n}\sigma_0 (i)^{\mu (i)+i}=\prod\limits_{i=1}^{n}\sigma_0(i)^i\cdot\prod\limits_{i=1}^{n}\sigma_0(i)^{\mu (i)}
\]
首先看第一部分的處理,我們記\(Sum(i)=\frac{n(n+1)}{2}\),那麽有:
具體一點的話就是,因為這裏是連乘,然後\(\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\)),所以:
\[
\begin{aligned}
&\prod\limits_{p>\sqrt n}\prod\limits_{k}(k+1)^W=\prod\limits_{p>\sqrt n}2^{W'}\&W'=Sum(\lfloor \frac{n}{p} \rfloor)\\end{aligned}
\]
這個部分的話,因為\(\lfloor \frac{n}{p}\rfloor\)的取值只有\(2\sqrt n\)個,所以我們可以直接用莫比烏斯反演中的那個常用技巧,分塊一下求就可以了
然後整個過程中我們需要預先處理出來的東西有兩個:約數個數和約數和
這兩個東西用min_25篩一下就好了
然而!
這樣我們求出來的是\(2\)的指數部分,這個數在計算的過程中如果不加取模操作的話會奇大無比,這就很尷尬
註意到題目中的模數是一個素數,那麽我們可以用一下由歐拉定理得到的一個結論:
\[
x^y\equiv x^{y\ mod\ \phi(p)}(mod\ p)
\]
也就是說我們可以對指數部分的計算\(\%\phi(p)\)(也就是\(p-1\)),並不會影響結果,問題解決ovo
然後看第二部分,這一個部分看上去會比較麻煩(實際上也確實有點麻煩)
把式子單獨列出來:
\[
\prod\limits_{i=1}^{n}\sigma_0(i)^{\mu(i)}
\]
首先來分析一下什麽樣的\(i\)才會被算進去
由於\(\mu(i)\)有一個很好的性質:當\(i\)有平方因子的時候\(\mu(i)=0\),所以我們接下來對\(i\)的討論可以多一個前提就是\(i\)分解質因數之後每一個素數的指數都是\(1\)
那麽沿用處理第一部分的時候的想法,我們將每個\(i\)分解質因數拆開來寫
用\(d(i)\)表示\(i\)的質因數個數,那麽可以得到:
\[
\prod\limits_{i=1}^{n}\sigma_0(i)^{\mu(i)}=\prod\limits_{i=1}^{n}2^{d(i)\cdot\mu(i)}
\]
有了上面的歐拉定理得到的結論支撐,我們就可以放心大膽地先求指數再求結果啦
用類似min_25篩中設狀態的方式:
記\(minp(i)\)表示\(i\)的最小質因子,\(d(i)\)表示\(i\)的質因數個數,\(P\)為素數列表
首先大力定義
\[
f(i,j)=\sum\limits_{k=1}^{i}[i\in Prime或minp(i)>P_j]\mu(i)\cdot d(i)
\]
然後為了方便轉移我們還需要定義一個輔助數組
\[
f1(i,j)=\sum\limits_{k=1}^{i}[i\in Prime或minp(i)>P_j]\mu(i)
\]
那麽我們可以得到轉移:
\[
\begin{aligned}
&f1(i,j)=f1(i,j+1)+(-1)*(f1(\lfloor \frac{n}{P_j}\rfloor,j+1)-f1(P_j,j+1))\&f(i,j)=f(i,j+1)+(-1)*(f(\lfloor \frac{n}{P_j}\rfloor,j+1)-f(P_j,j+1))+(-1)*(f1(\lfloor \frac{n}{P_j}\rfloor,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(\lfloor\frac{n}{P_j}\rfloor,j+1)\)和\(f1(\lfloor\frac{n}{P_j}\rfloor,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\n",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);
}
【51nod1965】奇怪的式子