狄利克雷卷积

积性函数

定义:

对于数论函数$f$,若对于任意互质的数$x,y$,满足$f(xy)=f(x)f(y)$,则$f$为一个积性函数。

事实上,我们见过的大部分数论函数都是积性函数,常见的如:

  • $\mu(x)$,莫比乌斯函数,在莫比乌斯反演有讨论过。
  • $\varphi(x)$,欧拉函数。
  • $d(x)$,表示$x$的约数个数。
  • $\sigma(x)$,约数和函数。
  • $\epsilon(x)$,狄利克雷卷积的原函数,即$\epsilon(x)=[x=0]$。
  • $id(x)$,定义$id(x)=x$。
  • $I(x)$,定义$I(x)$函数值恒为$1$,即$I(x)=1$。

这只是一小部分,其中有些函数的性质待会会分析。

狄利克雷卷积

定义:

对于数论函数$f$和$g$,定义它们的狄利克雷卷积为:

其中$(f*g)$表示$f$和$g$的狄利克雷卷积。

然后有一个很重要的性质,对于积性函数$f$和$g$,它们的狄利克雷卷积也是个积性函数。

证明很简单,设质数$p$和$n$互质,则:

然后展开:

由于$f,g$为积性函数,可以拆开提出来:

注意到第一项就是$(fg)(1)$,第二项为$(fg)(n)$,所以:

命题得证。

根据这个,可以有效的判断函数是否为积性函数。

狄利克雷卷积具有的一些性质:

  • 交换律:$fg=gf$。
  • 结合律:$f(gh)=(fg)h$。

证明很显然,这里就不赘述了。

然后考虑一些常见函数的狄利克雷卷积:

1.$\epsilon(x)$。

对于任意数论函数$f$,根据定义,拿一个卷上一个$\epsilon$,得到:

即$f*\epsilon=f$,

然后我们可以惊奇的发现,卷完了之后得到了他本身,所以说这个函数就是狄利克雷卷积的单位元。

2.$\mu(x)$。

首先我们知道这样一个式子:

写成狄利克雷卷积形式就是:

3.$\varphi(x)$。

这个函数有一些很妙的性质,比如说(众所周知):

证明如下:

设$f(n)=\sum_{d|n}\varphi(d)$,由上面提到的狄利克雷卷积性质可得,$f$为积性函数。

对于$n$的每一个质因数进行考虑,即:

因为$\varphi(x^p)=\varphi(x^{p-1})*p$,$\varphi(x)=x-1$可得:

由等比数列求和可得$f(x^p)=x^p$,又由$f$为积性函数可得$f(n)=n$,得证。

写成狄利克雷卷积的形式就是:

然后注意到$\mu*I=\epsilon$,于是尝试在等式两边分别卷上一个$\mu$:

写成一般形式就是:

这个式子也很常用的。

杜教筛

前置知识终于讲完了

对于积性函数$f$,现在考虑求$\sum_{i=1}^nf(i)$。

当然我知道你会线筛,但是有一些小清新(无良)出题人把数据开到$1e9$级别,你就做不了了,所以需要更快的算法。

设$S(n)=\sum_{i=1}^{n}f(i)$。

考虑先拿一个函数和$f$做卷积,先不考虑这个函数是啥,先设它为$g$,则:

对这个函数前缀和一下:

中间交换求和符号就不赘述了。

然后有一个很显然的式子:

然后把等号右边第一项换一下:

剩下的事情就很显然了,如果我们能凑出一个$g$,使得它和$f$的卷积的前缀和很好算,它也很好算,我们就可以先预处理出前$1e7$左右的数据,然后记忆化递归处理$S$。

至于$g$,因题目不同而不同,做多了也就有经验了。

注意上试是杜教筛的套路试子,如果不会推也没关系,记下来就好了。

这里还是举几个栗子吧:

1.$\sum_{i=1}^{n}\mu(i)$.

首先我们知道一个这样的东西:

对于$\epsilon$的前缀和,显然是$1$,

所以,令$g(n)=I(n)=1$,得到:

然后数论分块下就行了。

代码也很简单,记住一定要记忆化。

1
2
3
4
5
6
7
8
9
10
11
map<int,int > Mu;
int sum_mu(int n) {
if(n<maxn) return mu[n];
if(Mu[n]) return Mu[n];
int T=2,res=1;
while(T<=n) {
int pre=T;T=n/(n/T);
res=res-(T-pre+1)*sum_mu(n/T);T++;
}
return Mu[n]=res;
}

2.$\sum_{i=1}^{n}\varphi(i)*i$。

这次考虑一个难一点的。

我们的目的是要找一个$g$,使他们卷起来很好算。

那么先列出式子:

然后发现中间有个$i$的系数,可以考虑消去它,于是暂且令$g(n)=n$:

然后可以发现这个东西非常好算,因为:

所以:

总结

对于一个陌生的函数,如果要求前缀和,先判是不是积性函数,然后通过这个函数的性质进行分析凑$g$,也可以像栗2一样凑,实在不行就一个一个试。

当然有时候线筛也是很好的,或者可以枚举因数大力算,不要总是纠结于能不能杜教筛。

习题

[bzoj3944] Sum

[bzoj4176] Lucas的数论

[bzoj4916] 神犇和蒟蒻

[bzoj4652]|[Noi2016]循环之美

[bzoj3930] [CQOI2015]选数

题目大意

定义整数$a,b$,求所有满足一下条件的二元组$(i,j)$的$\rm lcm(i,j)$的和。

  • $1\leqslant i\leqslant a,1\leqslant j\leqslant b$。
  • $\forall n>1,n^2\not \mid \gcd(i,j)$。
阅读全文 »

Min_25筛简介

$\text{min_25}$筛是一种处理一类积性函数前缀和的算法。

其中这类函数$f(x)$要满足$\sum_{i=1}^{n}[i\in prime]\cdot f(i)$可以被$\sum_{i=1}^{n}[i\in prime]\cdot i^k$简单表示或者快速计算,其中$k$为较小的常数。

时间复杂度好像是$O(\frac{n^{0.75}}{\log n})$,不过据说被证伪了…也有人说是$O(n^{1-\epsilon})$,反正可以做$n=1e10$,貌似跑的比杜教筛快。

空间复杂度$O(\sqrt{n})$。

阅读全文 »

引入

关于斯特林数的定义可以上百度搜…这部分很简单就没写了。

本文中,$s_1(n,k)$表示第一类斯特林数,$s_2(n,k)$表示第二类。

这里不加说明的给出递推式:

卷积求第二类斯特林数

首先有一个组合意义很明显的式子:

就是说$k$个位置填$n$种颜色,我们可以分成若干组,每组颜色不同。

我们可以把后面的值为$0$的项补全:

对其进行二项式反演:

这个可以写成卷积形式,顺便把变量换下:

阶乘幂

阶乘幂的定义是这样的:

上升幂:$x^{\overline n}=\prod_{i=1}^{n}(x+i-1)$。

下降幂:$x^{\underline n}=\prod_{i=1}^{n}(x-i+1)$。

显然可以根据组合数的定义得到:

还有两个这样的式子:

把右边拆开就能证明了。

阶乘幂与斯特林数的转化

所以上面卷积的第一个式子可以写成这样:

还有一个第一类斯特林数和上升幂的关系

可以通过数学归纳来证明:

注意由于$s_1(n,0)=s_1(n,n+1)=0$,所以中间边界条件其实不需要在意。

我们把这个式子变一下可以得到一个第一类斯特林数和下降幂的关系

最上面那个普通幂转第二类斯特林数和上升幂的式子也可以类似的变一下:

总结一下一共四个转化的式子:

反转公式

首先搬出下降幂的公式:

套上面的公式变成上升幂:

套上升幂转第一类斯特林数的公式:

把$x$的幂提前,换一下求和符号:

由于这里我们把$x$看成未知量,其他的都是已知量,所以我们可以把左右当作多项式,那么对比系数可得:

这个叫做反转公式。

同理我们可以类似的得出第二个反转公式:

注意:反转公式$-1$的指数也可以写成$n-i$,稍加分析可以发现$m=n$时成立,$m\ne n$时有两种情况,一种不变,另一种会将答案取相反数,但是由于结果为$0$所以不影响。

斯特林反演

这里先给出反演的式子在加以证明:

考虑一般反演的套路,先写出一个$[i=n]$的形式:

在把和斯特林数以及$[m=n]$的式子套进去,也就是上面的反转公式(注意$-1$的指数):

证毕。

当然由于反转公式的对称性,所以互换$s_1,s_2$依然成立。

这其实就是一个小结论….我做题的时候碰到的,研究了下就写了篇博客。

(其实主要就是把wiki上一篇论文翻译简化了下)

结论是这样的:设$n$个点$s$个连通块的有标号森林个数为$f(n,s)$,则:

阅读全文 »

做到了一个$\rm cf$题需要用到这个东西我就去学了下…

以前觉得这玩意好毒瘤,不过其实也不难。

类欧几里得算法一般是用来快速计算以下几个函数的(当然有其他的作用但是一般不会出现,以后看到了在写吧):

其复杂度为$O(\log n)$。

计算F

我们分情况讨论,先考虑$a\geqslant c$或者$b\geqslant c$的情况。

