数学杂谈 #4

数论函数的级数

在此主要介绍狄利克雷级数和贝尔级数。

狄利克雷级数

狄利克雷级数是定义在任意数论函数上的一种级数,对于数论函数 \(f\) ,我们定义它的狄利克雷级数为:

\[F(z)=\sum_{k\ge1 }f(k)k^{-z}=\sum_{k\ge 1}\frac{f(k)}{k^{z}} \]

这种级数的重要意义在于,两个数论函数的狄利克雷级数的卷积恰好为函数狄利克雷卷积后的级数

也即:

\[(f*g)(n)=h(n)\Leftrightarrow F(z)G(z)=H(z) \]

有了这个工具我们就可以很方便地研究数论函数的狄利克雷卷积卷积相关的问题;

关于这个级数为什么定义成 \(k^{-z}\) 而非 \(k^z\) ,这个牵扯到黎曼 \(\zeta\) 函数,因为它长这个样子:

\[\zeta(z)=\sum_{k\ge 1}k^{-z}=\sum_{k\ge 1}\frac{1}{k^z} \]

可以发现 \(\zeta\) 就是 \(I\) 的狄利克雷级数。

也许定义成类似于 \(\zeta\) 的形式会有许多优美的性质。

贝尔级数

对于积性函数 \(f\) ,由于 \(f\) 的值实际上只和它在质数的幂处的取值相关,因而我们只需要研究这些地方 \(f\) 的性质。贝尔级数应运而生。

积性函数 \(f\) 对于质数 \(p\) 的贝尔级数定义为:

\[F_p(z)=\sum_{k\ge 0}f(p^k)z^k \]

类似于狄利克雷级数,贝尔级数在同一质数上时,两个积性函数的贝尔级数的卷积就是狄利克雷卷积后的贝尔级数

更重要的是,由于 \(F_p(z)\) 其实就是常规多项式的形式(通常而言常见数论函数在确定质数的幂的取值也仅与指数相关),因此我们可以使用研究普通生成函数的方法来研究它。

此外,贝尔级数的另一个优点是,由于指标变为了指数,因此在 OI 常见的数据范围中, \(k\) 的范围通常很小,这允许我们对贝尔级数进行常规的多项式计算而(通常)不会带来较大的复杂度。


对于常见的数论函数,可以得到它们的贝尔级数为:

  • \(\varepsilon(n)\Rightarrow \Epsilon_p(z)=1\)
  • \(I(n)\Rightarrow I_p(z)=\sum_{k\ge 0}z^k=\frac{1}{1-z}\)
  • \(id(n)\Rightarrow ID_p(z)=\sum_{k\ge 0}p^kz^k=\frac{1}{1-pz}\)
  • \(id^k(n)\Rightarrow (ID^k)_p(z)=\sum_{j\ge 0}(p^{k})^{j}z^j=\frac{1}{1-p^kz}\)
  • \(\varphi(n)\Rightarrow \Phi_p(n)=\sum_{k\ge 0}(p^k-p^{k-1})z^k=\frac{1-z}{1-pz}\)
  • \(\mu(n)\Rightarrow \Mu_p(n)=\sum_{k\ge 0}\mu(p^{k})z^k=1-z\)
  • \(\mu^2(n)\Rightarrow (\Mu^2)_p(n)=1+z\)

这里需要说明的是,数论函数的指数指的是自己和自己做点积的结果

从而一些重要的性质可以很容易地观察出来,比如:

  • \(\mu*I=\varepsilon\Leftrightarrow \forall p\in P,\Mu_p(z)\cdot I_p(z)=\Epsilon_p(z)\)
  • \(\varphi*I=id\Leftrightarrow \forall p\in P,\Phi_p(z)\cdot I_p(z)=ID_p(z)\)

其中 \(P\) 是质数的集合。

Powerful Numbers💣

对于积性函数的前缀和,我们已经有了许多具有优秀复杂度的筛法,例如杜教筛、 min_25 筛。在此,我们将介绍一种理论—— Powerful Numbers ,它允许我们快速计算一大类特殊的积性函数的前缀和。

Powerful Numbers 的简要内容

假如我们得到了一个积性函数 \(f(n)\) ,我们希望能快速求出以下式子:

\[\sum_{k=1}^n f(k) \]

在该理论中,我们可以寻找某一特殊的积性函数 \(g(n)\) ,使得在质数 \(p\) 处均有 \(f(p)=g(p)\) 。那么我们不难构造出另一个函数 \(h\) ,使得 \(f=g*h\) ,而 \(h=f/g\)

