数论
数论
数论基础
整除
整除的定义 : 设 a,b∈Z,a=0 ,如果 ∃q∈Z ,使得 b=aq, 那么就可以说 b 可被 a 整除, 记作 a∣b ,b 不被 a 整除记作 a∣b。
约数 : 若 a∣b , 则称 b 是 a 的倍数, a 是 b 的约数。
素数与合数
设整数 p=0,±1。如果 p 除了平凡约数外没有其他约数,那么称 p 为素数(不可约数)。
若整数 a=0,±1 且 a 不是素数,则称 a 为合数。
所有大于 3 的素数都可以表示为 6n±1 的形式。
算术基本定理
设正整数 a , 那么必有表示 :
a=p1k1p2k2⋯psks,p1<p2<⋯<ps
其中 pi(1≤i≤s) 为素数。
同余
同余的定义 : 设整数 m=0 。 若 m∣(a−b), 称 m 为模数, a 同余于 b 模 m, b 是 a 对模 m 的剩余。 记作 a≡b(modm) 。
同余有自反性, 对称性, 传递性的性质, 可以进行线性运算,
数论函数
数论函数指定义域为正整数的函数。数论函数也可以视作一个数列。
积性函数
若函数 f(n) 满足 f(1)=1 且 ∀x,y∈N∗,gcd(x,y)=1 都有 f(x)(y)=f(x)f(y) , 则 f(n) 为积性函数。
若函数 f(n) 满足 f(1)=1 且 ∀x,y∈N∗ 都有 f(x)(y)=f(x)f(y) , 则 f(n) 为完全积性函数。
若 f(x) 和 g(x) 均为积性函数, 则以下函数也为积性函数 :
h(x)=f(xp)h(x)=fp(x)h(x)=f(x)g(x)h(x)=d∣x∑f(d)g(dx)
比较常见的积性函数 :
d(x)=i∣n∑1σ(x)=i∣n∑iφ(x)=i=1∑x1[gcd(x,i)=1]μ(x)=⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧1(−1)k0x=1i=1∏kqi=1max{qi}>1ϵ(x)=[x=1]
素数
素数计数函数, 将 π(x) 记为 小于或等于 x 的素数的个数, 随着 x 的增大,有 π(x)∼lnxx。
费马小定理
定义 : 若 p 为 素数 ,gcd(a,p)=1 , 则 ap−1≡1(modp)。
证明 : 设一个质数 p, 然后取一个与 p 互质的数 a 。
首先构造一个序列 A=1,2,3,⋯,p−1, 有 :
i=1∏p−1Ai≡i=1∏p−1(Ai×a)(modp)
首先,
∵gcd(Ai,p)=1,gcd(Ai×a,p)=1
可以得到每一个 Ai×a 对应一个 Ai, 所以该性质成立。
所以进一步可得 :
i=1∏p−1Ai≡ap−1×i=1∏p−1Ai
即 :
ap−1≡1(modp)
Miller Rabin
Miller Rabin 是进阶的概率性素性测试, 主要是 Fermat 的优化。
朴素的 Fermat 只使用了费马小定理来检验 :
1 2 3 4 5 6 7 8 9 10 11
| bool fermat(int n) { if(n < 3)return n == 2; for(int i = 1; i <= T; i++) { int a = rand() % (n - 2) + 2; if(qpow(a, n - 1, n) != 1) return 0; } return 1; }
|
但是仍有满足费马小定理的合数, 也就是伪素数,这样会被特定的数据卡掉。
费马小定理的逆定理是不成立的。
上述的费马伪素数,也称为卡迈克尔数, 即对于合数 n ,如果对于所有正整数 a , a 和 n 互质, 都有同余式 an−1≡1(modn) 成立,则合数 n 为卡迈克尔数。
并且, 若 n 为卡迈克尔数, 则 m=2n−1 也是一个卡迈克尔数。
Miller Rabin 则使用了更多确定性算法, 来进行素性测试,比如二次探测定理。
二次探测定理 : 对于奇素数 p , 若 x2≡1(modp) , 则小于 p 的解只有两个, x1=1,x2=p−1。
将式子移项即可得证。
然后我们将其用在素性测试中,对于 an−1 ,可以拆分为 n−1=u×2t, 先求出 au ,然后依次做 t 次平方, 同时检验,当前值是否满足二次探测定理,如果出现不是 1 或 n−1 的解, 说明没有通过素性测试。
最后 Miller Rabin 如下
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
| bool miller_rabin(ll p) { if(p < 2)return 0; if(p == 2 || p == 3)return 1; ll d = p - 1, r = 0; while(!(d & 1))r++, d >>= 1; for(int k = 1; k <= 10; k++) { ll a = rand() % (p - 2) + 2; ll x = qpow(a, d, p); if(x == 1 || x == p - 1)continue; for(int i = 1; i <= r; i++) { x = (i128)x * x % p; if(x == p - 1)break; } if(x != p - 1)return 0; } return 1; }
|
最大公约数
一组整数的公约数,是指同时是这组数中每一个数的约数的数。±1 是任意一组整数的公约数。
一组整数的最大公约数,是指所有公约数里面最大的一个。
欧几里得算法
对于两个数 a,b ,有 gcd(a,b)=gcd(b,amodb)。
证明 :
假设 a=kb+c, 对于 a,b 的公因数 d ,因为 d∣a,d∣b, 所以 d∣c, 所以 a,b 的公因数均为 c 的约数, 反之 b,c 的公约数同时也是 a 的约数, 所以公约数集合相同,得证。
即可递归实现。
1 2 3 4 5
| int gcd(int a, int b) { if(!a || !b)return a + b; else return gcd(b, a % b); }
|
最小公倍数
对于两个数 a,b, 有 lcm(a,b)×gcd(a,b)=a×b, 由算术基本定理可得。
更相减损法
对于两个数 a,b, 有 gcd(a,b)=gcd(b,a−b), 证明同欧几里得算法。
拓展欧几里得算法
拓展欧几里得算法用于求解形如 ax+by=gcd(a,b) 方程的一组可行解,其推导过程如下,
∵gcd(a,b)=gcd(b,amodb)∴ax+by=gcd(a,b)=gcd(b,amodb)
考虑递归下去,当 b=0 时, 此时解为 x=1,y=0,然后对于当前解与上次迭代的关系为,
bx+(amodb)y=gcd(b,amodb)bx+(a−⌊ba⌋×b)y=gcd(a,b)
整理得,
ay+b(x−⌊ba⌋×y)=gcd(a,b)
递归求解即可,
1 2 3 4 5 6
| void exgcd(int a, int b, int &x, int &y) { if(!b){x = 1, y = 0; return;} exgcd(b, a % b, x, y); int z = x; x = y; y = z - a / b * y; }
|
求出的 x0,y0 为一组特解,通解为 x=x0+kgcd(a,b)b,y=y0−kgcd(a,b)a 。
数论分块
求解 ∑i=1n⌊in⌋, 复杂度 O(n),对于连续的一段 i, 其值可能是相同的,假设左边界为 l , 可得右边界为 ⌊⌊ln⌋n⌋ 。
证明 :
设 x=⌊in⌋,有 x≤in
所以对于每个 i 和其对应的 in, ⌊xn⌋≥⌊inn⌋, 即 i≤⌊⌊in⌋n⌋, 得证。
对于求解 ∑i=1n⌊in⌋
1 2 3 4 5
| for(int l = 1, r; l <= n; l = r + 1) { r = n / (n / l); sum = sum + (r - l + 1) * (n / l); }
|
欧拉函数
定义欧拉函数 φ(n) 表示 1 到 n 中与 n 互质的数的个数。
性质
假设将 n 唯一分解为 ∏i=1mpiki, 有 φ(n)=n∏i=1m(1−pi1)。
证明 :
n 个数中是 pi 倍数的个数有 pin 个, n−∑i=1mpin 后, 会减去重复的部分,还要加上,使用容斥原理即可得证。
欧拉函数为积性函数。
对于任意正整数 n, 有 n=∑d∣nφ(d) 。
证明 :
设 f(n)=∑d∣nφ(d), m⊥n,根据欧拉函数的积性函数的性质有:
f(n)×f(m)=i∣n∑φ(i)j∣m∑φ(j)=i∣n∑j∣m∑φ(i)φ(j)=i∣n∑j∣m∑φ(ij)=d∣nm∑φ(d)=f(nm)
将 n 质因数分解 p1c1×p2c2×p3c3×⋯×pkck。
所以有,
f(n)=f(p1c1)×f(p2c2)×⋯×f(pkck)
其中的每一项有:
f(pc)=i=0∑cφ(pi)=pc
所以得证:
f(n)=f(p1c1)×f(p2c2)×⋯×f(pkck)=Πi=1kpici=n=d∣n∑φ(d)
简单筛法
埃拉托斯特尼筛法
假如我们要筛出 小于或等于 n 的所有质数, 单纯做 n 次素性测试是达不到最优复杂度的。
首先考虑这样一种筛法, 对于任意一个大于 1 的正整数,其倍数一定是合数,然后枚举即可。
1 2 3 4 5 6
| for(int i = 2; i <= n; i++) { if(!vis[i])prime[++cnt] = i; for(int j = i; i * j <= n; j++) vis[i * j] = 1; }
|
时间复杂度为 O(nlnn)。
欧拉筛
埃氏筛的过程中有的合数被标记了不止一次, 还可以进一步优化,让每个数只被其最小质因子标记一次,也就是欧拉筛。
1 2 3 4 5 6 7 8 9
| for(int i = 2; i <= n; i++) { if(!vis[i])prime[++cnt] = i; for(int j = 1; j <= cnt && prime[j] * i <= n; j++) { vis[prime[j] * i] = 1; if(i % prime[j] == 0)break; } }
|
线性筛也可以用来筛积性函数,
筛欧拉函数,
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
| phi[1] = 1; for(int i = 2; i <= n; i++) { if(!vis[i]) { prime[++cnt] = i; phi[i] = i - 1; } for(int j = 1; j <= cnt && prime[j] * i <= n; 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] * phi[prime[j]]; } }
|
筛莫比乌斯函数,
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
| mu[1] = 1; for(int i = 2; i <= n; i++) { if(!vis[i]) { prime[++cnt] = i; mu[i] = -1; } for(int j = 1; j <= cnt && prime[j] * i <= n; j++) { vis[i * prime[j]] = 1; if(i % prime[j] == 0) { mu[i * prime[j]] = 0; break; } mu[i * prime[j]] = -mu[i]; } }
|
筛约数个数函数
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
| d[1] = 1; for(int i = 2; i <= n; i++) { if(!vis[i]) { prime[++cnt] = i; d[i] = 2; c[i] = 1; } for(int j = 1; j <= cnt && prime[j] * i <= n; j++) { vis[i * prime[j]] = 1; if(i % prime[j] == 0) { c[i * prime[j]] = c[i] + 1; d[i * prime[j]] = d[i] + d[i] / (c[i] + 1); break; } c[i * prime[j]] = 1; d[i * prime[j]] = d[i] * d[prime[j]]; } }
|
Pollard Rho
Pollard Rho 可以用来快速分解非平凡因数, 相比朴素的 O(n) 会快很多,但是是一个随机性算法, 实际运行效率不定。
随机选一个数来判断其是不是 n 因数, 显然概率低, 是不可行的,但是我们可以利用生日悖论,采用组合随机采样的方法,来提高正确率。
首先考虑最大公约数一定是某个数的约数, 也就是说 ∀k∈N∗,gcd(k,n)∣n ,随机选取一个 k 使得 gcd(k,n)>1 就可以求得一个约数。
比如选取一个序列 xi, 若有 gcd(∣xi−xj∣,n)>1, 则求得一个约数,但是这样比较的话,是 O(n2) 的,还需要降低时间复杂度。
我们可以构造一伪随机函数, f(x)=(x2+c)modn 来生成一个随机数序列, 其中 c=rand(), 为一个随机数。然后再随机取一个 x0 , 令 xi=f(xi−1), 生成随机序列, 取相邻两项求 gcd。
但是这个序列会出现循环节,需要进一步优化, 也就是需要找出循环节。
找循环节可以用一个简单的 Floyd 判圈算法, 也就是类似 ”龟兔赛跑“, 一个以另一个两倍的速度去遍历, 相遇即为找到环。
但是由于求 gcd 的时间复杂度为 O(logn) 的频繁调用的话会导致算法变慢,此时可以乘法累计减少求 gcd 的次数,显然如果 gcd(a,b)>1 , 则有 gcd(ac,b)>1, 又因为 gcd(ac,b)=gcd(acmodb,b), 所以我们可以把所有测试样本在模 n 的意义下乘起来, 然后使用倍增优化,每次在路径上取一段区间 [l,r], 每次取测试样本为 ∣xi−xl∣ 这样每次累计的样本个数就为 r−l+1 个,减少了取 gcd 的次数。 样本长度建议选为 127 。
下面是找一个数的最大质因子的方法,
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92
| #include <bits/stdc++.h>
using namespace std;
#define ll long long #define i128 __int128_t
ll n, Max;
ll gcd(ll a, ll b) { if(!a || !b)return a + b; else return gcd(b, a % b); }
ll qpow(ll a, ll b, ll mod) { ll t = 1; while(b != 0) { if(b & 1)t = (i128)t * a % mod; a = (i128)a * a % mod; b >>= 1; } return t; }
bool miller_rabin(ll p) { if(p < 2)return 0; if(p == 2 || p == 3)return 1; ll d = p - 1, r = 0; while(!(d & 1))r++, d >>= 1; for(int k = 1; k <= 10; k++) { ll a = rand() % (p - 2) + 2; ll x = qpow(a, d, p); if(x == 1 || x == p - 1)continue; for(int i = 1; i <= r; i++) { x = (i128)x * x % p; if(x == p - 1)break; } if(x != p - 1)return 0; } return 1; }
ll pollard_rho(ll x) { ll s = 0, t = 0, c = (ll)rand() % (x - 1) + 1; int step = 0, ed = 1; ll v = 1; for(ed = 1; ; ed <<= 1, s = t, v = 1) { for(step = 1; step <= ed; step++) { t = ((i128)t * t + c) % x; v = (i128)v * abs(t - s) % x; if(step % 127 == 0) { ll d = gcd(v, x); if(d > 1)return d; } } ll d = gcd(v, x); if(d > 1)return d; } }
void get_fact(ll x) { if(x <= Max || x < 2)return; if(miller_rabin(x)) return void(Max = max(Max, x)); ll p = x; while(p >= x)p = pollard_rho(x); while(x % p == 0)x /= p; get_fact(x); get_fact(p); }
int main() { srand(time(0)); int T; scanf("%d", &T); while(T--) { scanf("%lld", &n); Max = 0; get_fact(n); if(Max == n)printf("Prime\n"); else printf("%lld\n", Max); } return 0; }
|
裴蜀定理
设 a,b 是不全为零的整数, 则存在整数 x,y 使得 ax+by=gcd(a,b)。
对于方程 ax+by=c 若 gcd(a,b)∣c 则方程无整数解。
欧拉定理
若 gcd(a,m)=1 , 则 aφ(m)≡1(modm)
证明 :
跟费马小定理相似的证法, 构造一个与 m 互质的数列, 为 A1,A2,⋯Aφ(m),1≤Ai≤m, 然后选取一个与 m 互质的正整数 a 有,
i=1∏φ(m)Ai≡i=1∏φ(m)(Ai×a)(modp)
不难得出其是一一对应的, aφ(m)≡1(modm) 得证。
扩展欧拉定理 :
ab≡⎩⎪⎪⎪⎨⎪⎪⎪⎧abmodφ(m)ababmodφ(m)+φ(m)gcd(a,m)=1gcd(a,m)=1,b≤φ(m),(modm)gcd(a,m)=1,b≥φ(m)
乘法逆元
对于一个线性同余方程 ax≡1(modb) , 则 x 称为 amodb 的逆元, 记作 a−1 。
费马小定理
若求 a 在模 p 下的逆元,且 p 为质数,可以使用费马小定理。
∵ap−1≡1(modp)∴ap−2≡a1(modp)
扩展欧几里得算法
当 p 不为质数时, 就不可以使用费马小定理, 假设 a,m 互质。
∵ax≡1(modm)∴ax−my=1
直接用扩展欧几里得方程解出的解即为乘法逆元。
线性求逆元
上述方法单次求逆元复杂度都是 O(logn) , 存在线性递推的做法。
线性递推求逆元要求 p 为质数。
设 p=ki+r(r<i,1<i<p), 有,
ki+r≡0(modp)
两边同时乘 i−1,r−1, 得,
i−1r−1(ki+r)≡0(modp)
展开得,
i−1+kr−1≡0(modp)
即,
i−1≡−kr−1(modp)
根据设得 p 可以得知
i−1≡−⌊ip⌋×(pmodi)−1(modp)
又等价于
i−1≡(p−⌊ip⌋)×(pmodi)−1(modp)
此时即可递推, inv[1]=1 。
中国剩余定理
中国剩余定理可以用来求解如下形式得一元线性同余方程组。其中 n1,n2,⋯,nk 两两互质。
⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧x≡a1(modm)1x≡a2(modm)2⋮x≡an(modm)n
主要是用构造的方式,设 M=∏i=1nmi,Mi=miM , 解即为 :
i=1∑naiMiMi−1(modM)
其中 Mi−1 为模 mi 意义下的逆元。
扩展中国剩余定理
当模数不两两互质时,假设有两个方程为 x≡a1(modm)1,x≡a2(modm)2。
将其转化为不定方程为 x=m1p+a1=m2q+a2, 则有 m1p−m2q=a2−a1。
然后用扩展欧几里得算法求出一组可行解 p,q。
构成一个新的同余方程为 x≡b(modM), 其中 b=m1p+a1,M=lcm(m1,m2)。
威尔逊定理
对于素数 p 有 (p−1)!≡−1(modp)。
证明 :
首先 p=2 时显然成立,然后讨论奇素数的情况, 此时 1,2,⋯,p−1 均存在逆元, 且一一对应,那么只需要将一个数与其逆元一一对应即可,那么显然 (p−1)!modp=p−1,即 (p−1)!≡−1(modp) 。
BSGS
BSGS(baby-step giant-step) 是用来解决形如
ax≡b(modm)
的高次同余方程,其中 a⊥m
首先0≤x<m(就m个不同的得数),设x=Am−B,其中0≤A,B≤m,方程变为:
aAm−B≡b(modm)
由于 a⊥m移项可得:
aAm≡baB(modm)
由于a,b,m均已知可以考虑枚举B用map存储对应的值,然后再枚举A,检查map中是否有值对应
因为A,B均不大于m,所以时间复杂度为O(n)。(没算map)
exBSGS
问题仍然是去解决一个形如
ax≡b(modm)
的高次同余方程,但不保证a⊥m
考虑用扩展BSGS(exBSGS) 来解决这个问题
此时a和m不一定互质,我们可以将式子进一步化简为:
gcd(a,m)ax≡gcd(a,m)b(modgcd(a,m)m)
将左边提出一项
ax−1gcd(a,m)a≡gcd(a,m)b(modgcd(a,m)m)
接下来不断除以gcd(a,m)直到a与m互质,然后用BSGS求解
由于最后我们将式子递归执行了k层,此时式子变为
ax−kgcdk(a,m)ak≡gcdk(a,m)b(modgcdk(a,m)m)
最后利用BSGS求解出x后应当再加上k,x+k即为该方程的解。
卢卡斯定理
卢卡斯定理用于求解大组合数取模的问题,其中模数要求为质数。
对于质数 p ,有
(mn)modp=(⌊pm⌋⌊pn⌋)×(mmodpnmodp)modp
可以直接递归求解。
证明 :
考虑
(np)modp 的取值,注意到
(np)=n!(p−n)!p! ,分子的质因子分解中 p 的次数恰好为 1 ,因此只有当 n=0 或 n=p 的时候 n!(p−n)! 的质因子分解中含有 p ,因此
(np)modp=[n=0∨n=p] 。进而我们可以得出
(a+b)p=n=0∑p(np)anbp−n≡n=0∑p[n=0∨n=p]anbp−n≡ap+bp(modp)
注意过程中没有用到费马小定理,因此这一推导不仅适用于整数,亦适用于多项式。因此我们可以考虑二项式 fp(x)=(axn+bxm)pmodp 的结果
(axn+bxm)p≡apxpn+bpxpm≡axpn+bxpm≡f(xp)
考虑二项式 (1+x)nmodp ,那么
(mn) 就是求其在 xm 次项的取值。使用上述引理,我们可以得到
(1+x)n≡(1+x)p⌊pn⌋(1+x)nmodp≡(1+xp)⌊pn⌋(1+x)nmodp
注意前者只有在 p 的倍数位置才有取值,而后者最高次项为 nmodp≤p−1,因此这两部分的卷积在任何一个位置只有最多一种方式贡献取值,即在前者部分取 p 的倍数次项,后者部分取剩余项,即
(mn)modp=(⌊pm⌋⌊pn⌋)×(mmodpnmodp)modp
扩展卢卡斯定理
对于 p 不为质数的情况,就需要用到扩展卢卡斯定理。
对于 (mn)modM, 我们先将 M 唯一分解, M=∏i=1rpiki, 然后构造一个同余方程组。
⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧a1≡(mn)(modp1k1)a2≡(mn)(modp2k2)⋮ar≡(mn)(modprkr)
先求每个方程的解,然后用中国剩余定理合并即可。
对于求解 a≡(mn)(modpk),
∵(mn)=m!(n−m)!n!∴(mn)modpk=m!(n−m)!n!modpk
然后考虑分母如何计算,由于分母分子和模数不一定互质,考虑提出 n! 中的 p 的倍数。
n!=p⌊pn⌋×⌊pn⌋!×p∣i∏ni
前面的直接快速幂,⌊pn⌋! 直接递归计算。
后面不整除的数显然会出现循环节,也就是
p∣i∏ni=(i=1,p∣i∏pki)⌊pkn⌋(i=pkpkn,p∣i∏ni)
分成了循环节和余项。
递归处理把 p 的倍数处理完即可。
式子其实就是,
pym!pz(n−m)!pxn!px−y−zmodpk
对于 x−y−z, 其实就是 ∑i=1⌊pin⌋。
类欧几里得算法
类欧几里德算法是由洪华敦在 2016 年冬令营营员交流中提出的内容,其本质可以理解为,使用一个类似辗转相除法来做函数求和的过程。
问题1
设
f(a,b,c,n)=i=0∑n⌊cai+b⌋(1,1)
其中 a,b,c,n∈Z 。
一眼不可做,似乎像数论分块,但又好像不太行的样子。
a,b,c,n 之间的关系也没给,那就分情况吧
首先考虑第一种情况 a≥c,b≥c,(也是最简单的一种情况)
此时式子可以进一步化简为 :
f(a,b,c,n)=i=0∑n⌊cai+b⌋=i=0∑n⌊c(⌊ca⌋×c+amodc)i+⌊cb⌋×c+bmodc⌋=i=0∑n(⌊ca⌋i+⌊cb⌋+⌊c(amodc)i+bmodc⌋)=i=0∑n⌊ca⌋i+i=0∑n⌊cb⌋+i=0∑n⌊c(amodc)i+bmodc⌋=2n(n+1)⌊ca⌋+(n+1)⌊cb⌋+f(amodc,bmodc,c,n)(1,2)
此时我们把式子化为了一个含有 f 的递归式,接下来 a,b 必定会变得小于 c,还需讨论另一种情况
a<c,b<c 时,这时候也不能用上面的方法,🤔
利用枚举贡献的方法,将式子化为 :
f(a,b,c,n)=i=0∑nj=0∑⌊cai+b⌋−11(1,3)
似乎什么都没有做的样子。。
需要进一步化简:
i=0∑nj=0∑⌊cai+b⌋−11=i=0∑nj=0∑⌊can+b⌋−1[j<⌊cai+b⌋](1,4)
改了第二层求和的上界,好像有点头绪了
现在限制我们的就是后面的条件表达式了,也许可以将转化为某种方便计算的形式:
∵j<⌊cai+b⌋∴j≤cai+b−1∴jc≤ai+b−c∴jc<ai+b−c+1∴i>ajc+c−b−1∴i>⌊ajc+c−b−1⌋(1,5)
这样的话原来的式子就变为了:
f(a,b,c,n)=i=0∑nj=0∑⌊can+b⌋−1[i>⌊ajc+c−b−1⌋]=j=0∑⌊can+b⌋−1i=0∑n[i>⌊ajc+c−b−1⌋]=j=0∑⌊can+b⌋−1i=⌊ajc+c−b−1⌋+1∑n1=j=0∑⌊can+b⌋−1(n−⌊ajc+c−b−1⌋)=j=0∑⌊can+b⌋−1n−j=0∑⌊can+b⌋−1⌊ajc+c−b−1⌋=⌊can+b⌋n−f(c,c−b−1,a,⌊can+b⌋−1)(1,6)
将 ⌊can+b⌋ 设为 m 再整理一下,又得到一个递归式:
f(a,b,c,n)=mn−f(c,c−b−1,a,m−1)(1,7)
综上(1,2)(1,6):
f(a,b,c,n)=⎩⎪⎪⎪⎨⎪⎪⎪⎧2n(n+1)⌊ca⌋+(n+1)⌊cb⌋+f(amodc,bmodc,c,n)⌊can+b⌋n−f(c,c−b−1,a,⌊can+b⌋−1)a≥c∨b≥ca<c∧b<c(1,8)
可以在 O(logn) 时间内求出。
问题2
设
g(a,b,c,n)=i=0∑ni⌊cai+b⌋(2,1)
其中 a,b,c,n∈Z。
上个问题的升级版, 接着分情况,a≥c,b≥c 时:
g(a,b,c,n)=i=0∑ni⌊cai+b⌋=i=0∑ni⌊c(⌊ca⌋×c+amodc)i+⌊cb⌋×c+bmodc⌋=i=0∑ni(⌊ca⌋i+⌊cb⌋+⌊c(amodc)i+bmodc⌋)=i=0∑n⌊ca⌋i2+i=0∑n⌊cb⌋i+i=0∑n⌊c(amodc)i+bmodc⌋i=6n(n+1)(2n+1)⌊ca⌋+2n(n+1)⌊cb⌋+g(amodc,bmodc,c,n)(2,2)
看来我们完成了一半了(似乎不是)
a<c,b<c 时,枚举贡献将式子化为:
g(a,b,c,n)=i=0∑nj=0∑⌊cai+b⌋−1i(2,3)
改变上界:
i=0∑nj=0∑⌊cai+b⌋−1i=i=0∑nj=0∑⌊can+b⌋−1i[j<⌊cai+b⌋](2,4)
转化条件表达式:
g(a,b,c,n)=i=0∑nj=0∑⌊can+b⌋−1i[i>⌊ajc+c−b−1⌋]=j=0∑⌊can+b⌋−1i=0∑ni[i>⌊ajc+c−b−1⌋]=j=0∑⌊can+b⌋−1i=⌊ajc+c−b−1⌋+1∑ni(2,5)
设 k=⌊ajc+c−b−1⌋,m=⌊can+b⌋,接着化简:
g(a,b,c,n)=i=0∑m−1k+1∑ni=i=0∑m−12(n−k)(n+k+1)=i=0∑m−12n2+n−k2−k=21(mn2+mn−j=0∑m−1(⌊ajc+c−b−1⌋)2−i=0∑m−1⌊ajc+c−b−1⌋)(2,6)
。。。好像化的没问题啊,出现了一个没见过的平方,先定义为 h(a,b,c,n) 吧,整理一下得:
g(a,b,c,n)=⎩⎪⎪⎪⎨⎪⎪⎪⎧6n(n+1)(2n+1)⌊ca⌋+2n(n+1)⌊cb⌋+g(amodc,bmodc,c,n)21(⌊can+b⌋n(n+1)−h(c,c−b−1,a,⌊can+b⌋−1)−f(c,c−b−1,a,⌊can+b⌋−1))a≥c∨b≥ca<c∧b<c(2,7)
接下来解决带平方的问题。
问题3
设
h(a,b,c,n)=i=0∑n(⌊cai+b⌋)2(3,1)
其中 a,b,c,n∈Z。
上个问题还没有完全解决,需要解出 h(a,b,c,n),(不就多了个平方嘛,轻车熟路了,别打了别打了)
h(a,b,c,n)=i=0∑n(⌊cai+b⌋)2=i=0∑n(⌊c(⌊ca⌋×c+amodc)i+⌊cb⌋×c+bmodc⌋)2=i=0∑n(⌊ca⌋i+⌊cb⌋+⌊c(amodc)i+bmodc⌋)2(3,2)
手动多项式乘法!
原式变为:
h(a,b,c,n)=h(amodc,bmodc,c,n)+2⌊cb⌋f(amodc,bmodc,c,n)+2⌊ca⌋g(amodc,bmodc,c,n)+⌊ca⌋26n(n+1)(2n+1)+⌊cb⌋2(n+1)+⌊ca⌋⌊cb⌋n(n+1)(3,3)
接着考虑 a<c,b<c 的情况,设 k=⌊ajc+c−b−1⌋,m=⌊can+b⌋,
平方可以拆成:
n2=22n(n+1)−n=(2i=0∑ni)−n(3,3)
避免了出现 ∑×∑,数学技巧
h(a,b,c,n)=i=0∑n⌊cai+b⌋2=i=0∑n⎣⎢⎢⎡⎝⎜⎜⎛2j=1∑⌊cai+b⌋j⎠⎟⎟⎞−⌊cai+b⌋⎦⎥⎥⎤=⎝⎜⎜⎛2i=0∑nj=1∑⌊cai+b⌋j⎠⎟⎟⎞−f(a,b,c,n)(3,4)
接着化简前面的部分:
i=0∑nj=1∑⌊cai+b⌋j=i=0∑nj=0∑⌊cai+b⌋−1(j+1)=j=0∑m−1(j+1)i=0∑n[j<⌊cai+b⌋]=j=0∑m−1(j+1)i=0∑n[i>k]=j=0∑m−1(j+1)(n−k)=21nm(m+1)−j=0∑m−1(j+1)⌊ajc+c−b−1⌋=21nm(m+1)−g(c,c−b−1,a,m−1)−f(c,c−b−1,a,m−1)(3,5)
整理一下:
h(a,b,c,n)=⎩⎪⎪⎪⎨⎪⎪⎪⎧h(amodc,bmodc,c,n)+2⌊cb⌋f(amodc,bmodc,c,n)+2⌊ca⌋g(amodc,bmodc,c,n)+⌊ca⌋26n(n+1)(2n+1)+⌊cb⌋2(n+1)+⌊ca⌋⌊cb⌋n(n+1)n⌊can+b⌋(⌊can+b⌋+1)−2g(c,c−b−1,a,⌊can+b⌋−1)−2f(c,c−b−1,a,⌊can+b⌋−1)−f(a,b,c,n)a≥c∨b≥ca<c∧b<c(3,6)
时间复杂度均为 O(logn)!
狄利克雷卷积
狄利克雷卷积是定义在数论函数间的二元运算。
数论函数,是指定义域为 N ,值域为 C , 的一类函数,每个数论函数可以视为一个复数的序列。
定义式为:
(f∗g)(n)=d∣n∑f(d)g(dn)(d∈N)
也可以写为 :
(f∗g)(n)=d∣n∑f(dn)g(d)(d∈N)
同时由于对称性,又可以写为 :
(f∗g)(n)=xy=n∑f(x)g(y)(x,y∈N)
“+” 定义为了数论函数直接的直接相加, " ∗ ",其实就是乘。
常见的数论函数
单位函数
ε(n)=[n=1]
幂函数
Idk(n)=nk
约数函数
σ(n)=d∣n∑dk(d∈N)
欧拉函数
φ(n) 表示 1 n 中与 n 互质的整数的个数有多少个
φ(n)=np∣n∏(1−p1)(p∈P)
(容斥可证)
有趣的是上面所写的函数均为积性函数,其中单位函数和幂函数为完全积性函数。
相关定理
- 若 f , g 为积性函数,则 f∗g 也为积性函数
证明:
设gcd(a,b)=1∴(f∗g)(a)⋅(f∗g)(b)=i∣a∑f(i)g(ia)⋅j∣b∑f(j)g(jb)=i∣a∑j∣b∑f(i)g(ia)⋅f(j)g(jb)=d∣ab∑f(d)g(dab)=(f∗g)(ab)
用到了积性函数的性质
- f∗g=g∗f
证明 :
(f∗g)(n)=ij=n∑f(i)g(j)=ji=n∑f(j)g(i)=(g∗j)(n)
对称性。
- (f∗g)∗h=f∗(g∗h)
利用对称性的式子去拆开即可得证。
- f∗(g+h)=f∗g+f∗h
证明:
(f∗(g+h))(n)=ij=n∑f(i)(g+h)(j)=ij=n∑f(i)[g(j)+h(j)]=ij=n∑f(i)g(j)+f(i)h(j)=ij=n∑f(i)g(j)+ij=n∑f(i)h(j)=(f∗g+f∗h)(n)
特殊的卷积
Idk∗1=σk
证明:
(Idk∗1)(n)=d∣n∑Idk(d)1(dn)=d∣n∑Idk(d)=d∣n∑dk=σ(n)
同时可以得到一个应用广泛的式子 :
(f∗1)(n)=d∣n∑f(d)
φ∗1=Id
证明 :
∵(φ∗1)(n)=d∣n∑φ(d)将n拆分为∏pk∴(φ∗1)(n)=(φ∗1)(∏pk)=∏(φ∗1)(piki)=∏j=1∑kiφ(pij)=∏piki=n∴φ∗1=Id
1∗1=d
证明 :
(1∗1)(n)=d∣n∑1(d)1(dn)=d∣n∑1=d(n)
上述的运算加以结合可以得到更多结论 。
狄利克雷卷积逆
需要用到单位元,有 :
(f∗ε)(n)=d∣n∑ε(d)f(dn)=f(n)
定义 :
f∗f−1=ϵ
即为狄利克雷卷积逆
积性函数一定有狄利克雷逆,且它也是积性函数。
莫比乌斯反演
定义莫比乌斯函数为:
μ(x)=⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧1(−1)k0x=1i=1∏kqi=1max{qi}>1
推导
g=f∗1⟺f=g∗μ
也就是 :
f(n)=d∣n∑g(d)⟺g(n)=d∣n∑μ(d)f(dn)
莫比乌斯函数的性质:
d∣n∑nμ(d)=[n=1]
接下来证明莫比乌斯反演定理 :
f(n)=d∣n∑g(d)=d∣n∑g(dn)d∣n∑μ(d)f(dn)=d∣n∑μ(d)d1∣dn∑g(d1)d∣n∑d1∣dn∑μ(d)g(d1)=d1∣n∑d∣d1n∑μ(d)g(d1)=d1∣n∑g(d1)d∣d1n∑μ(d)=g(n)
上述仅为充分性证明,必要性证明逆推即可。
杜教筛
对于求一个数论函数的前缀和,杜教筛可以在低于线性时间的复杂度内求解。
对于数论函数 f , 要求计算 S(n)=∑i=1nf(i) 。
首先构造一个 S(n) 关于 S(⌊in⌋) 的递推式
对于任意一个数论函数 g ,必须满足
i=1∑nd∣i∑g(d)f(di)=i=1∑ng(i)S(⌊in⌋)⟺i=1∑n(f∗g)(i)=i=1∑ng(i)S(⌊in⌋)
简单的证明:
i=1∑nd∣i∑g(d)f(di)=i=1∑nj=1∑⌊in⌋g(i)f(j)=i=1∑ng(i)j=1∑⌊in⌋f(j)=i=1∑ng(i)S(⌊in⌋)
求出递推式
g(1)S(n)=i=1∑n(f∗g)(i)−i=2∑ng(i)S(⌊in⌋)
可以用数论分块对后半部分快速求出结果。
对于一些筛法的小技巧:
对于比较大的数据时,筛法在筛前半段时花费的时间显然是比较长的,这时我们可以直接线性筛一遍先记录下来,在求后半段时就可以省下大部分时间,称为根号分治。
Powerful Number 筛
Powerful Number 筛类似于杜教筛,一样可以解决积性函数的前缀和。
要求为 :
存在一个函数 g 满足 :
- g 是积性函数
- g 易求前缀和
- 对于质数 p , f(p)=g(p)
现在即可求 ∑i=1nf(i) 。
对于所有正整数 n, 记 n 的质因数分解为 n=∏i=1mpiei。 n 是 PN 当球仅当 ∀1≤i≤m,ei>1 。
所有 PN 都可以表示成 a2b3 的形式,由此可得 n 以内的 PN 至多有 O(n) 个。
对于求解 F(n)=∑i=1nf(i), 首先构造出一个积性函数 g ,设 G(n)=∑i=1ng(i) 。
然后再去构造一个积性函数 h ,满足 f=g∗h,对于素数 p, 有 f(p)=g(1)h(p)+g(p)h(1)=h(p)+g(p),即 h(p)=0。
根据 h(p)=0 和 h 是积性函数可以推出对于非 PN 的数 n 有 h(n)=0 ,也就是说 h 在非 PN 处值为零。
然后就是推导 F(n),
F(n)=i=1∑nf(i)=i=1∑nd∣i∑h(d)g(di)=d=1∑nh(d)i=1∑⌊dn⌋g(i)=d=1∑nh(d)G(⌊dn⌋)=d=1,d is PN∑nh(d)G(⌊dn⌋)
然后 O(n) 找出所有 PN 计算 h, h 的值其实只需要计算 h(pc) 的值即可。
对于计算 h(pc),有两种方法,一种是根据公式计算,直接找出 h(pc) 仅与 p,c 有关的计算公式,另一种是根据 f=g∗h 有 f(pc)=∑i=0c=g(pi)h(pc−i),移项即可得出, h(pc)=f(pc)−∑i=1cg(pi)h(pc−i) ,就可以直接枚举素数 p 再枚举指数 c 求解出所有 h(pc)。
Min_25 筛
Min_25 筛可以在 O(lognn43) 的时间复杂度下求出积性函数的前缀和。
其要求 f(p) 是一个关于 p 的项数较少的多项式或可以快速求值, f(pc) 可以快速求值。
记号
如无特别说明,本节中所有记为 p 的变量的取值集合均为全体质数。
- x/y:=⌊yx⌋
- isprime(n):=[∣{d:d∣n}∣=2] ,即 n 为质数时其值为 1 ,否则为 0。
- pk :全体质数中第 k 小的质数(如:p1=2,p2=3)。特别地,令 p0=1 。
- lpf(n):=[1<n]min{p:p∣n}+[1=n] ,即 n 的最小质因数。特别地,n=1 时,其值为 1。
- Fprime(n):=∑2≤p≤nf(p)
- Fk(n):=∑i=2n[pk≤lpf(i)]f(i)
解释
观察 Fk(n) 的定义,可以发现答案即为 F1(n)+f(1)=F1(n)+1 。
考虑如何求出 Fk(n) 。通过枚举每个 i 的最小质因子及其次数可以得到递推式:
Fk(n)=i=2∑n[pk≤lpf(i)]f(i)=k≤ipi2≤n∑c≥1pic≤n∑f(pic)([c>1]+Fi+1(n/pic))+k≤ipi≤n∑f(pi)=k≤ipi2≤n∑c≥1pic≤n∑f(pic)([c>1]+Fi+1(n/pic))+Fprime(n)−Fprime(pk−1)=k≤ipi2≤n∑c≥1pic+1≤n∑(f(pic)Fi+1(n/pic)+f(pic+1))+Fprime(n)−Fprime(pk−1)
最后一步推导基于这样一个事实:对于满足
pic≤n<pic+1 的 c,有 pic+1>n⟺n/pic<pi<pi+1,故 Fi+1(n/pic)=0。
其边界值即为 Fk(n)=0(pk>n)。
假设现在已经求出了所有的 Fprime(n),那么有两种方式可以求出所有的 Fk(n) :
直接按照递推式计算。
从大到小枚举 p 转移,仅当 p2<n 时转移增加值不为零,故按照递推式后缀和优化即可。
现在考虑如何计算 Fprime(n)。
观察求 Fk(n) 的过程,容易发现 Fprime 有且仅有 1,2,…,⌊n⌋,n/n,…,n/2,n 这 O(n) 处的点值是有用的。
一般情况下,f(p) 是一个关于 p 的低次多项式,可以表示为 f(p)=∑aipci 。
那么对于每个 pci ,其对 Fprime(n) 的贡献即为 ai∑2≤p≤npci 。
分开考虑每个 pci 的贡献,问题就转变为了:给定 n,s,g(p)=ps ,对所有的 m=n/i ,求 ∑p≤mg(p)。
Notice:g(p)=ps 是完全积性函数!
于是设 Gk(n):=∑i=1n[pk<lpf(i)∨isprime(i)]g(i),即埃筛第 k 轮筛完后剩下的数的 g 值之和。
对于一个合数 x,必定有 lpf(x)≤x,则Fprime=G⌊n⌋,故只需筛到 G⌊n⌋ 即可。
考虑 G 的边界值,显然为 G0(n)=∑i=2ng(i)。(还记得吗?特别约定了 p0=1)
对于转移,考虑埃筛的过程,分开讨论每部分的贡献,有:
对于 n<pk2 的部分,G 值不变,即 Gk(n)=Gk−1(n)。
对于 pk2≤n 的部分,被筛掉的数必有质因子 pk,即 −g(pk)Gk−1(n/pk)。
对于第二部分,由于 pk2≤n⟺pk≤n/pk,故会有 lpf(i)<pk 的 i 被减去。这部分应当加回来,即 g(pk)Gk−1(pk−1)。
则有:
Gk(n)=Gk−1(n)−[pk2≤n]g(pk)(Gk−1(n/pk)−Gk−1(pk−1))
贝尔级数
对于积性函数 f(n) ,定义其在质数 p 意义下的贝尔级数为 :
Fp(x)=i=0∑∞f(pi)xi
相当于是把狄利克雷卷积变成了一般多项式卷积。
两个数论函数做狄利克雷卷积卷积,其贝尔级数相乘。
常见的贝尔级数,
- e⇒1
- I⇒∑i=0∞xi=x1
- idk⇒∑i=0∞pikxi=1−pkx1
- μ⇒1−x
- d⇒∑i=0∞(i+1)xi=(1−x)21⇒I∗I=d
- σ⇒∑i=0∞xi∑j=0ipjk=(1−x)(1−pkx)1⇒I∗idk=σk
- φ⇒1+∑i=1∞pi(1−p1)xi=1−px1−x⇒φ∗I=id
- w(n)=2k(n), k(n) 为不同质因子个数, w⇒1+∑i=1∞2xi=1−x1+1⇒w=μ2∗I
常用在筛法中构造卷积。