世界总决赛选手带你玩转数论 3——同余方程原来如此简单

笔者

笔者曾获得 ICPC 2020 世界总决赛资格，ICPC 2020 亚洲区域总决赛第五名。

本次内容

本次主要针对一次同余方程和同余方程组展开，主要内容如下：

一次同余方程
一次同余方程组
同时补充讲解慢速乘

预告下一次，我们会针对二次同余和一些特殊形式的高次同余方程展开讲解。

一次同余方程

定理1

设

min Z^{ } ,a,bin Z,(a,m)=1

，则一次同余式

axequiv b(mod ; m)

恰有一个解，且其解为

xequiv ba^{varphi(m)-1}

。

定理2

设

m,nin Z^{ } ,a_1.a_2,cdots ,a_n,bin Z

，则一次同余式

a_1x_1 a_2x_2 cdots a_nx_nequiv b(mod ; m)

，有解的充分必要条件是

d|b

，其中

d=(a_1,a_2,cdots a_n,m)

，此时同余方程解的个数为

m^{n-1}d

。

一次同余方程组

一次同余方程组也被称为线性同余方程组。

由

个方程

xequiv x_i(mod ; m_i)

组成的方程组，求解

。

中国剩余定理

前提：

m_1,m_2,cdots ,m_n

两两互质。

令

M=prod_{i=1}^n m_i,M_i=frac{M}{m_i}

以及

M_i^{-1}

为

M_i ; mod ; m_i

的逆元，则有

xequiv sum_{i=1}^nx_iM_iM_i^{-1}(mod ; M)

。

证明

考虑

(x_iM_iM_i^{-1}); mod ; m_j

的结果：

当

i = j

时，

M_i^{-1}

为

M_i; mod ; m_i

的逆元，所以

(x_iM_iM_i^{-1})equiv x_i(mod ; m_i)

；当

ineq j

时，

M_i=frac{M}{m_i}

其中应该包含了

m_j

，则

(x_iM_iM_i^{-1})equiv 0(mod ; m_j)

。

所以

xequiv x_i (mod ; m_i)

显然成立。

一般情况的解法

所谓一般情况，其实就是考虑

m_1,m_2,cdots ,m_n

两两存在不互质的时候如何解。

考虑增量法来解线性同余方程组。即每次合并两个方程为一个方程，不断这样的往复操作，直到只剩下一个方程为止。

假设当前有两个方程

xequiv a(mod; b)\xequiv c(mod; d)

。

由第一个方程可以得到

x=bt_1 a

，其中

t_1

为任意整数，于是有

bt_1 aequiv c(mod; d)

，也就是

bt_1 a dt_2=c

，即

bt_1 dt_2=c-a

。

这个方程中有两个未知数

t_1,t_2

，我们可以考虑通过拓展 GCD 来求得一个任意解。

于是可以得到

xequiv bt_1 a(mod ; [b,d])

。

例题

https://codeforces.com/problemset/problem/338/D

Codeforces338D GCD Table

题目大意

给出一张

ntimes m

的表格，表格的第

行，第

列为

gcd(i,j)

，现在给出一个长度为

的数列，询问是否是矩阵中某一行连续的一部分。

分析

因为对于所有的

a_i

都有

gcd(x,y i-1)=a_i

，也就是

a_i|x

，那么

至少是

lcm(a_1,a_2,cdots ,a_k)

。

又因为

a_i|(y i-1)

，所以有

yequiv 1-i(mod ; a_i)

，于是我们可以通过解

个线性同余方程组求得一个最小的正整数

。

显然此时我们求得的一个可能解

x=x_0=lcm(a_1,a_2,cdots ,a_k),y=y_0

是最小的。

对于任意一个数

a_i=gcd(x,y)

，如果我们增大

或者

，显然必须满足

x=xcdot t_1

或者

y=y x_0t_2

，而显然，这样的变化只会让

gcd(x,y)

