积性函数小结

定义及前置芝士

数论函数：指定义域为正整数、定义域为复数的函数，在OI中这类函数的值域极大多数也为整数

积性函数：指对于数论函数(f(x))和任意一对互质整数(a,b),均有性质(f(ab)=f(a)f(b))。

莫比乌斯反演和狄利克雷卷积：很久以前自己写过的博客

在OI中，有一类经典的问题是求某个数论函数的值（或前缀和），它可能是像(mu,varphi)这样广为人所熟知的数论函数，也有可能是一些连完整的表达式都写不出来的函数，但是它们都有一个特点——它们都是积性函数。于是我们可以考虑一些特殊的方法将它们“筛”出来

线性筛求积性函数

积性函数的定义是基于互质的，而互质的一个特例就是其中一个数是质数，于是我们考虑在筛素数的时候来同时求出这个积性函数的值

具体的，我们记录(low[i])表示(i)的最小质因子的最大次幂的数使得它仍然是(i)的因数，这个根据线筛素数可以很方便的维护，接下来就是根据当前的数和筛它的质数是否互质来处理。

根据上面的叙述不难发现我们其实只是关心此积性函数的下列值:(f(1),f(p)(质数的初值),f(p^k)(若用来筛的素数与被筛的数不互质))，当然实际上很多时候(f(p^k))的值都是从(f(p^{k-1}))递推得到的

奉上一段伪代码

rep(i,2,N)
{
    if (!nopri[i]) {pri[++tot]=i;low[i]=i;f[i]=...(根据定义);}
    int j;
    for (j=1;j<=tot&&i*pri[j]<=N;j++)
    {
        nopri[i*pri[j]]=1;
        if (i%pri[j])
        {
            low[i*pri[j]]=low[i];
            f[i*pri[j]]=f[i]*f[pri[j]];
        }
        else
        {
            low[i*pri[j]]=low[i]*pri[j];
            f[i*pri[j]]=f[i/low[i]]*f[low[i]*pri[j]]\注意后者可能需要根据定义求得
            break;
        }
    }
}

(未尝试编译，可能会有问题，欢迎指出！）

杜教筛

线性筛的时间复杂度显然是(O(n))的，在许多问题中已经是十分优秀的时间了

然而事实是总有毒瘤出题人写什么(nleq 10^9)甚至更大，于是我们需要一些跟更优秀的做法

dls为我们提供了一种在(O(n^{frac{2}{3}}))的时间内求积性函数前缀和的方法

对于一个积性函数(f)，找到两个方便求前缀和的积性函数(g,h)使得(f*g=h)，于是

[egin{aligned} sum_{i=1}^nf(i)=&sum_{i=1}^n(h(i)-sum_{d|i,d eq i}(f(d)g(frac{i}{d}))\ =&H(n)-sum_{i=1}^nsum_{d|i,d eq i}(f(d)g(frac{i}{d}))\ =&H(n)-sum_{d=1}^nf(d)sum_{i=2}^{lfloor frac{n}{d} floor}g(i)\ =&H(n)-sum_{i=2}^ng(i)sum_{d=1}^{lfloorfrac{n}{i} floor}f(i)\ =&H(n)-sum_{i=2}^ng(i)F(lfloorfrac{n}{i} floor) end{aligned} ]

注意上面的推导过程中用大写字母表示了对应积性函数的前缀和

于是我们可以考虑用线筛处理一部分的(f)，然后再运用记忆化搜索和数论分块求解最后的答案，经证明，当(N=n^frac{2}{3})时时间最优，为(O(n^{frac{2}{3}}))

例：求(mu)，(varphi)的前缀和，(nleq 10^9)

考虑杜教筛，(mu*I=epsilon),(varphi*I=id_1)，套用上面的推导过程即可

因为懒得手写hash于是就用了unordered_map

#include<iostream>
#include<string.h>
#include<string>
#include<stdio.h>
#include<algorithm>
#include<math.h>
#include<vector>
#include<queue>
#include<map>
#include<unordered_map>
using namespace std;
const int maxd=1000000007,N=5000000;
const double pi=acos(-1.0);
typedef long long ll;
int n,mu[5005000],pri[1001000],cnt=0;
ll phi[5005000];
bool ispri[5005000];
unordered_map<int,int> ansmu;
unordered_map<int,ll> ansphi;

int read()
{
    int x=0,f=1;char ch=getchar();
    while ((ch<'0') || (ch>'9')) {if (ch=='-') f=-1;ch=getchar();}
    while ((ch>='0') && (ch<='9')) {x=x*10+(ch-'0');ch=getchar();}
    return x*f;
}

void init()
{
    mu[1]=1;phi[1]=1;
    int i,j;
    for (i=2;i<=N;i++)
    {
        if (!ispri[i]) {pri[++cnt]=i;mu[i]=-1;phi[i]=i-1;}
        for (j=1;j<=cnt && i*pri[j]<=N;j++)
        {
            ispri[i*pri[j]]=1;
            if (i%pri[j]==0) {phi[i*pri[j]]=phi[i]*pri[j];break;}
            else {mu[i*pri[j]]-=mu[i];phi[i*pri[j]]=phi[pri[j]]*phi[i];}
        }
    }
    for (i=1;i<=N;i++) {phi[i]+=phi[i-1];mu[i]+=mu[i-1];}
}

int get_mu(int n)
{
    if (n<=N) return mu[n];
    if (ansmu[n]) return ansmu[n];
    int ans=1;
    unsigned int l,r;
    for (l=2;l<=n;l=r+1)
    {
        r=n/(n/l);
        ans-=(r-l+1)*get_mu(n/l);
    }
    ansmu[n]=ans;
    return ans;
}

ll get_phi(int n)
{
    if (n<=N) return phi[n];
    if (ansphi[n]) return ansphi[n];
    ll ans=((1ll+n)*n)/2ll;
    unsigned int l,r;
    for (l=2;l<=n;l=r+1)
    {
        r=n/(n/l);
        ans-=1ll*(r-l+1)*get_phi(n/l);
    }
    ansphi[n]=ans;
    return ans;
}

void work()
{
    int i,T=read();
    while (T--)
    {
        int n=read();
        ll ans1=get_phi(n);
        int ans2=get_mu(n);
        printf("%lld %d
",ans1,ans2);
    }
}

int main()
{
    init();
    work();
    return 0;
}

(min25)筛

杜教筛求积性函数的前缀和已经能做到优秀的(O(n^{frac{2}{3}}))了，但是还是有一个很大的缺陷：我们必须找到两个便于求前缀和的积性函数(g,h)使得它们与所求的积性函数(f)有(f*g=h)的关系，这在一些常规的积性函数中比较常见，但是在某些毒瘤出题人所给出的式子中却无法找到（有的函数甚至都不会给你完整的定义式），看起来这个时候我们就不得不回到最开始的线筛了，有没有更优秀的做法呢？

(min25)提出了一种积性函数筛法可以在(O(frac{n^{0.75}}{log n}))的时间内求出积性函数的前缀和（但是我并不会证这个复杂度，另一种广为流传的时间复杂度为(O(n^{1-epsilon}))

(min25)主要有以下两个步骤
（1）对(forall x=lfloorfrac{n}{i} floor)，求(sum_{i=1}^x[iin P]i^k)，其中大写(P)表示素数集合，(k)是一个常数，其取值往往会根据题中所给函数的形式（主要是(f(p^k))）而确定，有时候(k)的取值会多于一个

我们记(g(n,j)=sum_{i=1}^n[iin P or min(i)geq p_j]i^k)，其中(min(i))表示(i)的最小质因子，(p_j)表示第(j)个素数，很明显我们在这一步最终要求的是(forall x=lfloorfrac{n}{i} floor,g(x,|P|))的值

如何求(g)的值呢？我们考虑一下它的实际意义，让我们回到一个古老的筛素数方法——埃氏筛法，它考虑的是选取当前还未被筛掉的最小的数，认定它是一个质数，并用它筛掉剩余数中是该质数的倍数的数。我们发现(g(n,j))就是在用(j)个素数进行了埃筛之后留下的数的(f)值

根据当前(j)的值进行分类讨论
i) 若(p_jleq sqrt{n})，那么我们在这一次中会筛掉最小质因子恰好是(p_j)的数，也就是在除掉一个(p_j)后剩下的数的最小质因子必然大于等于(p_j)，也就是(g(lfloorfrac{n}{p_j} floor,j-1)-sum_{i=1}^{j-1}f(p_i))

ii)若(p_j>sqrt{n})，那么这一次筛掉的数至少为(p_j^2)，而这是一个大于(n)的数，所以这一次筛不会改变留下来的数，即为(g(n,j-1))

综上所述，有

[g(n,j)= egin{cases} g(n,j-1) & p_j>sqrt{n}\ g(n,j-1)-p_j^k[g(lfloorfrac{n}{p_j} floor,j-1)-g(p_j-1,j-1)]& p_jleq sqrt{n} end{cases} ]

所以我们可以直接筛出(sqrt{n})以内的素数，然后递推得到(g(n,|P|))即可

（2）求(S(n,j)=sum_{i=1}^n[min(i)geq p_j]f(i))
那么很明显(sum_{i=1}^nf(i)=S(n,1)+f(1))
求(S(n,j))的话考虑两部分，一部分是质数的时候，直接减去(<p_j)的质数即可，这一部分是(g(n,|P|)-sum_{i=1}^{j-1}f(p_i))

接下来考虑合数，注意到满足条件的数一定可以写成(p_k^e*x)的形式（其中(min(x)geq p_k)），于是我们可以枚举(k)和(e)，分开计算(x)中是否还有(p_k)即可
写在一起就是

[S(n,j)=g(n,|P|)-sum_{i=1}^{j-1}f(p_i)+sum_{k=j}^{p_k^2leq n}sum_{e=1}^{p^{e+1}leq n}(S(frac{n}{p_k^e},k+1)f(p_k^e)+f(p_k^{e+1})) ]

边界条件就是(nleq 0)或者是(p_j>n)时(S(n,j)=0)

例题：DIVCNT2

据说可以大力杜教筛。。。感觉智商太低不大会

去spoj官网看一下的话发现这个题数据单位还是比较资瓷min25筛的，于是就写了一发min25

首先积性函数很容易就能证了，其次就有(f(1)=1,f(p)=3,f(p^k)=2k+1),由于此处(f(p))与(f(p^k))均与(p)无关，于是我们直接求(sum_{i=1}^{lfloorfrac{n}{j} floor}[iin P])即可

献上常数巨大的代码

#pragma GCC optimize(2)
#include<iostream>
#include<string.h>
#include<string>
#include<stdio.h>
#include<algorithm>
#include<vector>
#include<math.h>
#include<queue>
#include<set>
#include<map>
using namespace std;
typedef long long ll;
typedef long double db;
typedef unsigned long long ull;
const int N=2000000;
const db pi=acos(-1.0);
#define lowbit(x) (x)&(-x)
#define sqr(x) (x)*(x)
#define rep(i,a,b) for (register int i=a;i<=b;i++)
#define per(i,a,b) for (register int i=a;i>=b;i--)
#define fir first
#define sec second
#define mp(a,b) make_pair(a,b)
#define pb(a) push_back(a)
#define maxd 998244353
#define eps 1e-8
int m,ptot=0,id1[N+10],id2[N+10],tot=0,pri[N+10];
ll n,w[N+10];
bool nopri[N+10];
ull g[N+10];

ll read()
{
    ll x=0;int f=1;char ch=getchar();
    while ((ch<'0') || (ch>'9')) {if (ch=='-') f=-1;ch=getchar();}
    while ((ch>='0') && (ch<='9')) {x=x*10+(ch-'0');ch=getchar();}
    return x*f;
}

void sieve()
{
    rep(i,2,N)
    {
        if (!nopri[i]) pri[++ptot]=i;
        int j;
        for (j=1;j<=ptot && i*pri[j]<=N;j++)
        {
            nopri[i*pri[j]]=1;
            if (i%pri[j]==0) break;
        }
    }
}

void calc(ll n)
{
    m=sqrt(n);tot=0;
    ll l,r;
    for (l=1;l<=n;l=r+1)
    {
        r=n/(n/l);w[++tot]=(n/l);
        if (w[tot]<=m) id1[w[tot]]=tot;else id2[n/w[tot]]=tot;
        g[tot]=w[tot]-1;
    }
    int i,j;
    for (i=1;i<=ptot && 1ll*pri[i]*pri[i]<=n;i++)
    {
        for (j=1;j<=tot && 1ll*pri[i]*pri[i]<=w[j];j++)
        {
            int k;
            if (w[j]/pri[i]<=m) k=id1[w[j]/pri[i]];else k=id2[n/(w[j]/pri[i])];
            g[j]=g[j]-g[k]+(i-1);
        }
    }
}

ull s(ll lim,int cnt)
{
    if ((lim<=1) || (pri[cnt]>lim)) return 0;
    ull ans=0;
    if (lim<=m) ans=3ull*(g[id1[lim]]-(cnt-1));else ans=3ull*(g[id2[n/lim]]-(cnt-1));
    int i,j;
    for (i=cnt;i<=ptot&&1ll*pri[i]*pri[i]<=lim;i++)
    {
        ll pro=pri[i];
        for (j=1;pro*pri[i]<=lim;j++,pro=pro*pri[i])
        {
            ans+=s(lim/pro,i+1)*(j*2+1)+j*2+3;
        }
    }
    return ans;
}

signed main()
{
    sieve();
    int T=read();
    while (T--)
    {
        n=read();
        calc(n);
        printf("%llu
",s(n,1)+1);
    }
    return 0;
}