扩展欧几里得算法及其应用

扩展欧几里得算法不仅可以算出a,b的最大公约数 ，还可以给出不定方程ax+by = gcd(a,b)的一组解。这里a,b的正负不影响答案
对于任意的不定方程ax+by = d。若d%gcd(a,b) 不为0，则该不定方程无解。
证明
解不定方程：ax1 + by1 = gcd(a,b)
由欧几里得算法得 gcd(a,b) = gcd(b,a%b)
gcd(b,a%b) = bx2 + (a%b)y2
a%b = a - a/b
将a%b带入得
bx2 + (a - a/b)y2 = ay2 + b(x2 - a/b)
将 该式与原来的式子比较得：
ay2 + b(x2 - a/b) = ax1 + by1
x1 = y2
y1 = x2 - (a/b)y2
当b = 0时（递归的边界），x = 1,y = 0,回溯回利用上式求得x1,y1
当我们求出一个解时，有时需要求出通解或最小整数解。
ax + by = d
ax + b * y + n * a * (b/d) - n * (a / d ) b = d
(x+n*(b/d)) * a + (y-n*(a/d)) * b = d
所以通解为 x = x1 + n * (b/d)，y = y1 - n * (a/d)
所以要求最小整数解就很简单了，要求x的最小整数解，只要模b/d即可。若等式右边是gcd的倍数时，求出结果后先乘上倍数，再取模（右边不管是什么，取模都成立），当b/d是负数时，求最小整数解直接取绝对值就好。

//扩展欧几里得，返回gcd(a,b),并且可以得到不定方程 ax + by = gcd(a,b)的一组解 
//通解为x = x1 + n * (b/d)，y = y1 - n * (a/d)
//若等号右边为t,则通解为x= x1 * t /d + n * (b/d) 
//求最小整数解时直接取模（b/d）即可。
int ext_gcd(int a,int b,int &x,int &y)
{
	int d = a;
	if( !b )
	{
		x = 1;
		y = 0;
	}else
	{
		d = ext_gcd(b,a%b,y,x);
		y -= a / b * x;
	}
	return d;
} 

利用欧几里得算法计算逆元
x为b的逆元，则bx = 1(mod c) = 1 + cy
bx + cy = 1
若gcd(b,c) = 1,则不定方程有解，b的逆元存在，利用扩展欧几里得计算出x即可

void ext_gcd(int a,int b,int &x,int &y)
{
	if( !b )
	{
		x = 1;
		y = 0;
	}else
	{
		ext_gcd(b,a%b,y,x);
		y -= a / b * x;
	}
} 

int inv(int t,int p)  //求t关于p的逆元
{
	int x,y;
	ext_gcd(t,p,x,y);
	x = ( x % p + p ) % p;
	return x; 
}

求1~n中的数的逆元
利用逆元的线性递推关系就可以在O(n)的复杂度下求出1~n中的数关于p的逆元
推导：
设 p = k * i + r;(k = p / i , r = p % i)
k * i + r = 0(mod p)
两边同时乘以i^-1, r^-1
k * r^-1 + i^-1 = 0 (mod p)
i^-1 = -k * r^-1
将k与r带入得
i^-1 = -(p/i) * (p % i)^-1 mod p
由于p%i小于i，所以我们在计算i^ -1时，(p%i)^-1已经被算出来了
由于需要保证i^ -1>0，所以我们在右边加上p
i^-1 = p - (p/i) * (p % i)^-1 mod p

#include <cstdio>
using namespace std;
typedef long long ll;

ll inv[3000005];

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

求1~n阶乘的逆元
对于阶乘的逆元也是有递推公式的，f[i+1] * inv[i+1] =1 = f[i] * (i+1*inv[i+1]) = f[i] * inv[i]，所以inv[i] = inv[i+1] * (i+1)。

f[0]=1;
for(int i=2;i<=MAXN;i++)fac[i]=fac[i-1]*i%mod;
inv[MAXN]=q_pow(fac[MAXN],mod-2);
for(int i=MAXN-1;i>=0;i--)inv[i]=inv[i+1]*(i+1)%mod;