注意到根据整除的性质我们可以把式子写成这样:

注意到这里递归了,而且递归后一定是一个$a<c,b<c$的情况,那么此时用上面的展开就没用了。

注意到这个和式实质上是求一条直线下的整点个数,一开始我们是枚举$x$,然后算$y$,我们现在先枚举$y$,也就是横着的一条,可以得到:

推一下后面那个不等式,因为是大于等于号所以整除可以去掉:

那么式子就可以化简了:

观察下递归时$a,c$的变化:$(a,c)\to (a\bmod c,a)\to (a,a\bmod c)$。

这就是欧几里得算法的递归过程,所以最终$a$一定会递归到$0$,此时式子可以直接计算,递归次数为$O(\log n)$。

计算G

和上面思路相同,当$a\geqslant c$或$b\geqslant c$时:

否则:

其中$f$我们已经会算了,这里还用到了$h$,所以我们在推推$h$。

计算H

当$a\geqslant c$或$b\geqslant c$时,思路一样,但是式子有点复杂注意别搞错:

其中函数后面参数都是$(a\bmod c,b\bmod c,c,n)$,我省略没写了。

第二种情况这时候就有点难算了,因为平方无法处理,这里有一个小技巧可以把平方拆开:

那么就很简单了,一个整点$(x,y)$会被统计$y$的贡献,式子写出来就是:

代码

如果按照上面写三个函数互相调用显然是不优的,因为一个状态有可能被算很多次,当然可以采取记忆化,但是这样会多一个查找的$\log $。

有一个小技巧就是,我们可以写一个$\rm dfs$,返回三个值,分别是$f,g,h$,注意到这三个函数递归的时候参数都是一样的,所以这样只会递归$O(\log n)$次,并且每次递归都是$O(1)$算。

所以总复杂度$O(\log n)$。

代码长这样:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
#include<bits/stdc++.h>
using namespace std;

void read(int &x) {
x=0;int f=1;char ch=getchar();
for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}

void print(int x) {
if(x<0) putchar('-'),x=-x;
if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('\n');}

#define lf double
#define ll long long

#define pii pair<int,int >
#define vec vector<int >

#define pb push_back
#define mp make_pair
#define fr first
#define sc second

#define FOR(i,l,r) for(int i=l,i##_r=r;i<=i##_r;i++)

const int maxn = 1e6+10;
const int inf = 1e9;
const lf eps = 1e-8;
const int mod = 998244353;
const int inv2 = 499122177;
const int inv6 = 166374059;

struct data {int f,g,h;};

data calc(int a,int b,int c,int n) {
data ans;
int s1=1ll*n*(n+1)%mod*inv2%mod,s2=1ll*n*(n+1)%mod*(n*2+1)%mod*inv6%mod;
if(!a) {
ans.f=1ll*(n+1)*(b/c)%mod;
ans.g=1ll*s1*(b/c)%mod;
ans.h=1ll*(b/c)*(b/c)%mod*(n+1)%mod;
return ans;
}
if(a>=c||b>=c) {
data res=calc(a%c,b%c,c,n);
ans.f=(0ll+res.f+1ll*s1*(a/c)%mod+1ll*(n+1)*(b/c)%mod)%mod;
ans.g=(0ll+res.g+1ll*s2*(a/c)%mod+1ll*s1*(b/c)%mod)%mod;
ans.h=(0ll+res.h+1ll*(a/c)*res.g*2%mod+1ll*(b/c)*res.f*2%mod)%mod;
ans.h=(0ll+ans.h+1ll*s2*(a/c)%mod*(a/c)%mod+1ll*n*(n+1)%mod*(b/c)%mod*(a/c)%mod)%mod;
ans.h=(ans.h+1ll*(n+1)*(b/c)%mod*(b/c)%mod)%mod;
return ans;
}int m=(1ll*a*n+b)/c;
data res=calc(c,c-b-1,a,m-1);
ans.f=(1ll*m*n%mod-res.f+mod)%mod;
ans.g=(1ll*m*n%mod*(n+1)%mod-res.h-res.f)%mod;ans.g=(1ll*ans.g*inv2%mod+mod)%mod;
ans.h=(0ll+1ll*m*n%mod*(m+1)%mod-2ll*res.g-2ll*res.f-ans.f)%mod;ans.h=(ans.h+mod)%mod;
return ans;
}

int main() {
int t;read(t);
while(t--) {
int a,b,c,n;read(n),read(a),read(b),read(c);
data ans=calc(a,b,c,n);
printf("%d %d %d\n",ans.f,ans.h,ans.g);
}
return 0;
}

看起来很复杂其实就是对着式子抄了一遍,提交至洛谷模板可$\rm AC$。