杜教筛
前置技能树:积性函数
就是对于函数(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])就可以了。
这样就可以跑得飞快,不用太费心卡常。
/*
@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;
}