同余方程与离散对数

同余方程：

一、解线性同余方程组的方法：

1、方程组只有两个的情况：

$$\begin{gathered} x\equiv r_0(mod{m}_0) \\ x\equiv r_1(mod{m}_1) \end{gathered}$$

​ 这两个方程等价于：

$$\begin{gathered} x+k_0\times m_0=r_0 \\ x+k_1\times m_1=r_1 \end{gathered}$$
​ 若方程组有解，则不定方程 $k_0\times m_0-k_1\times m_1=r_0-r_1$ 有解（这个方程有无解可以通过裴蜀定理来判断），且将该不定方程的任意一组整数解带入到原方程组中都将使原先的两个方程同时成立。

​ 我们假设 $k_{00}$,$k_{10}$ 是方程 $k_0\times m_0-k_1\times m_1=r_0-r_1$ 的一组特解，则方程的通解为：

$$\begin{gathered} k_0=k_{00}+t\times \frac{m_1}{gcd(m_0,m_1)} \\ {k}_1=k_{10}+t\times \frac{m_0}{gcd(m_0,m_1)} \end{gathered}$$
（两式中t取值相同且为整数）

​ 将 $k_0=k_{00}+t\times \frac{m_1}{gcd(m_0,m_1)}$ 带入到式子 $x+k_0\times m_0=r_0$ 中,得:
$$x+k_{00}\times m_0+t\times \frac{m_0\times m_1}{gcd(m_0,m_1)}=r_0$$
​ 也就是
$$x+t\times [m_0,m_1]=r_0-k_{00}\times m_0$$
​ 所以
$$x\equiv (r_0-k_{00}\times m_0)(mod[m_0,m_1])$$
​ 所以，原方程组的解与方程 $x\equiv (r_0-k_{00}\times m_0)(mod[m_0,m_1])$ 相同，即我们可以将原来的两个方程合并为一个方程，这样，我们就只要求解一个同余方程就可以了，用拓展欧几里得算法可以快速解决。

2、方程组含有多个方程的情况：

不断将两个方程合并为一个方程即可：

code:

#include<iostream>
#include<cmath>
#include<vector>
#define int long long
#define pii pair<int,int> 
using namespace std;

int exgcd(int a,int b,int &x,int &y)
{
    if(b==0)
    {
        x=1,y=0;
        return a;
    }
    int d = exgcd(b,a%b,y,x);
    y -= x * (a/b);
    return d;
}

int x,y;

pii solve(vector<pii> p)
{
    int r0=0,m0=1;
    int r1,m1;
    int m;
    int d,t;
    for(int i=0;i<p.size();i++)
    {
        r1 = p[i].first,m1=p[i].second;
        m = r0-r1;
        d=exgcd(m0,m1,x,y);
        if(m%d!=0)return {0,0};
        x*=(m/d);
        x = (x%(m1/d)+m1/d)%(m1/d);
        r0 = r0 - x*m0;
        m0 = (m1/d*m0);
        r0 = (r0%m0+m0)%m0;
       // cout<<r0<<" "<<a0<<endl;
    }
    return {r0,m0};
}

vector<pii> p;

signed main()
{
    int k;
    cin>>k;
    for(int i=1;i<=k;i++)
    {
        int l,r;
        cin>>l>>r;
        p.push_back({r,l});
    }
    pii ans = solve(p);
    if(ans.second ==0)cout<<-1;
    else cout<<ans.first;
}

3、$x$ 的系数不为1的情况 $a_i\times x\equiv r_i(mod{m}_i)$

方法一：

​ 我们考虑将每一个方程 $a_i\times x\equiv r_i(mod{m}_i)$ 转换成和它等价的 $x\equiv r(modm)$，转换之后，就转换成了系数为1的情况，再套用上面的代码就可以了。

那么，如何转换呢？