现在我们可以再考虑一下目标式:

\[\begin{aligned} \sum_{k=1}^nf(k) &=\sum_{k=1}^n\sum_{d|k}g(d)h(\frac{k}{d})\\ &=\sum_{d=1}^n\sum_{j=1}^{\lfloor\frac n d\rfloor}h(d)g(j)\\ &=\sum_{d=1}^nh(d)\left(\sum_{j=1}^{\lfloor\frac n d\rfloor} g(j)\right) \end{aligned} \]

注意到 \(h\) 的特殊性:对于 \(p\in P\) ,可以发现:

\[h(p)=\frac{f(p)-g(p)h(1)}{g(1)}=0 \]

而因为 \(f,g\) 都是积性函数,因而 \(h\) 必然也是积性函数。因此,对于任意的 \(n\) ,如果其有某一个质因子的质数为 1 ,那么 \(h(n)=0\) 。所以,真正有用的 \(h\) 其实并不多(我们将稍后分析这样的数的数量)。

考虑对于每一个质数 \(p\) ,计算 \(h(p^k)\) 的值。由于 \(f,g\) 都是已知的,所以我们可以得到 \(F_p(z),G_p(z)\) ,因而可以暴力地算出 \(H_p(z)\) (当然也可以尝试直接算出 \(H_p(z)\) 的封闭形式)。

因此,只要我们可以快速计算 \(g\) 在所有 \(\lfloor\frac n d\rfloor\) 处的值,我们就可以直接搜索有用的 \(h\) ,进而快速算出 \(f\) 的前缀和了。

Powerful Numbers 的时间复杂度

  1. 求出 \(g\) 的前缀和,这一部分我们可以考虑为 \(O(n^{\frac{2}{3}})\) (可能为 \(O(\frac{n^{\frac{3}{4}}}{\log n})\))。

  2. 由于当 \(p>\sqrt n\) 的时候有用的指数最高为 1 ,因此我们只需要考虑 \(\le \sqrt n\) 的指数。暴力做除法的复杂度是 \(O(\log_p^2n)\) ,而有用的质数数量大约为 \(O(\frac{\sqrt n}{\log n})\) ,因而可以估计为 \(O(\sqrt n\log n)\)

  3. 我们将令 \(h(n)\) 有值的 \(n\) 称为 Powerful Numbers 。不难发现,所有的 Powerful Numbers 一定可以被表示为 \(a^2b^3(a,b\in \mathbb N_+)\) 的形式。那么在 \(n\) 范围内的 Powerful Numbers 的数量为:

    \[\sum_{k=1}^{\lfloor\sqrt n\rfloor}\left\lfloor\sqrt[3]{\frac{n}{k^2}}\right\rfloor \]

    使用积分近似得到这个算式结果是 \(O(\sqrt n)\) 的:

    \[\int_0^{n^{\frac 1 2}} n^{\frac 1 3}\cdot x^{-\frac{2}{3}} \mathop{}\!\mathrm dx=n^{\frac 1 3}\cdot \left.(3x^{\frac 1 3})\right|^{n^{\frac 1 2}}_0 \]

    同时也不难知道它的下界的确是 \(\lfloor\sqrt n\rfloor\) ,因此单纯搜索的复杂度为 \(O(\sqrt n)\)

总的来说, Powerful Numbers 的主要复杂度其实就是求 \(g\) 的复杂度。一般我们可以将它认为是 \(O(n^{\frac 2 3})\) 或者 \(O(\frac{n^{\frac 3 4}}{\log n})\) 的,取决于你的筛法实现。

一道例题

题目链接: 【模板】Min_25筛

究竟有多少份 AC 提交是真的 min_25 ?

这里有 \(f(p)=p(p-1)\) ,一个符合要求的 \(g\)\(g(n)=(id\cdot \varphi)(n)\) ,在质数处的值与 \(f\) 相同。

求解 \(g\) 的前缀和可以利用性质:

\[g*id=\sum_{d|n}d\varphi (d)\cdot \frac{n}{d}=n^2=id^2 \]

其实到这里差不多就可以直接算了,不过为了取得更高效的算法,我们可以进行几步推导。

对于质数 \(p\) ,考虑 \(f\)\(g\) 的贝尔级数:

\[\begin{aligned} F_p(z)&=\sum_{k\ge 0}(p^{2k}-p^k)z^k=\frac{1}{1-p^2z}-\frac{1}{1-pz}+1\\ G_p(z)&=\frac{(ID^2)_k(z)}{ID_k(z)}=\frac{1-pz}{1-p^2z}\\ \end{aligned} \]

