數(shù)論二

1、《算法藝術(shù)與信息學(xué)競賽》標(biāo)準(zhǔn)課件,數(shù)論(二): 經(jīng)典問題和算法劉汝佳,問題1: 大整數(shù)取余,給一個n位的大整數(shù)P和正整數(shù)m求P mod m的值輸入int nint digit[]：digit[0]為P第首位，digit[n-1]末位輸出int d：P mod m的值,分析,公式一: (a*b) mod m = (a mod m) * (b mod m) mod m公式二: (a+b) mod m = ((a mod m

2、) + (b mod m) mod m注意: 在兩個公式中, 都需要在最后對m取模,分析,由公式, 可以由前i-1位的余數(shù)d[i-1]推出第i位的余數(shù)d[i] = (d[i-1]*10+digit[i]) mod m時間復(fù)雜度: O(n)每個d[i]都只使用一次，空間復(fù)雜度O(1),d = 0;for(i = 0; i < n; ++i) d = (d * 10 + digit[i]) % m,問題2: 模取冪,給

3、出a, n, p, 求an mod p的值輸入int aint nint p輸出int d,分析,算法一: 利用公式an mod p = (an-1 mod p) * a mod p時間復(fù)雜度: O(n)算法二: 把n寫成二進(jìn)制n = (nknk-1…n0)2設(shè)di為a的2i次方模p的結(jié)果, 則an mod p = (dk^nk) * (dk-1^nk-1) * …即把nk為1的那些dk乘起來時間復(fù)雜度: O(

4、logn),分析,設(shè)si為數(shù)(nknk-1…ni), 則nk = sk & 1從右往左遞推dk = (dk-1 * dk-1) mod psk = sk-1 >> 1由于每個數(shù)只被用一次, 空間是O(1)的問題: dk-1*dk-1可能會overflow! (n=105時已經(jīng)會) 解決方法: 使用更大的整數(shù)范圍,分析,下面的數(shù)據(jù)類型LL取決于編譯器. gcc使用long long, 而VC++使用__int

5、64對于其他操作符, 把*d%p替換掉即可,LL ans = 1, d = a % p;do{ if(n & 1) ans = ans * d % p; d = (d * d) % p;}while(n >>= 1);,問題3: 求1~n的所有素數(shù),求1~n的所有素數(shù)輸入int n輸出bool isPrime[]：數(shù)i(1<=i<=n)是否為素數(shù)int prime[]：第i

6、(1<=i<=π(n+1))個素數(shù)int primeCount：素數(shù)總數(shù),分析,假設(shè)要求1~100的素數(shù)2是素數(shù), 刪除2*2, 2*3, 2*4, …, 2*50第一個沒被刪除的是3, 刪除3*3, 3*4, 3*5,…,3*33第一個沒被刪除的是5, 刪除5*5, 5*6, … 5*20得到素數(shù)p時, 需要刪除p*p, p*(p+1), … p*[n/p], 運(yùn)算量為[n/p]-p, 其中p不超過n1/2(想一

7、想, 為什么),Eratosthenes的篩子,primeCount = 0;for(i = 2; i < n; i++) isPrime[i] = true;for(i = 2; i < n; i++) if(isPrime[i]){ primes[primeCount++] = i;if(i <= n/i) for(j = i*i; j <= n

8、; j += i) isPrime[j] = false; },分析,問題：i*j可能超界！解決辦法： 改用除法判斷，并預(yù)先判斷是否有i2>n. 更好的辦法是遞推求i2, 利用(i+1)2-i2=2i+1, 另設(shè)一個變量k=i2100, 1000, 10000…109內(nèi)的素數(shù)個數(shù)分別為：25, 168, 1229, 9592, 78498, 664579, 5761455, 50847534用篩法一

9、般最多用來計算107內(nèi)的素數(shù),優(yōu)化,枚舉過程也可以優(yōu)化一下優(yōu)化一：除了2以外偶數(shù)都不是素數(shù)，所以每次i可以增加2優(yōu)化二：除了2、3以外，素數(shù)p除以6的余數(shù)只能是1或5，所以可以修改為每次交替增加2, 4時間優(yōu)化并不明顯, 但空間分別縮小為原來的1/2和1/3,分析,時間復(fù)雜度顯然為篩的時間復(fù)雜度+O(n)篩的過程不超過n/2+n/3+n/5+…由公式得: 篩的過程是nlnlnn的, 總O(nloglogn)其中B1為M

10、ertens常數(shù)0.2614972128…,問題4: 素數(shù)判定,給定正整數(shù)p判斷p是否為素數(shù)輸入int p輸出bool isPrime,算法一,樸素的枚舉法, 枚舉到n1/2為止任務(wù): 統(tǒng)計1~n的素數(shù)個數(shù)n=105時瞬間, n=106時需要數(shù)秒,for(i = 2; i * i <= p; i++) if(p % i == 0) return false;return true;,優(yōu)化,可以利用前面所說的

11、模6優(yōu)化, 速度為原來的若干倍. 優(yōu)化后106次一秒之內(nèi)出解, 107需要約10秒,if(p==2||p==3) return true;if(p%2==0||p%3==0) return false;for(i=5,k=4; i*i<=p; i+=(k=6-k)) if(p%i==0) return false;return true;,優(yōu)化,也可以改成只用素數(shù)試除, 速度變快但速度并不是很明顯(n=107時約兩倍

12、). 其中primes初始化為2和3, 隨著循環(huán)測試的進(jìn)行不斷更新,for(i=0; primes[i]*primes[i]<= p; i++) if(p % primes[i] == 0) return false;return true;,Miller-Rabin測試,對于奇數(shù)n, 記n=2r*s+1, 其中s為奇數(shù)隨機(jī)選a(1<=a<=n-1), n通過測試的條件是as≡1(mod n), 或者

13、存在0<=j<=r-1使得a(2^j)*s≡-1(mod n)素數(shù)對于所有a通過測試, 合數(shù)通過測試的概率不超過1/4注意: 先要判斷此數(shù)本身是否為a的倍數(shù),Miller-Rabin測試,能通過a=2的最小合數(shù)是2047=23*89能通過a=2, 3的最小合數(shù)是1373653能通過a=2, 3, 5的最小合數(shù)是能通過a=2, 3, 5, 7的, 2.5*1013以內(nèi)唯一的合數(shù)是3215031751,實(shí)現(xiàn)與實(shí)驗,計算

14、出k = as mod n是首先判斷k是否為1或者-1, 然后j每次加1, 相當(dāng)于每次t平方, 而不需要每次都調(diào)用冪取模的函數(shù)重算一次任務(wù)一: 統(tǒng)計1~n的素數(shù)個數(shù)預(yù)先判斷是否能被2, 3, 5, 7整除(而不是在每個測試中做), n=107減小到14秒, n=108約2.5分鐘任務(wù)二: 測試109~109+106的數(shù)答案: 共48155個素數(shù). 比較試除法和Miller-Rabin測試的速度,用a測試n的代碼,int r =

15、0, s = n - 1, j;if(!(n%a)) return false;while(!(s&1)){ s >>= 1; r++; }LL k = pow_mod(a, s, n);if(k == 1) return true;for(j = 0; j < r; j++, k = k * k % n) if(k == n - 1) return true;return false;,M

16、iller-Rabin測試的代碼,以下代碼適合于大于7, 小于109的數(shù)(如果要擴(kuò)展到2.5*1013, 則需要單獨(dú)判斷3215031751, 并避免LL乘法溢出),if(!miller_rabin(n, 2)) return false;if(!miller_rabin(n, 3)) return false;if(!miller_rabin(n, 5)) return false;if(!miller_rabin(n, 7))

17、 return false;return true;,問題5: 最大公約數(shù),給a, b, 求a和b的最大公約數(shù)gcd(a, b)輸入int aint b輸出int d,分析,方法一: 先分解素因數(shù), 然后求最大公約數(shù)方法二: (Euclid算法)利用公式gcd(a, b)=gcd(b, a mod b), 時間復(fù)雜度為O(logb)方法三: (二進(jìn)制算法) gcd(a,a)=a, a>b時a和b均為偶數(shù), gcd

18、(a,b)=2*gcd(a/2,b/2)a為偶數(shù), b為奇數(shù), gcd(a,b)=gcd(a/2,b)如果a和b均為奇數(shù), gcd(a,b)=gcd(a-b,b)不需要除法, 只需減法和右移, 適合大整數(shù),Euclid算法,測試: 求1<=a, b<=n的所有g(shù)cd之和n=1000時為4449880n=5000時為135637352需要的次數(shù)滿足[Lamé]:即steps<=4.785lgN

19、 + 1.6723Lamé定理: 算法的最壞情況為計算gcd(Fn, Fn-1)Heilbronn定理: 平均次數(shù)12ln2/π2lgn=0.843lgn,Euclid算法,遞歸形式和迭代形式(效率基本相當(dāng)),int gcd(int a, int b){ return (b? gcd(b, a%b) : a); },int gcd(int a, int b){ int t; while(b){ t

20、= a % b; a = b; b = t; } return a; },二進(jìn)制算法,int t = 1, c, d;while(a != b){ if(a >= 1; c = 1; }else c=0; if(!(b&1)){ b >>= 1; d = 1; }else d=0; if(c && d) t <<= 1; else if(!

21、c && !d) a -= b;}return t * a;,問題6: 線性不定方程,給出整數(shù)a和b和c, 求出一組整數(shù)x, y, 滿足ax + by = c輸入int a, int b, int c輸出int x, int y,分析,方程: ax+by=c設(shè)gcd(a,b)=d, 則如果c不是d的倍數(shù), 一定無解, 因為方程左邊為d的倍數(shù)否則可以令c=c’*d, 當(dāng)求出方程ax+by=d的任意解(x0

22、, y0)后, (c’x0, c’y0)就是原方程的解, 因ax+by=ac’x0 + bc’y0 =c’(ax0+by0)=c’d=c問題轉(zhuǎn)化為求方程ax+by=d的一組解,分析,ax+by=d的一組解可以用以下算法求出:,gcd(int a, int b, int& d, int& x, int& y){ if(!b){ d = a; x = 1; y = 0; } else{

23、 gcd(b, a%b, d, y, x); y -= x*(a/b);},證明,邊界: b = 0時a*1+b*0 = a成立, 否則遞歸調(diào)用后, y’和x’滿足by’+a%b*x’=d賦值后y = y’-a/b*x’, 計算ax+by即可注意a/b是整除運(yùn)算, 它等于(a-a%b)/b, 故ax+by = ax’+by’- ((a-a%b)/b*x’)*b= ax’+by’- (a-a%b)*x’

24、= by’+a%b*y’ = d,分析,如何通過這組解求出所有解?如果有兩組解(x1, y1)和(x2, y2), 有ax1+by1=d, ax2+by2=d因此, a(x1-x2)=b(y2-y1), 消去公因子d得a’(x1-x2)=b’(y2-y1)因為a’和b’互素, x1-x2必須是b’的整數(shù)倍. 設(shè)x1-x2+kb’, 則由等式得y2-y1=ka’, 因此x1-x2必須取整數(shù)值，且可取任意整數(shù)值,問題7: 模線性

25、方程,給出正整數(shù)a, b, n求解模線性方程ax≡b(mod n)輸入:Int a, int b, int n輸出Int d: 解剩余系的個數(shù)int sol[]: 對0<=i<d, x≡sol[i](mod n)是解,分析,首先明確一點(diǎn), 如果x是解(0<=x<n), 則對于任意整數(shù)k, x+kn也是解, 所以解應(yīng)表示成一些剩余類x≡xiax≡b(mod n)等價于存在整數(shù)y, 使得ax-ny=b

26、這是一個線性同余方程, 首先判斷d=(a,n)是不是b的約數(shù), 如果是, 等價于方程a’x-n’y=b’, 相當(dāng)于求解a’x≡b’(mod n’), (a’,n’)=1,分析,這個方程有兩種解法直接解ax-ny=d, 再兩邊同時乘以b/d. 注意這只是一個解, 一共應(yīng)該有d個, 間隔為n/d(因為對于模n/d來說解是唯一的)在模n’的縮系中a是存在逆a-1的. 因此解為x≡a’-1b’(mod n’), 同樣擴(kuò)展得d個解.在縮

27、系中求逆也有兩種方法求a’x-n’y=1的解, 和方法一等價利用歐拉定理a-1≡aφ(n)-1(mod n),分析,特別地, (a, n)=1時, 求a-1(mod n)程序如下:,int inverse(int a, int n){ int d, x, y; gcd(a, n, d, x, y); return x;},分析,方法一的程序如下,gcd(a, n, d, x, y);if(

28、b%d) d = 0; // no answerelse{ sol[0] = x * (b/d) % n; for(i = 1; i < d; i++) sol[i] = (sol[i-1] + n/d) % n;},分析,由剛才的分析, 方程ax≡b(mod n)完全被轉(zhuǎn)化為了x≡x0 (mod m0), 其中x0=sol[0], m0=n/(a, n)下面考慮由形如x≡xi (mod m

29、i)的方程所組成的方程組,問題8: 模線性方程組,有n個方程組成一個方程組x≡ai (mod mi)其中所有模mi兩兩互素求任意一個解x輸入int n, int a[], int m[]輸出int x,分析,令M=m1m2…mn中國剩余定理指出, 這樣的解在[0,M-1]的范圍內(nèi)是唯一的, 設(shè)它為x, 則任意解都能寫成x+kM的形式. 關(guān)鍵是要求出x記wi=M/mi，則(Mi,mi)=1, 因此存在pi,qi使得wi

30、pi+miqi=1,分析,關(guān)鍵等式: wipi+miqi=1記ei=wipi，則j=i時ei≡1(mod mj) (等式兩邊取模mj)j≠i時ei≡0(mod mj) (ei中含有因子mj)因此e1a1+e2a2+…+enan就是一個解調(diào)整到范圍[0,M-1]中程序?qū)崿F(xiàn): 一個一個方程考慮, 每次求出wi, pi, ei, 并累加eixi到結(jié)果中,分析,注意調(diào)用gcd時不要把x覆蓋掉了,M=1; for(i = 0; i

31、< n; i++) M*=m[i];for(i = 0; i < n; i++){ w = M / m[i]; gcd(m[i], w, dummy, dummy, y); x = (x + y*w*a[i]) % M;}return (n+x%M)%M;,分析,順便說一下, 對于任意的模m, 可以分解成素數(shù)冪乘積的形式, 則問題問題完全轉(zhuǎn)化為本問題(需要先判斷無解),問題9: 求歐拉函數(shù)φ,給

