利用powerful number求積性函數前綴和

很久沒更博客了，先水一篇再說。其實這個作法應該算是杜教篩的一個拓展。

powerful number的定義是每一個質因子次數都 $\geq 2$ 的數。首先，$\leq n$ 的powerful number個數是 $O(\sqrt{n})$ 的，這是由於全部powerful number顯然能夠表示成 $a^2b^3$，因此個數不超過 $\sum_{i=1}^{\sqrt{n}} (n/i^2)^{1/3}$，積分積一下就算出來了。求全部 $\leq n$ 的powerful number只要暴搜質因子分解式便可。

例題1  pe63?

有一個積性函數 $F$ 知足對於全部質數 $p$，$F(p^q)=p~(q \geq 1)$，求 $F$ 的前綴和。

咱們發現有一個跟它長得很像的積性函數 $G$！$G(x)=x$，咱們會求 $G$ 的前綴和！而且對於全部質數 $p$，$G(p)=F(p)=p$。

咱們求出 $H=F/G$，其中除法指的是狄利克雷除法，即狄利克雷卷積的逆運算。$H$ 也是一個積性函數，那麼由 $F(p^q)=\sum_{i=0}^q G(p^i)H(p^{q-i})$ 和 $H(1)=1$ 不難發現對於全部質數 $p$，$H(p)=0$。

咱們欲求的是 $\sum_{i=1}^n F(i)$，因爲 $F=H*G$（乘法爲狄利克雷卷積），那麼有 $\sum_{i=1}^n F(i)=\sum_{ij \leq n} H(i)G(j)=\sum_{i=1}^n H(i) \sum_{j=1}^{n/i} G(j)$。

因爲 $H(p)=0$，全部 $H(i) \neq 0$ 的位置顯然都是powerful number，咱們只需枚舉全部powerful number，算出對應的 $H$ 便可。

例題2  pe48?

求知足對質數 $p$，$F(p^d)=p^{d-[d \bmod p]}$ 的積性函數 $F$ 的前綴和。

$F(p)=1$。同上構造 $G(x)=1$ 便可。

例題3  loj6053

求知足對質數 $p$，$F(p^c)=p \oplus c$ 的積性函數 $F$ 的前綴和。

對 $p \neq 2$，$F(p)=p-1$。構造 $G=\varphi$ 便可，注意要特殊處理一下 $p=2$ 的情形。求歐拉函數的前綴和能夠杜教篩，這裏再也不贅述。 跑過min25篩是不可能的，這輩子都不可能跑過的

這裏給出loj6053的代碼：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MOD=1e9+7;
#define SZ 10000099
bool np[SZ];
int ph[SZ],ps[SZ/10],pn;
void shai()
{
    ph[1]=1;
    for(int i=2;i<SZ;++i)
    {
        if(!np[i]) ps[++pn]=i,ph[i]=i-1;
        for(int j=1;j<=pn&&i*ps[j]<SZ;++j)
        {
            np[i*ps[j]]=1;
            if(i%ps[j]==0)
            {
                ph[i*ps[j]]=ph[i]*ps[j];
                break;
            }
            ph[i*ps[j]]=ph[i]*(ps[j]-1);
        }
    }
    for(int i=1;i<SZ;++i)
        ph[i]=(ph[i-1]+ph[i])%MOD;
}
ll n,u,s[1005];
ll S2(ll a)
{
    ll b=a+1;
    if(a&1) b>>=1; else a>>=1;
    return (a%MOD)*(b%MOD)%MOD;
}
int h[100099][66],d[100099];
ll ans=0;
void dfs(ll x,ll v,int w)
{
    ans=(ans+v*((n/x<SZ)?ph[n/x]:s[n/(n/x)]))%MOD;
    if(w>1&&x>n/ps[w]/ps[w]) return;
    for(int s=w;s<=pn;++s)
    {
        if(s>1&&x*ps[s]*ps[s]>n) break;
        ll y=x*ps[s];
        for(int j=1;y<=n;++j,y*=ps[s])
        {
            if(d[s]<j)
            {
                ++d[s];
                ll F=ps[s]^j,G=ps[s]-1;
                for(int k=1;k<=j;++k)
                    F=(F-G%MOD*h[s][j-k])%MOD,G*=ps[s];
                h[s][j]=F;
            }
            if(!h[s][j]) continue;
            dfs(y,v*h[s][j]%MOD,s+1);
        }
    }
}
int main()
{
    for(int i=0;i<=100000;++i)
        h[i][0]=1;
    shai(); cin>>n;
    u=1; while(n/u>=SZ) ++u;
    for(int i=u;i>=1;--i)
    {
        //s[i]=phi(n/i)
        ll t=n/i,a=2,b,p; s[i]=S2(t);
        for(;a<=t;a=b+1)
            p=t/a,b=t/p,
            s[i]=(s[i]-(b-a+1)%MOD*((p<SZ)?ph[p]:s[b*i]))%MOD;
    }
    dfs(1,1,1);
    ans=(ans%MOD+MOD)%MOD;
    cout<<ans<<"\n";
}