杜教筛
杜教筛
前置技能树：积性函数
就是对于函数\(f(x)\)
对于任意两个互质整数\(a,b\),如果有\(f(a)·f(b)=f(ab)\)
则\(f\)为积性函数
如果对任意\(a,b\)成立，\(f\)为完全积性函数。
前置技能树：狄利克雷卷积
狄利克雷卷积是一种运算定义。
\(f*g=\sum\limits_{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)=\sum\limits_{i=1}^{n}f(i)\)。
我们构造两个积性函数\(h,g\)，使得\(h=g*f\)。
\(\sum\limits_{i=1}^n h(i)=\sum\limits_{i=1}^n(g*f)(i)\)
\(=\sum\limits_{i=1}^n \sum\limits_{d|i}g(d)*f(\frac nd)\)
考虑套路，把\(d\)提出来，枚举倍数。
\(=\sum\limits_{d=1}^n g(d)\sum\limits_{i=1}^{\lfloor\frac nd\rfloor}f(i)\)
\(=\sum\limits^{n}_{d=1}g\left( d\right) S\left( \lfloor \dfrac {n}{d}\rfloor \right)\)
我们发现这个里面已经有我们要求的\(S(n)\)了，于是单独拿出来。
\(=g(1)S(n)+\sum\limits_{d=2}^n S\left(\lfloor \dfrac {n}{d}\rfloor \right)\)
我们移一下项，\(\Large g(1)S(n)=\sum\limits_{i=1}^nh(i)-\sum\limits_{d=2}^nS(\lfloor\frac nd \rfloor)\)
这就是杜教筛的套路柿子，于是大功告成，我们就沿着这个柿子递归就可以了
...吗？
还没完...
我们发现如果直接递归，复杂度是\(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-\sum\limits_{d=2}^nS(\lfloor\frac nd \rfloor)\)
例题\(2\)：\(f=\varphi\)
仍然构造\(g=I\)
\(\varphi * I=id\)
所以\(S(n)=\frac{n(n+1)}2-\sum\limits_{d=2}^nS(\lfloor\frac nd \rfloor)\)
卡常小技巧
我们平时写杜教筛的时候，筛出的结果显然要存下来。
但是这样要手写哈希(麻烦)或者使用\(STL\)的\(unorded\_map\)与\(\_\_gnu\_pbds::cc\_hash\_table\)(慢)
有没有两全其美的办法呢？
有！
我们发现求解的总是n的约数。
于是我们开一个两倍的\(dp\)数组。
当\(x\le \sqrt n\)时
返回\(dp[x]\)
否则返回\(dp[\sqrt n+n/x]\)就可以了。
这样就可以跑得飞快，不用太费心卡常。
/*
@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\n",getphi(n),getmu(n));
}
return 0;
}
浙公网安备 33010602011771号