32、整數(shù)n, 求φ(n)輸入int n輸出Int phi,分析,對于素數(shù)的冪pk, 和pk不互素的數(shù)一定有因子p, 即p, 2p, 3p, …pk-1*p, 共pk-1個. 因此φ(pk)=pk-pk-1=pk(1-1/p)因為φ是積性函數(shù), 因此有公式需要得到素因數(shù)分解式. 注意每得到一個素因子p, 需要約去p的所有指數(shù),分析,注意乘法總在除法之后, 因此不會overflow,phi = n;for(i=2,j=4;j1)

33、 phi = phi / n * (n-1);,問題10: 求1~n的歐拉函數(shù),給正整數(shù)n, 求1~n每個數(shù)的歐拉函數(shù)值輸入Int n輸出Int phi[],分析,n的每個素因子p對n的phi值都有(1-1/p)的貢獻(xiàn), 因此可以借助素數(shù)篩法實(shí)現(xiàn), 但是必須從p開始而不是p*p開始, 因此更接近nlnlnn不需要額外的isPrime標(biāo)志, 直接判斷目前的phi. 如果為0, 說明為素數(shù)先檢查j的phi函數(shù)值. 如果為0, 說

34、明從來沒有被篩過, 賦初值j.篩的時候順便把phi[j]乘以(1-1/p), 注意先除再乘避免overflow,分析,可以順便求出素數(shù)表, 只需要在if(!phi[i])時先進(jìn)行primes[primeCount++] = i即可.,phi[1] = 1;for(i = 2; i <= n; i++) if(!phi[i]) for(j=i; j<=n; j+=i){ if(!ph

35、i[j]) phi[j] = j; phi[j] = phi[j] / i * (i-1); },問題11: 離散對數(shù),給a, b, n, 求方程ax≡b(mod n)的一個解輸入Int a, int b, int n輸出Int x,分析,離散對數(shù)是一個困難的問題這里只考慮n為素數(shù)的情況. 由歐拉定理, 只需要檢查x=0,1,2,…,n-1的情況, 時間復(fù)雜度為O(n)由于n為素數(shù),

