zoukankan      html  css  js  c++  java
  • 基本数论算法

    dalao博客,至少很好看。。

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

    最大公因数

    普通gcd

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

    1 inline int gcd(int x,int y)
    2 {
    3     return y == 0 ? x : gcd(y, x % y)
    4 }
    View Code

     

    二进制优化gcd

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

    1 inline int bsgcd(int x, int y)
    2 {
    3     if(x == y) return x;
    4     if(x < y) x ^= y ^= x ^= y;
    5     if(!(x & 1)) //x偶 y偶 gcd(x,y)=2*gcd(x/2,y/2),x偶 y奇 gcd(x,y)=gcd(x/2,y)
    6          return (!(y & 1)) ? bsgcd(x >> 1, y >> 1) << 1 : bsgcd(x >> 1, y);
    7      //x奇 y偶 gcd(x, y)=gcd(x,y/2)
    8     return (!(y & 1)) ? bsgcd(x, y >> 1) : bsgcd(y, x - y); 
    9 }
    View Code

    素数判定

    埃拉托色尼筛法

    O(nloglogn)

    不够优越,所以不学。

    线性筛

    模板题

    O(n)

     1 const int MAXN = 10000001;
     2 
     3 bool notpri[MAXN];
     4 int n, m, cnt, prime[MAXN];
     5 
     6 inline void get_prime()
     7 {
     8     int i, j;
     9     notpri[0] = notpri[1] = 1; // 1 表示不是质数 
    10     for(i = 2; i <= n; i++)
    11     {
    12         if(!notpri[i]) prime[++cnt] = i;
    13         for(j = 1; i * prime[j] <= n; j++)
    14         {
    15             notpri[i * prime[j]] = 1;
    16             if(!(i % prime[j])) break;
    17         }
    18     }
    19 }
    View Code

    Miller-Rabin 大素数判定

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

    模板题

    O(log32n)

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

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

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

    还有一个定理,能提高Miller测试的效率:

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

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

     1 #include <cstdio>
     2 #include <ctime>
     3 #include <cstdlib>
     4 #define LL long long
     5 
     6 int n, ans;
     7 const int S = 5;
     8 
     9 //快速计算 a * b % p 
    10 inline LL mod_mul(LL a, LL b, LL p)
    11 {
    12     LL ret = 0;
    13     while(b)
    14     {
    15         if(b & 1) ret = (ret + a) % p;
    16         a = (a << 1) % p;
    17         b >>= 1;
    18     }
    19     return ret;
    20 }
    21 
    22 //快速计算 a ^ b % p 
    23 inline LL mod_exp(LL a, LL b, LL p)
    24 {
    25     LL ret = 1;
    26     while(b)
    27     {
    28         if(b & 1) ret = mod_mul(ret, a, p);
    29         a = mod_mul(a, a, p);
    30         b >>= 1;
    31     }
    32     return ret;
    33 }
    34 
    35 inline bool miller_rabin(LL p)
    36 {
    37     if(p == 2 || p == 3 || p == 5 || p == 7 || p == 11) return 1;
    38     if(p == 1 || !(p % 2) || !(p % 3) || !(p % 5) || !(p % 7) || !(p % 11)) return 0;
    39     LL x, pre, u;
    40     int i, j, k = 0;
    41     u = p - 1; // 要求 u = p - 1 , (x ^ u) % p
    42     while(!(u & 1)) k++, u >>= 1; // 如果 u 为偶数则右移,k 记录位数 
    43     srand((LL)time(0));
    44     for(i = 0; i < S; i++)
    45     {
    46         x = rand() % (p - 2) + 2; // 在 [2, p) 中选取随机数
    47         if(!(x % p)) continue; // 如果不是互质
    48         x = mod_exp(x, u, p); // 先计算 (x ^ u) % p
    49         pre = x;
    50         for(j = 0; j < k; j++) // 把移位补
    51         {
    52             x = mod_mul(x, x, p); // 运用二次探测 
    53             if(x == 1 && pre != 1 && pre != p - 1) return 0;
    54             // 如果 (x ^ 2) % p == 1 但是 x != 1, x != p - 1, 就不是素数 
    55             pre = x;
    56         }
    57         if(x != 1) return 0; // 费马小定理
    58     }
    59     return 1;
    60 }
    61 
    62 int main()
    63 {
    64     LL x;
    65     while(~scanf("%d", &n))
    66     {
    67         ans = 0;
    68         while(n--)
    69         {
    70             scanf("%lld", &x);
    71             if(miller_rabin(x)) ans++;
    72         }
    73         printf("%d
    ", ans);
    74     }
    75     return 0;
    76 } 
    View Code

    拓展欧几里得

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

    正常版

    1 inline int exgcd(int a, int b, int &x, int &y) {
    2   if (!b) {x = 1, y = 0; return a;}
    3   int r = exgcd(b, a % b, x, y);
    4   int t = x, x = y, y = t - a / b * y;
    5   return r;
    6 }
    View Code

    极简版

    1 inline int exgcd(int a, int b, int &x, int &y)
    2 {
    3     if(!b){x = 1, y = 0; return a;}
    4     int r = exgcd(b, a % b, y, x);
    5     y -= a / b * x;
    6     return r;
    7 }
    View Code

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

    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'

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

    快速幂

    递归版

    1 #define LL long long
    2 inline LL pow(LL a, LL b, LL p)
    3 {
    4     if(!b) return 1;
    5     LL ans = pow(a, b >> 1, p);
    6     ans = ans * ans % p;
    7     if(b & 1) ans = a * ans % p;
    8     return ans;
    9 }
    View Code

     

    非递归版

     1 #define LL long long
     2 inline LL pow(LL a, LL b, LL p)
     3 {
     4     LL ans = 1;
     5     while(b)
     6     {
     7         if(b & 1) ans = ans * a % p;
     8         a = a * a % p;
     9         b >>= 1;
    10     }
    11     return ans;
    12 }
    View Code

    膜意义下的逆元

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

    逆元的求解方法:

    1.快速幂

    如果p是质数,那么 a 的逆元为 ap-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;
    }
    

      

  • 相关阅读:
    How do I change a .txt file to a .c file?
    [CQOI2007]余数求和
    CSP-J总结&题解
    【CSP游记S】
    [LuoguP1462]通往奥格瑞玛的道路
    归并排序——逆序对
    [NOIP 2011]选择客栈
    [二分图初步]【模板】二分图匹配,匈牙利算法
    [NOIP 2018]旅行
    黑魔法师之门 (magician)-并查集
  • 原文地址:https://www.cnblogs.com/zhenghaotian/p/6832636.html
Copyright © 2011-2022 走看看