基本数论算法

dalao博客，至少很好看。。

因为本人数论实在渣渣，但是考试确是得考的，只好尽早学，尽早掌握。

最大公因数

普通gcd

O(log(min(a,b)))

inline int gcd(int x,int y)
{
    return y == 0 ? x : gcd(y, x % y)
}

 

二进制优化gcd

O(log(min(a,b))) 没有%常数减小

inline int bsgcd(int x, int y)
{
    if(x == y) return x;
    if(x < y) x ^= y ^= x ^= y;
    if(!(x & 1)) //x偶 y偶 gcd(x,y)=2*gcd(x/2,y/2)，x偶 y奇 gcd(x,y)=gcd(x/2,y)
         return (!(y & 1)) ? bsgcd(x >> 1, y >> 1) << 1 : bsgcd(x >> 1, y);
     //x奇 y偶 gcd(x, y)=gcd(x,y/2)
    return (!(y & 1)) ? bsgcd(x, y >> 1) : bsgcd(y, x - y); 
}

素数判定

埃拉托色尼筛法

O(nloglogn)

不够优越，所以不学。

线性筛

模板题

O(n)

const int MAXN = 10000001;

bool notpri[MAXN];
int n, m, cnt, prime[MAXN];

inline void get_prime()
{
    int i, j;
    notpri[0] = notpri[1] = 1; // 1 表示不是质数 
    for(i = 2; i <= n; i++)
    {
        if(!notpri[i]) prime[++cnt] = i;
        for(j = 1; i * prime[j] <= n; j++)
        {
            notpri[i * prime[j]] = 1;
            if(!(i % prime[j])) break;
        }
    }
}

Miller-Rabin 大素数判定

http://www.cnblogs.com/vongang/archive/2012/03/15/2398626.html

模板题

O(log₃₂n)

费马小定理：对于素数p和任意整数a，有a^p ≡ a(mod p)（同余）。反过来，满足a^p ≡ a(mod p)，p也几乎一定是素数。

伪素数：如果p是一个正整数，如果存在和p互素的正整数a满足 a^p-1 ≡ 1(mod p)，我们说p是基于a的伪素数。如果一个数是伪素数，那么它几乎肯定是素数。

Miller-Rabin测试：不断选取不超过p-1的基a(s次)，计算是否每次都有a^p-1 ≡ 1(mod p)，若每次都成立则p是素数，否则为合数。

还有一个定理，能提高Miller测试的效率：

二次探测定理：如果p是奇素数，则 x² ≡ 1(mod p)的解为 x = 1 || x = p - 1(mod p);

这个算法的正确率据说是3/4，所以多测几次。。

#include <cstdio>
#include <ctime>
#include <cstdlib>
#define LL long long

int n, ans;
const int S = 5;

//快速计算 a * b % p 
inline LL mod_mul(LL a, LL b, LL p)
{
    LL ret = 0;
    while(b)
    {
        if(b & 1) ret = (ret + a) % p;
        a = (a << 1) % p;
        b >>= 1;
    }
    return ret;
}

//快速计算 a ^ b % p 
inline LL mod_exp(LL a, LL b, LL p)
{
    LL ret = 1;
    while(b)
    {
        if(b & 1) ret = mod_mul(ret, a, p);
        a = mod_mul(a, a, p);
        b >>= 1;
    }
    return ret;
}

inline bool miller_rabin(LL p)
{
    if(p == 2 || p == 3 || p == 5 || p == 7 || p == 11) return 1;
    if(p == 1 || !(p % 2) || !(p % 3) || !(p % 5) || !(p % 7) || !(p % 11)) return 0;
    LL x, pre, u;
    int i, j, k = 0;
    u = p - 1; // 要求 u = p - 1 , (x ^ u) % p
    while(!(u & 1)) k++, u >>= 1; // 如果 u 为偶数则右移，k 记录位数 
    srand((LL)time(0));
    for(i = 0; i < S; i++)
    {
        x = rand() % (p - 2) + 2; // 在 [2, p) 中选取随机数
        if(!(x % p)) continue; // 如果不是互质
        x = mod_exp(x, u, p); // 先计算 (x ^ u) % p
        pre = x;
        for(j = 0; j < k; j++) // 把移位补
        {
            x = mod_mul(x, x, p); // 运用二次探测 
            if(x == 1 && pre != 1 && pre != p - 1) return 0;
            // 如果 (x ^ 2) % p == 1 但是 x != 1, x != p - 1, 就不是素数 
            pre = x;
        }
        if(x != 1) return 0; // 费马小定理
    }
    return 1;
}

