【挑战程序设计竞赛】- 2.6 数论基础 (扩展欧几里得、素数筛法、快速幂)

2.6 数论基础 (扩展欧几里得、素数筛法、快速幂)

2.6.1辗转相除法

1. 求最大公约数

例题1：给定平面上两个格点P1(x1,y1), P2(x2,y2)。线段P1P2上，（除了P1 P2）有几个格点？（格点即整数点）
$-e^9<=x_1,x_2,y_1,y_2<=e^9$

【思路】：把(x1,y1)挪到(0,0)，P2(x2-x1,y2-y1)，实际上是找gcd(|x2-x1|,|y2-y1|)-1。
这里-1是因为算的是中间端点不是线段，线段是不用-1。

int gcd(int a, int b){
	return b==0?a:gcd(b,a%b);
}

辗转相除法复杂度：$O(log(max(a,b)))$

2. 扩展欧几里得

例题2：双六

有向前向后无限延伸的格子（…-3、-2、-1、0、1、2、3…）。
每个格子上都是整数。
其中0号格子是起点，1号格子是终点。
骰子上只有a \ -a \ b \ -b 四个数，问四个数各投几次能到达终点？无解输出1。
1<=a,b<=1e9

【思路】：实质是求解ax+by = 1。

若gcd(a,b)!=1, 设gcd(a,b)=t, a = k1t, b=k2t，那么ax+by = (ak1+bk2)t = 1显然不可能。

因此gcd(a,b) = 1。化简为$ax+by=gcd(a,b)$。

若已经求得

$bx+(a\%b)y=gcd(a,b)$

又因为 $a\%b=a-(a/b)\ast b$，因此

$bx+(a-(a/b)\ast b)y=gcd(a,b)$

化简为：$ay+b(x-a/b\ast y)=gcd(a,b)$

因此当前新的x’和y‘为：$x^{\prime }=y$, $y^{\prime }=x-(a/b)y=x-(a/b)x^{\prime }$

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

2.6.2 素数

1. 素数判定

例题1 判定n是不是素数，1 <= n <=1e9

【思路】素数，检查2到$\sqrt{n}$是否有约束

复杂度，每次$O(\sqrt{n})$

【容易bug的地方】
判断1的特殊情况！
因数分解要把另一个也加上
整数分解要把因子除尽。

2. 埃氏筛法

例题2 请问n以内有几个素数。1 <= n <=1e6

【思路】将2-n所有整数写下来，把所有2的倍数划去、3的倍数划去…

当前没画掉的最小数是素数，总是划去素数的倍数。

int prime[NMAX]; //存素数
bool isprime[NMAX];
int calprime(int n){
	int p = 0;
	for(int i=2; i<=n; i++) isprime[i]=true;
    for(int i=2; i<=n; i++){
        if(isprime[i]){
            prime[p++]=i;
            for(int j=2*i; j<=n; j+=i) isprime[i]=false;
        }
    }
    return p;//素数个数
}

3. 区间筛法（埃氏筛法+素数判定）

例题3 给定整数a和b，问区间 [a,b) 内有多少素数。
1<=a<b<=1e12, b-a<=1e6

【思路】b（及小于b）的最大质因子必然小于等于$\sqrt{b}$。

做2~$\sqrt{b}$ 的质数表，将a~b范围内的倍数筛去。

typedef long long ll;
const int nmax=1e6+30;
bool isprime_ab[NMAX]; //0 ~ b-a
bool isprime[NMAX]; //2 ~ sqrt(b)
int calprime(ll a, ll b){
	int p = 0;
	for(ll i=2; i*i<=b; i++) 
        isprime[i]=true;
    for(ll i=2; i*i<=b; i++){
        if(isprime[i]){
			//筛2~sqrt(b)
            for(ll j=2*i; j*j<=b; j+=i) 
                isprime[i]=false;
            ll st = max(2LL, (a+i-1)/i)*i;
            for(ll j = st; j<=b; j+=i)
                isprime_ab[j-a]=false;
        }
    }
    int cnt = 0;
    for(ll i=a; i<=b; i++){
        if(isprime_ab[i-a]) cnt++;
    }
    return cnt;//素数个数
}

2.6.3 取模

为了不让其溢出整数（2147483647）。

a是负数是a%m也是负数，记得用a%m+m 给他纠正到0~m-1。

加减乘都非常自由，但是除法不好模。

2.6.4 快速幂计算

例题1： uva10006 carmichael number

C数：对于任意的1<x<n，都有x^n=x(mod)n成立的合数。

对于给定的n，判断其是否是C数。2<n<=65000

【思路】朴素计算要$O(n^2)$。快速幂。

$x^n=(x^2{)}^2$，想到把n化成若干2的幂次和。

$x^n=x^{2^{k_1}}{x}^{2^{k_2}}{x}^{2^{k_3}}...$

因此依次求2的幂次即可。

$x^{22}=x^{2^4}\ast x^{2^2}\ast x^{2^1}$，也是因为22写成二进制be like：10110

快速幂板子：

typedef long long ll;
ll ppow(ll x, ll n, ll mod){
    ll res = 1;
    while(n){
        if(n & 1) res = res*x % mod;
        x = x*x%mod;
        n>>=1;
    }
    return res;
}