• 杜教筛


    杜教筛

    前置技能树:积性函数

    就是对于函数(f(x))

    对于任意两个互质整数(a,b),如果有(f(a)·f(b)=f(ab))

    (f)为积性函数

    如果对任意(a,b)成立,(f)为完全积性函数。

    前置技能树:狄利克雷卷积

    狄利克雷卷积是一种运算定义。

    (f*g=sumlimits_{d|n}f(d)g(frac nd))

    其显然满足交换律。

    前言

    一般我们求积性函数有优秀的(O(n))欧拉筛

    但是实际运用中,我们往往需要得出积性函数的前缀和来进行运算,而且总是有一些毒瘤出题人把数据出到(1e10)之类,此时线性筛就不够用了。

    为了解决这个问题,我们就需要一种新的筛法—杜教筛。

    杜教筛是一种筛法,能够以(O(n^{frac23}))的时间复杂度求积性函数的前缀和

    此外好像还有(min\_25)筛,复杂度为(O(n^{frac{3/4}{log_n}}))

    但是学不动了(Orz)

    具体推导

    这个东西都是套路。。。

    如果不想看公式了其实翻到下面加粗加大地方背个板子也挺不错

    (f(n))为你要筛的函数,(S(n)=sumlimits_{i=1}^{n}f(i))

    我们构造两个积性函数(h,g),使得(h=g*f)

    (sumlimits_{i=1}^n h(i)=sumlimits_{i=1}^n(g*f)(i))

    (=sumlimits_{i=1}^n sumlimits_{d|i}g(d)*f(frac nd))

    考虑套路,把(d)提出来,枚举倍数。

    (=sumlimits_{d=1}^n g(d)sumlimits_{i=1}^{lfloorfrac nd floor}f(i))

    (=sumlimits^{n}_{d=1}gleft( d ight) Sleft( lfloor dfrac {n}{d} floor ight))

    我们发现这个里面已经有我们要求的(S(n))了,于是单独拿出来。

    (=g(1)S(n)+sumlimits_{d=2}^n Sleft(lfloor dfrac {n}{d} floor ight))

    我们移一下项,(Large g(1)S(n)=sumlimits_{i=1}^nh(i)-sumlimits_{d=2}^nS(lfloorfrac nd floor))

    这就是杜教筛的套路柿子,于是大功告成,我们就沿着这个柿子递归就可以了

    ...吗?

    还没完...

    我们发现如果直接递归,复杂度是(O(n^{frac 34}))的,非常不优秀。

    于是我们就先预处理出前(n^{frac 23})的数字,这样递归的复杂度也是(O(n^{frac 23}))的了,(复杂度我也不会证)

    于是我们就成功的在(O(n^{frac 23}))的复杂度内求出了(f)的前缀和。

    n=1e10~11都可以做

    使用方法

    我们发现性能瓶颈主要在(h)函数的前缀和上,所以我们需要构造一个(g)使得(h)的前缀和能在(O(1))时间复杂度上解决。

    例题(1)(f=mu)

    构造(g=I)

    (mu * I = epsilon)

    (epsilon)的前缀和显然很菜,就是1

    所以(S(n)=1-sumlimits_{d=2}^nS(lfloorfrac nd floor))

    例题(2)(f=varphi)

    仍然构造(g=I)

    (varphi * I=id)

    所以(S(n)=frac{n(n+1)}2-sumlimits_{d=2}^nS(lfloorfrac nd floor))

    卡常小技巧

    我们平时写杜教筛的时候,筛出的结果显然要存下来。

    但是这样要手写哈希(麻烦)或者使用(STL)(unorded\_map)(\_\_gnu\_pbds::cc\_hash\_table)(慢)

    有没有两全其美的办法呢?

    有!

    我们发现求解的总是n的约数。

    于是我们开一个两倍的(dp)数组。

    (xle sqrt n)

    返回(dp[x])

    否则返回(dp[sqrt n+n/x])就可以了。

    这样就可以跑得飞快,不用太费心卡常。

    洛谷P4213 【模板】杜教筛

    /*
    @Date    : 2019-08-02 19:52:32
    @Author  : Adscn (adscn@qq.com)
    @Link    : https://www.cnblogs.com/LLCSBlog
    */
    #include<bits/stdc++.h>
    using namespace std;
    #define IL inline
    #define RG register
    #define gi getint()
    #define gc getchar()
    #define File(a) freopen(a".in","r",stdin);freopen(a".out","w",stdout)
    IL int getint()
    {
    	RG int xi=0;
    	RG char ch=gc;
    	bool f=0;
    	while(ch<'0'||ch>'9')ch=='-'?f=1:f,ch=gc;
    	while(ch>='0'&&ch<='9')xi=(xi<<1)+(xi<<3)+ch-48,ch=gc;
    	return f?-xi:xi;
    }
    template<typename T>
    IL void pi(T k,char ch=0)
    {
    	if(k<0)k=-k,putchar('-');
    	if(k>=10)pi(k/10,0);
    	putchar(k%10+'0');
    	if(ch)putchar(ch);
    }
    const int MAXN=5000000;
    const int N=INT_MAX;
    const int SQRN=sqrt(N)+7;
    typedef long long ll;
    int n,sqrn;
    bool npr[MAXN+7];
    int pr[MAXN+7],cnt;
    ll mu[MAXN+7];
    unsigned long long phi[MAXN+7];
    ll dpp[SQRN<<1],dpm[SQRN<<1];
    //map<int,ll>dpp,dpm;
    void init()
    {
    	npr[1]=1;
    	mu[1]=phi[1]=1;
    	for(int i=2;i<=MAXN;++i)
    	{
    		if(!npr[i])pr[++cnt]=i,mu[i]=-1,phi[i]=i-1;
    		for(int j=1;j<=cnt&&1ll*pr[j]*i<=MAXN;++j)
    		{
    			npr[i*pr[j]]=1;
    			if(i%pr[j]==0){mu[i*pr[j]]=0,phi[i*pr[j]]=phi[i]*pr[j];break;}
    			mu[i*pr[j]]=-mu[i],phi[i*pr[j]]=phi[i]*(pr[j]-1);
    		}
    	}
    	for(int i=1;i<=MAXN;++i)mu[i]+=mu[i-1],phi[i]+=phi[i-1];
    }
    inline ll getphi(int x)
    {
    	if(x<=MAXN)return phi[x];
    	if(x<=SQRN&&dpp[x])return dpp[x];
    	else if(x>SQRN&&dpp[N/x+SQRN])return dpp[N/x+SQRN];
    	unsigned long long ans=1ull*x*(x+1ull)/2ull;
    	for(int l=2,r;r<INT_MAX&&l<=x;l=r+1)r=x/(x/l),ans-=1ull*(r-l+1)*getphi(x/l);
    	if(x<=SQRN)dpp[x]=ans;
    	else dpp[N/x+SQRN]=ans;
    	return ans;
    }
    inline ll getmu(int x)
    {
    	if(x<=MAXN)return mu[x];
    	if(x<=SQRN&&dpm[x])return dpm[x];
    	else if(x>SQRN&&dpm[N/x+SQRN])return dpm[N/x+SQRN];
    	ll ans=1;
    	for(int l=2,r;r<INT_MAX&&l<=x;l=r+1)r=x/(x/l),ans-=1ll*(r-l+1)*getmu(x/l);
    	if(x<=SQRN)dpm[x]=ans;
    	else dpm[N/x+SQRN]=ans;
    	return ans;
    }
    int main(void)
    {
    	init();
    	int T=gi;
    	while(T--)
    	{
    		n=gi;
    		memset(dpp,0,sizeof dpp);
    		memset(dpm,0,sizeof dpm);
    		printf("%lld %lld
    ",getphi(n),getmu(n));
    	}
    	return 0;
    }
    
  • 相关阅读:
    html5shiv.js-让IE浏览器支持HTML5标准
    CSS2系列:外边距合并问题(margincollapse)
    HTML5:离线存储(缓存机制)-IndexDB
    CSS3系列:流式(弹性)布局(flex布局)
    Sublime Text 3 常用插件以及安装方法(转)
    后台配置参数写在文件上
    20160414
    2016413
    20160412
    网页设计素材网站
  • 原文地址:https://www.cnblogs.com/LLCSBlog/p/11305973.html
Copyright © 2020-2023  润新知