最近有点颓废啊,写篇blog振作一下…(不过没图的数论blog真是不对我胃口)
emmm …首先介绍一下这是一篇关于数论中较为重要 (主要可能经常要用到) 的一个知识分支—逆元
相信你听到数论之后可能鼠标就想往网页上的 X 键上移了,但是还是劝你看一看, 总能有点收获的吧
(反正我觉得自己讲的相对网上其他 blog 还算是蛮清楚的了)
解释一下逆元:
首先同学们小学的时候都学过除法吧(这个问题有点奇怪,滑稽)?
那么你一定知道,“除以一个数等同于乘以这个数的倒数”对吧?
所以逆元是什么东西呢? 首先这里有个式子: (a/b) %p ,这个式子的答案怎么求?
没错,暴力求是一种方法,但是当 b 非常大的时候呢 ? 这个时候就要用到逆元了
这里有一句口诀 “除以一个数再取模等同于乘以这个数的逆元再取模 ”(自创)
什么意思? 设 inv[b] 是 b 的逆元, 那么 (a/b) %p = (a*inv[b]) %p
这下就你能大概理解逆元是个什么东西(以及它大概是用来干什么的)了吧 ?
那么逆元怎么求 ? 首先这里得纠正一点, 我们不能直接说一个数的逆元是多少,
应该这么说: 一个数 x 在模 p 的条件下的逆元是多少.(这点概念还是蛮重要的)
其次,我们不难得知一个数的逆元有多个,但是我们只需要求得一个数的最小正整数逆元就行了
另外,一个数 x 在模 p 的条件下不一定有逆元, x 关于 p 的逆元存在 当且仅当 x 和 p 互质
这里有一个推导: (设 a 为 x 的逆元, b为任意整数)
x * a ≡ 1 (mod p) = 将p连入式子=> x*a = 1 -b*p => x*a+b*p=1
那么我们就得到了一组新方程: x*a+b*p=1
若x 和 p不互质, 则 x和p 存在公约数 d=gcd(x,p)>1 ,
提取出d, 得到: d(x/d*a+p/d*b) = 1
移项得到: (x/d * a+ p/d *b) = 1/d
易知x , p 能整除 d, 所以括号内定为整数,
又因为d > 1 ,等式右边必为真分数, 等式无解
证毕
那么说了半天,逆元到底怎么求呢 ? 别着急,往下看—-
求逆元有好几种方法, 但是在这里我们只讨论两种方法(一个偷懒用的 ,一个不容易翻车的),
不过在此之前我们先来讲讲费马小定理 (别害怕,这玩意儿知道怎么用就好,不会叫你推的)
这是费马小定理的原始状态: a^p ≡ a (mod p) ,其中a , p 互质 [我也没有查过这玩意儿怎么来的,不过会用就好]
然后, 原式 =移项 => a^p - a ≡0 (mod p) =提取 a => a(a^(p-1) - 1) ≡0 (mod p)
=两边同除以 a => a^(p-1)-1 ≡0 (mod p) =再次移项 => a^(p-1) ≡1 (mod p)
于是 费马小定理的结论就产生了(注意,这个结论很重要),求逆元时经常要用(主要偷懒求逆元时必用)
好了,现在我们再来考虑逆元的求法
一、 蒙哥马利快速幂模(偷懒求逆元大法):
首先听到这么高大上(主要是长)的名字,相信你的内心一定是崩溃的,仿佛见到了高斯定理或者欧拉定理一样(别着急啊,扩欧等下才讲)
但其实蒙哥马利快速模这个算法并没有那么高大上,它的局限性很大,只有在 p 是质数的情况下才可以使用
首先我们设 inv_a 是 a 的逆元 那么根据定义, inv_a * a ≡1 (mod p)
再根据 费马小定理 a^(p-1) ≡1 , 易得 inv_a * a ≡a^(p-1) (mod p)
移项,得: inv_a ≡ a^(p-2)
于是我们得到了快速幂模算法的一个前提条件: inv_a ≡ a^(p-2) (mod p)
下面附上代码(自己即兴手打):
[这里是位运算处理的方法,当然也有递归求的方法,那个代码短一些,不过这个代码我打的顺手]
inline ll quick_pow(ll x,rint p){
ll res=1;
while(p){
if(p&1) res=(res*x)%mod;
x=(x*x)%mod, p>>=1;
}
return res;
}
inline inv(ll a){
inv_a=quick_pow(a,mod-2);
return inv_a;
}
下面我们不着急讲第二种算法,先来讲讲这种快速模算法的应用:
求 C(n,r) % p 的值 , 其中 p 是大质数(如1e9+7) ,那么最后的答案是: 原式 = n! * (r! * (n-r)!)^(p-2) %p
推导过程:
C(n,r) (mod p) = n! / ( r! * (n-r)!) (mod p) =多乘一个 1 => n! / (r! * (n-r)!) * 1 (mod p) ----1式
∵ p是质数,再根据费马小定理可得: (r! * (n-r)!)^(p-1) ≡1 (mod p) ----2式
再将2式带入1式中得到: C(n,r) ≡n! /(r! * (n-r)!) * (r! *(n-r)!)^(p-1) (mod p)
指数相消得到 C(n,r) ≡n! * (r! * (n-r)!)^(p-2) (mod p)
或者我们也可以换元一下, 令 tmp = r! * (n-r)! ,那么:
C(n,r)%p = ( n! / tmp )%p = ( n! / tmp ) * 1 %p
= ( n! * tmp^-1 ) * tmp^(p-1) %p
= ( n! * tmp^(p-2) )%p = 上述的式子
代码实现:
//by Judge
const ll mod=; //1e9+7或者其他题目给定的大质数
inline ll quick_pow(ll x,rint p){
ll res=1;
while(p){
if(p&1) res=(res*x)%mod;
x=(x*x)%mod, p>>=1;
}
return res;
}
inline C_mod(rint n,rint m){
ll a=1,b=1;
for(rint i=2;i<=n<<1;++i)
a=(a*i)%mod;
for(rint i=2;i<=m;++i)
b=(b*i)%mod;
for(rint i=2;i<=n-m;++i)
b=(b*i)%mod;
return (a*quick_pow(b,mod-2))%mod;
}
然后是上文提到过的扩欧
二、扩展欧几里得算法(一般我们都讲扩欧) :
这个东西讲起来真是烦,一开始我还一直以为扩欧就是在求逆元…其实扩欧是求同余方程的解…
首先如果你不是萌新的话(如果是个老司机),那么你应该看到过ex_gcd这样(恶心)的字眼,
这玩意儿其实就是扩欧…(ex 表示扩展, gcd 最大公约数是可以用欧拉求的嘛)
不废话了,进入正题, 首先扩欧公式如下: ax+by = (a,b) [注意这里是等于]
[(a,b) 代表 gcd(a,b), a 和 b 是已知数, 要求的是 x 和 y 的一组解]
那么很容易可以知道, 当 a 和 b 互质的情况下,原来的等式就变成了: ax+by = 1
然后我们怎么用这玩意儿求逆元呢? 我们先考虑求逆元的式子: a * inv[a] ≡1 (mod p)
其中inv[a] 同上表示 a 的逆元, 现在我们令 x = inv[a] , b = p;
那么同余的式子就成了: ax ≡1 (mod b) , 然后我们把 b 加入式子中,就成了这样:
ax+by = 1 (这里我就不慢慢推了,看博客也是要动动脑子的嘛 )其中 y 和 x 一样,也是未知数
emmmm还是讲讲, ax ≡1 (mod b) => ax = by+1 (然后就把同余去掉了,变成了等式)
移项=> ax+by = 1 (由于 y 是不确定的,而且我们求的主要是 x 所以 y 的正负性没什么影响)
这样你应该就能看懂了, 扩欧公式中求解的 x 就是 a 的逆元(一切好像都非常的显然了)
ok, 我还是没有讲 扩欧中的 x 和 y 怎么解 (这个才是真正麻烦的地方,解释起来很麻烦) 对吧?
哎, 式子先列出来再说, 我再斟酌斟酌…
ax+by =(a,b) =1 , bX+(a%b)Y =(b,a%b) =1 (注意大小写,并且显然这里要满足 a b 互质的条件)
=> ………… => x=Y , y=X-(a/b)*Y (其中 X 和 Y 已知)
好了,自己推导吧! (试图想混过去的样子) 咳咳… 首先 ax+by = bX+(a%b)Y 对吧? (两式合并)
然后我们把 a%b 转化一下,变成: ax+by = bX+(a-b*(a/b))Y (显然这里a/b是取整的)
=> ax+by = bX+aY-b(a/b)Y =提取公因式 a,b => a(x-Y) = b(X-(a/b)Y-y)
那么很显然, 我们令 x = Y , 令 y = X-(a/b)Y 就能轻松的构造出一组 x 和 y 的解了
然而现在又出现了一个问题: X 和 Y 怎么求? 其实我们同上的方法求出 X 和 Y 就行了:
bX+(a%b)Y=(a%b)X’ + (b%(a%b))Y’ 如果你还不能理解, 那么我们令 p = b, q = a%b ,
那么原式 => pX+qY=qX’+(p%q)Y’ 这不就和之前的式子差不多了吗?
现在你应该很容易就能看出来这是在用递归的方法求解原式的解
但是递归总是有有边界的,不能永远递归下去对吧(非常显然)
那么这个递归式子的边界是什么呢? emmmm 提示一下, 这个边界和 gcd 的边界有点关系
对, 边界就是 b 等于 0 的时候, 我们不能再递归下去了,那么原式就成了:
1x+0y=1 [a此时为 (a,b) ,b为 0 ] 这应该很好解吧? 我们令x=1,y=0 就行了(y=0是为了…方便嘛)
然后递归求解同上,下面附上代码(也是即兴手打的)
[下面常规版的,偷懒版和这个其实本质上没有区别,只是简化了一下代码,那么这里就不附上了,有兴趣的同学可以自行百度]
#define ll long long
const int mod=;
void ex_gcd(ll a,ll b,ll &x,ll &y){
if(!b){ x=1,y=0; return ; }
ex_gcd(b,a%b,x,y);
ll t=x; x=y,y=t-(a/b)*y;
}
inline ll inv(ll a){
ll inv_a,y;
ex_gcd(a,mod,inv_a,y);
return inv_a;
}
另外这里补充一下线性推逆元的方法 (inv[i] 表示 i 的逆元 , fac[i] 表示 i 的阶乘 , p为模数)
三、线性推逆元
1.求 1! ~ n! 的逆元:(因为我觉的阶乘的逆元反倒比较好推所以就先讲)
公式: inv[i] ≡inv[i+1]*(i+1) (mod p)
推导过程:
fac[i] * inv[i] ≡1 (mod p)
fac[i+1] * inv[i+1] ≡1 (mod p) => fac[i]* (i+1) * inv[i+1] ≡1(mod p)
由 同余的除法原理可得 : inv[i] ≡inv[i+1] * (i+1)
推导完毕
2.求 1 ~ n 的逆元(这个还稍微麻烦些)
公式:inv[i] ≡inv[p%i] * (- p/i) (mod p)
推导过程:
令 s = p/i , t = p%i , 则有: s*i + t = p (显然)
然后 s*i + t ≡ 0 (mod p)
移项得 t ≡ -s*i (mod p)
同除以 t * i 得 t / (t*i) ≡ -s*i / (t*i) (mod p)
将除法转化为乘法 => inv[i] ≡ -s * inv[t] (mod p)
于是将 s=p/i , t=p%i带入就可以得到: inv[i] ≡ inv[p%i] * (-p/i) (mod p)
推导完毕
四、Lucas定理(求大组合数取模)
再讲讲lucas定理这个东西(扩展lucas就不讲了,因为不大会…咳咳,然后也不怎么会用到吧)
基本公式: C(n,m) ≡ C(n/p,m/p)*C(n%p,m%p) (mod p)
(也就是: C(n,m)%p=C(n/p,m/p)*C(n%p,m%p)%p )
适用范围:n 和 m 非常大,而模数 p 比较小的情况 (偷懒)
运用:
将组合数中的 n 和 m 不断除以 p , 同时用除 p 的余数做组合数,累计入答案。
即,不断调用基本公式,递归求解,直至 m 等于 0 .
用上述例子就是:令C(n/p,m/p) ≡ C( (n/p)/p,(m/p)/p ) * C( (n/p)%p,(m/p)%p ) (mod p)
然后把 C(n/p,m/p) 求得的解 带入基本公式, 求出 C(n,m)%p 的值
证明及推导:
以下证明推导源自
Lucas定理——推导及证明
要证:C(n,m)%p=C(n/p,m/p)*C(n%p,m%p)%p
即证:C(n,m)=C(n/p,m/p)*C(n%p,m%p)
证明条件:已知p是素数,n、m、p为整数。
1.由二项式定理得:
,且
2.由费马小定理得:
(1)式 :(x + 1) ^p ≡ x+1 (mod p)
[ 由x ^p ≡ x (mod p) , x 用 (x+1) 代替得到]
(2)式 :x ^p + 1 ≡ x+1 (mod p)
[ 由x ^p ≡ x (mod p) , 等式两边同时+1得到]
合并 1 、2 两式,得到:
(3)式 :(x + 1) ^p ≡ x ^p + 1 (mod p)
3.由 n = n/p *p + n%p . 则:
(x+1) ^p ≡ (x+1) ^(n/p *p) * (x+1)^(n%p) (mod p)
带入三式: (x+1)^p ≡ (x^p +1)^(n/p)*(x+1)^(n%p) (mod p)
将上式由二项式公式可转化为:(可以把式子写一下再看)
∑(n,z=0) C(n,z)* x^z ≡ ∑(n/p,i=0) C( n/p,i )* x^(p * i )* ∑(n%p,j=0) C(n%p,j)* x^j (mod p)
任意一个Z,必存在一个
i , j 满足:
x ^ z = x ^ (p* i) * x^ j (即满足:n = n/p *p + n%p),
当且仅当:
i=z/p,j=z%p 时成立
此时,
C(n,m)=C(n/p,m/p)*C(n%p,m%p) 成立
证毕
代码(同上随意手打):
define ll long long
ll mod;
inline ll quick_pow(ll x,ll p){
ll ans=1;
while(p){
if(p&1) ans=ans*x%p;
x=x*x%p, b>>=1;
}
return ans;
}
inline ll C(ll n,ll m){
ll cn=1,cm=1,res=1;
if(nreturn 0;
if(n==m) return 1;
if(m>n-m) m=n-m;
for(ll i=0;imod;
cm=(cm*(m-i))%mod;
}
res=(cn*quick_pow(cm,mod-2))%mod;
return res;
}
ll Lucas(ll n,ll m){
ll ans=1;
while(n && m && ans){
ans=(ans*C(n%mod,m%mod))%mod;
n/=mod, m/=mod;
}
return ans;
}
如果你看完这些不过瘾,给你几个网站用来生啃代码:
四法求逆元
逆元的几种求法
最后的最后,逆元的这些东西…随便学学就好,到时候用得到(也可能用不到)
于是Judge的课堂又这样水过去了一节…
欢迎下次光临!
_(:з」∠)_