【杜教筛】初次见面，请多关照 (公式推导+题集)
常见积性函数
- \(\mu\)
- \(φ\)
- \(d(n)\)：\(n\)的约数个数
- \(σ(n)\)：\(n\)的约数和
- \(ϵ(n)\)：单位元函数，\(e(n)=[n=1]\) （完全积性）
- \(I(n)\)：\(I(n)=1\)（完全积性）
- \(id(n)\)：\(id(n)=n\) （完全积性）
\[\sum_{d|n}φ(d)=d
\]
\[\sum_{d|n}\mu(d)=[n=1]
\]
\[\sum_{d|n}φ(d)\times d=\frac{n\times φ(n)+[n=1]}{2}
\]
公式推导
求$$\sum_{i=1}^nf(i)$$
狄利克雷卷积
定义两个数论函数\(f,g\)的卷积为\(h\)
\[(f*g)(n)=h(n) \Rightarrow\ h(n)=\sum_{d|n}f(d)\times g(\frac{n}{d})=\sum_{d|n}g(d)\times f(\frac{n}{d})
\]
卷积满足乘法交换/结合律/分配律 (这不是重点😓)
\(f*g=g*f\)
\((f*g)*h=f*(g*h)\)
\((f+g)*h=f*h+g*h\)
杜教筛
杜教筛是一种以低于线性筛时间复杂度筛积性函数前缀和的筛法
时间复杂度\(O(n^{\frac{2}{3}})\) 不会证😔
将所求问题定义成函数，令\(S(n)=\sum_{i=1}^nf(i)\)
先推\(\sum_{i=1}^nh(i)\)
\[\sum_{i=1}^nh(i)
\]
\[=\sum_{i=1}^n\sum_{d|i}g(d)\times f(\frac{i}{d})
\]
\[=\sum_{d=1}^ng(d) \sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(i)
\]
\[=\sum_{d=1}^ng(d)S(\lfloor\frac{n}{d}\rfloor)
\]
将\(d=1\)的一项单独提出来\(\Downarrow\)
\[=g(1)\times S(n)+\sum_{d=2}^ng(d)S(\lfloor\frac{n}{d}\rfloor)
\]
\[g(1)\times S(n)=\sum_{i=1}^nh(i)-\sum_{d=2}^ng(d)S(\lfloor\frac{n}{d}\rfloor)
\]
一般套路情况下\(g(1)=1\)
实现
因为\(g,h\)都是我们自己选择构造出来的，通过上面化出的最后式子来看
我们希望\(h(i)\)的前缀和是好求的，换言之就是很快就能求出来的
然后按\(S(\lfloor\frac{n}{d}\rfloor)\)可以对后面式子进行分块处理
至于怎么选择\(g\)和\(h\)，说实话只能观察式子，或者结合以前的做题经验套路
常见卷积
Ⅰ\(S(n)=\sum_{i=1}^n\mu(i)\)
因为知道在莫比乌斯反演有个关于莫比乌斯函数的公式
\[\sum_{d|n}\mu(d)=[n=1]
\]
很容易想到\(ϵ=\mu*I\)
\[ϵ(n)=\sum_{d|n}\mu(d)\times I(\frac{n}{d})=\sum_{d|n}\mu(d)=[n=1]
\]
Ⅱ \(s(n)=\sum_{i=1}^nφ(i)\)
上述提到过一个公式
\[n=\sum_{d|n}φ(d)
\]
所以我们想到\(id=φ*I\)
\[id(n)=\sum_{d|n}φ(d)\times I(\frac{n}{d})=\sum_{d|n}φ(d)=n
\]
Ⅲ \(S(n)=\sum_{i=1}^nφ(i)\times i\)
乍一看似乎配不出什么鸟玩意儿
尝试从狄利克雷卷积入手
\[h=f*g=\sum_{d|n}f(d)\times g(\frac{n}{d})=\sum_{d|n}φ(d)\times d\times g(\frac{n}{d})
\]
我们想着如果能把\(f\)顺带的多余的\(d\)除掉，尝试给\(g\)配成\(id\)
\[h=\sum_{d|n}φ(d)\times d\times \frac{n}{d}=\sum_{d|n}\{φ(d)\times n\}=n\sum_{d|n}φ(d)=n^2
\]
\(i^2\)的前缀和，欸好巧哦，由公式可以快速算出，所以成功配出来了\(g,h\)
\[\sum_{i=1}^ni^2=\frac{n\times (n+1)\times (2n+1)}{6} \]\[\sum_{i=1}^ni^3=(\sum_{i=1}^ni)^2=[\frac{n\times (n+1)}{2}]^2=\frac{n^2(n+1)^2}{4} \]
习题集
Sum
- solution
前面提过了\(id=φ*I,ϵ=\mu*I\)，直接带公式即可\(id,ϵ\)就是\(h\)，\(φ,\mu\)就是\(f\)，\(I\)是\(g\)
\[g(1)\times S(n)=\sum_{i=1}^nh(i)-\sum_{d=2}^ng(d)S(\lfloor\frac{n}{d}\rfloor)
\]
- 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
\(A\)很简单，\(i^2\)的\(\mu\)一定是\(0\)，所以\(A\)一定是\(1\)
考虑\(B\)，这里要用到欧拉函数的一个性质
如果\(i\% j=0\)，则\(φ(i\times j)=φ(i)\times j\)
所以\(φ(i^2)=φ(i)\times 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
取模的操作代码里加入就可以了，下面式子中就省略不写了
\[\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}ij
\]
\[=\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
\]
前面杜教筛，后面数论分块
- 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;
}
else
phi[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;
}