36、只要a不為0, 一定存在逆a-1首先檢查前m項, 即a0, a1, a2, … am-1模n的值是否為b, 并把a(bǔ)i mod n保存在e[i]里求出am的逆a-m,分析,下面檢查am, am+1, … a2m-1如果有解, 相當(dāng)于存在i使得ei*am≡b(mod n)兩邊左乘v = a-m, 得到ei≡v*b(mod n)只需檢查是否存在e[i]等于給定值事先給e排序, 每次可二分查找每處理一輪, 查找內(nèi)容需要乘以v,

37、分析,需要O(m)的時間計算e數(shù)組, O(n/m*logm)的時間進(jìn)行剩下的工作, 取m=n1/2, 則時間復(fù)雜度為O(n1/2logn)可以用sort+二分查找, 也可以用map下面的程序使用了一個map x, 其中x[j]代表滿足ei=j的最小i(可能有多個i, 如a=1時),m = (int)ceil(sqrt(n));v = inverse(pow_mod(a, m, n), n);e = 1; x[1] = 0;fo

38、r(i = 1; i < m; i++){ e = e * (LL)a % n; if(!x.count(e)) x[e] = i; }for(i = 0; i < m; i++){ if(x.count(b)) return i*m+x[b]; b=b * (LL)v % n; },問題12: 模平方根,給整數(shù)a, n求滿足x2≡a(mod n)的所有根輸入Int a, int