GCD、LCM、拓展欧几里得、逆元

这个数论很早就知道，一直没学，昨天好不容易学了出来，很煎熬，怕以后忘了，把自己的想法写出来。
先帖一下GCD、LCM的模板：

int gcd(int a,int b){
	return  b == 0 ? a : gcd(b, a % b);
}
int lcm(int a,int b){
	return a  / gcd(a,b) * b;// 公式:a / gcd(a,b) * b / gcd(a,b) * gcd(a,b);
}

--------拓展欧几里得算法
主要是求二元一次方程的通解：ax + by = n 有整数解的条件是 n 可以整除 gcd(a,b)。如果有解的话（x0,y0），那么通解就是 x = x0 + (b/gcd)t ; y = yo - (a/gcd)t;
解：假设c为x的解的最小间距，此时d为y的解的间距，所以x=x0+ct，y=y0-dt（x0，y0为一组特解，t为任意整数）

  带入方程得：a*x0+a*c*t+b*y0-b*d*t=n,因为a*x0+b*y0=n，所以a*c*t-b*d*t=0，t不等于0时，a*c=b*d

  因为a,b,c,d都为正整数，所以用最小的c，d，使得等式成立，ac，bd就应该等于a，b的最小公倍数a*b/gcd(a,b),

  所以c=b/gcd(a,b)，d就等于a/gcd(a,b)。

欧几里得算法能求出ax + by = gcd(a,b)的一个解
贴欧几里得算法模板：

void extend_gcd(int a, int b, int &x, int &y){
	if(b == 0){
		x = 1,y = 0;//当b为0时，a = gcd(a,b) ，所以 a * 1 + b * 0 = gcd(a,b)
		return ;
	}
	extend_gcd(b, a % b, x, y);
	int tmp = x;//这里推导出了由上一层到这一层x与y的公式，证明可在B站看到
	x = y;
	y = tmp - (a / b) * y;
}

上面先求出了ax + by = gcd(a,b)的解，可利用它求出ax + by = n的通解：
1、判断ax + by = n 是否有整数解，有解的条件是n可以整除gcd(a,b)
2、用拓展欧几里得算法求ax + by = gcd(a,b)的一个解（x0,y0）
3、在ax0 + by0 = gcd(a,b)两边同时乘以n/gcd(a,b),得：
ax0 * n/gcd(a,b) + by0 * n/gcd(a,b) = n;
4、对照ax + by = n，得到它得一个解为(x1,y1):
x1 = x0 * n/gcd(a,b)
y1 = y0 * n/gcd(a,b)
通解为 : x = x1 + b/gcd(a,b) * t; y = y1 + a/gcd(a,b)
求正数最小解为 ：(b/gcd(a,b) + x1 %(b/gcd(a,b)) % b/gcd(a,b)
关于通解得证明方式可以借鉴关于poj2115的题解
关于这题可拿poj1061练手

#include <iostream>
#include <stdio.h>
using namespace std;
typedef long long ll;
ll abs(ll a){
    return a >= 0 ? a : -a;
}
ll gcd(ll a,ll b){
    return b == 0 ? a : gcd(b,a % b);
}
void extend_gcd(ll a, ll b, ll &x, ll &y){
    if(b == 0){
        x = 1, y = 0;
        return;
    }
    extend_gcd(b,a % b,x,y);
    ll tmp = x;
    x = y;
    y = tmp - (a / b) * y;
}
int main()
{
    ll x,y,m,n,l;
    while(~scanf("%lld%lld%lld%lld%lld",&x,&y,&m,&n,&l)){
    ll a = n - m, b = l, c = x - y;
    ll d = gcd(a,b);
    if(c % d != 0) printf("Impossible\n");
    else{
        ll p,q;
        extend_gcd(a,b,p,q);

        p = 1ll * p * c / d;
        q = 1ll * q * c / d;
        //printf("%lld %lld\n",p,q);
        long long t = ((b/d) + p % (b/d)) % (b/d);
       printf("%lld\n",t);
    }
    }
    return 0;
}

------同余与逆元
如果 a % m = b % m 则可以记作a ≡ b (mod m)
若有一个方程 ax ≡ 1(mod m)，则意味着 ax % m = 1，这个式子可以转换成 ax + my = 1，求出(x0,y0),则ax0 + my0 = 1,实则为ax0 % m = 1;则这种情况下称x为a模m的逆；
逆元的意义是:取模没有除法运算，所以将除变成乘倒数
设b的逆元为k，则bk % m = 1；
（a/b）mod m = (a / b) % m * ((bk) % m) = (a / b * bk) % m = ak % m
求逆元的模板：
1、拓展欧几里得

int mod_inverse(int a, int m){
	int x, y;
	extend_gcd(a,m,x,y);
	return (m + x % m) % m;
]

2、费马小定理 (注意条件)
我们知道，当p为素数，并且gcd(a,p)=1时，我们有ap−1≡1(modp)。那么我们就有ap−2×a≡1(modp)。所以逆元就是ap−2了。这里使用快速幂来计算

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

3、线性求逆元：
直接出公式，证明网上好多:A[i]=-(p/i)*A[p%i];注意，一般会爆int,看情况加1ll

inv[0]=inv[1]=1;  
    for(int i=2;i<N;i++)
        inv[i]=((MOD-MOD/i)*inv[MOD%i])%MOD;  

4、阶乘求逆元:
如果我们需要求0!到n!的逆元，对于每个元素都求一遍，就显得有点慢。（虽然exPower的时间快到可以认为是小常数。）
前面我们说了，逆元就可一看做是求倒数。那么不就有1/(n+1)!×(n+1)=1/n!。

int inv( int b, int p ) {
    int a, k;
    exPower( b, p, a, k );
    if( a < 0 ) a += p;
    return a;
}
void init( int n ) {
    Fact[ 0 ] = 1;
    for( int i = 1; i <= n; ++i ) Fact[ i ] = Fact[ i - 1 ] * i % Mod;
    INV[ n ] = inv( Fact[ n ], Mod );
    for( int i = n - 1; i >= 0; --i ) INV[ i ] = INV[ i + 1 ] * ( i + 1 ) % Mod;
    return;
}

拿hdu5976练手：

#include <iostream>
using namespace std;
const long long model = 1e9 + 7;
const int maxn = 1e5 + 50;
int sum[maxn],ff[maxn],inv[maxn];
int main(){

    inv[1] = 1;ff[1] =1;
    for(int i = 2; i < maxn; i++){
        sum[i] = sum[i - 1] + i;
        ff[i] = 1ll * ff[i - 1] * i % model;
        inv[i]= (model-model/i)*inv[model%i]%model;
    }
        int T; scanf("%d",&T);
        while(T--){
        int num; scanf("%d",&num);
        if(num <= 4){
            printf("%d\n",num);
            continue;
        }

        int l = 2, r = maxn, mid = (l + r)>>1;
        while(l  + 1 < r){
            if(sum[mid] > num) r = mid;
            else l = mid;
            mid = (l + r)>>1;
        }
        int k = num - sum[l];
        long long ans = 0;

        if(k == l) ans =1ll *  ff[l] * inv[2] % model * (l + 2) % model;
        else ans = 1ll * ff[l] * inv[l - k + 1] % model * (l + 1) % model;


        printf("%lld\n",ans);
    }
}