数论代码总结

欧几里得:求a,b的最大公约数。

1
2
3
4
int gcd(int a,int b){
    if(b==0)return a;
    return gcd(b,a%b);
}

扩展欧几里得:(求a,b的最大公约数和不定方程ax+by=d(=gcd(a,b))的一对解)
裴蜀定理
对于任意整数 a , b,存在无穷多组整数对 (x,y) 满足不定方程 ax+by=d ,其中 d=gcd(a,b)
在求 d=gcd(a,b)的同时,可以求出关于x,y的不定方程ax+by=d的一组整数解。
考虑递归计算:假设已经算出了(b,a \bmod b)的一对解(x_0,y_0)满足bx_0+(a \bmod b)y_0=d
可以得到bx_0+(a – b \left \lfloor \frac{a}{b} \right \rfloor)y_0=d
整理可得b(x_0-\left \lfloor \frac{a}{b} \right \rfloor y_0) +ay_0=d
ay_0+b(x_0-\left \lfloor \frac{a}{b} \right \rfloor y_0)=d
x=y_0,y=x_0-\left \lfloor \frac{a}{b} \right \rfloor y_0
假设先行swap(x_0,y_0),则有x=x_0,y=y_0-\left \lfloor \frac{a}{b} \right \rfloor x_0

1
2
3
4
5
6
7
8
9
int exgcd(int a,int b,int &x,int &y){
    if(b==0){
        x=1;y=0;
        return a;
    }
    int d=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return d;
}

乘法逆元
方法1:利用扩展欧几里得求逆元,参见:https://renjikai.com/luogu-p1082/
适用条件:被求数a和模数b互质。
方法2:利用费马小定理(快速幂)求逆元,参见:https://renjikai.com/luogu-p1313/
适用条件:模数为质数。(保证了被求数a和模数b互质。)
方法3:不用记复杂的线性递推公式,就可以用O(N)复杂度推出1到n的所有乘法逆元(模意义下,在后面不再赘述)。除此之外,在此过程中求出1!n!的阶乘及其逆元。
前提结论:易知\prod_{i=1}^{n}\ inv(i)=inv(n!)
乘法逆元定义:inv(i)\cdot i\equiv 1
令数组A_i=i!,容易求出数组A。令数组B_i=inv(i!),那么数组B应该怎么求呢?
用方法1、2容易求出inv(A_n),即inv(n!),即求出了B_n
由前提结论知,B数组的递推公式为B_i=B_{i-1} \cdot inv(i),两边同乘i,则有B_i \cdot i =B_{i-1} \cdot inv(i)\cdot i,由乘法逆元定义有B_{i-1} = B_i \cdot i,只需回推,即可求出由i!的乘法逆元。
B_n=\prod_{i=1}^{n}\ inv(i)=inv(n)\cdot inv(n-1) \cdot inv(n-2) … \cdot inv(1)A_{n-1}=(n-1)!=(n-1)\cdot(n-2)…\cdot(1),因此两式相乘则有inv(n)=B_n\cdot A_{n-1},即求出了1到n的所有乘法逆元,是线性复杂度O(N),可能常数大一点,但是好记啊。
C++伪代码如下:(不一定完全正确,请注意识别)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
#define MOD
#define ll long long

ll quickpow(ll a, ll b);
ll fInv(ll x) {
    return quickpow(x, MOD - 2) % MOD;
}
ll n;
ll fact[n + 2];
ll invFact[n + 2];
ll inv[n + 2];
void gen() {
    fact[1] = 1;
    for (ll i = 2; i <= n; i++) {
        fact[i] = (fact[i - 1] * i) % MOD;
    }
    invFact[n] = fInv(fact[n]);
    for (ll i = n - 1; i >= 1; i--) {
        invFact[i] = invFact[i + 1] * (i + 1) % MOD;
    }
    inv[1] = 1;
    for (ll i = 2; i <= n; i++) {
        inv[i] = invFact[i] * fact[i - 1] % MOD;
    }
}

发表评论

您的电子邮箱地址不会被公开。 必填项已用*标注