数论与数论函数
全文参考oi-wiki.org
基础知识
快速幂
ll qpow(ll a,ll b,ll mod){
ll res=1;
a=a%mod;
while(b){
if(b&1)res=res*a%mod;
a=a*a%mod;
b>>=1;
}
return res;
}
可以分块预处理任意进制的快速幂,或者“光速幂”。
快速乘法
用于 a ∗ b a*b a∗b爆 l o n g l o n g long long longlong的情况
ll qmul(ll a,ll b,ll mod){
ll res=0;
while(b){
if(b&1)res=(res+a)%mod;
a=(a+a)%mod;
b>>=1;
}
return res;
}
还有更快的方法,传说中的 O ( 1 ) O(1) O(1)快速乘:
ll qmul(ll x,ll y,ll mod){
return (x*y-(ll)((long double)x/mod*y)*mod+mod)%mod;
}
整除及其性质
素数
1、素数的个数:
素数计数函数:小于等于
x
x
x的素数的个数,用
π
(
x
)
\pi(x)
π(x)表示,随着
x
x
x的增大,有近似结论:
π
(
x
)
∼
x
ln
(
x
)
\pi(x)\sim\frac{x}{\ln(x)}
π(x)∼ln(x)x
2、素数判定:
暴力判定:
O
(
n
)
O(\sqrt {n})
O(n)
素性测试:
O
(
k
log
3
n
)
O(k \log^3n)
O(klog3n)(详见匡斌模板)
最大公约数
欧几里得算法
ll gcd(ll a,ll b){
return b==0?a:gcd(b,a%b);
}
多个数的最大公约数
g c d ( a , b , c ) = g c d ( g c d ( a , b ) , c ) gcd(a,b,c)=gcd(gcd(a,b),c) gcd(a,b,c)=gcd(gcd(a,b),c)
最小公倍数
l c m ( a , b ) = a ∗ b / g c d ( a , b ) lcm(a,b)=a*b/gcd(a,b) lcm(a,b)=a∗b/gcd(a,b)
多个数的最小公倍数
l
c
m
(
a
,
b
,
c
)
=
l
c
m
(
l
c
m
(
a
,
b
)
,
c
)
lcm(a,b,c)=lcm(lcm(a,b),c)
lcm(a,b,c)=lcm(lcm(a,b),c)
谨记
l
c
m
(
a
,
b
,
c
)
≠
a
∗
b
∗
c
/
g
c
d
(
a
,
b
,
c
)
lcm(a,b,c)\neq a*b*c/gcd(a,b,c)
lcm(a,b,c)=a∗b∗c/gcd(a,b,c)
扩展欧几里得定理
常用于求解 a x + b y = g c d ( a , b ) ax+by=gcd(a,b) ax+by=gcd(a,b)
int Exgcd(int a, int b, int &x, int &y) {
if (!b) {
x = 1;
y = 0;
return a;
}
int d = Exgcd(b, a % b, x, y);
int t = x;
x = y;
y = t - (a / b) * y;
return d;
}
函数返回值为 g c d ( a , b ) gcd(a,b) gcd(a,b),如果 a x + b y = c ax+by=c ax+by=c中 c c c不是 g c d ( a , b ) gcd(a,b) gcd(a,b)的倍数,则方程无解。
欧拉函数
φ
(
n
)
\varphi(n)
φ(n)表示小于等于
n
n
n和
n
n
n互质的数的个数。
特别地,令
φ
(
1
)
=
1
\varphi(1)=1
φ(1)=1
设
n
=
p
1
k
1
p
2
k
2
…
p
s
k
s
n=p_1^{k_1}p_2^{k_2}…p_s^{k_s}
n=p1k1p2k2…psks,那么
φ
(
n
)
=
n
∗
∏
i
=
1
s
p
i
−
1
p
i
\varphi(n)=n*\prod_{i=1}^s{\frac{p_i-1}{p_i}}
φ(n)=n∗∏i=1spipi−1
若
n
=
p
k
n=p^k
n=pk,
p
p
p是质数,那么
φ
(
n
)
=
p
k
−
1
∗
(
p
−
1
)
\varphi(n)=p^{k-1}*(p-1)
φ(n)=pk−1∗(p−1)
欧拉定理
若
g
c
d
(
a
,
m
)
=
1
gcd(a,m)=1
gcd(a,m)=1,则
a
φ
(
m
)
≡
1
(
m
o
d
m
)
a^{\varphi(m)}\equiv 1\pmod{m}
aφ(m)≡1(modm)
当m为素数时即为费马小定理。
扩展欧拉定理
a
b
≡
{
a
b
m
o
d
φ
(
p
)
,
gcd
(
a
,
p
)
=
1
a
b
,
gcd
(
a
,
p
)
≠
1
,
b
<
φ
(
p
)
a
b
m
o
d
φ
(
p
)
+
φ
(
p
)
,
gcd
(
a
,
p
)
≠
1
,
b
≥
φ
(
p
)
(
m
o
d
p
)
a^b\equiv\begin{cases} a^{b\bmod{\varphi(p)}},\,&\gcd(a,\,p)=1\\ a^b,& \gcd(a,p)\ne1,b<\varphi(p)\\ a^{b\bmod{\varphi(p)+\varphi(p)}},&\gcd(a,p)\ne1,b\ge\varphi(p)\end{cases}\pmod{p}
ab≡⎩⎪⎨⎪⎧abmodφ(p),ab,abmodφ(p)+φ(p),gcd(a,p)=1gcd(a,p)=1,b<φ(p)gcd(a,p)=1,b≥φ(p)(modp)
这个定理常用来欧拉降幂
同余与同余方程
剩余系
完全剩余系:对于
n
n
n,就是集合
{
0
,
1
,
…
,
n
−
1
}
\{0,1,…,n-1\}
{0,1,…,n−1}。
简化剩余系:又叫既约剩余系,对于
n
n
n,是完全剩余系中与
n
n
n互质的元素构成的子集。
裴蜀定理
设
a
,
b
a,b
a,b是不全为零的整数,则存在整数
x
,
y
x,y
x,y,使得
a
x
+
b
y
=
gcd
(
a
,
b
)
ax+by=\gcd(a,b)
ax+by=gcd(a,b)
这个定理可以扩展到
n
n
n个数的情况
乘法逆元
如果一个线性同余方程 a x ≡ 1 ( m o d b ) ax\equiv 1\pmod b ax≡1(modb),则 x x x称为 a m o d b a \bmod b amodb 的逆元,记作 a − 1 a^{-1} a−1.
扩展欧几里得法
void exgcd(int a, int b, int& x, int& y) {
if (b == 0) {
x = 1, y = 0;
return;
}
exgcd(b, a % b, y, x);
y -= a / b * x;
}
费马小定理
a − 1 = a p − 2 m o d p a^{-1}=a^{p-2}\bmod p a−1=ap−2modp,可以用快速幂求 a p − 2 a^{p-2} ap−2
线性求逆元
求
1
,
2
,
…
,
n
1,2,…,n
1,2,…,n中每个数关于
p
p
p的逆元
首先,有
1
−
1
≡
1
(
m
o
d
p
)
1^{-1}\equiv 1\pmod p
1−1≡1(modp)
设
p
=
k
i
+
j
,
j
<
i
,
1
<
i
<
p
p=ki+j,j<i,1<i<p
p=ki+j,j<i,1<i<p
有
k
i
+
j
≡
0
(
m
o
d
p
)
ki+j\equiv 0\pmod p
ki+j≡0(modp)
两边同时乘
i
−
1
,
j
−
1
i^{-1},j^{-1}
i−1,j−1:
k
j
−
1
+
i
−
1
≡
0
(
m
o
d
p
)
kj^{-1}+i^{-1}\equiv 0\pmod p
kj−1+i−1≡0(modp)
i
−
1
≡
−
k
j
−
1
(
m
o
d
p
)
i^{-1}\equiv -kj^{-1}\pmod p
i−1≡−kj−1(modp)
i
−
1
≡
−
(
p
i
)
(
p
m
o
d
i
)
−
1
(
m
o
d
p
)
i^{-1}\equiv -(\frac{p}{i})(p\bmod i)^{-1}\pmod p
i−1≡−(ip)(pmodi)−1(modp)
inv[1] = 1;
for (int i = 2; i <= n; ++i)
inv[i] = (long long)(p - p / i) * inv[p % i] % p;
线性求任意 n n n个数的逆元
如果要求任意给定
n
n
n个数的逆元(
1
≤
a
i
<
p
1\le a_i<p
1≤ai<p)。
首先计算
n
n
n个数的前缀积,记为
s
n
s_n
sn,求得其逆元,记为
s
v
n
sv_n
svn,当把它乘上
a
n
a_n
an时,就会得到
a
1
a_1
a1到
a
n
−
1
a_{n-1}
an−1的积的逆元,记为
s
v
n
−
1
sv_{n-1}
svn−1,同理可得所有的
s
v
i
sv_i
svi,于是
a
i
−
1
a_i^{-1}
ai−1就可以用
s
i
−
1
∗
s
v
i
s_{i-1}*sv_i
si−1∗svi得到。
线性同余方程
形如
a
x
≡
b
(
m
o
d
c
)
ax\equiv b\pmod c
ax≡b(modc)的方程被称为线性同余方程,
其解和
a
x
+
c
y
=
b
ax+cy=b
ax+cy=b是等价的
我们可以先解出
a
x
+
b
y
=
gcd
(
a
,
b
)
ax+by=\gcd(a,b)
ax+by=gcd(a,b)的一组
x
0
,
y
0
x_0,y_0
x0,y0,然后两边同时除以
gcd
(
a
,
b
)
\gcd(a,b)
gcd(a,b),在乘
c
c
c。就找到了方程的一个解。
若
gcd
(
a
,
b
)
=
1
\gcd(a,b)=1
gcd(a,b)=1,且
x
0
,
y
0
x_0,y_0
x0,y0为方程
a
x
+
b
y
=
c
ax+by=c
ax+by=c的一组解,则该方程的任意解可表示为:
x
=
x
0
+
b
t
,
y
=
y
0
−
a
t
x=x_0+bt,y=y_0-at
x=x0+bt,y=y0−at,对任意整数t都成立。
int ex_gcd(int a, int b, int& x, int& y) {
if (b == 0) {
x = 1;
y = 0;
return a;
}
int d = ex_gcd(b, a % b, x, y);
int temp = x;
x = y;
y = temp - a / b * y;
return d;
}
bool liEu(int a, int b, int c, int& x, int& y) {
int d = ex_gcd(a, b, x, y);
if (c % d != 0) return 0;
int k = c / d;
x *= k;
y *= k;
return 1;
}
中国剩余定理
一般这玩意都是用扩展中国剩余定理(模数不互质)
#include<iostream>
#include<vector>
#include<algorithm>
#include<queue>
#include<cstring>
#include<cstdio>
using namespace std;
typedef long long lt;
lt read()
{
lt f=1,x=0;
char ss=getchar();
while(ss<'0'||ss>'9'){if(ss=='-')f=-1;ss=getchar();}
while(ss>='0'&&ss<='9'){x=x*10+ss-'0';ss=getchar();}
return f*x;
}
const int maxn=100010;
int n;
lt ai[maxn],bi[maxn];
lt mul(lt a,lt b,lt mod)
{
lt res=0;
while(b>0)
{
if(b&1) res=(res+a)%mod;
a=(a+a)%mod;
b>>=1;
}
return res;
}
lt exgcd(lt a,lt b,lt &x,lt &y)
{
if(b==0){x=1;y=0;return a;}
lt gcd=exgcd(b,a%b,x,y);
lt tp=x;
x=y; y=tp-a/b*y;
return gcd;
}
lt excrt()
{
lt x,y,k;
lt M=bi[1],ans=ai[1];//第一个方程的解特判
for(int i=2;i<=n;i++)
{
lt a=M,b=bi[i],c=(ai[i]-ans%b+b)%b;//ax≡c(mod b)
lt gcd=exgcd(a,b,x,y),bg=b/gcd;
if(c%gcd!=0) return -1; //判断是否无解,然而这题其实不用
x=mul(x,c/gcd,bg);
ans+=x*M;//更新前k个方程组的答案
M*=bg;//M为前k个m的lcm
ans=(ans%M+M)%M;
}
return (ans%M+M)%M;
}
int main()
{
n=read();
for(int i=1;i<=n;++i)
bi[i]=read(),ai[i]=read();//读入时看好先读哪个
printf("%lld",excrt());
return 0;
}
BSGS
常用来在
O
(
p
)
O(\sqrt{p})
O(p)的时间复杂度内求解
a
x
≡
b
m
o
d
p
a^x\equiv b\bmod p
ax≡bmodp,其中
a
a
a和
p
p
p互质,方程的解
x
x
x满足
0
≤
x
<
p
0\le x< p
0≤x<p。
令
x
=
A
⌈
p
⌉
−
B
x=A\lceil\sqrt{p}\rceil-B
x=A⌈p⌉−B,其中
0
≤
A
,
B
≤
⌈
p
⌉
0\le A,B\le \lceil\sqrt{p}\rceil
0≤A,B≤⌈p⌉,则有
a
A
⌈
p
⌉
−
B
≡
b
(
m
o
d
p
)
a^{A\lceil\sqrt{p}\rceil-B}\equiv b\pmod p
aA⌈p⌉−B≡b(modp),则有
a
A
⌈
p
⌉
≡
b
a
B
a^{A\lceil\sqrt{p}\rceil}\equiv ba^B
aA⌈p⌉≡baB.
首先,计算并记录所有
b
a
B
ba^B
baB的取值,再枚举
A
A
A,计算
a
A
⌈
p
⌉
a^{A\lceil\sqrt{p}\rceil}
aA⌈p⌉是否有
b
a
B
ba^B
baB与之相等,从而可以得到所有的
x
x
x,
x
=
A
⌈
p
⌉
−
B
x=A\lceil\sqrt{p}\rceil-B
x=A⌈p⌉−B
待补充……
原根
阶
若
(
a
,
m
)
=
1
(a,m)=1
(a,m)=1,使得
a
l
≡
1
(
m
o
d
m
)
a^l\equiv 1\pmod m
al≡1(modm)成立最小的
l
l
l,称为
a
a
a关于模
m
m
m的阶,记为
o
r
d
m
a
ord_ma
ordma。
若
o
r
d
m
a
=
l
ord_ma=l
ordma=l,则
o
r
d
m
a
t
=
l
(
t
,
l
)
ord_ma^t=\frac{l}{(t,l)}
ordmat=(t,l)l
由欧拉定理,设
o
r
d
m
a
=
l
ord_ma=l
ordma=l,则
a
n
≡
1
(
m
o
d
m
)
a^n\equiv 1\pmod m
an≡1(modm)当且仅当
l
∣
n
l|n
l∣n,特别地,
l
∣
φ
(
m
)
l|\varphi(m)
l∣φ(m).
性质:
1、设
p
p
p是素数,
o
r
d
p
a
=
l
ord_pa=l
ordpa=l,那么有且仅有
φ
(
l
)
\varphi(l)
φ(l)个关于模
p
p
p的阶为
l
l
l且两两互不同余的数
2、设
o
r
d
m
a
=
l
ord_ma=l
ordma=l,则
1
,
a
,
a
2
,
…
,
a
l
−
1
1,a,a^2,…,a^{l-1}
1,a,a2,…,al−1关于模
m
m
m两两互不同余
3、设
p
p
p是素数,
l
∣
p
−
1
l|p-1
l∣p−1,则存在
φ
(
l
)
\varphi(l)
φ(l)个关于模
p
p
p的阶为
l
l
l且两两互不同余的数。
4、若
m
=
p
1
a
1
p
2
a
2
…
p
k
a
k
m=p_1^{a_1}p_2^{a_2}…p_k^{a_k}
m=p1a1p2a2…pkak,则
o
r
d
m
a
=
[
o
r
d
p
1
a
1
,
o
r
d
p
2
a
2
,
…
,
o
r
d
p
k
a
k
]
ord_ma=[ord_{p_1}^{a_1},ord_{p_2}^{a_2},…,ord_{p_k}^{a_k}]
ordma=[ordp1a1,ordp2a2,…,ordpkak]
原根
(
g
,
m
)
=
1
,
(g,m)=1,
(g,m)=1,若
o
r
d
m
g
=
φ
(
m
)
ord_mg=\varphi(m)
ordmg=φ(m),则称
g
g
g为
m
m
m的一个原根。
g
g
g为
m
m
m的一个原根当且仅当{
g
,
g
2
,
…
,
g
φ
(
m
)
g,g^2,…,g^{\varphi(m)}
g,g2,…,gφ(m)}构成模
m
m
m的一个既约剩余系。
是否有原根
若 m m m有原根,则 m m m一定是: 2 , 4 , p a , 2 p a 2,4,p^a,2p^a 2,4,pa,2pa,其中 p p p是奇素数, a a a为正整数。
求所有原根
设 g g g为 m m m的一个原根,则集合 S = { g s ∣ 1 ≤ s ≤ φ ( m ) , ( s , φ ( m ) ) = 1 } S=\{g^s|1\le s\le \varphi(m),(s,\varphi(m))=1\} S={gs∣1≤s≤φ(m),(s,φ(m))=1}给出 m m m的全部原根。
求一个原根
对
p
−
1
p-1
p−1进行质因数分解得到不同的质因子
d
1
,
d
2
,
…
,
d
m
d_1,d_2,…,d_m
d1,d2,…,dm,对于任意的
1
<
a
<
p
1<a<p
1<a<p,要判定
a
a
a是否是模
p
p
p的原根,只需要检验
a
p
−
1
d
1
,
a
p
−
1
d
2
,
…
,
a
p
−
1
d
m
a^{\frac{p-1}{d_1}},a^{\frac{p-1}{d_2}},…,a^{\frac{p-1}{d_m}}
ad1p−1,ad2p−1,…,admp−1这
m
m
m个数中是否存在一个数在模
p
p
p意义下与1同余。若存在,则
a
a
a不是
p
p
p的原根;若不存在,则
a
a
a是
p
p
p的原根。
或者说:
(
g
,
m
)
=
1
(g,m)=1
(g,m)=1,设
p
1
,
p
2
,
…
,
p
k
p_1,p_2,…,p_k
p1,p2,…,pk是
φ
(
m
)
\varphi(m)
φ(m)的所有不同的素因子,则
g
g
g是
m
m
m的原根,当且仅当对任意
1
≤
i
≤
k
1\le i\le k
1≤i≤k,都有
g
φ
(
m
)
p
i
≢
1
(
m
o
d
m
)
g^{\frac{\varphi(m)}{p_i}}\not\equiv 1\pmod m
gpiφ(m)≡1(modm)
注:998244353 的原根是3;1e9+7的原根是5
下面给出暴力求原根的代码:(python版本)
def check(p):
fac = []
n = p - 1
i = 2
while i * i <= n:
if n % i == 0:
fac.append(i)
while n % i == 0:
n = n // i
i += 1
fac.append(n)
i = 2
while i < p:
flag = True
for j in fac:
if pow(i, (p - 1) // j, p) == 1:
flag = False
break
if flag:
return i
i += 1
n = int(input())
print(check(n))
卢卡斯定理(大组合数取模)
Lucas定理内容如下:
C
n
m
m
o
d
p
=
C
⌊
n
p
⌋
⌊
m
p
⌋
⋅
C
n
m
o
d
p
m
m
o
d
p
m
o
d
p
C_n^m\bmod p=C_{\lfloor\frac{n}{p}\rfloor}^{\lfloor\frac{m}{p}\rfloor}\cdot C_{n\bmod p}^{m\bmod p}\bmod p
Cnmmodp=C⌊pn⌋⌊pm⌋⋅Cnmodpmmodpmodp
long long Lucas(long long n, long long m, long long p) {
if (m == 0) return 1;
return (C(n % p, m % p, p) * Lucas(n / p, m / p, p)) % p;
}
exLucas定理
当模数
p
p
p不是素数时,就需要用到exLucas定理,待补充。。
下面是代码,求
C
n
m
m
o
d
p
C_n^m\bmod p
Cnmmodp
#include<bits/stdc++.h>
#define ll long long
using namespace std;
#ifndef Fading
inline char gc(){
static char now[1<<16],*S,*T;
if (T==S){T=(S=now)+fread(now,1,1<<16,stdin);if (T==S) return EOF;}
return *S++;
}
#endif
#ifdef Fading
#define gc getchar
#endif
void exgcd(ll a,ll b,ll &x,ll &y){
if (!b) return (void)(x=1,y=0);
exgcd(b,a%b,x,y);
ll tmp=x;x=y;y=tmp-a/b*y;
}
ll gcd(ll a,ll b){
if (b==0) return a;
return gcd(b,a%b);
}
inline ll INV(ll a,ll p){
ll x,y;
exgcd(a,p,x,y);
return (x+p)%p;
}
inline ll lcm(ll a,ll b){
return a/gcd(a,b)*b;
}
inline ll mabs(ll x){
return (x>0?x:-x);
}
inline ll fast_mul(ll a,ll b,ll p){
ll t=0;a%=p;b%=p;
while (b){
if (b&1LL) t=(t+a)%p;
b>>=1LL;a=(a+a)%p;
}
return t;
}
inline ll fast_pow(ll a,ll b,ll p){
ll t=1;a%=p;
while (b){
if (b&1LL) t=(t*a)%p;
b>>=1LL;a=(a*a)%p;
}
return t;
}
inline ll read(){
ll x=0,f=1;char ch=gc();
while (!isdigit(ch)) {if (ch=='-') f=-1;ch=gc();}
while (isdigit(ch)) x=x*10+ch-'0',ch=gc();
return x*f;
}
inline ll F(ll n,ll P,ll PK){
if (n==0) return 1;
ll rou=1;//循环节
ll rem=1;//余项
for (ll i=1;i<=PK;i++){
if (i%P) rou=rou*i%PK;
}
rou=fast_pow(rou,n/PK,PK);
for (ll i=PK*(n/PK);i<=n;i++){
if (i%P) rem=rem*(i%PK)%PK;
}
return F(n/P,P,PK)*rou%PK*rem%PK;
}
inline ll G(ll n,ll P){
if (n<P) return 0;
return G(n/P,P)+(n/P);
}
inline ll C_PK(ll n,ll m,ll P,ll PK){
ll fz=F(n,P,PK),fm1=INV(F(m,P,PK),PK),fm2=INV(F(n-m,P,PK),PK);
ll mi=fast_pow(P,G(n,P)-G(m,P)-G(n-m,P),PK);
return fz*fm1%PK*fm2%PK*mi%PK;
}
ll A[1001],B[1001];
//x=B(mod A)
inline ll exLucas(ll n,ll m,ll P){
ll ljc=P,tot=0;
for (ll tmp=2;tmp*tmp<=P;tmp++){
if (!(ljc%tmp)){
ll PK=1;
while (!(ljc%tmp)){
PK*=tmp;ljc/=tmp;
}
A[++tot]=PK;B[tot]=C_PK(n,m,tmp,PK);
}
}
if (ljc!=1){
A[++tot]=ljc;B[tot]=C_PK(n,m,ljc,ljc);
}
ll ans=0;
for (ll i=1;i<=tot;i++){
ll M=P/A[i],T=INV(M,A[i]);
ans=(ans+B[i]*M%P*T%P)%P;
}
return ans;
}
signed main(){
ll n=read(),m=read(),P=read();
printf("%lld\n",exLucas(n,m,P));
return 0;
}
数论函数
线性筛法
线性筛法求素数:
void Prime(){
memset(vis,false,sizeof(vis));
memset(prime,0,sizeof(prime));
for(int i=2;i<=maxn;i++){
if(!vis[i]){
prime[++cnt]=i;
}
for(int j=1;j<=cnt&&i*prime[j]<=maxn;j++){
vis[i*prime[j]]=true;
if(i%prime[j]==0){
break;
}
}
}
}
线性筛法求欧拉函数:
void euler(){
m = 0;
memset(v, 0, sizeof(v));
for(int i = 2; i < maxn; i++){
if(v[i] == 0){
v[i] = i;
p[m++] = i;
phi[i] = i - 1;
}
for(int j = 0; j < m; j++){
if(p[j] > v[i] || p[j] > maxn / i)
break;
v[i * p[j]] = p[j];
phi[i * p[j]] = phi[i] * (i % p[j] ? p[j] - 1 : p[j]);
}
}
}
线性筛法求莫比乌斯函数:
void pre() {
mu[1] = 1;
for (int i = 2; i <= 1e7; ++i) {
if (!v[i]) mu[i] = -1, p[++tot] = i;
for (int j = 1; j <= tot && i <= 1e7 / p[j]; ++j) {
v[i * p[j]] = 1;
if (i % p[j] == 0) {
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = -mu[i];
}
}
线性筛法求约数个数:
void pre() {
d[1] = 1;
for (int i = 1; i <= n; ++i) {
if (!v[i]) v[i] = 1, p[++tot] = i, d[i] = 2, num[i] = 1;
for (int j = 1; j <= tot && i <= n / p[j]; ++j) {
v[p[j] * i] = 1;
if (i % p[j] == 0) {
num[i * p[j]] = num[i] + 1;
d[i * p[j]] = d[i] / num[i * p[j]] * (num[i * p[j]] + 1);
break;
} else {
num[i * p[j]] = 1;
d[i * p[j]] = d[i] * 2;
}
}
}
}
线性筛法求约数和:
void pre() {
g[1] = f[1] = 1;
for (int i = 2; i <= n; ++i) {
if (!v[i]) v[i] = 1, p[++tot] = i, g[i] = i + 1, f[i] = i + 1;
for (int j = 1; j <= tot && i <= n / p[j]; ++j) {
v[p[j] * i] = 1;
if (i % p[j] == 0) {
g[i * p[j]] = g[i] * p[j] + 1;
f[i * p[j]] = f[i] / g[i] * g[i * p[j]];
break;
} else {
f[i * p[j]] = f[i] * f[p[j]];
g[i * p[j]] = 1 + p[j];
}
}
}
for (int i = 1; i <= n; ++i) f[i] = (f[i - 1] + f[i]) % Mod;
}
非常见积性函数的线性筛法
这里主要讨论求两个积性函数狄利克雷卷积的情况,这个函数显然也是积性函数。
设
l
o
w
(
i
)
low(i)
low(i)表示
p
1
a
1
p_1^{a_1}
p1a1,对于线性筛法最为关键的
i
%
p
j
=
0
i\%p_j=0
i%pj=0;
有了
l
o
w
(
i
)
low(i)
low(i),我们可以分两种情况来讨论。
1、若
l
o
w
(
i
)
=
i
low(i)=i
low(i)=i,此时
i
i
i一定是某个素数的幂次的形式,(否则会被
b
r
e
a
k
break
break),
此时只要我们可以快速地从
f
(
p
k
)
f(p^k)
f(pk)来更新
f
(
p
k
+
1
)
f(p^{k+1})
f(pk+1)
2、若
l
o
w
(
i
)
!
=
i
{low(i)}!=i
low(i)!=i,此时
i
/
l
o
w
(
i
)
i/low(i)
i/low(i)一定与
l
o
w
(
i
)
∗
p
j
low(i)*p_j
low(i)∗pj是互质的,
我们可以直接利用积性函数的性质去更新。
下面是伪代码:
vis[1] = low[1] = 1; H[1] = 初始化
for(int i = 2; i <= N; i++) {
if(!vis[i]) prime[++tot] = i, mu[i] = -1, H[i] = 质数的情况, low[i] = i;
for(int j = 1; j <= tot && i * prime[j] <= N; j++) {
vis[i * prime[j]] = 1;
if(!(i % prime[j])) {
low[i * prime[j]] = (low[i] * prime[j]);
if(low[i] == i) H[i * prime[j]] = 特殊判断;
else H[i * prime[j]] = H[i / low[i]] * H[prime[j] * low[i]];
break;
}
H[i * prime[j]] = H[i] * H[prime[j]];
low[i * prime[j]] = prime[j];
}
}
莫比乌斯反演
数论分块
引理1
∀ a , b , c ∈ Z , ⌊ a b c ⌋ = ⌊ ⌊ a b ⌋ c ⌋ \forall a,b,c\in Z,\lfloor \frac{a}{bc} \rfloor=\lfloor\frac{\lfloor\frac{a}{b}\rfloor}{c}\rfloor ∀a,b,c∈Z,⌊bca⌋=⌊c⌊ba⌋⌋
引理2
∀
n
∈
N
,
∣
{
⌊
n
d
∣
d
∈
N
⌋
}
∣
≤
⌊
2
n
⌋
\forall n\in N,\vert\{\lfloor\frac{n}{d}|d\in N \rfloor\}\vert \leq \lfloor 2\sqrt{n} \rfloor
∀n∈N,∣{⌊dn∣d∈N⌋}∣≤⌊2n⌋
∣
V
∣
表
示
集
合
V
的
元
素
个
数
\vert V \vert表示集合V的元素个数
∣V∣表示集合V的元素个数
数论分块
对于含有
⌊
n
i
⌋
\lfloor \frac{n}{i} \rfloor
⌊in⌋的求和式子(n为常数)
对于任意一个
i
(
i
≤
n
)
i(i \leq n)
i(i≤n),我们需要找到一个最大的
j
(
i
≤
j
≤
n
)
j(i \leq j \leq n)
j(i≤j≤n),使得
⌊
n
i
⌋
=
⌊
n
j
⌋
.
\lfloor\frac{n}{i}\rfloor=\lfloor \frac{n}{j} \rfloor.
⌊in⌋=⌊jn⌋.则有
j
=
⌊
n
⌊
n
i
⌋
⌋
j=\lfloor \frac{n}{\lfloor\frac{n}{i}\rfloor}\rfloor
j=⌊⌊in⌋n⌋.
所以我们可以每次以
[
i
,
j
]
[ i,j]
[i,j]为一块,分块求和即可。
积性函数
定义:
若 g c d ( x , y ) = 1 且 f ( x y ) = f ( x ) f ( y ) , 则 f ( n ) 为 积 性 函 数 若gcd(x,y)=1且f(xy)=f(x)f(y),则f(n)为积性函数 若gcd(x,y)=1且f(xy)=f(x)f(y),则f(n)为积性函数
性质:
若
f
(
x
)
和
g
(
x
)
均
为
积
性
函
数
,
则
以
下
函
数
也
为
积
性
函
数
:
若f(x)和g(x)均为积性函数,则以下函数也为积性函数:
若f(x)和g(x)均为积性函数,则以下函数也为积性函数:
h
(
x
)
=
f
(
x
p
)
\quad\quad\quad\quad h(x)=f(x^p)
h(x)=f(xp)
h
(
x
)
=
f
p
(
x
)
\quad\quad\quad\quad h(x)=f^p(x)
h(x)=fp(x)
h
(
x
)
=
f
(
x
)
g
(
x
)
\quad\quad\quad\quad h(x)=f(x)g(x)
h(x)=f(x)g(x)
h
(
x
)
=
∑
d
∣
x
f
(
d
)
g
(
x
d
)
\quad\quad\quad\quad h(x)=\sum\limits_{d|x}f(d)g(\frac{x}{d})
h(x)=d∣x∑f(d)g(dx)
常见积性函数:
1.
幺
元
函
数
:
ϵ
(
n
)
=
[
n
=
1
]
1.幺元函数: \epsilon(n)=[n=1]
1.幺元函数:ϵ(n)=[n=1]
2.
常
函
数
:
o
n
e
(
x
)
=
1
2.常函数:one(x)=1
2.常函数:one(x)=1
3.
标
号
函
数
:
i
d
(
x
)
=
x
3.标号函数:id(x)=x
3.标号函数:id(x)=x
4.
欧
拉
函
数
:
φ
(
n
)
=
∑
i
=
1
n
[
g
c
d
(
i
,
n
)
=
1
]
4.欧拉函数:\varphi(n)=\sum_{i=1}^n[gcd(i,n)=1]
4.欧拉函数:φ(n)=∑i=1n[gcd(i,n)=1]
5.
除
数
函
数
:
5.除数函数:
5.除数函数:
σ
(
k
,
x
)
=
∑
d
∣
x
d
k
\quad\quad\quad \sigma(k,x)=\sum_{d|x}d^k
σ(k,x)=∑d∣xdk
当
k
=
1
时
,
表
示
x
的
因
子
之
和
,
简
记
为
σ
(
n
)
\quad当k=1时,表示x的因子之和,简记为\sigma(n)
当k=1时,表示x的因子之和,简记为σ(n)
当
k
=
0
时
,
表
示
x
的
因
子
个
数
,
简
记
为
d
(
n
)
\quad当k=0时,表示x的因子个数,简记为d(n)
当k=0时,表示x的因子个数,简记为d(n)
\quad
k省略时默认为1
6.
莫
比
乌
斯
函
数
:
μ
(
n
)
=
{
1
n
=
1
0
∃
d
:
d
2
∣
n
(
−
1
)
ω
(
n
)
o
t
h
e
r
w
i
s
e
6.莫比乌斯函数:\mu(n)=\begin{cases}1&n=1\\0&\exists d:d^2|n\\(-1)^{\omega(n)}&otherwise\end{cases}
6.莫比乌斯函数:μ(n)=⎩⎪⎨⎪⎧10(−1)ω(n)n=1∃d:d2∣notherwise其中
ω
(
n
)
表
示
n
的
本
质
不
同
的
质
因
子
个
数
,
是
一
个
加
性
函
数
\omega(n)表示n的本质不同的质因子个数,是一个加性函数
ω(n)表示n的本质不同的质因子个数,是一个加性函数
Dirichlet 卷积
定义
两 个 数 论 函 数 f , g 的 D i r i c h l e t 卷 积 为 ( f ∗ g ) ( n ) = ∑ d ∣ n f ( d ) g ( n d ) 两个数论函数f,g的Dirichlet卷积为(f*g)(n)=\sum\limits_{d|n}f(d)g(\frac{n}{d}) 两个数论函数f,g的Dirichlet卷积为(f∗g)(n)=d∣n∑f(d)g(dn)
性质
D
i
r
i
c
h
l
e
t
卷
积
满
足
交
换
律
、
结
合
律
和
分
配
率
.
Dirichlet卷积满足交换律、结合律和分配率.
Dirichlet卷积满足交换律、结合律和分配率.
交
换
律
:
f
∗
g
=
g
∗
f
;
交换律:f*g=g*f;
交换律:f∗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;
ε
为
D
i
r
i
c
h
l
e
t
卷
积
中
的
单
位
元
(
任
何
函
数
卷
ε
都
是
其
本
身
)
\varepsilon 为Dirichlet卷积中的单位元(任何函数卷\varepsilon都是其本身)
ε为Dirichlet卷积中的单位元(任何函数卷ε都是其本身)
例子
ε
=
μ
∗
1
⇔
ε
(
n
)
=
∑
d
∣
n
μ
(
d
)
\varepsilon =\mu*1\Leftrightarrow\varepsilon(n)=\sum\limits_{d|n}\mu(d)
ε=μ∗1⇔ε(n)=d∣n∑μ(d)
d
=
1
∗
1
⇔
d
(
n
)
=
∑
d
∣
n
1
d =1*1 \Leftrightarrow d(n)=\sum\limits_{d|n}1
d=1∗1⇔d(n)=d∣n∑1
σ
=
d
∗
1
⇔
σ
(
n
)
=
∑
d
∣
n
d
\sigma=d*1\Leftrightarrow\sigma(n)=\sum\limits_{d|n}d
σ=d∗1⇔σ(n)=d∣n∑d
φ
=
μ
∗
i
d
⇔
φ
(
n
)
=
∑
d
∣
n
d
⋅
μ
(
n
d
)
\varphi =\mu*id\Leftrightarrow\varphi(n)=\sum\limits_{d|n}d\cdot \mu(\frac{n}{d})
φ=μ∗id⇔φ(n)=d∣n∑d⋅μ(dn)
i
d
=
φ
∗
1
⇔
n
=
∑
d
∣
n
φ
(
d
)
id=\varphi*1\Leftrightarrow n=\sum\limits_{d|n}\varphi(d)
id=φ∗1⇔n=d∣n∑φ(d)
g
=
f
∗
μ
⇔
f
=
g
∗
1
g=f*\mu\Leftrightarrow f=g*1
g=f∗μ⇔f=g∗1
莫比乌斯反演
μ
(
n
)
=
{
1
n
=
1
0
n
含
有
平
方
因
子
(
−
1
)
k
k
为
n
的
本
质
不
同
质
因
子
个
数
\mu(n)=\begin{cases}1&n=1\\0&n含有平方因子\\(-1)^k&k为n的本质不同质因子个数\end{cases}
μ(n)=⎩⎪⎨⎪⎧10(−1)kn=1n含有平方因子k为n的本质不同质因子个数
反演重要结论:
ε
(
n
)
=
∑
d
∣
n
μ
(
d
)
=
[
n
=
1
]
\varepsilon(n)=\sum\limits_{d|n}\mu(d)=[n=1]
ε(n)=d∣n∑μ(d)=[n=1]
[
g
c
d
(
i
,
j
)
=
1
]
⇔
∑
d
∣
g
c
d
(
i
,
j
)
μ
(
d
)
[gcd(i,j)=1] \Leftrightarrow \sum\limits_{d|gcd(i,j)}\mu(d)
[gcd(i,j)=1]⇔d∣gcd(i,j)∑μ(d)
void getMu() {
mu[1] = 1;
for (int i = 2; i <= n; ++i) {
if (!flg[i]) p[++tot] = i, mu[i] = -1;
for (int j = 1; j <= tot && i * p[j] <= n; ++j) {
flg[i * p[j]] = 1;
if (i % p[j] == 0) {
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = -mu[i];
}
}
}
莫比乌斯反演有两种形式:
一、
F(n)=
∑
d
∣
n
f
(
d
)
⇒
f
(
n
)
=
∑
d
∣
n
μ
(
d
)
F
(
n
d
)
\sum\limits_{d|n}f(d)\Rightarrow f(n)=\sum\limits_{d|n}\mu(d)F(\frac{n}{d})
d∣n∑f(d)⇒f(n)=d∣n∑μ(d)F(dn)
二、
F(n)=
∑
n
∣
d
f
(
d
)
⇒
f
(
n
)
=
∑
n
∣
d
μ
(
d
n
)
F
(
d
)
\sum\limits_{n|d}f(d)\Rightarrow f(n)=\sum\limits_{n|d}\mu(\frac{d}{n})F(d)
n∣d∑f(d)⇒f(n)=n∣d∑μ(nd)F(d)
例题:
1、Luogu 2303
题意:求
∑
i
=
1
N
g
c
d
(
i
,
N
)
,
(
1
≤
N
≤
2
32
)
\sum\limits_{i=1}^{N}gcd(i,N),(1\leq N\leq2^{32})
i=1∑Ngcd(i,N),(1≤N≤232)
解法:
∑
i
=
1
N
g
c
d
(
i
,
N
)
=
∑
d
∣
N
d
∑
i
=
1
N
[
g
c
d
(
i
,
N
)
=
d
]
=
∑
d
∣
N
d
∑
i
=
1
⌊
N
d
⌋
[
g
c
d
(
i
,
⌊
N
d
⌋
)
=
1
]
=
∑
d
∣
N
d
φ
(
⌊
N
d
⌋
)
\sum\limits_{i=1}^Ngcd(i,N)=\sum\limits_{d|N}d\sum\limits_{i=1}^N[gcd(i,N)=d]=\sum\limits_{d|N}d\sum\limits_{i=1}^{\lfloor \frac{N}{d}\rfloor}[gcd(i,\lfloor \frac{N}{d}\rfloor)=1]=\sum\limits_{d|N}d\varphi(\lfloor \frac{N}{d}\rfloor)
i=1∑Ngcd(i,N)=d∣N∑di=1∑N[gcd(i,N)=d]=d∣N∑di=1∑⌊dN⌋[gcd(i,⌊dN⌋)=1]=d∣N∑dφ(⌊dN⌋)
2、求
∑
i
=
1
n
∑
j
=
1
m
[
g
c
d
(
i
,
j
)
=
1
]
\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[gcd(i,j)=1]
i=1∑nj=1∑m[gcd(i,j)=1]
解法:
∑
i
=
1
n
∑
j
=
1
m
[
g
c
d
(
i
,
j
)
=
1
]
=
∑
i
=
1
n
∑
j
=
1
m
∑
d
∣
g
c
d
(
i
,
j
)
μ
(
d
)
=
∑
d
=
1
m
i
n
(
n
,
m
)
μ
(
d
)
⌊
n
d
⌋
⌊
m
d
⌋
\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[gcd(i,j)=1]=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\sum\limits_{d|gcd(i,j)}\mu(d)=\sum\limits_{d=1}^{min(n,m)}\mu(d)\lfloor \frac{n}{d}\rfloor\lfloor \frac{m}{d}\rfloor
i=1∑nj=1∑m[gcd(i,j)=1]=i=1∑nj=1∑md∣gcd(i,j)∑μ(d)=d=1∑min(n,m)μ(d)⌊dn⌋⌊dm⌋
杜教筛
问题:求积性函数
f
(
i
)
f(i)
f(i)的前缀和
∑
i
=
1
n
f
(
i
)
\sum\limits_{i=1}^nf(i)
i=1∑nf(i)
构造两个积性函数
h
,
g
h,g
h,g,使得
h
=
f
∗
g
h=f*g
h=f∗g
∑
i
=
1
n
h
(
i
)
=
∑
i
=
1
n
∑
d
∣
i
g
(
d
)
f
(
i
d
)
\sum\limits_{i=1}^n h(i)=\sum\limits_{i=1}^n\sum\limits_{d|i}g(d)f(\frac{i}{d})
i=1∑nh(i)=i=1∑nd∣i∑g(d)f(di)
=
∑
d
=
1
n
g
(
d
)
∑
i
=
1
⌊
n
d
⌋
f
(
i
)
\quad\quad\quad=\sum\limits_{d=1}^{n}g(d)\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(i)
=d=1∑ng(d)i=1∑⌊dn⌋f(i)
∑
i
=
1
n
h
(
i
)
=
∑
i
=
1
n
g
(
d
)
S
(
⌊
n
d
⌋
)
\sum\limits_{i=1}^{n}h(i)=\sum\limits_{i=1}^n g(d)S(\lfloor\frac{n}{d}\rfloor)
i=1∑nh(i)=i=1∑ng(d)S(⌊dn⌋)
∑
i
=
1
n
h
(
i
)
=
g
(
1
)
S
(
n
)
+
∑
d
=
2
n
g
(
d
)
S
(
⌊
n
d
⌋
)
\sum\limits_{i=1}^{n}h(i)=g(1)S(n)+\sum\limits_{d=2}^{n}g(d)S(\lfloor\frac{n}{d}\rfloor)
i=1∑nh(i)=g(1)S(n)+d=2∑ng(d)S(⌊dn⌋)
g
(
1
)
S
(
n
)
=
∑
i
=
1
n
h
(
i
)
−
∑
d
=
2
n
g
(
d
)
S
(
⌊
n
d
⌋
)
g(1)S(n)=\sum\limits_{i=1}^{n}h(i)-\sum\limits_{d=2}^{n}g(d)S(\lfloor\frac{n}{d}\rfloor)
g(1)S(n)=i=1∑nh(i)−d=2∑ng(d)S(⌊dn⌋)
构
造
出
的
h
(
i
)
要
使
得
∑
i
=
1
n
h
(
i
)
能
较
快
的
求
出
来
,
构造出的h(i)要使得\sum\limits_{i=1}^{n}h(i)能较快的求出来,
构造出的h(i)要使得i=1∑nh(i)能较快的求出来,
先
用
线
性
筛
把
前
m
项
筛
出
来
,
当
m
=
n
2
3
时
,
复
杂
度
为
O
(
n
2
3
)
;
先用线性筛把前m项筛出来,当m=n^{\frac{2}{3}}时,复杂度为O(n^{\frac{2}{3}});
先用线性筛把前m项筛出来,当m=n32时,复杂度为O(n32);
当
m
=
n
1
2
时
,
复
杂
度
为
O
(
n
3
4
)
;
当m=n^{\frac{1}{2}}时,复杂度为O(n^{\frac{3}{4}});
当m=n21时,复杂度为O(n43);
ll GetSum(int n){//f的前缀和
if(n<=1000000)return sum[n];//线性筛预处理
if(mp[n]!=0)return mp[n];//记忆化
ll ans = f_g_sum(n)%mod; //h的前缀和
for(ll l = 2, r; l <= n; l = r + 1) { //整除分块
r = (n / (n / l));
ans = (ans - (g_sum(r) - g_sum(l - 1)) * GetSum(n / l) + mod)%mod;
//g_sum是g的前缀和,GetSum递归求
}
mp[n]=ans;
return ans;
}
Min_25筛
前提:一个积性函数
f
(
i
)
f(i)
f(i),且
f
(
p
)
f(p)
f(p)是一个关于
p
p
p的项数较少的多项式或可以快速求值;即
f
(
p
c
)
f(p^c)
f(pc)可以快速求值。
用法:在
O
(
n
3
4
log
n
)
O(\frac{n^{\frac{3}{4}}}{\log n})
O(lognn43)的时间复杂度下求
∑
i
=
1
n
f
(
i
)
\sum\limits_{i=1}^nf(i)
i=1∑nf(i)
一些定义:
p
k
p_k
pk:第
k
k
k个质数,特别地,令
p
0
=
1
p_0=1
p0=1。
l
p
f
(
n
)
:
lpf(n):
lpf(n):为
n
n
n的最小质因子,特别地,令
l
p
f
(
1
)
=
1
lpf(1)=1
lpf(1)=1
F
p
r
i
m
e
(
n
)
:
=
∑
2
≤
p
≤
n
f
(
p
)
F_{prime}(n):=\sum_{2\le p\le n}f(p)
Fprime(n):=∑2≤p≤nf(p)
F
k
(
n
)
:
=
∑
i
=
2
n
[
p
k
≤
l
p
f
(
i
)
]
f
(
i
)
F_k(n):=\sum\limits_{i=2}^{n}[p_k\le lpf(i)]f(i)
Fk(n):=i=2∑n[pk≤lpf(i)]f(i)
根据定义:可以发现
∑
i
=
1
n
f
(
i
)
=
F
1
(
n
)
+
f
(
1
)
\sum\limits_{i=1}^nf(i)=F_1(n)+f(1)
i=1∑nf(i)=F1(n)+f(1)
F
k
(
n
)
=
∑
i
=
2
n
[
p
k
≤
l
p
f
(
i
)
]
f
(
i
)
①
=
∑
k
≤
i
,
p
i
≤
n
f
(
p
i
)
+
∑
k
≤
i
,
p
i
2
≤
n
∑
c
≥
1
,
p
i
c
≤
n
f
(
p
i
c
)
(
[
c
>
1
]
+
F
i
+
1
(
n
p
i
c
)
)
②
=
F
p
r
i
m
e
(
n
)
−
F
p
r
i
m
e
(
p
k
−
1
)
+
∑
k
≤
i
,
p
i
2
≤
n
∑
c
≥
1
,
p
i
c
≤
n
f
(
p
i
c
)
(
[
c
>
1
]
+
F
i
+
1
(
n
p
i
c
)
)
③
=
F
p
r
i
m
e
(
n
)
−
F
p
r
i
m
e
(
p
k
−
1
)
+
∑
k
≤
i
,
p
i
2
≤
n
∑
c
≥
1
,
p
i
c
+
1
≤
n
(
f
(
p
i
c
)
F
i
+
1
(
n
p
i
c
)
+
f
(
p
i
c
+
1
)
)
④
\begin{aligned} F_k(n) &=\sum\limits_{i=2}^{n}[p_k\le lpf(i)]f(i)&①\\ &=\sum\limits_{k\le i,p_i\le n}f(p_i)+\sum\limits_{k\le i,p_i^2\le n}\sum\limits_{c\ge 1,p_i^c\le n}f(p_i^c)([c>1]+F_{i+1}(\frac{n}{p_i^c}))& ②\\ &=F_{prime}(n)-F_{prime}(p_{k-1})+\sum\limits_{k\le i,p_i^2\le n}\sum\limits_{c\ge 1,p_i^c\le n}f(p_i^c)([c>1]+F_{i+1}(\frac{n}{p_i^c}))& ③\\ &=F_{prime}(n)-F_{prime}(p_{k-1})+\sum\limits_{k\le i,p_i^2\le n}\sum\limits_{c\ge 1,p_i^{c+1}\le n}(f(p_i^c)F_{i+1}(\frac{n}{p_i^c})+f(p_i^{c+1}))& ④\end{aligned}
Fk(n)=i=2∑n[pk≤lpf(i)]f(i)=k≤i,pi≤n∑f(pi)+k≤i,pi2≤n∑c≥1,pic≤n∑f(pic)([c>1]+Fi+1(picn))=Fprime(n)−Fprime(pk−1)+k≤i,pi2≤n∑c≥1,pic≤n∑f(pic)([c>1]+Fi+1(picn))=Fprime(n)−Fprime(pk−1)+k≤i,pi2≤n∑c≥1,pic+1≤n∑(f(pic)Fi+1(picn)+f(pic+1))①②③④
②
②
②推导:枚举素因子
p
i
p_i
pi,前面部分不解释了,后面部分是
p
i
p_i
pi作为最小素因子时的合数,考虑合数中
p
i
p_i
pi出现了
c
c
c次,
F
i
+
1
(
n
p
i
c
)
F_{i+1}(\frac{n}{p_i^c})
Fi+1(picn)即是剩余部分,因为剩余部分最小的质因子也要大于
p
i
p_i
pi,所以从
p
i
+
1
p_{i+1}
pi+1开始。
③
③
③推导是傻子问题,不解释了。
④
④
④推导:对于满足
p
i
c
≤
n
<
p
i
c
+
1
p_i^c\le n<p_i^{c+1}
pic≤n<pic+1的
c
c
c,有
p
i
c
+
1
>
n
⟺
n
p
i
c
<
p
i
<
p
i
+
1
p_i^{c+1}>n\iff \frac{n}{p_i^c}<p_i<p_{i+1}
pic+1>n⟺picn<pi<pi+1,所以
F
i
+
1
(
n
p
i
c
)
=
0
F_{i+1}(\frac{n}{p_i^c})=0
Fi+1(picn)=0;
然后又有
∑
c
≥
1
,
p
i
c
+
1
≤
n
f
(
p
i
c
+
1
)
=
∑
c
≥
1
,
p
i
c
≤
n
f
(
p
i
c
)
[
c
>
1
]
\sum\limits_{c\ge 1,p_i^{c+1}\le n}f(p_i^{c+1})=\sum\limits_{c\ge 1,p_i^c\le n}f(p_i^c)[c>1]
c≥1,pic+1≤n∑f(pic+1)=c≥1,pic≤n∑f(pic)[c>1],其边界值为
F
k
(
n
)
=
0
(
p
k
>
n
)
F_k(n)=0(p_k>n)
Fk(n)=0(pk>n)。
那假设现在已经求出了所有的
F
p
r
i
m
e
(
n
)
F_{prime}(n)
Fprime(n),那么有两种方式可以求出所有的
F
k
(
n
)
F_k(n)
Fk(n):
1.直接按照递推式计算。(一般数据不是特别大的话用这个)
2.从大到小枚举
p
p
p的转移,仅当
p
2
<
n
p^2<n
p2<n时转移增加值不为0,按照递推式后缀和优化即可。
如何计算
F
p
r
i
m
e
(
n
)
F_{prime}(n)
Fprime(n)呢?
容易发现
F
p
r
i
m
e
F_{prime}
Fprime有且仅有
1
,
2
,
…
,
⌊
n
⌋
,
n
n
,
…
,
n
2
,
n
1,2,…,\lfloor\sqrt{n}\rfloor,\frac{n}{\sqrt{n}},…,\frac{n}{2},n
1,2,…,⌊n⌋,nn,…,2n,n处的值是有用的。
Min_25可用的前提是
f
(
p
)
f(p)
f(p)是一个关于
p
p
p的低次多项式,可以表示为
f
(
p
)
=
∑
a
i
p
c
i
f(p)=\sum a_ip^{c_i}
f(p)=∑aipci,那么每个
p
c
i
p_{c_i}
pci对
F
p
r
i
m
e
(
n
)
F_{prime}(n)
Fprime(n)的贡献是
a
i
∑
2
≤
p
≤
n
p
c
i
a_i\sum_{2\le p\le n}p^{c_i}
ai∑2≤p≤npci。
问题就转变成了:给定
n
,
s
,
g
(
p
)
=
p
s
n,s,g(p)=p^s
n,s,g(p)=ps
构造一个新的函数
G
k
(
n
)
=
∑
i
=
1
n
[
p
k
<
l
p
f
(
i
)
∣
i
s
p
r
i
m
e
(
i
)
]
g
(
i
)
G_k(n)=\sum\limits_{i=1}^{n}[p_k<lpf(i)|isprime(i)]g(i)
Gk(n)=i=1∑n[pk<lpf(i)∣isprime(i)]g(i)。
对于一个合数
x
x
x,必定有
l
p
f
(
x
)
≤
x
lpf(x)\le \sqrt{x}
lpf(x)≤x,则
F
p
r
i
m
e
(
n
)
=
G
⌊
n
⌋
F_{prime}(n)=G_{\lfloor\sqrt{n}\rfloor}
Fprime(n)=G⌊n⌋(n),因为不存在
l
p
f
(
i
)
>
p
⌊
n
⌋
lpf(i)>p_{\lfloor\sqrt{n}\rfloor}
lpf(i)>p⌊n⌋。
考虑G的边界值,有
G
0
(
n
)
=
∑
i
=
2
n
g
(
i
)
,
p
0
=
1
G_0(n)=\sum\limits_{i=2}^{n}g(i),p_0=1
G0(n)=i=2∑ng(i),p0=1。
对于转移,有:
对于
n
<
p
k
2
n<p_k^2
n<pk2的部分,
G
k
(
n
)
=
G
k
−
1
(
n
)
G_k(n)=G_{k-1}(n)
Gk(n)=Gk−1(n)。
对于
p
k
2
≤
n
p_k^2\le n
pk2≤n的部分,被筛掉的数必有质因子
p
k
p_k
pk,即
−
g
(
p
k
)
G
k
−
1
(
n
p
k
)
-g(p_k)G_{k-1}(\frac{n}{p_k})
−g(pk)Gk−1(pkn);但是会有
l
p
f
(
i
)
<
p
k
lpf(i)<p_k
lpf(i)<pk的
i
i
i被减去,应该加回来,即
g
(
p
k
)
G
k
−
1
(
p
k
−
1
)
g(p_k)G_{k-1}(p_{k-1})
g(pk)Gk−1(pk−1)。
所以,有:
G
k
(
n
)
=
G
k
−
1
(
n
)
−
[
p
k
2
≤
n
]
g
(
p
k
)
(
G
k
−
1
(
n
p
k
)
−
G
k
−
1
(
p
k
−
1
)
)
G_k(n)=G_{k-1}(n)-[p_k^2\le n]g(p_k)(G_{k-1}(\frac{n}{p_k})-G_{k-1}(p_{k-1}))
Gk(n)=Gk−1(n)−[pk2≤n]g(pk)(Gk−1(pkn)−Gk−1(pk−1))
代码实现方案:
对于
F
p
r
i
m
e
(
n
)
F_{prime}(n)
Fprime(n),直接按照递推式实现;
对于
p
k
2
≤
n
p_k^2\le n
pk2≤n,可以用线性筛预处理
s
k
=
F
p
r
i
m
e
(
p
k
)
s_k=F_{prime}(p_k)
sk=Fprime(pk)来代替
F
k
F_k
Fk递推式中的
F
p
r
i
m
e
(
p
k
−
1
)
F_{prime}(p_{k-1})
Fprime(pk−1)。
G
G
G递推式中的
G
k
−
1
(
p
k
−
1
)
=
∑
i
=
1
k
−
1
g
(
p
i
)
G_{k-1}(p_{k-1})=\sum_{i=1}^{k-1}g(p_i)
Gk−1(pk−1)=∑i=1k−1g(pi)也可以预处理。
使用Min_25时,首先要明确:
如何线性筛出前
n
\sqrt{n}
n个
f
f
f值;
f
(
p
)
f(p)
f(p)的多项式表示;
如何快速求出
f
(
p
c
)
f(p^c)
f(pc)。
然后按照顺序实现:
1、筛出
[
1
,
n
]
[1,\sqrt{n}]
[1,n]内的质数与前
n
\sqrt{n}
n个
f
f
f值;
2、对
f
(
p
)
f(p)
f(p)多项式表示中的每一项筛出对应的
G
G
G,合并得到
F
p
r
i
m
e
F_{prime}
Fprime的所有
O
(
n
)
O(\sqrt{n})
O(n)个有用点值;
3、按照
F
k
F_k
Fk的递推式实现递归,求出
F
1
(
n
)
F_1(n)
F1(n)。
以下模板以dengtesla大佬的模板进行地修改(完全是因为码风不同,所以改的),大佬模板链接dengtesla大佬的Min_25模板
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int maxn = 2000;
const int N = 710000;
const int mod = 1e9+7;
int b[maxn+1],c[maxn+1][maxn+1],Inv[maxn+1];
ll sqr,n; /// sqr为sqrt(n)
ll w[N],id1[N],id2[N];
///w[i]记录了形如n/k的第i大的取值是多少
int tot; ///记录对于要筛的n,sqrt(n)以内质数的个数
int isp[N],p[N];
ll zh[N][3]; ///zh[i][k]记录(p[1])^k + (p[2])^k + ... + (p[i])^k
ll g[N][3];
///这个k要根据题目改范围,就是这个3要改成多项式的最高次幂
ll qpow(ll a,ll b){
ll ans=1;
a=a%mod;
while(b){
if(b&1)
ans=ans*a%mod;
a=a*a%mod;
b>>=1;
}
return ans;
}
ll sigma_f(ll n,int k){ ///得到∑i^k, i:1~n
if(k==0) return n;
n++;
n%=mod;
ll tmp = n;
ll ans=0;
for (int i=1;i<=k+1;i++){
ans += 1LL*c[k+1][i]*b[k+1-i]%mod*n%mod;
ans %= mod;
n *= tmp%mod;
n %= mod;
}
ans = ((ans*Inv[k+1])%mod+mod)%mod;
return ans;
}
void get_p(int n,int w){
tot = 0;
memset(isp,1,sizeof(isp));
isp[0] = 0;
isp[1] = 0;
for(int i=2;i<=n;i++){
if(isp[i]){
p[++tot] = i;
ll wait = 1;
for(int j=0;j<=w;j++){
zh[tot][j] = zh[tot-1][j] + wait;
zh[tot][j] %= mod;
wait *= i;
}
}
for(int j=1;p[j]*i<=n&&j<=i;j++){
isp[i*p[j]] = 0;
if(i%p[j]==0) break;
}
}
}
void get_g(ll n,int t){
///对每个x=n/i,求出∑[i是质数](i^t) (i from 1 to x)。每个对应的值存储在g[x][t]中
int m = 0;
ll i=1,r;
while(i<=n){
ll len = n/i;
r = n/len;
if(len<=sqr) id1[len] = ++m;
else id2[r] = ++m;
for(int ww=0;ww<=t;ww++){
g[m][ww] = sigma_f(len,ww)-1;
g[m][ww] %= mod;
g[m][ww] += mod;
g[m][ww] %= mod;
}
w[m] = len; ///w[i]记录了形如n/k的第i大的取值是多少
i = r+1;
}
for(int i=1;i<=tot;i++){
for(int j=1;j<=m;j++){
if(1LL*p[i]*p[i]>w[j]) break;
else{
int op;
if(w[j]/p[i]<=sqr) op = id1[w[j]/p[i]];
else op = id2[n/(w[j]/p[i])];
for(int ww=0;ww<=t;ww++){
g[j][ww] = g[j][ww] - qpow(p[i],ww)*((g[op][ww]-zh[i-1][ww])%mod);
g[j][ww] %= mod;
g[j][ww] += mod;
g[j][ww] %= mod;
}
}
}
}
}
inline ll get_value(int wz,int k){
ll w = (g[wz][2]+2*g[wz][1]-g[wz][0])-(zh[k-1][2]+2*zh[k-1][1]-zh[k-1][0]);
w %= mod;
w += mod;
w %= mod;
///比如f(x) = x - 1,
///ll w = (g[wz][1]-g[wz][0])-(zh[k-1][1]-zh[k-1][0]);
return w; ///自己填写f(x)的表达式(在质数时)
///比如f(x) = x^2 + 2*x - 1,
///就写(g[wz][2]+2*g[wz][1]-g[wz][0])-(zh[k-1][2]+2*zh[k-1][1]-zh[k-1][0])
}
ll f(ll p,ll k){ ///计算f(p^k)处的值
if(p==1) return f(1)%mod;//f(1)的值
else if(k==1) return ...%mod;
return ...%mod; ///自己填写
}
ll get_s(ll x,int k){
if(x<=1||p[k]>x) return 0;
int wz;
if(x<=sqr) wz = id1[x];
else wz = id2[n/x];
ll ans = get_value(wz,k);
//if(k==1) ans += 2;
for(int i=k;i<=tot&&1LL*p[i]*p[i]<=x;++i){
for(ll l=p[i],e=1;l*p[i]<=x;l=l*p[i],++e){
ans = ans + (get_s(x/l,i+1)*f(p[i],e)%mod)%mod+f(p[i],e+1);
ans %= mod;
}
}
ans += mod;
ans %= mod;
return ans;
}
void init(){
c[0][0]=1;
for (int i=1;i<maxn;i++){//求组合数
for (int j=1;j<=i;j++)
c[i][j]=(c[i-1][j-1]+c[i-1][j]) % mod;
c[i][0]=1;
}
Inv[1]=1;//线性求逆元
for (int i=2;i<maxn;i++)
Inv[i]=1LL*Inv[mod % i] * (mod-mod/i) % mod;
b[0]=1;
for (int i=1;i<maxn;i++){//伯努利数
b[i]=0;
for (int k=0;k<i;k++)
b[i]=(b[i]+1LL*c[i+1][k]*b[k] % mod) % mod;
b[i]=(1LL*b[i]*(-Inv[i+1]) % mod+mod)%mod;
}
}
void solve(ll n){
sqr = sqrt(n);
get_p(sqr,2);
get_g(n,2);
ll ans = get_s(n,1)+f(1,1);
cout << ans << endl;
}
int main(){
init();
while(cin >> n){
solve(n);
}
}
代码中的一些说明:
求
∑
i
=
1
n
i
k
\sum\limits_{i=1}^{n}i^k
i=1∑nik时用了伯努利数,即生成函数。
有
∑
i
=
1
n
i
k
=
∑
j
=
0
k
(
−
1
)
j
C
k
+
1
j
B
j
n
k
+
1
−
j
k
+
1
\sum\limits_{i=1}^{n}i^k=\frac{\sum\limits_{j=0}^{k}(-1)^jC_{k+1}^{j}B_jn^{k+1-j}}{k+1}
i=1∑nik=k+1j=0∑k(−1)jCk+1jBjnk+1−j,其中
B
j
B_j
Bj是伯努利数。
B
n
=
[
n
=
0
]
−
∑
k
=
0
n
−
1
C
n
+
1
k
B
k
n
+
1
B_n=[n=0]-\sum\limits_{k=0}^{n-1}C_{n+1}^{k}\frac{B_k}{n+1}
Bn=[n=0]−k=0∑n−1Cn+1kn+1Bk,且
B
0
=
1
B_0=1
B0=1。