转自https://www.cnblogs.com/bolderic/p/7406204.html
一些常用的简单数论模板以及书中的定理
组合数取模
1:N M < 1000, 杨辉三角双循环,边加边取模(代码未添加取模)
1 c[1][1] = c[1][0] = 1; 2 for(int i = 2; i <= 50; i++){ 3 c[i][0] = 1; 4 for(int j = 1; j <= i; j++){ 5 c[i][j] = (c[i - 1][j - 1] + c[i - 1][j]) % mod;
6 }
7 }
2. C(n, k) = C(n, k - 1) * (n - k + 1) / k,对2~maxn的逆元打表后,利用这个公式可以O(n)时间内求出C(n, k)
1 void getcom(int x){ 2 c[0] = c1[0] = 1; c[1] = x; c1[1] = n - x; 3 for(int i = 2; i <= x; i++) 4 // inv[i] 指i的逆元,可以使用扩展gcd打表 5 c[i] = c[i - 1] % mod * (x - i + 1) % mod * inv[i] % mod; 6 }
GCD求最大公约数
1 //欧几里得算法,gcd() 2 int gcd(int a, int b){ 3 return b == 0 ? a: gcd(b, a%b); 4 }
扩展欧几里得算法求ax + by = c 的解以及判断是否有解 -- 当c为gcd(a,b)的倍数
1 void gcd(int a, int b, int &d, int& x, int& y) { 2 if(!b) { d = a; x = 1; y = 0;} 3 else{ 4 gcd(b, a%b, d, y, x); y -=x*(a/b); 5 } 6 }
Eratosthenes's sieve 埃氏筛选法求素数
1 memset(vis, 0, sizeof vis); 2 for(int i = 2; i <= n; i++){ 3 for(int j = i*2; j <=n; j+=i){ 4 vis[j] = 1; 5 } 6 }
欧拉函数:φ(n) = 小于n的正整数中与n互质的数的数目,φ(1) = 1.
φ(n) = n * (1 - 1 / p1) * ( 1 - 1 / p2) …… (1 - 1 / pk). pi为n的素因子
O(√n)求n的 欧拉函数值:
1 int euler(int n){ 2 int m = (int)sqrt(n + 0.5); 3 int ans = n; 4 for(int i = 2; i <= m; i++){ 5 if(n % i == 0){ 6 ans = ans / i * (i - 1); 7 while(n % i == 0) n / i; 8 } 9 } 10 if(n > 1) ans = ans / n * (n - 1); 11 return ans; 12 }
筛选法求1~n欧拉函数值 -- 与埃氏筛选法一样,如果一个数j是i的倍数,那么其欧拉函数值就phi[j] = phi[j] / i * (i - 1) 并标记当前值,同时这也保证了i一定是素数 。当然初始值phi[i] = i
1 memset(phi, 0, sizeof(phi)); 2 phi[1] = 1; 3 for(int i = 2; i <= n; i++){ 4 if(!phi[i]){ 5 for(int j = i; j <= n; j += i){ 6 if(!phi[j]) phi[j] = j; 7 phi[j] = phi[j] / i * (i - 1); 8 } 9 } 10 }
欧拉公式:若gcd(a, m) = 1, 则 aφ(m) ≡ 1 (mod m).
φ函数公式: a): 如果p是素数且k >= 1, 则 φ(pk) = pk - pk-1.
b): 如果 gcd(m, n) = 1, 则 φ(mn) = φ(m)φ(n).
指数降幂:f (1) = 1, f (n) = nf(n - 1). 求 f(n) mod m.
n <= 5 直接计算即可
n > 5: z = f(n - 1) mod φ(m), f(n) mod m = nφ(m) + z mod m.
然后使用快速幂:
1 ll n, m; 2 ll pow_mod(ll a, ll b, int mod){ 3 if(b == 0) return 1; 4 ll tmp = pow_mod(a, b / 2, mod); 5 ll ans = tmp * tmp % mod; 6 if(b & 1) ans = a % mod * ans % mod; 7 return ans; 8 } 9 int euler(int x){ 10 int y = (int)sqrt(x + 0.5); 11 int ans = x; 12 for(int i = 2; i <= y; i++){ 13 if(x % i == 0){ 14 ans = ans / i * (i - 1); 15 while(x % i == 0) x /= i; 16 } 17 } 18 if(x > 1) ans = ans / x * (x - 1); 19 return ans; 20 } 21 ll cacu(ll x, ll mod){ 22 if(mod == 1) return 0; 23 if(x == 1) return 1 % mod; 24 else if(x == 2) return 2 % mod; 25 else if(x == 3) return 9 % mod; 26 else if(x == 4) return (1 << 18) % mod; 27 else if(x == 5) return pow_mod(5, 262144, mod); 28 else{ 29 int phi = euler(mod); 30 ll z = cacu(x - 1, phi); 31 ll ans = pow_mod(x, phi + z, mod); 32 return ans; 33 } 34 }
快速幂
1 //快速幂 2 int pow(int a, int b){ 3 if(b == 0) 4 return 1; 5 if(b & 1) { 6 return a * pow(a, b - 1); 7 } 8 else { 9 int t = pow(a, b/2); 10 return t * t; 11 } 12 } 13 //快速幂,迭代 14 int pow(int a, int b) { 15 int result = 1; 16 int base = a; 17 while(b) { 18 if(b & 1) { 19 result *= base; 20 } 21 base *= base; 22 b >> 1; 23 } 24 }
快速幂取模
1 //幂取模 2 int pow_mod(int a, int n, int m){ 3 if(n == 0) return 1; 4 int x = pow_mod(a, n/2, m); 5 long long ans = (long long)x * x % m; 6 if(n%2 == 1) ans = ans * a % m; 7 return (int) ans; 8 } 9 //快速求 a^b % c,要避免中间结果溢出 10 int powmod(int a, int b, int c) { 11 int result = 1; 12 int base = a%c; 13 while(b) { 14 if(b&1){ 15 result = (result * base)%c; 16 } 17 base = (base * base) %c; 18 b >>= 1; 19 } 20 return result; 21 }
合数的拉宾-米勒测试(SPOJ PON):设n是奇素数,记n-1=2kq,q为奇数。对不被n整除的某个a,如果下述两个条件都成立,则n是合数。
a): aq !≡ 1 (mod n), 即a的q次方不与1模n同余
b): 对所有i = 0,1,2,……,k-1, a2^i * q !≡ -1 (mod n)
代码:
1 VI primes; 2 bool vis[maxn]; 3 4 ll mul_mod(ll a, ll b, ll m){ 5 ll ans = 0, base = a; 6 while(b){ 7 if(b & 1) ans = (ans + base) % m; 8 base %= m; 9 base = base * 2 % m; 10 b >>= 1; 11 } 12 return ans; 13 } 14 ll pow_mod(ll a, ll b, ll m){ 15 ll ans = 1, base = a % m; 16 while(b){ 17 if(b & 1) ans = mul_mod(ans, base, m); 18 base %= m; 19 base = mul_mod(base, base, m); 20 b >>= 1; 21 } 22 return ans % m; 23 } 24 void dowork(){ 25 memset(vis, 0, sizeof(vis)); 26 for(int i = 2; i < maxn / 2; i++) 27 for(int j = i * 2; j > 0 && j < maxn; j += i) 28 vis[j] = 1; 29 for(int i = 2; i < 1e5; i++) 30 if(!vis[i]) primes.push_back(i); 31 } 32 bool check(ll x){ 33 if(x == 2) return true; 34 if(x % 2 == 0) return false; 35 ll q = x - 1; 36 int k = 0; 37 while(!(q & 1)){ 38 k++; q >>= 1; 39 } 40 bool flag = true; 41 for(int i = 0; i < 100; i++){ 42 ll a = primes[i]; 43 if(a == x) return true; 44 if(pow_mod(a, q, x) == 1) 45 continue; 46 bool none = true; 47 for(int j = 0; j < k; j++){ 48 ll qi = 1 << j; 49 ll tmp = pow_mod(a, qi * q, x); 50 if(tmp == x - 1){ 51 none = false; 52 break; 53 } 54 else if(tmp == 1) 55 break; 56 } 57 if(none) { 58 flag = false; 59 break; 60 } 61 } 62 return flag; 63 }
a模p的阶:a模p的阶 ep(a) = (使得 ae ≡ 1 (mod p) 的最小指数 e >= 1)
原根 :具有最高次 ep(g) = p - 1 的数g成为模p的原根。
阶整除性质:设a是不被素数p整除的整数,假设an ≡ 1 (mod p), 则阶 ep(a) 整除n. 特别的,ep(a)总整除 p - 1.
判断x是否是素数p的原根:参考:http://www.apfloat.org/prim.html
若x不是原根,则 ep(x) 一定整除 p - 1。基于此,我们可以将p - 1质因子分解,然后对每个质因子 f 判断是否 x(p - 1) / f ≡ 1 (mod p),如果都不是,则x为模p的原根.
SPOJ模板题:http://www.spoj.com/problems/PROOT/
威尔逊定理:当且仅当p为素数时:( p -1 )! ≡ -1 ( mod p )
费马小定理 -- 假如p是质数,且gcd(a,p)=1,那么 a(p-1)≡1(mod p)
逆元 -- 对于正整数a和m,如果有ax ≡ 1 % m,那么把这个同余方程中 x 的最小正整数解叫做a 模 m 的逆元. 可以使用扩展欧几里得求逆元。
中国剩余定理:设m与n是整数,gcd(m, n) = 1, b 与 c是任意整数,则同余式组 x ≡ b (mod m) 与 x ≡ c (mod n) 恰有一解 0 <= x < mn.
欧拉完全数定理:如果n是偶完全数,则n形如 n = 2p-1(2p - 1), 其中2p - 1是梅森素数。
σ函数公式:σ(n) = n 的所有因数之和(包括1和n).
a): 如果p是素数,k >= 1,则 σ(pk) = 1 + p + p2 + …… + pk = (pk+1 - 1) / (p - 1)
b): 如果gcd(m, n) = 1, 则 σ(mn) = σ(m)σ(n)
等比数列公式
斐波那契数列通式: