有点复杂的同余算法

目录

• 单变量同余方程组
  • 膜数两两互质：CRT
  • 膜数不两两互质（挺多人称这个算法为ExCRT）
• 组合数取模
  • 模数为质数
  • 模数不为质数（经常被称为 ExLucas 的算法）

---

单变量同余方程组

问题就是如何解这个方程：

[egin{cases} x equiv a_1 (mod m_1) \ x equiv a_2 (mod m_2) \ vdots \ x equiv a_n (mod m_n) \ end{cases} ]

膜数两两互质：CRT

思想时对每个单独的方程 (xequiv a_i(mod m_i)) 构造一个特殊的 “元”， 称为 (o_i)， 使得这个 (o_i) 在膜 (m_i) 意义下是 (a_i)， 且对于 (forall j eq i)， (o_i equiv 0(mod m_j))， 最后 (x) 就是 (sum_i o_i) 了。

这个 (o_i) 其实挺好构造， 设 (M = prod_i m_i)， 那么设 (z_i = M/m_i)， 再设 (z_i) 在膜 (m_i) 意义下的逆元是 (v_i)， 那么 (o_i) 就等于 (a_i*z_i*v_i) 了。

最后就是 (x) 在过程中的取模问题： 满足了多少个方程， 就用恒等式合并的方法， 对 (m) 做个 (lcm)， 最后模数就是这个 (lcm)。（这是因为若干个同余恒等式的合并是可以拆分成原来的这些恒等式的）

代码大概长这样

int n;
long long ri[MN], riri[MN];//ri意为right， riri意为right的right， 分别方程的 a 和 m
long long CRT() {
	long long M=1ll, res=0ll;
    for(int i=1;i<=n;++i) M *= riri[i];
    for(int i=1;i<=n;++i) {
        long long Mi = M/riri[i];
        long long iMi = inv(Mi, riri[i]);
        res += Mi*iMi%M * ri[i] % M;
        res %= M;
	}
    res = (res%M+M)%M;
    return res;
}

膜数不两两互质（挺多人称这个算法为ExCRT）

考虑做增量（以同余恒等式的合并与拆分为原理）， 对于前 (k) 个方程， 算出在膜 (LCM_{i=1}^k m_i) 意义下满足前 (k) 个方程的一个解 (x_k)， 然后用这个解找出 (x_{k+1})。

首先通解是： (x_k + q*LCM_{i=1}^k m_i)， 其中 (q) 是整数。

那么可以解这个方程： (x_k+ q*LCM_{i=1}^k m_i equiv a_{k+1} quad(mod m_{k+1})) 。（注意这个方程的变量是 (q)）

然后再给 (x_k) 做增量变成 (x_{k+1}) ， 当然是在膜 (LCM_{i=1}^{k+1} m_i) 意义下， 就行了。

代码大概长这样：

int n;
int ri[MN], riri[MN];
long long ExCRT() {
	long long X=ri[1], LCM=riri[1];
    for(int i=2;i<=n;++i) {
        long long c = ((a[i]-x)%m[i]+m[i])%m[i];
        long long x,y,d;
        exgcd(x,y,LCM,m[i],d); // 这里若 d 不整除 c 就无解
        x = slow_mul(x, c/d, m[i]/d); // 由于 c 已经被修正成正数所以不用担心 slow_mul 这个函数会出问题
        X += x*LCM;
        LCM = LCM/d * m[i];
        X = ((X%LCM)+LCM)%LCM;
	}
}

组合数取模

模数为质数

[inom{n}{m} = inom{lfloor n/p floor}{lfloor m/p floor}inom{nmod p}{mmod p} ag{$mod p$} ]

不会证弃了弃了。

模数不为质数（经常被称为 ExLucas 的算法）

首先将模数分解为若干个质数的幂次， 然后用 CRT 合并。

那么对于 (inom{n}{m} mod p^k) ：

首先设 (num_p(x)) 表示 (lfloordfrac xp floor)， 意即 (x!) 中 (p) 因子的个数。

然后 (inom{n}{m} mod p^k) 就等于

[dfrac{dfrac{n!}{p^{num_p(n)}}}{dfrac{m!}{p^{num_p(m)}} dfrac{(n-m)!}{p^{num_p(n-m)}}} *p^{num_p(n)-num_p(m)-num_p(n-m)} mod p^k ]

这是为了使得除法有逆元。

插一句， (num_p(n)-num_p(m)-num_p(n-m) ge 0)， 具体看具体数学吧……我不会， 第三章我跳了。

求解 (dfrac{n!}{p^{num_p(n)}} mod p^k) ：

[egin{aligned} n! & = 1*2*cdots*n \ & = (p*2p*cdots*lfloor n/p floor p) * 1 * cdots * unknown \ & = p^{lfloor n/p floor} * (1*2*cdots*lfloor n/p floor) * 1 * cdots * unnkown end{aligned} ]

去掉 (p^{lfloor n/p floor}) 就是 (dfrac{n!}{p^{num_p(n)}})， 这样， 就可以把问题缩小。（不过后面小问题的 p 不能去掉了）

对于后面那一串 (1*cdots *unknown)， 由于 (xequiv x+p^k(mod p^k))， 可以循环节算。

代码长这样：

#define li long long
li ksm(li a, li b, li p) {
	li res = 1ll;
	for(;b;b>>=1, a=(a*a)%p)
		if(b&1) res = (res*a)%p;
	return res%p;
}
void exgcd(li a, li b, li &x, li &y) {
	!b ? x=1,y=0 : (exgcd(b,a%b,y,x), y-=x*(a/b));
}
li inv(li a, li p) {
	// ax + py = 1;
	li x, y;
	exgcd(a, p, x, y);
	x = ((x%p+p)%p);
	return x;
}

li fac(li n, li p, li pk) {
	if(!n) return 1ll;
	li res = 1ll;
	for(li i=0; i<pk; ++i) if(i%p) res = res*i % pk;
	res = ksm(res, n/pk, pk);
	for(li i=0; i<=n%pk; ++i) if(i%p) res = res*i % pk;
	return res * fac(n/p,p,pk) % pk;
}

li C(li n, li m, li p, li pk) {
	if(!m) return 0ll;
	li f1 = fac(n,p,pk), f2=fac(m,p,pk), f3=fac(n-m,p,pk), cnt = 0ll;
	for(li i=n;i;i/=p) cnt += i/p;
	for(li i=m;i;i/=p) cnt -= i/p;
	for(li i=n-m;i;i/=p) cnt -= i/p;
	return f1 * inv(f2,pk) % pk * inv(f3,pk) % pk * ksm(p,cnt,pk) % pk;
}

int tot;
li ri[10001], riri[10001];
li CRT() {
	li M = 1ll, res = 0ll;
	for(int i=1; i<=tot; ++i) M *= riri[i];
	for(int i=1; i<=tot; ++i) {
		int nd = M/riri[i];
		res += ri[i] * nd % M * inv(nd,riri[i]) % M;
		res %= M;
	}
	return ((res%M+M)%M);
}
li exlucas(li n, li m, li p) {
	for(li i=2; i<=sqrt(p); ++i) {
		li md = 1ll;
		while(p%i==0) p/=i, md*=i;
		if(md>1) ri[++tot] = C(n,m,i,md), riri[tot] = md;
	}
	if(p>1) ri[++tot] = C(n,m,p,p), riri[tot] = p;
	return CRT();
}