1. 程式人生 > >LOJ.6053.簡單的函式(Min_25篩)

LOJ.6053.簡單的函式(Min_25篩)

題目連結
Min_25篩見這裡:


\(Description\)

給定\(n\),求積性函式\(f(p^c)=p\oplus c\)的字首和。\(\oplus\)表示異或運算。
\(n\leq 10^{10}\)

\(Solution\)

所求積性函式為\(f(p^c)=p\oplus c,\quad p\in Prime\)
先考慮質數的貢獻。因為除\(2\)以外的質數\(p\)都是奇數,所以\(f(p)=p-1\)\(f(2)\)就是\(2+1=3\)了。

不妨先把\(f(2)\)也看做\(f(2)=p-1\)
這樣\(f(p)=p-1\),但還不是積性函式,但是我們可以拆成兩個積性函式的和:\(f(p)=g(p)-h(p)\)

,其中\(g(p)=p,\ h(p)=[p是質數]\)。分別計算這兩個函式的字首和。

然後,首先計算初值\(g(n,0)\),把所有合數看做質數,那麼\(g(n),h(n)\)的字首和也就是初值分別是\(\frac{n(n+1)}{2}-1\)\(n-1\)(不考慮\(f(1)\))。

然後計算\(g(x,|P|)\)。在外層列舉\(j\)把第二維滾動掉。
另外每次從\(\frac{n}{P_j}\)轉移,都是整除,結果只有\(O(sqrt(n))\)種(好像是\(2sqrt(n)\)?不知道為什麼...),所以可以先離散化。
然後套式子得到\(g(x,|P|)\)

然後套式子計算\(S(n,1)\)

。當\(j=1\)時給結果加個\(2\),因為\(f(P_1)=f(2)\)是按\(2-1=1\)計算的,應該是\(3\)
這裡直接遞迴算就行了,而且不需要記憶化。

//1214ms    4.12M
#include <cmath>
#include <cstdio>
#include <algorithm>
#define mod 1000000007
#define Mod(x) x>=mod&&(x-=mod)
typedef long long LL;
const int N=2e5+5;

int Sqr,cnt,P[N>>2],sp[N],id1[N],id2[N],g[N],h[N];
LL n,w[N];
bool notP[N];

void Init(int n)
{
    notP[1]=1;
    for(int i=2; i<=n; ++i)
    {
        if(!notP[i]) P[++cnt]=i, sp[cnt]=sp[cnt-1]+i, Mod(sp[cnt]);
        for(int j=1; j<=cnt && i*P[j]<=n; ++j)
        {
            notP[i*P[j]]=1;
            if(!(i%P[j])) break;
        }
    }
}
int S(LL x,int y)
{
    if(x<=1||(y!=cnt+1&&P[y]>x)) return 0;//注意y=cnt+1時也需要計算!
    int k=x<=Sqr?id1[x]:id2[n/x];
    LL res=g[k]-sp[y-1]-h[k]+y-1;//g-h
    if(y==1) res+=2;//f(2)還是就放到裡面算吧 否則要判下n<2。
    for(int i=y; 1ll*P[i]*P[i]<=x; ++i)
    {
        LL p=P[i],p1=p,p2=p*p;
        for(int e=1; p2<=x; ++e,p1=p2,p2*=p)
            res+=1ll*(p^e)*S(x/p1,i+1)%mod+(p^(e+1));
        res%=mod;
    }
    return res%mod;
}

 main()
{
    scanf("%lld",&n); Init(Sqr=sqrt(n+0.5));
    int m=0;
    for(LL i=1,j; i<=n; i=j+1)
    {
        w[++m]=n/i, j=n/w[m];
        if(w[m]<=Sqr) id1[w[m]]=m;
        else id2[j]=m;
        g[m]=w[m]&1?w[m]%mod*((w[m]+1>>1)%mod)%mod-1:(w[m]>>1)%mod*(w[m]%mod+1)%mod-1;
        h[m]=(w[m]-1)%mod;
    }
    P[cnt+1]=1e9, w[m+1]=-1;
    for(int j=1; j<=cnt; ++j)
    {
        int pj=P[j]; LL lim=1ll*pj*pj;
        for(int i=1; lim<=w[i]; ++i)
        {
            int k=w[i]/pj<=Sqr?id1[w[i]/pj]:id2[n/(w[i]/pj)];//n/(w/pj)! id1[x]=x,但id2[x]的編號是對[n/x]的。
            (g[i]-=1ll*pj*(g[k]-sp[j-1])%mod)%=mod;//g[k]-sp[j-1]有可能是負的,如果+mod會爆int!
            h[i]+=mod-h[k]+j-1, Mod(h[i]);
        }
    }
    printf("%d\n",((S(n,1)+1)%mod+mod)%mod);

    return 0;
}