int main()
{
    LL x;
    while(~scanf("%d", &n))
    {
        ans = 0;
        while(n--)
        {
            scanf("%lld", &x);
            if(miller_rabin(x)) ans++;
        }
        printf("%d
", ans);
    }
    return 0;
} 

拓展欧几里得

能求出 ax + by = gcd(a,b) 的一组解，还能顺带求出 gcd

正常版

inline int exgcd(int a, int b, int &x, int &y) {
  if (!b) {x = 1, y = 0; return a;}
  int r = exgcd(b, a % b, x, y);
  int t = x, x = y, y = t - a / b * y;
  return r;
}

极简版

inline int exgcd(int a, int b, int &x, int &y)
{
    if(!b){x = 1, y = 0; return a;}
    int r = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return r;
}

证明——其实就是推一下式子，

ax+by=gcd(a,b) ①

bx'+(a%b)y'=gcd(b,a%b) ②

已知x'和y'求x和y

∵gcd(a,b)=gcd(b,a%b)

∴ax+by=bx'+(a%b)y'

  ax+by=bx'+(a-a/b*b)y'

  ax+by=bx'+ay'-a/b*by'

  ax+by=ay'+b(x'-a/b*y')

∴x=y'，y=x'-a/b*y'

证毕，但要注意除法都是下取整的

快速幂

递归版

#define LL long long
inline LL pow(LL a, LL b, LL p)
{
    if(!b) return 1;
    LL ans = pow(a, b >> 1, p);
    ans = ans * ans % p;
    if(b & 1) ans = a * ans % p;
    return ans;
}

 

非递归版

#define LL long long
inline LL pow(LL a, LL b, LL p)
{
    LL ans = 1;
    while(b)
    {
        if(b & 1) ans = ans * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return ans;
}

膜意义下的逆元

若 a*x≡1(mod b)，且a和b互质，则称x为a的逆元，记作a^-1

逆元的求解方法：

1.快速幂

如果p是质数，那么 a 的逆元为 a^p-2，好像是通过什么费马小定理证明的，记住就行

2.扩展欧几里得

定义式可以转化为a*x+b*y==1，所以自然可以利用扩展欧几里得求啦！

3.线性递推

首先1^-1≡1(mod p)，

假设p=k*i+r，r<i，1<i<p，那么k*i+r≡0(mod p)

等式两边同时乘上i^-1，r^-1，式子就变成k*r^-1+i^-1≡0(mod p)

i^-1≡-k*r^-1

i^-1≡-(p/i)*(p mod i)^-1

这样便可以线性递推求解

#include <cstdio>
#define N 3000001

int n, p;
int inv[N] = {0, 1};

int main()
{
	int i;
	scanf("%d %d", &n, &p);
	puts("1");
	for(i = 2; i <= n; i++)
	{
		inv[i] = 1ll * -(p / i) * inv[p % i] % p;
		printf("%d
", (inv[i] + p) % p);
	}
	return 0;
}

卢卡斯定理——用于解决组合数的取膜问题

卢卡斯定理

C(n,m)%p=C(n%p,m%p)*C(n/p,m/p)%p

至于证明。。算了吧，背过得了

因为p不大

所以C(n%p,m%p)可直接求出，

但是如果n%p<m%p时返回0

C(n/p,m/p)可继续递归求解

当m==0时返回1

#include <cstdio>
#define N 200001
#define LL long long

int T;
int F[N], inv[N];

inline int C(int n, int m, int p)
{
	if(m > n) return 0;
	return (LL)F[n] * inv[m] % p * inv[n - m] % p;
}

inline int Lucas(int n, int m, int p)
{
	if(!m) return 1;
	return (LL)C(n % p, m % p, p) * Lucas(n / p, m / p, p) % p;
}

int main()
{
	int i, n, m, p;
	scanf("%d", &T);
	while(T--)
	{
		scanf("%d %d %d", &n, &m, &p);
		n += m;
		F[0] = F[1] = inv[0] = inv[1] = 1;
		for(i = 2; i <= p; i++)
			inv[i] = -(LL)(p / i) * inv[p % i] % p;
		for(i = 2; i <= p; i++)
			F[i] = (LL)F[i - 1] * i % p,
			inv[i] = (LL)inv[i] * inv[i - 1] % p;
		printf("%d
", (Lucas(n, m, p) + p) % p);
	}
	return 0;
}