从exgcd到exCRT

Posted 2021-02-09 captain1

从最基础的开始。
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 - lfloorfrac{a}{b} floor * b)
(ax_1 + by_1 = bx_2 + ay_2 - blfloorfrac{a}{b} floor y_2)
合并同类项，得到(ax_1 + by_1 = b(x_2 - lfloorfrac{a}{b} floor y_2) + ay_2)
所以得到(x_1 = y_2,y_1 = x_2 - lfloorfrac{a}{b} floor 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 eq i))
想要满足后面那项比较容易，我们直接取(prod_{j eq 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(‘
‘)

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
",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}^nlceilfrac{a_i}{c_i} ceil),特判掉就好了。
回到这个问题上，但是这个同余方程和我们熟悉的形式不一样……我们也不能直接把(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(‘
‘)

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
",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
");continue;}
      ll ans = excrt();
      printf("%lld
",ans);
   }
   return 0;
}