不变或者变大，于是我们知道，如果

x=x_0=lcm(a_1,a_2,cdots ,a_k),y=y_0

是不合法的，显然通过变化得到的

x,y

也不会合法。

上面推导得到

x,y

的解显然只满足必要性，所以我们必须要对结果进行验证。

在合并线性方程的过程中，需要用到慢速乘。

【补充】关于慢速乘

如果

x,y,ple 10^{18}

，显然

x,y,p

均可以用 long long 存储，但是如果要计算

xymod p

的结果，因为中间过程会超过 long long，而不能直接计算。

在这种时候我们可以采用快速幂类似的方法得到慢速乘，即把乘法拆成加法来做，原理类似。

代码语言：javascript复制

LL mul(LL a,LL b,LL m){
    LL ret=0;
    while (b){
        if (b&1){
            ret =a;
            if (ret>=m) ret-=m;
        }
        a =a;
        if (a>=m) a-=m;
        b>>=1;
    }
    return ret;
}

还有一种更加神奇的方法是用 long double，同时通过分类讨论避免可能存在的精度问题：

代码语言：javascript复制

inline LL mul(LL x,LL y,LL p)
{
    x%=p;y%=p;
    if (p<=1000000000) return x*y%p;
    if (p<=1000000000000LL) return (((x*(y>>20)%p)<<20) x*(y&((1<<20)-1)))%p;
    LL d=(LL)floor(x*(long double)y/p 0.5);
    LL res=(x*y-d*p)%p;if (res<0) res =p;
    return res;
}

下面是例题的代码：

代码语言：javascript复制

#include<bits/stdc  .h>
#define LL long long
using namespace std;

inline LL mul(LL x,LL y,LL p)
{
    x%=p;y%=p;
    if (p<=1000000000) return x*y%p;
    if (p<=1000000000000LL) return (((x*(y>>20)%p)<<20) x*(y&((1<<20)-1)))%p;
    LL d=(LL)floor(x*(long double)y/p 0.5);
    LL res=(x*y-d*p)%p;if (res<0) res =p;
    return res;
}

LL ex_gcd(LL a, LL b, LL &x, LL &y) {//ax by=gcd(a,b)
    if (b == 0) { x = 1; y = 0; return a; }
    LL ret = ex_gcd(b, a % b, y, x);
    y -= a / b * x;
    return ret;
}

LL GCD(LL x,LL y)
{
    LL r=x%y;
    while(r) x=y,y=r,r=x%y;
    return y;
}

LL n,m,a[10010];
int k;

bool check(LL x,LL y)
{
    for (int i=1;i<=k;i  )
    {
        LL z=GCD(x,y);
        //cout<<x<<" "<<y<<" "<<a[i]<<endl;
        if (z!=a[i]) return false;
        y  ;
    }
    return true;
}

int main()
{
    scanf("%lld%lld%d",&n,&m,&k);
    scanf("%lld",&a[1]);
    LL x=a[1],y,l=a[1];
    for (int i=2;i<=k;i  )
    {
        scanf("%lld",&a[i]);
        x=GCD(l,a[i]);
        if (l/x>n/a[i]) {puts("NO");return 0;}
        l=l/x*a[i];
    }
    LL A=0,B=a[1];
    for (int i=2;i<=k;i  )
    {
        LL z=GCD(B,a[i]),C=1-i,D=a[i];
        //cout<<A<<" "<<B<<" "<<C<<" "<<D<<endl;
        if ((C-A)%z!=0) {puts("NO");return 0;}
        ex_gcd(B,D,x,y);
        x*=(C-A)/z;y=B;
        B=B/z*D;A=(mul(y,x,B) A)%B;
    }
    y=(A%B B)%B;
    if (y==0) y =B;
    if (y k-1>m) {puts("NO");return 0;}
    if (check(l,y)) puts("YES");
    else puts("NO");
    return 0;
}

海外加速

0 人点赞