我们原来的方程就等价于 ：
$$a_i\times x+m_i\times y=r_i$$
如果 $gcd(a_i,m_i)|r_i$，则该方程有解，那么此时我们我们可以给方程两边同时除以 $gcd(a_i,m_i)$，我们记 $d=gcd(a_i,m_i)$则 :
$$\frac{a_i}{d}\times x+\frac{m_i}{d}\times y=\frac{r_i}{d}$$
这个式子等价于：
$$\frac{a_i}{d}\times x\equiv \frac{r_i}{d}(mod\frac{m_i}{d})$$
此时
$$gcd(\frac{a_i}{d},\frac{m_i}{d})=1$$
所以我们可以通过逆元将其再次转换成下面这个式子：
$$x\equiv \frac{r_i}{d}\times (\frac{a_i}{d}{)}^{-1}(mod\frac{m_i}{d})$$
code 1:(蒟蒻的码）

#include<bits/stdc++.h>
#define int long long
#define pii pair<int,int>

using namespace std;

int exgcd(int a,int b,int &x,int &y)
{
    if(b==0)
    {
        x=1,y=0;
        return a;
    }
    int d = exgcd(b,a%b,y,x);
    y -= x * (a/b);
    return d;
}
pii get_pii(vector<pii> p)
{
    int x,y;
    int r0=0,m0=1;
    int r1,m1;
    int m;
    int d,t;
    for(int i=0;i<p.size();i++)
    {
        r1 = p[i].first,m1=p[i].second;
        m = r0-r1;
        d=exgcd(m0,m1,x,y);
        if(m%d!=0)return {-1,0};//pii.second=0为无解的标致。
        x*=(m/d);
        x = (x%(m1/d)+m1/d)%(m1/d);
        r0 = r0 - x*m0;
        m0 = (m1/d*m0);
        r0 = (r0%m0+m0)%m0;
        // cout<<r0<<" "<<a0<<endl;
    }
    return {r0,m0};
}

pii solve(vector<int> a,vector<int> r,vector<int> m)
{
    pii res;
    int x,y;
    vector<pii> p;
    for(int i=0;i<a.size();i++)
    {
        int d = __gcd(a[i],m[i]);
        if(r[i]%d)return {-1,0};//pii.second = 0 为无解的标志.
        a[i]/=d,r[i]/=d,m[i]/=d;
        exgcd(a[i],m[i],x,y);
        x = (x % m[i] + m[i]) % m[i];
        r[i] = r[i] * x ;
        p.push_back({r[i],m[i]});
    }
    return get_pii(p);
}

signed main()
{
    vector<int> a,r,m;
    int n;
    cin>>n;
    for(int i=1;i<=n;i++)
    {
        int sr,sm;
        cin>>sm>>sr;
        a.push_back(1);
        r.push_back(sr);
        m.push_back(sm);
    }
    cout<<solve(a,r,m).first;
    return 0;
}

方法二：

上面我们提出了一种将两个同余方程合并为一个的方法，现在，我们介绍一种更为通用的方法：

对于任意一个线性同余方程 $a_{i}x\equiv r_i(mod{m}_i)$，若有解，则解一定可以表示为 $x\equiv r(modm)$，这样求解线性同余方程组的问题就转换成了求解 $r$ ，$m$ 的问题。

假设我们现在已经有 $x\equiv r(modm)$,再加入一个方程 $a_1\times x\equiv r_1(mod{m}_1)$，那么如何进行合并呢？

首先：
$$x\equiv r(modm)等价于x=r+t\times m$$
将其代入到 $a_1\times x\equiv r_1(mod{m}_1)$ 中去,得：
$$a_1\times m\times t\equiv r_1-r\times a_1(mod{m}_1)$$
我们令
$$A=a_1\times m,R=r_1-r\times a_1$$
则式子 $a_1\times m\times t\equiv r_1-r\times a_1(mod{m}_1)$就等价于：
$$A\times t=R+k\times m_1(其中k为任意整数)$$
令
$$gcd(A,m_1)=d$$
则让两边同时除以 $d$，得：
$$\frac{A}{d}\times t=\frac{R}{d}+k\times \frac{m_1}{d}$$
这个式子就等价于 ：
$$\frac{A}{d}\times t\equiv \frac{R}{d}(mod\frac{m_1}{d})$$
由于 $gcd(\frac{A}{d},\frac{m_1}{d})$ 必然等于 $1$，所以 $\frac{A}{d}$ 模 $\frac{m_1}{d}$ 的逆元存在，所以再对式子变形得：
$$t\equiv \frac{R}{d}\times (\frac{A}{d}{)}^{-1}(mod\frac{m_1}{d})$$
也就是说
$$t=\frac{R}{d}\times (\frac{A}{d}{)}^{-1}+k\times \frac{m_1}{d}（其中k为整数），$$
代入到 $x$ 的表达式中得：
$$x=r+m\times \frac{R}{d}\times (\frac{A}{d}{)}^{-1}+k\times (m\times \frac{m_1}{d})$$
等价于：
$$x\equiv r+m\times \frac{R}{d}\times (\frac{A}{d}{)}^{-1}(modm\times \frac{m_1}{d})$$
这样，我们就完成了合并，对于方程组，我们可以一直合并下去，所有方程合并完之后就得到了方程组的解。

code2:(大佬的码，来自挑战程序设计竞赛)

#include<bits/stdc++.h>
#define int long long
#define pii pair<int,int>
using namespace std;

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-=x*(a/b);
    return d;
}

int get_inv_mod(int a,int m)//保证a,m互质
{
    int x,y;
    exgcd(a,m,x,y);
    return (x%m+m)%m;
}

pii get_ans(vector<int> &a,vector<int> &r,vector<int> &m)
{
    int ans_r=0,ans_m=1;
    for(int i=0;i<a.size();i++)
    {
        int A=a[i]*ans_m%m[i],R=((r[i]-a[i]*ans_r)%m[i]+m[i])%m[i];
        int d= __gcd(m[i],A);
        if(R%d)return {-1,0};
        A/=d,R/=d,m[i]/=d;
        int t = R * get_inv_mod(A,m[i])%m[i];
        ans_r += ans_m * t ;
        ans_m *= m[i];
        ans_r %= ans_m;
    }
    return {ans_r%ans_m,ans_m};
}

signed main()
{
    vector<int> a,r,m;
    int n;
    cin>>n;
    for(int i=1;i<=n;i++)
    {
        int sm,sr;
        cin>>sm>>sr;
        a.push_back(1),r.push_back(sr),m.push_back(sm);
    }
    cout<<get_ans(a,r,m).first;
    return 0;
}

二、离散对数：解 $a^t\equiv b(modp)$ 类型的同余方程（ $t$ 为未知量）：

1、$gcd(a,p)=1$ 的情况（ $bsgs$ 算法）：

• 首先，很显然，如果有解，根据费马小定理和欧拉定理我们可以知道解一定在区间 $\left[0,\phi \left(p\right)-1\right]$ 上。

• 那么，我们考虑将一个更大的区间分块。

• 比如我们将区间 $\left[1,p\right]$ 分成 $k=\lfloor \sqrt{p}+1\rfloor$ 块，每块的长度为 $k$ ,那么实际上我们是将区间 $\left[1,k\times k\right]$ 分成了k块。

