LG4213 【【模板】杜教筛(Sum)】

sum\(\mu\)求法


\[S(n)=\sum_{i=1}^n \mu(i)\]
回顾公式
\[\sum_{d|n}\mu(d)=[n=1]\]
对\(n\)求和
\[\sum_{i=1}^n\sum_{d|i}\mu(d)=1\]
换一种求和
\[\sum_{i=1}^n\sum_{d=1}^{\lfloor n/i \rfloor}\mu(d)=1\]
拆成两部分
\[\sum_{i=1}^n\mu(i)=1-\sum_{i=2}^n\sum_{d=1}^{\lfloor n/i \rfloor}\mu(d)\]
注意到右边遇到的都是子问题
\[S(n)=1-\sum_{i=2}^nS(\lfloor\frac{n}{i}\rfloor)\]
注意到连除有结合律
\[\lfloor\frac{\lfloor\frac{a}{b}\rfloor}{c}\rfloor=\lfloor\frac{a}{bc}\rfloor\]

时间复杂度

本质上只需要计算\(S(\lfloor\frac{n}{i}\rfloor)(1\leq i\leq n)\)。

大家都知道\(\lfloor\frac{n}{i}\rfloor\)本质不同的结果只有\(2\sqrt{n}\)个,整除分块的套路。

在已知\(S(\lfloor\frac{n}{i}\rfloor)(2\leq i\leq n)\)的情况下,计算\(S(n)\)的复杂度为\(O(\sqrt{n})\)。

所以总时间复杂度为
\[\sum_{i=1}^{\sqrt{n}}\sqrt{\frac{n}{i}}+\sum_{i=1}^{\sqrt{n}}\sqrt{i}\]
把求和换成积分的技巧,原式等于
\[\int_{0}^{\sqrt{n}}\sqrt{\frac{n}{x}}dx+\int_{0}^{\sqrt{n}}\sqrt{x}dx\]
不定积分的结果为
\[\int\sqrt{\frac{n}{x}}dx=2\sqrt{nx}\]
\[\int\sqrt{x}dx=\frac{2}{3}x^{\frac{3}{2}}\]
牛顿-莱布尼茨定理求出定积分\(=\frac{8}{3}n^{\frac{3}{4}}=O(n^{\frac{3}{4}})\)

引入参数\(B\),假设\(1,2,\dots,\frac{n}{B}\)用线性筛,其余部分用原办法,总复杂度为
\[\sum_{i=1}^B\sqrt{\frac{n}{i}}+\frac{n}{B}=2\sqrt{nB}+\frac{n}{B}\]
当\(B=O(n^{\frac{1}{3}})\)时,两部分相等,取到最小值,时间复杂度为\(O(n^{\frac{2}{3}})\)。

sum\(\varphi\)求法

同样的技巧可以用在欧拉函数上
\[\sum_{d|n}\varphi(n)=n\]
只是右侧求和要换成\(n(n+1)/2\)。

代码

借鉴@lazyzhong的代码,不过我都是用上述办法算的两个函数,没有嵌套

#include<cstdlib>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<ctime>
#include<iostream>
#include<string>
#include<vector>
#include<list>
#include<deque>
#include<stack>
#include<queue>
#include<map>
#include<set>
#include<algorithm>
#include<complex>
#define rg register
#pragma GCC optimize ("O3")
using namespace std;
template<class T> inline T read(T&x){
    T data=0;
    int w=1;
    char ch=getchar();
    while(!isdigit(ch))
    {
        if(ch=='-')
            w=-1;
        ch=getchar();
    }
    while(isdigit(ch))
        data=10*data+ch-'0',ch=getchar();
    return x=data*w;
}
typedef long long ll;
const int INF=0x7fffffff;

const int MAXN=1700000; // 2147483647^{2/3}=1664510.6449518746191262425357001
int prime[MAXN],pcnt; // use prime as well as isprime
int mu[MAXN];
ll phi[MAXN];

inline void sieve()
{
    fill(prime,prime+MAXN,1);
    mu[1]=1,phi[1]=1,pcnt=0;
    for(rg int i=2;i<MAXN;++i)
    {
        if(prime[i])
        {
            prime[++pcnt]=i,mu[i]=-1,phi[i]=i-1;
        }
        for(rg int j=1;j<=pcnt&&i*prime[j]<MAXN;++j)
        {
            prime[i*prime[j]]=0;
            if(i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
            mu[i*prime[j]]=-mu[i];
            phi[i*prime[j]]=phi[i]*(prime[j]-1);
        }
    }
    for(rg int i=2;i<MAXN;++i)
        mu[i]+=mu[i-1],phi[i]+=phi[i-1];
}

map<int,ll>ans_mu;

ll cal_mu(int n)
{
    if(n<MAXN)
        return mu[n];
    if(ans_mu.count(n))
        return ans_mu[n];
    ll ans=1; // first positive 1
    for(rg int i=2,j;i<=n;i=j+1) // i should be long long because n max=INF,then i can be -1
    {
        j=n/(n/i),ans-=(j-i+1)*cal_mu(n/i);
    }
    return ans_mu[n]=ans;
}

map<int,ll>ans_phi;

ll cal_phi(int n)
{
    if(n<MAXN)
        return phi[n];
    if(ans_phi.count(n))
        return ans_phi[n];
    ll ans=(ll)n*(n+1)/2;
    for(rg int i=2,j;i<=n;i=j+1)
    {
        j=n/(n/i),ans-=(j-i+1)*cal_phi(n/i);
    }
    return ans_phi[n]=ans;
}

int main()
{
    sieve();
    int T;
    read(T);
    while(T--)
    {
        int n;
        read(n);
        printf("%lld %lld\n",cal_phi(n),cal_mu(n));
    }
    return 0;
}

Hint

AC记录

上一篇:一些你可能需要的okhttp实现


下一篇:transform:rotate3d/tranlate3d