迪利克雷生成函数

迪利克雷生成函数

定义

对于数列 ({f}),定义 (f) 的 Dirichlet 生成函数 (( m DGF)) 为:

[f(z)=sum_{n=1}^{infty} frac{f_n}{n^z} ]

不难发现 ((fcdot g)(z)) 为:

[(fcdot g)(z)=sum_{n=1} frac{sum_{d|n}f_dg_{frac{n}{d}}}{n^z} ]

  • 考虑到 (a^z imes b^z=(ab)^z)

对于 Dirichlet 卷积意义下的单位元 (epsilon),其 DGF 为 (epsilon(z)=1)

同时,我们定义的 (I(n)=1),其对应的 PGF 即为 (zeta(z)=sum_{n=1}^{infty} frac{1}{n^z}),即 Zeta 函数


运算


( imes)

[(fcdot g)(z)=sum_{n=1} frac{sum_{d|n}f_dg_{frac{n}{d}}}{n^z} ]


( m Inv)

等价于迪利克雷卷积逆。

(f_1=1,f_{1}^{-1}=1)

(sum_{d|n}f_df^{-1}_{frac{n}{d}}=0 o f_n^{-1}=-sum_{d|n,d e 1} f_df^{-1}_{frac{n}{d}})

可以在 (mathcal O(nlog n)) 的时间复杂度解决。


(ln & exp)

例题:

Dirichlet k 次根,即对于 (g,k) 计算 (f^{k}=g)

要求一个 (log)(即倍增卷积无法通过)保证 (g_1=1)

(nle 10^6,kle 10^9),答案对 (10^9+7) 取模。


(g=ln f),那么 (g'=frac{f'}{f})

于是

[g=int f' imes f^{-1} ]

考虑求导:

[egin{aligned} &f(z)'=igg(sum_n frac{f_n}{n^z}igg)' \&=sum_n f_n igg(Big(frac{1}{n}Big)^zigg)' end{aligned}]

注意到 ((n^x)'=exp(xln n)=n^xln n)

于是有:

[egin{aligned} &f(z)'=sum_n frac{f_n}{n^z}ln (frac{1}{n}) \&f(z)'=-sum_n frac{f_n}{n^z}ln n end{aligned}]

于是对于级数的变换即 (f_n'=-f_nln n)

问题来了,(ln n) 在有理数域没有定义。

不过我们知道最后会求一次积分,所以找一个在求导意义下可以替换他的即可。

(c(n)) 表示 (n) 的质因子个数,那么 (c(nm)=c(n)+c(m))

这样就可以解决 (ln) 的问题了。

另一边,我们考虑积分的问题,显然积分即导数的逆运算,所以有 (f_n'=-frac{f_n}{c(n)})

接下来考虑 (exp) 了。

(g=exp(f) o g'=gf')

于是有 (g_n'=sum_{d|n} g_df'_{frac{n}{d}}),故 (-c(n)g_n=f_1c(1) imes g_n+sum_{d e n,d|n} g_df_{frac{n}{d}}')

(g_n=-frac{1}{c(n)}sum_{d e n,d|n} g_df_{frac{n}{d}}')

对于 (n=1)(g_1=1(f_1=1))

(Code:)

#include<bits/stdc++.h>
using namespace std ;
#define Next( i, x ) for( register int i = head[x]; i; i = e[i].next )
#define rep( i, s, t ) for( register int i = (s); i <= (t); ++ i )
#define drep( i, s, t ) for( register int i = (t); i >= (s); -- i )
#define re register
#define int long long
int gi() {
	char cc = getchar() ; int cn = 0, flus = 1 ;
	while( cc < '0' || cc > '9' ) {  if( cc == '-' ) flus = - flus ; cc = getchar() ; }
	while( cc >= '0' && cc <= '9' )  cn = cn * 10 + cc - '0', cc = getchar() ;
	return cn * flus ;
}
const int P = 998244353 ; 
const int N = 1e6 + 5 ; 
int n, K, f[N], g[N], ig[N], c[N], inv[N] ;
int top, isp[N], p[N / 10] ; 
int fpow(int x, int k) {
	int ans = 1, base = x ;
	while(k) {
		if(k & 1) ans = 1ll * ans * base % P ;
		base = 1ll * base * base % P, k >>= 1 ;
	} return ans ;
}
void init(int x) {
	rep( i, 2, x ) {
		if( !isp[i] ) p[++ top] = i, c[i] = 1 ; 
		for(re int j = 1; j <= top && p[j] * i <= x; ++ j) {
			isp[i * p[j]] = 1, c[i * p[j]] = c[i] + 1 ;
			if( i % p[j] == 0 ) break ;
		}
	} 
	rep( i, 1, x ) inv[i] = fpow(c[i], P - 2) ; 
}
int h[N] ; 
void fmul(int *a, int *b) {
	memset( h, 0, sizeof(h) ) ;
	rep( i, 1, n ) for(re int j = i; j <= n; j += i)
		h[j] = (h[j] + a[i] * b[j / i]) % P ; 
	rep( i, 1, n ) a[i] = h[i], h[i] = 0 ;
}
void fdev(int *a) {
	rep( i, 1, n ) a[i] = P - a[i] * c[i] % P ; 
}
void fint(int *a) {
	rep( i, 1, n ) a[i] = P - a[i] * inv[i] % P ; 
}
void fexp(int *a) {
	h[1] = 1 ; 
	rep( i, 1, n ) {
		if( i ^ 1 ) h[i] = P - inv[i] * h[i] % P ; 
		for(re int j = i * 2; j <= n; j += i)
			h[j] = (h[j] + h[i] * a[j / i]) % P ; 
	} rep( i, 1, n ) a[i] = h[i], h[i] = 0 ;
}
void fdiv(int *a) {
	h[1] = 1 ; 
	rep( i, 1, n ) {
		if( i ^ 1 ) h[i] = P - h[i] % P ; 
		for(re int j = i * 2; j <= n; j += i)
			h[j] = (h[j] + h[i] * a[j / i]) % P ; 
	} rep( i, 1, n ) a[i] = h[i], h[i] = 0 ;
}
signed main()
{
	n = gi(), K = gi(), init(n), K = fpow( K, P - 2 ) ; 
	rep( i, 1, n ) ig[i] = g[i] = gi() ; 
	fdev(ig), fdiv(g), fmul(g, ig), fint(g) ;
	rep( i, 1, n ) g[i] = g[i] * K % P ;
	fdev(g), fexp(g) ;
	rep( i, 1, n ) printf("%lld ", g[i] ) ; 
	return 0 ;
}
原文地址:https://www.cnblogs.com/Soulist/p/13741395.html