利用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) eq 0$ 的位置显然都是powerful number，我们只需枚举所有powerful number，算出对应的 $H$ 即可。

例题2  pe48?

求满足对质数 $p$，$F(p^d)=p^{d-[d mod p]}$ 的积性函数 $F$ 的前缀和。

$F(p)=1$。同上构造 $G(x)=1$ 即可。

例题3  loj6053

求满足对质数 $p$，$F(p^c)=p oplus c$ 的积性函数 $F$ 的前缀和。

对 $p eq 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<<"
";
}