从exgcd到exCRT
从最基础的开始。
1.gcd
这个不用说了吧……\(gcd(a,b) = gcd(b,a\%b)\),这个很显然。
2.exgcd
这玩意可以用来求形如\(ax+by = gcd(a,b)\)的不定方程的一组特解。
首先来证明一下为什么一定是有解的。
因为我们是像上面的gcd一样递归解决问题的,所以当\(b = 0\)时,我们返回a,此时方程必然有一个特解\(x = 1,y = 0\)成立。
我们假设现在已经求出了一组解\(x_1,y_1\),我们要求下一组解\(x_2,y_2\)
有\(ax_1 + by_1 = gcd(a,b) = gcd(b,a\%b) = bx_2 + (a\%b)y_2\)
\(a\%b = a - \lfloor\frac{a}{b}\rfloor * b\)
\(ax_1 + by_1 = bx_2 + ay_2 - b\lfloor\frac{a}{b}\rfloor y_2\)
合并同类项,得到\(ax_1 + by_1 = b(x_2 - \lfloor\frac{a}{b}\rfloor y_2) + ay_2\)
所以得到\(x_1 = y_2,y_1 = x_2 - \lfloor\frac{a}{b}\rfloor y_2\)
这样推下去一定可以得到方程的一组解,所以形如\(ax+by = gcd(a,b)\)的方程是一定有解的。
同理我们可以推出,对于不定方程\(ax+by = c\),此方程有解的充要条件是\(gcd(a,b) | c\),我们可以先求出\(ax+by = gcd(a,b)\)的特解,然后把解同时乘以\(\frac{c}{gcd(a,b)}\)就得到了解。
这玩意还可以用来求逆元(前提是这个数和模数必须互质),对于式子\(ax \equiv 1 (mod\ p)\),改写成\(ax - py = 1\),求不定方程解即可。\((gcd(a,p) = 1)\),必然有解。
实现方法如下。
int exgcd(int a,int b,int &x,int &y)
{
if(!b){x = 1,y = 0;return a;}
int d = exgcd(b,a%b,y,x);
y -= a / b * x;
return d;
}
3.CRT
中国剩余定理(CRT),是用于解同余方程的一种方法。
同余方程就是给定多个形如\(x \equiv a_i (mod\ b_i)\)的方程,其中\(b_i\)两两互质。求x的最小正整数解。
方法的思想个人认为其实是构造。
就是对于一个\(a_i\),我们构造一个数\(G_i\),使得\(G_i\)满足\(G_i \equiv a_i(mod\ b_i),G_i \equiv 0(mod\ b_j),(j\neq\ i)\)
想要满足后面那项比较容易,我们直接取\(\prod_{j\neq\ i}b_i\)即可。但是如何让他同时满足前一项?
这个也是比较简单的,对于\(G_i\),我们先求出其在\(mod\ b_i\)意义下的逆元\(inv_i\),那么\(G_i * inv_i * a_i\)即为我们要构造的答案。求逆元的时候用exgcd即可。
最后我们求出所有\(b_i\)的\(lcm\),因为两两互质其实就是\(\prod_{i=1}^nb_i\),把所有的\(G_i\)加在一起。对lcm取模即为答案。这个还是很好理解的。
有一道板子题[TJOI2009]猜数字代码实现如下。
#include<bits/stdc++.h>
#define rep(i,a,n) for(int i = a;i <= n;i++)
#define per(i,n,a) for(int i = n;i >= a;i--)
#define enter putchar('\n')
using namespace std;
typedef long long ll;
const int M = 500005;
const int INF = 1000000009;
const double eps = 1e-7;
ll read()
{
ll ans = 0,op = 1;char ch = getchar();
while(ch < '0' || ch > '9') {if(ch == '-') op = -1;ch = getchar();}
while(ch >= '0' && ch <= '9') ans = ans * 10 + ch - '0',ch = getchar();
return ans * op;
}
ll k,a[15],b[15],L = 1,ans;
ll exgcd(ll a,ll b,ll &x,ll &y)
{
if(!b){x = 1,y = 0;return a;}
ll d = exgcd(b,a%b,y,x);
y -= a / b * x;
return d;
}
ll mul(ll a,ll b,ll t)
{
ll cur = a * b - (ll)((long double)a / t * b + eps) * t;
return (cur % t + t) % t;
}
ll inc(ll a,ll b,ll t){return (a+b) % t;}
int main()
{
k = read();
rep(i,1,k) a[i] = read();
rep(i,1,k) b[i] = read(),L *= b[i];
rep(i,1,k)
{
ll x,y,cur = L / b[i];
exgcd(cur,b[i],x,y),x = (x % b[i] + b[i]) % b[i];
ans = inc(ans,mul(mul(cur,x,L),a[i],L),L);
}
printf("%lld\n",ans);
return 0;
}
4.exCRT
对于上面的问题我们解决的很顺利,因为我们保证了\(b_i\)是两两互质的,这也就使得所有构造出的\(G_i\)都是有逆元的。
但是如果\(b_i\)不是两两互质的,或者有时候有太多的同余方程,使得\(\prod{i=1}^nb_i\)根本无法计算,那这时候怎么解决问题呢?
于是就有了拓展中国剩余定理(exCRT).
首先,我们假设现在已经求出了前k-1个方程的一个解x,同时也知道\(M =LCM_{i=1}^{k-1} b_i\),那么前k-1个方程的通解就是\(x + tM\)
我们要求第k个方程的解,那么其实也就是求一个整数t,使得\(x + tM \equiv a_k(mod\ b_k)\)
这个式子是可以用\(exgcd\)求解的,如果无解那么就说明整个同余方程也是无解的。
否则的话前面k个方程的通解就是\(x + tM\),并且M要更新为\(LCM_{i=1}^k b_i\)
那么我们来看一道例题。[NOI2018]屠龙勇士
首先,每次攻击所用的剑是可以预先得知的,使用set处理一下即可。
那么我们的任务就变成了求解\(c_ix \equiv a_i (mod\ p_i)\)这样一个同余方程的最小正整数解。
这里有一个问题,就是只有当所有\(a_i < p_i\)的时候才成立,否则就会出现一个问题就是,你在攻击龙的时候可能没有把龙的血量打成负数,但是成为了\(p_i\)的倍数,这样在这个同余方程的计算中龙也是死了,但是你的解肯定是不对的。不过对于任意一个没有保证\(a_i < p_i\)的点,都有所有的\(p_i = 1\),那么答案很显然就是\(max_{i=1}^n\lceil\frac{a_i}{c_i}\rceil\),特判掉就好了。
回到这个问题上,但是这个同余方程和我们熟悉的形式不一样……我们也不能直接把\(c_i\)消掉 ,因为\(c_i\)不一定有\(mod\ p_i\)意义下的逆元。
那我们把这个方程组转化一下。
对于一个\(c_ix \equiv a_i (mod\ p_i)\),我们可以把它改写成为\(c_ix - p_iy = a_i\)的一个形式。这个式子是可以用exgcd求解的。我们只要求出这个式子的一组特解\(sx\),那么就有了\(x \equiv sx (mod\ gcd(c_i,p_i))\),于是同余方程就转化为我们熟悉的样子了。
之后只要用\(exCRT\)求解即可。注意本题需要使用龟速乘或者那个神奇的乘法,否则中间会爆longlong。
#include<bits/stdc++.h>
#define rep(i,a,n) for(int i = a;i <= n;i++)
#define per(i,n,a) for(int i = n;i >= a;i--)
#define enter putchar('\n')
using namespace std;
typedef long long ll;
const int M = 500005;
const int INF = 1000000009;
const double eps = 1e-7;
ll read()
{
ll ans = 0,op = 1;char ch = getchar();
while(ch < '0' || ch > '9') {if(ch == '-') op = -1;ch = getchar();}
while(ch >= '0' && ch <= '9') ans = ans * 10 + ch - '0',ch = getchar();
return ans * op;
}
ll T,n,m,a[M],p[M],c[M],g[M],s[M],bonus[M],att[M];
bool flag;
multiset<ll> q;
multiset<ll> :: iterator it;
ll inc(ll a,ll b,ll t){return (a + b) % t;}
ll mul(ll a,ll b,ll t)
{
ll tmp = (a * b - (ll)((long double)a / t * b + 1.0e-8) * t);
return (tmp % t + t) % t;
}
ll gcd(ll a,ll b){return b ? gcd(b,a%b) : a;}
ll exgcd(ll a,ll b,ll &x,ll &y)
{
if(!b){x = 1,y = 0;return a;}
ll d = exgcd(b,a%b,y,x);
y -= a / b * x;
return d;
}
void solve1()
{
ll ans = 0;
rep(i,1,n)
{
ll cur = (a[i] % c[i]) ? a[i] / c[i] + 1 : a[i] / c[i];
ans = max(ans,cur);
}
printf("%lld\n",ans);
}
void init()
{
n = read(),m = read(),q.clear(),flag = 0;
rep(i,1,n) a[i] = read();
rep(i,1,n) p[i] = read();
rep(i,1,n) bonus[i] = read();
rep(i,1,m) att[i] = read(),q.insert(att[i]);
rep(i,1,n)
{
ll k = *(q.begin());
if(a[i] < k) c[i] = k,q.erase(q.begin());
else it = q.upper_bound(a[i]),--it,c[i] = *(it),q.erase(it);
q.insert(bonus[i]);
}
rep(i,1,n) if(a[i] > p[i]) {solve1(),flag = 1;break;}
}
ll excrt()
{
ll N = p[1],ans = s[1],x,y;
rep(i,2,n)
{
ll a = N,b = p[i],c = (s[i] - ans % b + b) % b;
ll G = exgcd(a,b,x,y),b1 = b / G;
if(c % G) {return -1;}
x = mul(x,c/G,b1),ans += x * N,N *= b1;
ans = (ans % N + N) % N;
}
return (ans % N + N) % N;
}
int main()
{
T = read();
while(T--)
{
init();if(flag) continue;
rep(i,1,n)
{
ll x,y,G = gcd(c[i],p[i]);
if(a[i] % G) {flag = 1;break;}
p[i] /= G,exgcd(c[i] / G,p[i],x,y);
x = (x % p[i] + p[i]) % p[i],s[i] = mul(x,a[i] / G,p[i]);
}
if(flag){printf("-1\n");continue;}
ll ans = excrt();
printf("%lld\n",ans);
}
return 0;
}