一定要注意级数的常数项。如果这里不补 1 就会得到错误的结果!

进而可以直接得到 \(H_p(z)\)

\[\begin{aligned} H_p(z) &=\frac{F_p(z)}{G_p(z)}\\ &=\frac{\frac 1{1-p^2z}-\frac{1}{1-pz}+1}{\frac{1-pz}{1-p^2z}}\\ &=\frac{1}{1-pz}-\frac{1-p^2z}{(1-pz)^2}+\frac{1-p^2z}{1-pz} \end{aligned} \]

可以得到 \(h(p^k)=[z^k]H_p(z)=(k-1)(p^{k+1}-p^k)\)

所以就可以直接计算 \(h\) 的值了。亲测会比 min_25 筛本身跑得快。

以下是代码:

#include <cmath>
#include <cstdio>

#define rep( i, a, b ) for( int i = (a) ; i <= (b) ; i ++ )
#define per( i, a, b ) for( int i = (a) ; i >= (b) ; i -- )

typedef long long LL;

const int mod = 1e9 + 7, inv6 = 166666668, inv2 = 500000004;
const int MAXLOG = 55, MAXN = 5e6 + 5;

template<typename _T>
void read( _T &x )
{
	x = 0; char s = getchar(); int f = 1;
	while( s < '0' || '9' < s ) { f = 1; if( s == '-' ) f = -1; s = getchar(); }
	while( '0' <= s && s <= '9' ) { x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar(); }
	x *= f;
}

template<typename _T>
void write( _T x )
{
	if( x < 0 ) putchar( '-' ), x = - x;
	if( 9 < x ) write( x / 10 );
	putchar( x % 10 + '0' );
}

int val1[MAXN], val2[MAXN];

LL pw2[MAXN];
int phi[MAXN];
int prime[MAXN], pn;
bool isPrime[MAXN];

int ans = 0;
LL N, rt;

inline int Mul( int x, int v ) { return 1ll * x * v % mod; }
inline int Sub( int x, int v ) { return ( x -= v ) < 0 ? x + mod : x; }
inline int Add( int x, int v ) { return ( x += v ) >= mod ? x - mod : x; }

void EulerSieve( const int n )
{
	rep( i, 2, n )
	{
		if( ! isPrime[i] ) phi[prime[++ pn] = i] = i - 1;
		for( int j = 1 ; j <= pn && 1ll * i * prime[j] <= n ; j ++ )
		{
			isPrime[i * prime[j]] = true;
			if( ! ( i % prime[j] ) )
			{
				phi[i * prime[j]] = phi[i] * prime[j];
				break;
			}
			phi[i * prime[j]] = phi[i] * ( prime[j] - 1 );
		}
	}
	val1[1] = 1;
	rep( i, 2, n ) 
		val1[i] = Add( val1[i - 1], Mul( phi[i], i ) );
}

void Init()
{
	rt = pow( N, 0.667 );
	EulerSieve( rt );
	rep( k, 1, pn )
		pw2[k] = 1ll * prime[k] * prime[k];
}

int GetG( const LL n )
{
	if( n <= rt ) return val1[n];
	if( val2[N / n] ) return val2[N / n];
	int ret = n % mod * ( ( n + 1 ) % mod ) % mod * ( ( n * 2 + 1 ) % mod ) % mod * inv6 % mod;
	for( LL l = 2, r ; l <= n ; l = r + 1 )
	{
		r = n / ( n / l );
		ret = Sub( ret, Mul( ( l + r ) % mod * ( r - l + 1 ) % mod * inv2 % mod, GetG( n / l ) ) );
	}
	return val2[N / n] = ret;
}

void S( const int idx, const int h, const LL n )
{
	LL lim = N / n;
	if( idx > pn || pw2[idx] > lim )
	{
		if( lim <= rt ) ans = Add( ans, Mul( h, val1[lim] ) );
		else ans = Add( ans, Mul( h, val2[N / lim] ) );
		return ;
	}
	S( idx + 1, h, n );
	LL p = pw2[idx];
	for( int e = 2 ; p * n <= N ; p *= prime[idx], e ++ )
		S( idx + 1, p % mod * ( prime[idx] - 1 ) % mod * ( e - 1 ) % mod * h % mod, n * p );
}

signed main()
{
	read( N ), Init(); 
	GetG( N ); 
	S( 1, 1, 1 );
	write( ans ), putchar( '\n' );
	return 0;
}
posted @ 2021-04-13 22:18  crashed  阅读(245)  评论(0编辑  收藏  举报