• 注意，区间 $\left[1,k\times k\right]$ 包含了区间 $\left[0,\phi \left(p\right)-1\right]$ ,所以我们的解也一定会落在我们分块的区间上(对于 $0$ 我们进行特判）。

• 我们把第 $x$ 块的取值范围表示出来就是 $t=k\times x-y,y\in \left[0,k-1\right]$ 。

• 如果在该分块上有解，也就是说存在 $y\in \left[0,k-1\right]$ 使得 $a^{k\times x-y}\equiv b(modp)$。

• 对上面这个式子变形得： $a^{k\times x}\equiv b\times a^y(modp)$，再次变形得 $(a^k{)}^x\equiv b\times a^y(modp)$。

• 对于每个 $x $，上面这个式子中得 $y$ 都能取到$[0,k-1]$,所以我们可以将 $b\times a^y\mathrm{mod}p$ 的结果预处理出来，并将结果与对应的 $y$ 存起来，由于我们一般情况下只要最小整数解，即使得 $t=k\times x-y$ 最小，所以我们对于同一个结果，只要存储它对应的最大的 $y$ 就可以。

• 预处理之后，我们只需要从小到大再遍历一下 $x$ 找出第一个存在解的 $x$ 就可以了，如果没有找到，则无解。

  code :

  //将一个大于p的区间（1~p)分成k块，这里我们设置k = sqrt(p) +1 块，每块长度为 k ,
  // 那么这k块就包含了 1~k*k 这个区间的所有数，
  // 第 x 块所包含的范围是 k*x  + y ,其中 x 属于 [1,k] , y 属于 [0,k-1];
  // 那么 a^t = b (mod p ) 就等价于 a^(kx-y) = b (mod p);
  // 就会等价于 a^(kx) = b * a ^ (y) (mod p);
  // 这样的话，我们可以用一个哈希表将 a^y %p 的结果与 y 对应起来，如果出现冲突，只存最大的y ,因为我们要使得 kx - y 最小。
  // 然后再从小到大枚举x,如果存在一个y,使得a^kx = b * a ^ y (mod p) ,则输出 k*x - y 即可。
  #include<algorithm>
  #include<unordered_map>
  #include<cmath>
  #include<iostream>
  #define int long long
  using namespace std;
  const int inf = -1e9;
  int bsgs(int a, int b, int p)
  {
      if (1 % p == b % p) return 0;
      int k = sqrt(p) + 1;
      unordered_map<int, int> hash;
      for (int i = 0, j = b % p; i < k; i ++ )
      {
          hash[j] = i;
          j = j * a % p;
      }
      int ak = 1;
      for (int i = 0; i < k; i ++ ) ak =ak * a % p;
  
      for (int i = 1, j = ak; i <= k; i ++ )
      {
          if (hash.count(j)) return i * k - hash[j];
          j = j * ak % p;
      }
      return inf;
  }
  
  signed main() {
      int a, b, p;
      while (cin >> a >> p >> b, a || b || p)
      {
          int res = bsgs(a, b, p);
          if (res < 0)cout << "No Solution\n";
          else cout << res << "\n";
      }
  }
  

2、$gcd(a,p)\ne 1$的情况（ $exbsgs$ 算法）：

• 我们考虑将其转换成 $gcd(a,p)=1$ 的情况。

• 原式等价于 $a^t+k\times p=b$，设 $d=gcd(a,p)$ 则当 $d\mid b$ 时，原方程可能有解，否则一定无解。

• 当 $d\mid b$ 时，我们对式子两边同时除以 $d$ 得：$\frac{a}{d}\times a^{t-1}+k\times \frac{p}{d}=\frac{b}{d}$，即 $a^{t-1}\equiv \frac{b}{d}\times (\frac{a}{d}{)}^{-1}(mod\frac{p}{d})$。

• 更新 $t=t-1$，$b=\frac{b}{d}\times (\frac{a}{d}{)}^{-1}$ 和$p=\frac{p}{d}$，一直这样递归操作下去，直到可以判定为无解或者 $gcd(a,p)=1$ 时停止，若无解，则返回无解的标志，若 $gcd(a,p)=1$ 则返回 $bsgs(a,b,p)$ 的值。

  code:

  #include<algorithm>
  #include<unordered_map>
  #include<cmath>
  #include<iostream>
  #define int long long
  using namespace std;
  const int inf = -1e9;
  
  int exgcd(int a,int b,int &x,int &y)
  {
      if(b==0)
      {
          x=1,y=0;
          return a;
      }
      int d = exgcd(b,a%b,y,x);
      y -= a/b * x;
      return d;
  }
  
  int bsgs(int a, int b, int p)
  {
      if (1 % p == b % p) return 0;
      int k = (int)sqrt(p) + 1;
      unordered_map<int, int> hash;
      for (int i = 0, j = b % p; i < k; i ++ )
      {
          hash[j] = i;
          j = j * a % p;
      }
      int ak = 1;
      for (int i = 0; i < k; i ++ ) ak =ak * a % p;
  
      for (int i = 1, j = ak; i <= k; i ++ )
      {
          if (hash.count(j)) return i * k - hash[j];
          j = j * ak % p;
      }
      return inf;
  }
  
  int exbsgs(int a,int b,int p)
  {
      b = (b % p + p) % p;
      if (1 % p == b % p)return 0;
      int d = __gcd(a, p);
      if (d > 1)
      {
          if (b % d)return inf;
          b /= d, p /= d;
          int x, y;
          int t = exgcd(a / d, p, x, y);
          x = (x % (p / t) + p / t) % (p / t);
          return exbsgs(a, b * x, p) + 1;
      }
      else return bsgs(a, b, p);
  }
  
  signed main() {
      int a, b, p;
      while (cin >> a >> p >> b, a || b || p)
      {
          int res = exbsgs(a, b, p);
          if (res < 0)cout << "No Solution\n";
          else cout << res << "\n";
      }
  }