筛积性函数的前缀和
- 常见积性函数
- 公式推导
- 狄利克雷卷积
- 杜教筛
- 实现
- 常见卷积
- 习题集
- Sum
- 神犇和蒟蒻
- 简单的数学题
常见积性函数
- μ\muμ
- φφφ
- d(n)d(n)d(n):nnn的约数个数
- σ(n)σ(n)σ(n):nnn的约数和
- ϵ(n)ϵ(n)ϵ(n):单位元函数,e(n)=[n=1]e(n)=[n=1]e(n)=[n=1] (完全积性)
- I(n)I(n)I(n):I(n)=1I(n)=1I(n)=1(完全积性)
- id(n)id(n)id(n):id(n)=nid(n)=nid(n)=n (完全积性)
∑d∣nφ(d)=n\sum_{d|n}φ(d)=nd∣n∑φ(d)=n∑d∣nμ(d)=[n=1]\sum_{d|n}\mu(d)=[n=1]d∣n∑μ(d)=[n=1]∑d∣nφ(d)×d=n×φ(n)+[n=1]2\sum_{d|n}φ(d)\times d=\frac{n\times φ(n)+[n=1]}{2}d∣n∑φ(d)×d=2n×φ(n)+[n=1]
公式推导
求∑i=1nf(i)\sum_{i=1}^nf(i)i=1∑nf(i)
狄利克雷卷积
定义两个数论函数f,gf,gf,g的狄利克雷卷积为hhh
(f∗g)(n)=h(n)⇔h(n)=∑d∣nf(d)×g(nd)=∑d∣ng(d)×f(nd)(f*g)(n)=h(n) \Leftrightarrow\ h(n)=\sum_{d|n}f(d)\times g(\frac{n}{d})=\sum_{d|n}g(d)\times f(\frac{n}{d})(f∗g)(n)=h(n)⇔ h(n)=d∣n∑f(d)×g(dn)=d∣n∑g(d)×f(dn)
如果 f,gf,gf,g 是积性函数则 hhh 函数也是一个积性函数。
证明:令 n=x∗y∧(x,y)=1n=x*y\wedge(x,y)=1n=x∗y∧(x,y)=1。
h(n)=∑d∣nf(d)g(nd)⇔h(xy)=∑d1∣x,d2∣yf(d1d2)g(xd1yd2)=∑d1∣x,d2∣yf(d1)f(d2)g(xd1)g(yd2)h(n)=\sum_{d|n}f(d)g(\frac n d)\Leftrightarrow h(xy)=\sum_{d_1|x,d_2|y}f(d_1d_2)g(\frac x {d_1}\frac y {d_2})=\sum_{d_1|x,d_2|y}f(d_1)f(d_2)g(\frac x{d_1})g(\frac y{d_2})h(n)=d∣n∑f(d)g(dn)⇔h(xy)=d1∣x,d2∣y∑f(d1d2)g(d1xd2y)=d1∣x,d2∣y∑f(d1)f(d2)g(d1x)g(d2y)=∑d1∣xf(d1)g(xd1)∑d2∣yf(d2)g(yd2)=h(x)h(y)=\sum_{d_1|x}f(d_1)g(\frac x{d_1})\sum_{d_2|y}f(d_2)g(\frac{y}{d_2})=h(x)h(y)=d1∣x∑f(d1)g(d1x)d2∣y∑f(d2)g(d2y)=h(x)h(y)
卷积满足乘法交换/结合律/分配律 (这不是重点😓)
f∗g=g∗ff*g=g*ff∗g=g∗f
(f∗g)∗h=f∗(g∗h)(f*g)*h=f*(g*h)(f∗g)∗h=f∗(g∗h)
(f+g)∗h=f∗h+g∗h(f+g)*h=f*h+g*h(f+g)∗h=f∗h+g∗h
杜教筛
杜教筛是一种以低于线性筛时间复杂度筛积性函数前缀和的筛法。
将所求问题定义成函数,令S(n)=∑i=1nf(i)S(n)=\sum_{i=1}^nf(i)S(n)=∑i=1nf(i)
先推∑i=1nh(i)\sum_{i=1}^nh(i)∑i=1nh(i)
∑i=1nh(i)\sum_{i=1}^nh(i)i=1∑nh(i)=∑i=1n∑d∣ig(d)×f(id)=\sum_{i=1}^n\sum_{d|i}g(d)\times f(\frac{i}{d})=i=1∑nd∣i∑g(d)×f(di)=∑d=1ng(d)∑i=1⌊nd⌋f(i)=\sum_{d=1}^ng(d) \sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(i)=d=1∑ng(d)i=1∑⌊dn⌋f(i)=∑d=1ng(d)S(⌊nd⌋)=\sum_{d=1}^ng(d)S(\lfloor\frac{n}{d}\rfloor)=d=1∑ng(d)S(⌊dn⌋)
将d=1d=1d=1的一项单独提出来⇓\Downarrow⇓
=g(1)×S(n)+∑d=2ng(d)S(⌊nd⌋)=g(1)\times S(n)+\sum_{d=2}^ng(d)S(\lfloor\frac{n}{d}\rfloor)=g(1)×S(n)+d=2∑ng(d)S(⌊dn⌋)g(1)×S(n)=∑i=1nh(i)−∑d=2ng(d)S(⌊nd⌋)g(1)\times S(n)=\sum_{i=1}^nh(i)-\sum_{d=2}^ng(d)S(\lfloor\frac{n}{d}\rfloor)g(1)×S(n)=i=1∑nh(i)−d=2∑ng(d)S(⌊dn⌋)
一般套路情况下g(1)=1g(1)=1g(1)=1
时间复杂度是 O(n23)O(n^\frac{2}{3})O(n32),大概就是 SSS 函数可以记忆化递归,另外设定一个阀值 MMM 预处理 S(1∼M)S(1\sim M)S(1∼M) 内的信息。
当 MMM 大约取 n23n^\frac{2}{3}n32 最优,让预处理和后面记忆化的时间复杂度差不多。
实现
因为g,hg,hg,h都是我们自己选择构造出来的,通过上面化出的最后式子来看
我们希望h(i)h(i)h(i)的前缀和是好求的,换言之就是很快就能求出来的
然后按S(⌊nd⌋)S(\lfloor\frac{n}{d}\rfloor)S(⌊dn⌋)可以对后面式子进行分块处理
至于怎么选择ggg和hhh,说实话只能观察式子,或者结合以前的做题经验套路
常见卷积
ⅠS(n)=∑i=1nμ(i)S(n)=\sum_{i=1}^n\mu(i)S(n)=∑i=1nμ(i)
因为知道在莫比乌斯反演有个关于莫比乌斯函数的公式
∑d∣nμ(d)=[n=1]\sum_{d|n}\mu(d)=[n=1]d∣n∑μ(d)=[n=1]
很容易想到ϵ=μ∗Iϵ=\mu*Iϵ=μ∗I
ϵ(n)=∑d∣nμ(d)×I(nd)=∑d∣nμ(d)=[n=1]ϵ(n)=\sum_{d|n}\mu(d)\times I(\frac{n}{d})=\sum_{d|n}\mu(d)=[n=1]ϵ(n)=d∣n∑μ(d)×I(dn)=d∣n∑μ(d)=[n=1]
Ⅱ s(n)=∑i=1nφ(i)s(n)=\sum_{i=1}^nφ(i)s(n)=∑i=1nφ(i)
上述提到过一个公式
n=∑d∣nφ(d)n=\sum_{d|n}φ(d)n=d∣n∑φ(d)
所以我们想到id=φ∗Iid=φ*Iid=φ∗I
id(n)=∑d∣nφ(d)×I(nd)=∑d∣nφ(d)=nid(n)=\sum_{d|n}φ(d)\times I(\frac{n}{d})=\sum_{d|n}φ(d)=nid(n)=d∣n∑φ(d)×I(dn)=d∣n∑φ(d)=n
Ⅲ S(n)=∑i=1nφ(i)×iS(n)=\sum_{i=1}^nφ(i)\times iS(n)=∑i=1nφ(i)×i
乍一看似乎配不出什么鸟玩意儿
尝试从狄利克雷卷积入手
h=f∗g=∑d∣nf(d)×g(nd)=∑d∣nφ(d)×d×g(nd)h=f*g=\sum_{d|n}f(d)\times g(\frac{n}{d})=\sum_{d|n}φ(d)\times d\times g(\frac{n}{d})h=f∗g=d∣n∑f(d)×g(dn)=d∣n∑φ(d)×d×g(dn)
我们想着如果能把fff顺带的多余的ddd除掉,尝试给ggg配成ididid
h=∑d∣nφ(d)×d×nd=∑d∣n{φ(d)×n}=n∑d∣nφ(d)=n2h=\sum_{d|n}φ(d)\times d\times \frac{n}{d}=\sum_{d|n}\{φ(d)\times n\}=n\sum_{d|n}φ(d)=n^2h=d∣n∑φ(d)×d×dn=d∣n∑{φ(d)×n}=nd∣n∑φ(d)=n2
i2i^2i2的前缀和,欸好巧哦,由公式可以快速算出,所以成功配出来了g,hg,hg,h
∑i=1ni2=n×(n+1)×(2n+1)6\sum_{i=1}^ni^2=\frac{n\times (n+1)\times (2n+1)}{6}i=1∑ni2=6n×(n+1)×(2n+1)
∑i=1ni3=(∑i=1ni)2=[n×(n+1)2]2=n2(n+1)24\sum_{i=1}^ni^3=(\sum_{i=1}^ni)^2=[\frac{n\times (n+1)}{2}]^2=\frac{n^2(n+1)^2}{4}i=1∑ni3=(i=1∑ni)2=[2n×(n+1)]2=4n2(n+1)2
习题集
Sum
…
- solution
前面提过了id=φ∗I,ϵ=μ∗Iid=φ*I,ϵ=\mu*Iid=φ∗I,ϵ=μ∗I,直接带公式即可id,ϵid,ϵid,ϵ就是hhh,φ,μφ,\muφ,μ就是fff,III是ggg
g(1)×S(n)=∑i=1nh(i)−∑d=2ng(d)S(⌊nd⌋)g(1)\times S(n)=\sum_{i=1}^nh(i)-\sum_{d=2}^ng(d)S(\lfloor\frac{n}{d}\rfloor)g(1)×S(n)=i=1∑nh(i)−d=2∑ng(d)S(⌊dn⌋) - code
#include <cstdio>
#include <map>
using namespace std;
#define maxn 2000000
#define int long long
map < int, int > mp1, mp2;
int cnt;
bool vis[maxn + 5];
int phi[maxn + 5], mu[maxn + 5], prime[maxn];
int sum1[maxn + 5], sum2[maxn + 5];void init() {mu[1] = phi[1] = 1;for( int i = 2;i <= maxn;i ++ ) {if( ! vis[i] ) prime[++ cnt] = i, mu[i] = -1, phi[i] = i - 1;for( int j = 1;j <= cnt && 1ll * i * prime[j] <= maxn;j ++ ) {vis[i * prime[j]] = 1;if( i % prime[j] == 0 ) {mu[i * prime[j]] = 0;phi[i * prime[j]] = phi[i] * prime[j];break;}mu[i * prime[j]] = -mu[i];phi[i * prime[j]] = phi[i] * ( prime[j] - 1 );}}for( int i = 1;i <= maxn;i ++ ) sum1[i] = sum1[i - 1] + phi[i], sum2[i] = sum2[i - 1] + mu[i];
}int find_phi( int n ) {if( n <= maxn ) return sum1[n];if( mp1[n] ) return mp1[n];int ans = n * ( n + 1 ) / 2;int r;for( int i = 2;i <= n;i = r + 1 ) {r = n / ( n / i );ans -= ( r - i + 1 ) * find_phi( n / i );}return mp1[n] = ans;
}int find_mu( int n ) {if( n <= maxn ) return sum2[n];if( mp2[n] ) return mp2[n];int ans = 1;int r;for( int i = 2;i <= n;i = r + 1 ) {r = n / ( n / i );ans -= ( r - i + 1 ) * find_mu( n / i );}return mp2[n] = ans;
}signed main() {init();int T, n;scanf( "%lld", &T );while( T -- ) {scanf( "%lld", &n );printf( "%lld %lld\n", find_phi( n ), find_mu( n ) );}return 0;
}
神犇和蒟蒻
…
- solution
AAA很简单,i2i^2i2的μ\muμ一定是000,所以AAA一定是111
考虑BBB,这里要用到欧拉函数的一个性质
如果i%j=0i\% j=0i%j=0,则φ(i×j)=φ(i)×jφ(i\times j)=φ(i)\times jφ(i×j)=φ(i)×j
所以φ(i2)=φ(i)×iφ(i^2)=φ(i)\times iφ(i2)=φ(i)×i,就神奇的转化为了上面的应用公式三 - code
#include <cstdio>
#include <map>
using namespace std;
#define int long long
#define mod 1000000007
#define maxn 1000000
map < int, int > mp;
int inv6, inv2, cnt;
bool vis[maxn + 5];
int prime[maxn], sum[maxn + 5], phi[maxn + 5];void init() {phi[1] = 1;for( int i = 2;i <= maxn;i ++ ) {if( ! vis[i] ) prime[++ cnt] = i, phi[i] = i - 1;for( int j = 1;j <= cnt && i * prime[j] <= maxn;j ++ ) {vis[i * prime[j]] = 1;if( i % prime[j] == 0 ) {phi[i * prime[j]] = phi[i] * prime[j];break;}phi[i * prime[j]] = phi[i] * ( prime[j] - 1 );}}for( int i = 1;i <= maxn;i ++ )sum[i] = ( sum[i - 1] + i * phi[i] ) % mod;
}int find_phi( int n ) {if( n <= maxn ) return sum[n];if( mp[n] ) return mp[n];int ans = n * ( n + 1 ) % mod * ( 2 * n + 1 ) % mod * inv6 % mod;for( int l = 2, r;l <= n;l = r + 1 ) {r = n / ( n / l );ans -= ( l + r ) * ( r - l + 1 ) % mod * find_phi( n / l ) % mod * inv2 % mod;if( ans < 0 ) ans += mod; }return mp[n] = ans;
}int qkpow( int x, int y ) {int ans = 1;while( y ) {if( y & 1 ) ans = ans * x % mod;x = x * x % mod;y >>= 1;}return ans;
}signed main() {init();int n;inv2 = qkpow( 2, mod - 2 );inv6 = qkpow( 6, mod - 2 );scanf( "%lld", &n );printf( "1\n%lld\n", find_phi( n ) );return 0;
}
简单的数学题
…
-
solution
取模的操作代码里加入就可以了,下面式子中就省略不写了
∑i=1n∑j=1ni×j×gcd(i,j)=∑i=1n∑j=1ni×j∑d∣i,d∣jφ(d)=∑d=1nφ(d)∑d∣i∑d∣jij\sum_{i=1}^n\sum_{j=1}^ni\times j\times gcd(i,j)=\sum_{i=1}^n\sum_{j=1}^ni\times j\sum_{d|i,d|j}φ(d)=\sum_{d=1}^nφ(d)\sum_{d|i}\sum_{d|j}iji=1∑nj=1∑ni×j×gcd(i,j)=i=1∑nj=1∑ni×jd∣i,d∣j∑φ(d)=d=1∑nφ(d)d∣i∑d∣j∑ij
=∑d=1nφ(d)×d2∑i=1⌊nd⌋i∑j=1⌊nd⌋j=∑d=1nφ(d)×d2(∑i=1⌊nd⌋i)2=\sum_{d=1}^nφ(d)\times d^2\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}i\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}j=\sum_{d=1}^nφ(d)\times d^2(\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}i)^2=d=1∑nφ(d)×d2i=1∑⌊dn⌋ij=1∑⌊dn⌋j=d=1∑nφ(d)×d2(i=1∑⌊dn⌋i)2
前面杜教筛,后面数论分块
-
code
#include <cstdio>
#include <map>
using namespace std;
#define int long long
#define maxn 5000000
map < int, int > mp;
int mod, cnt, inv;
bool vis[maxn + 5];
int prime[maxn], sum[maxn + 5], phi[maxn + 5];void init() {phi[1] = 1;for( int i = 2;i <= maxn;i ++ ) {if( ! vis[i] ) prime[++ cnt] = i, phi[i] = i - 1;for( int j = 1;j <= cnt && i * prime[j] <= maxn;j ++ ) {vis[i * prime[j]] = 1;if( i % prime[j] == 0 ) {phi[i * prime[j]] = phi[i] * prime[j] % mod;break;}elsephi[i * prime[j]] = phi[i] * ( prime[j] - 1 ) % mod;}}for( int i = 1;i <= maxn;i ++ )sum[i] = ( sum[i - 1] + phi[i] * i % mod * i % mod ) % mod;
}int calc( int n ) {n %= mod;return n * ( n + 1 ) / 2 % mod;
}int js( int n ) {n %= mod;return n * ( n + 1 ) % mod * ( 2 * n + 1 ) % mod * inv % mod;
}int solve( int n ) {if( n <= maxn ) return sum[n];if( mp[n] ) return mp[n];int ans = calc( n ) * calc( n ) % mod, r;for( int i = 2;i <= n;i = r + 1 ) {r = n / ( n / i );ans = ( ans - ( js( r ) - js( i - 1 ) ) * solve( n / i ) ) % mod;}return mp[n] = ans;
}int qkpow( int x, int y ) {int ans = 1;while( y ) {if( y & 1 ) ans = ans * x % mod;x = x * x % mod;y >>= 1;}return ans;
}signed main() {int n;scanf( "%lld %lld", &mod, &n );init();inv = qkpow( 6, mod - 2 );int r, ans = 0;for( int l = 1;l <= n;l = r + 1 ) {r = n / ( n / l );ans = ( ans + calc( n / l ) * calc( n / l ) % mod * ( solve( r ) - solve( l - 1 ) ) + mod ) % mod;}printf( "%lld", ( ans + mod ) % mod );return 0;
}