幻灯片 1

Size: px

Start display at page:

Download "幻灯片 1"

右蓝
5 years ago
Views:

1 北京大学暑期课 ACM/ICPC 竞赛训练北京大学信息学院郭炜课程网页 :

2 信息科学技术学院程序设计与算法 ACM/ICPC 中的数学郭炜 / 林舒 2

3 一些数论基本定理 ( a + b ) mod c = ((a mod c) + (b mod c)) mod c ( a * b ) mod c = ((a mod c) * (b mod c)) mod c 3

4 一些数论基本定理消去律 : 若 gcd(c,p) = 1, 则 ac bc mod p => a b mod p 4

5 一些数论基本定理消去律证明 : ac bc mod p => ac = bc + kp => c(a-b) = kp c 和 p 没有除 1 以外的公因子 => 1) c 能整除 k 或 2) a = b 如果 2 不成立, 则 c kp c 和 p 没有公因子 => c k, 所以 k = ck' => c(a-b)=kp 可以表示为 c(a-b) =ck'p 因此 a-b = k'p, 得出 a b mod p 5

6 快速幂 int Pow(int a,int b) { // 快速求 a^b, 复杂度 log(b) if(b == 0) return 1; if(b & 1) { //b 是奇数 return a * Pow(a,b-1); } else { int t = Pow(a,b/2); return t * t; } } 6

7 快速幂 int Pow(int a,int b) { // 快速求 a^b, 复杂度 log(b) int result = 1; int base = a; while(b) { if( b & 1) result *= base; base *= base; b >>= 1; } return result; } 7

8 快速幂取模 int PowMod(int a,int b,int c) {// 快速求 a^b % c, 要避免计算中间结果溢出 int result = 1; int base = a%c; while(b) { if( b & 1) result = (result * base)%c; base = (base * base) % c; b >>= 1; } return result; } 8

9 等比数列二分求和取模 S n = a+a a n 要求 S n mod p 如果用公式算, 可能溢出, 因此用二分法求 1) 若 n 是偶数 S n = a+...+a n/2 + a n/2+1 + a n/ a n/2+n/2 =(a+...+a n/2 ) + a n/2 (a+...+a n/2 ) =S n/2 + a n/2 S n/2 =(1+a n/2 )S n/2 9

10 等比数列二分求和取模 2) 若 n 是奇数 S n = a+...+a (n-1)/2 + a (n-1)/ a (n-1)/2+(n-1)/2 + a (n-1)/2+(n-1)/2 + 1 =S (n-1)/2 + a (n-1)/2 (a+...+a (n-1)/2 )+a n =(1+a (n-1)/2 )S (n-1)/2 +a n 10

11 等比数列二分求和取模 int PowSumMod(int a,int n,int p) {// return (a+ a^ a^n) Mod p; if( n == 1) return a%p; if( n %2 == 0) return (1+PowMod(a,n/2,p))*PowSumMod(a,n/2,p) % p; else return ((1+PowMod(a,(n-1)/2,p)) * PowSumMod(a,(n-1)/2,p) + PowMod(a,n,p)) % p; } 11

12 POJ3233 Matrix Power Series 矩阵快速幂 + 等比数列二分求和取模给一个 n n 的整数矩阵 A 和正整数 k,m, 令 S = A + A 2 + A A k 求 S mod m (S 的每个都 mod m) n (n 30), k (k 10 9 ) and m (m < 10 4 ). 12

13 矩阵快速幂取模 struct Matrix { T a[32][32]; int r; // 行列数 Mat(int rr):r(rr) { } void MakeI() { // 变为单位矩阵 memset(a,0,sizeof(a)); } for(int i = 0;i < r; ++i) a[i][i] = 1; }; 13

14 矩阵快速幂取模 Matrix Pow(const Matrix & m,int k,int p) { // 求 m k mod p int r = m.r; Matrix result(r); result.makei(); //MakeI 是将 result 变为单位矩阵 Matrix base = m; while(k) { } if( k & 1) result = Product(result,base,p); //result*base mod p k >>= 1; base = Product(base,base,p); return result; } 14

15 欧几里得算法求最大公约数 int Gcd(int a,int b) { if( b == 0) return a; return Gcd(b,a%b); } 最小公倍数 :lcm(a,b) = a*b/gcd(a,b) 15

16 扩展欧几里得算法 ax+by=gcd(a,b) 有整数解 (x,y) <=> bx+(a%b)y = gcd(a,b) 有整数解 (x1,y1), 且 x = y1, y = x1-[a/b]*y1 因此, 可以在求 gcd(a,b) 的同时, 对 ax+by=gcd(a,b) 求解 16

17 扩展欧几里得算法 int GcdEx(int a,int b,int &x,int & y) // 求 ax+by=gcd(a,b) 的整数解, 返回 gcd(a,b) { if( b == 0) { x = 1; y = 0; return a; } int x1,y1; int gcd = GcdEx(b,a%b,x1,y1); x = y1; y = x1 - a/b * y1; return gcd; } 17

18 扩展欧几里得算法 ax+by=c 有解的充要条件是 gcd(a,b) c 设 d = gcd(a,b), k = c/d, ax+by = d 的解是 (x1,y1) 则 ax+by = c 的解集是 : x = k*x1 + t*(b/d) y = k*y1 t*(a/d) t 为任意整数 18

19 扩展欧几里得算法求 ax c(mod b) 求 ax c(mod b) 等价于求 ax+by = c 设 d = gcd(a,b), k = c/d, ax+by = d 的解是 (x1,y1) 则 ax c(mod b) 的解集是 : x = k*x1 + t*(b/d) t 为任意整数其中模 b 不同的解共有 d 个, 为 : x = k*x1 + t*(b/d) t=0,1,..d-1 19

20 扩展欧几里得算法求 ax c(mod b) 求 ax c(mod b) 的最小非负整数解 : 令 ans = k * x1; s = b/d; 则 x = ans + t*s t 为任意整数则最小非负整数解是 :(ans%s + s)%s 20

21 例题 : 百练 2115 C Looooops 对下面的循环 : for (variable = A; variable!= B; variable += C) statement; A,B,C 和 variable 都是 K(K<=32) 比特的无符号数,+= 运算的结果也是 K 比特的无符号数, 给定 A,B,C 和 K, 问循环执行多少次 21

22 例题 : 百练 2115 C Looooops 对下面的循环 : for (variable = A; variable!= B; variable += C) statement; A,B,C 和 variable 都是 K(K<=32) 比特的无符号数,+= 运算的结果也是 K 比特的无符号数, 给定 A,B,C 和 K, 问循环执行多少次 ( 可能无穷多次 ) 实际上就是求 A + CX B (mod 2 K ) 的最小非负整数解 CX (B-A) (mod 2 K ) 22

23 中国剩余定理孙子定理, 韩信点兵, 隔墙算, 鬼谷算, 大衍求一术... " 物不知数 " 问题 :" 今有物不知其数, 三三数之剩二, 五五数之剩三, 七七数之剩二, 问物几何? 答曰 :' 二十三.' 术曰 : 三三数之剩二, 置一百四十, 五五数之剩三, 置六十三, 七七数之剩二, 置三十, 并之, 得二百三十三, 以二百一十减之, 即得. 凡三三数之剩一, 则置七十, 五五数之剩一, 则置二十一, 七七数之剩一, 则置十五, 即得." -- 孙子算经 23

24 中国剩余定理设 n 组数 (ai, bi), 其中 bi 两两互素求 x 使得 x = a1 mod b1 x = a2 mod b2... x = an mod bn 24

25 中国剩余定理给定两两互质的正整数 n 1,n 2,...,n k, 要求找到最小的正整数 x, 满足方程组 x a i (mod n i ) (i=1,2...k) 算法步骤 : 令 n=n 1 n 2...n k, m i =n/n i 显然 gcd(m i,n i )=1, 利用扩展欧几里德算法计算出 x i 满足 m i x i 1(mod n i ) x = (a 1 x 1 m 1 +a 2 x 2 m a k x k m k ) mod n 此方程组任意两个解模 n 同余, 因此 x 就是最小的

26 中国剩余定理此方程组任意两个解模 n 同余证 : 设有两个解 x1,x0 则 : x1 a i (mod n i ) (i=1,2...k) x0 a i (mod n i ) (i=1,2...k) => n i (x1-x0) (i=1,2...k) n i (x1-x0) 且 n i 两两互质 => n (x1-x0) => x1 x0 (mod n)

27 中国剩余定理的一般情况给定正整数 n 1,n 2,...,n k ( 未必两两互质 ), 要求找到 x, 满足 x a i (mod n i ) (i=1,2...k) x a 1 ( mod n 1 ) x a 2 ( mod n 2 ) => x + u*n 1 = a1 且 x - v*n 2 = a 2 => n 1 *u + n 2 *v = (a 1 -a 2 ) 此关于 u,v 的方程, 当且仅当 gcd(n 1,n 2 ) (a 1 -a 2 ) 时有解

28 中国剩余定理的一般情况设用扩展欧几里得算法求得 n 1 *u + n 2 *v = (a 1 -a 2 ) 根为 (u 0,v 0 ) 则 : 则此方程解集为 : u = u 0 + t*(n 2 /gcd(n 1,n 2 )) v = v 0 t*(n 1 /gcd(n 1,n 2 )) t 为任意整数 x + u*n 1 = a 1 => x 的解集为 :a 1 u 0 n 1 tn 1 n 2 /gcd(n 1,n 2 ) 即 :x = a 1 u 0 n 1 t*lcm(n 1,n 2 ) t 为任意整数

29 中国剩余定理的一般情况 x = a 1 u 0 n 1 t*lcm(n 1,n 2 ) t 为任意整数等价于将 x a 1 ( mod n 1 ) x a 2 ( mod n 2 ) 合并为 : x a 1 u 0 n 1 (mod lcm(n 1,n 2 )) 再继续两个方程并为一个, 最后就能求得 x 的解空间

30 线性筛法求素数朴素筛法求 n 以内的所有质数初始时容器内为 2 到 n 的所有数取出最小的数 p,p 一定是质数, 删去除 p 外的所有 p 的倍数重复上述步骤直到向后找不到没被删掉的数缺陷 : 一个数可能被重复删去多次, 影响效率

31 线性筛法求素数改进 : 对每个素数 p 考虑所有 i, 若 i 的最小素因子 >=p, 则将 i*p 去掉 i = q 1 *q 2 *...q n q i 是素数, q 1 >= p i*p = p*q 1 *q 2 *...q n i*p 只会被删掉一次, 只在考察 p 的时候被删, 不会在考察 q 1 q 2 q n 的时候被删

32 #include <iostream> #include <vector> int main() { using namespace std; int n; } cin >> n; // 求 n 以内素数 vector<int> prime; vector<bool> isprime(n+1); for(int i = 1;i <= n; ++i) isprime[i] = true; for(int i = 2; i <= n; ++i) { if( isprime[i]) // 处理到 i 时它还没被删掉, 则 i 为素数 prime.push_back(i); for(int j = 0;j < prime.size() ; ++j) { if( i*prime[j] <= n) isprime[i*prime[j]] = false; else break; if( i % prime[j] == 0)// prime[j] 是 i 的最小素因子 break; } } for(int i = 0;i < prime.size(); ++i) cout << prime[i] << endl; return 0;

33 欧拉函数和欧拉定理欧拉函数 : φ(n) = 小于 n 且和 n 互质的正整数 ( 包括 1) 的个数 (n 为正整数 ) 完全余数集合 : Zn = { 小于 n 且和 n 互质的数 } Zn =φ(n) 对于素数 p,φ(p) = p -1

34 欧拉函数和欧拉定理两个不同素数 p,q, n=p*q, 则证 : φ(n) = (p -1) * (q -1) Zn = {1, 2, 3,..., n - 1} {p, 2p,..., (q - 1) * p} {q, 2q,..., (p - 1) * q} φ(n) = (n - 1) - (q - 1) - (p - 1) = (p * q p) (q-1) = p(q-1) (q-1) = (p -1) * (q -1) =φ(p) * φ(q)

35 欧拉函数和欧拉定理消去律 : 如果 gcd(c,p) = 1, 则 ac bc mod p a b mod p 欧拉定理 : 互质的正整数 a 和 n, 有 a φ(n) 1 mod n 证 : (1) 令 Zn = {x 1, x 2,..., x φ(n) }, S = {ax 1 mod n, ax 2 mod n,..., ax φ(n) mod n}, 则可证 Zn = S 1 a 与 n 互质且 x i 与 n 互质 => ax i 与 n 互质 => ax i mod n Zn 2 若 i j, 则 x i x j, 且由 a, n 互质可得 ax i mod n ax j mod n ( 否则根据消去率,x i =x j ) (2) a φ(n) x 1 x 2...x φ(n) mod n (ax 1 )(ax 2 )...(ax φ(n) ) mod n (ax 1 mod n)(ax 2 mod n)...(ax φ(n) mod n) mod n x 1 x 2... x φ(n) mod n 因为 x i (1 i φ(n)) 与 n 互质所以 a φ(n) 1 mod n ( 消去律 )

36 费马小定理若正整数 a 与素数 p 互质, 则有 a p mod p φ(p) = p -1, 代入欧拉定理即证

37 欧拉函数的公式 (1) p 为素数, 正整数 n = p k, 则证 : φ(n) = p k - p k-1 小于 p k 的正整数个数为 p k - 1 个, 其中和 p k 不互质的正整数有 { p,2p,3p..., (p k-1-1)p} 共计 p k 个所以 φ(n) = p k (p k-1-1) = p k - p k - 1

38 欧拉函数的公式 (2) p,q 是两个互质的正整数,n=pq 则 : φ(n) = φ(p)φ(q) 证 : 略 (3) 当 b 是质数,a%b=0, 则 : φ(ab)=φ(a)b 证 : 设 a=kb n 且 gcd(k,b)=1, 则 φ(a)=φ(k)φ(b n ),φ(k)=φ(a)/φ(b n ) φ(ab) = φ(kb n+1 )= φ(k)φ(b n+1 )= (φ(a)/φ(b n ))*φ(b n+1 ) = φ(a)*(φ(b n+1 ) /φ(b n )) = φ(a)b

39 欧拉函数的公式 (4) 对任意 n: φ(n) = n(1-1/p 1 )(1-1/p 2 )...(1-1/p n ) n = p 1 a1 p 2 a2... p n an (p i 为素数 ). 证 : φ(n)= φ(p 1 a1 )φ(p 2 a2 )... φ(p n an ) =(p 1-1)p 1 (a1-1) (p 2-1)p 2 (a2-1)...(p n -1) p n (an-1) =(p 1 a1 p 2 a2... p n an )(p 1-1)(p 2-1)...(p n -1)/(p 1 p 2... p n ) =n(1-1/p 1 )(1-1/p 2 )...(1-1/p n ) //=> 对 n>2 φ(n) 为偶数

40 欧拉函数的公式递推式 : 质数 p 满足 p x, 若 p 2 x, 则 φ(x)=φ(x/p)*p, 若 p 2 x 不成立则 φ(x)=φ(x/p)*(p-1) (1) 若 p 2 x 不成立, 因 x/p 和 p 互质, 故 : φ(x)=φ(x/p)*φ(p)=φ(x/p)*(p-1) (2) 若 p 2 x, 因 p (x/p) 且 p 是质数, 故 φ(x) = φ(x)=φ(x/p)*p ( 当 b 是质数,a%b=0, 则 : φ(ab)=φ(a)b)

41 用递推式求欧拉函数 int euler(int n) {// 递推式 : 质数 p 满足 p x, 若 p 2 x, 则 φ(x)=φ(x/p)*p, 若 p 2 x 不成立则 φ(x)=φ(x/p)*(p-1) int ret=1; for(int i=2;i*i<=n;i++){ // 只要考虑 sqrt(n) 以内的素数 if(n%i==0){//i 若是合数则不可能满足此条件 n/=i,ret*=i-1; // 最后总会剩一个 i,n/i 不整除 i while(n%i==0){ n/=i,ret*=i; } } } if(n>1) ret*=n-1;// 到此所有 <= sqrt( 最初的 n) 的因子都已经去掉如果还剩一个因子, 该因子一定 > sqrt( 最初的 n), 则该因子只有一个 return ret; }

42 POJ3101 Astronomy 题目大意 : 有 n 个行星, 它们的轨道是同一平面上的同向同心圆, 且它们始终做匀速圆周运动, 周期 t i 已知所有卫星都处于过圆心的某条直线上的现象, 称为卫星平行求相邻两次卫星平行现象的间隔时间, 用分数表示

43 POJ3101 Astronomy 算法思路 : 所有卫星平行 <=> 任意两个卫星平行 <=> 相邻两个卫星平行 ( 卫星平行具有传递性 ) 两个卫星 i,j 平行的最小时间间隔 t 为它们的运行圈数差半圈, 即 t/t 1 -t/t 2 = 0.5, 故 t = 0.5/(1/t i -1/t j ) 写出相邻两个卫星平行的时间间隔 d i =b i /a i, 则问题转化为求这 n-1 个分数的 " 最小公倍数 "(b i a i 互质 ) 分母 p=gcd(a 1,a 2,...,a (n-1) ), 分子 q=lcm(b 1,b 2,...,b (n-1) ), 约分即得最终答案

44 POJ2142 The Balance 题目大意 : 现有质量为 a 和 b 的砝码, 数量不限要求在天平上称出质量为 d 的物品, 天平左右均可放砝码求一种可行方案, 要求 : 放置砝码数量尽可能少 ; 数量相同时, 总质量尽可能少

45 POJ2142 The Balance 问题转化 : 求 ax+by=d 的一组整数解 (x,y), 要求 x + y 尽可能小, 若相等, 则 a x +b y 尽可能小 (x<0, 表示砝码和物体放在同一侧 ) 先求出不定方程的一组特解 (x 0,y 0 ), 令 m=gcd(a,b),a'=a/m,b'=b/m, 则通解为 x=x 0 +b't,y=y 0 -a't (a',b'>0,t 为整数 ) x + y 最小时, 要么 x 是最小正解, 要么 y 是最小正解比较一下即可

46 扩展欧几里得算法 ax+by=c 有解的充要条件是 gcd(a,b) c 设 d = gcd(a,b), k = c/d, ax+by = d 的解是 (x1,y1) 则 ax+by = c 的解集是 : x = k*x1 + t*(b/d) y = k*y1 t*(a/d) t 为任意整数 46

47 POJ2689 Prime Distance 给定区间 [L,U],L 和 U 可以很大, 但区间长度不超过 10 6 (1<=L< U<=2,147,483,647) 求这个区间中最近和最远的两对素数

48 POJ2689 Prime Distance 直接试除? 耗费时间太长直接筛法? 耗费空间太大区间长度不大 -> 筛法 + 试除先用筛法求出不大于 sqrt(u) 的所有素数, 然后用这些素数一一试除

49 POJ2478 Farey Sequence 求所有分母不大于 n 的既约真分数个数

50 POJ2478 Farey Sequence 分母为 x 的既约真分数有 φ(x) 个分母不大于 n 的既约真分数个数为 φ(1)+φ(2)+...+φ(n)

51 POJ3696 The Luckiest number 定义 : 只含有数字 8 的数为幸运数给定正整数 L(1 L 2,000,000,000)., 求 L 的所有倍数中最小的幸运数的位数

52 The Luckiest number 设最终答案为 x, 则 x 满足 (10 x -1)/9*8 0(mod L) 化简上述式子 : =>(10 x -1)*8 0(mod 9L)=>10 x -1 0(mod 9L/gcd(9L,8)) =>10 x 1(mod 9L/gcd(9L,8)) 令 m=9l/gcd(9l,8), 若 gcd(10,m)>1, 显然无解若 gcd(10,m)=1, 由欧拉定理 :10 φ(m) 1(mod m) 不过, 我们要求的是最小解, 而 φ(m) 只是一个可行解可以证明, 最小解一定是 φ(m) 的一个约数枚举 φ(m) 所有约数, 用快速幂验证即可 10 k 1 (mod n) 则成立的最小 k 是 ϕ (n) 的约数这个性质竞赛中经常用到

RSA 图为 RSA 公开密钥算法的发明人, 从左到右 Ron Rivest, Adi Shamir, Leonard Adleman. 照片摄于 1978 年裴士辉 QQ:

RSA 图为 RSA 公开密钥算法的发明人, 从左到右 Ron Rivest, Adi Shamir, Leonard Adleman. 照片摄于 1978 年裴士辉 QQ:168159305 RSA 的数论基础质数 ( 素数 ) (prime number) 一个大于 1 的自然数, 除了 1 和它本身以外不再有其他的因数, 那么这个数为素数 ; 否则称为合数最小的质数是目前为止, 人们未找到一个公式可求出所有质数