同余方程总结

zach0914 · 2021-01-31 19:51:25 · 个人记录

~有关逆元的一些常用结论：我的博客~

一、同余方程

ax \equiv c(mod\ b)

大致把它转化成：

ax+by=c

进而->

ax+by=d\ (d=(a,\ b))

跑扩展欧几里得算法。

重点在于后续处理：

首先，如果c不是d的倍数，那么一定要判断一下是否无解；
其次，对于每一个x来讲，扩欧跑出来的是d下的值，所以转化为原来的式子必须要乘上c/d才行；
对于x的通解留意：x+k*(b/d)所以，对于x取模的时候，一般我们习惯上将mod=|b/d|进行正数取模；

例题：同余方程

很板子吧：：

而且已经保证有解了。

例题：青蛙的约会

这题也是版子题，相比于上一道题，这道题需要判断是否无解，取模取正。

例题：扩展中国剩余定理

核心做法：类似于中国剩余定理，可以考虑构造。

不妨设处理到了第k个数，之前的答案为：x+k*lcm(lcm=[a_1,a_2,…,a_{k-1}])，那么对于当前的答案来讲我们只需要满足的k使得当前第k个方程成立。

用扩展欧几里得算法（线性同余方程的一般解法）。

要注意一些要点：

合并解的时候，新的答案（不妨设为x）应对刚刚计算的lcm取模，而x乘的是k-1轮的lcm。而在求解线性同余方程的时候，应对另一个数取模。

跑n次exgcd就可以解决出最终的答案。

不过，题目明显提示你可能会溢出的情况。

为什么会出现溢出的情况，我乘完一个数就立即取模呐！！

事实上，如果两个1e18范围内的数取模，在中间的乘法会溢出来。

所以，我们得使用“龟速乘”解决此类问题（方法同快速幂）。

const int N = 1e5 + 5;

typedef long long LL;

LL mul(LL x, LL y, LL mod)
{
    LL res = 0;
    x %= mod, y %= mod;
    bool mark = false;
    if(y < 0) mark = true, y = -y;
    while(y)
    {
        if(y & 1) res = (res + x) % mod;
        y >>= 1;
        x = (x + x) % mod;
    }
    if(!mark) return res;
    return -res;
}
LL gcd(LL x, LL y)
{
    if(!y) return x;
    return gcd(y, x % y);
}
LL lcm(LL x, LL y)
{
    return x / gcd(x, y) * y;
}

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;
}
int n;
LL a[N], b[N];
int main()
{
    scanf("%d", &n);
    FOR(i, 1, n) 
    {
        scanf("%lld %lld", &a[i], &b[i]);
        b[i] %= a[i];
    }
    LL ans = b[1], c, d, x, y, L = a[1], mod;
    FOR(i, 2, n)
    {
        c = b[i] - ans;
        d = exgcd(L, a[i], x, y);
        mod = abs(a[i] / d);
        c = c / d;
        x = mul(x, c, mod);
        x = (x % mod + mod) % mod;
        LL tmp = L;
        L = lcm(L, a[i]);
        ans = (ans + mul(x, tmp, L)) % L;
    }
    printf("%lld\n", ans);
    return 0;
}

例题：[NOI2002] 荒岛野人

这道题看起来好像有一个时间复杂度很险的做法（居然是正解）：

枚举m，对于每一个m判断是否满足要求，而n很小，枚举所有的约束条件即可。

……
typedef long long LL;

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

int n, C[N], P[N], L[N];

int main()
{
    scanf("%d", &n);
    int lim = 0;
    FOR(i, 1, n) 
    {
        scanf("%d %d %d", &C[i], &P[i], &L[i]);
        lim = max(lim, C[i]);
    }

    int a, c, d;//a = pi-pj, c = cj-ci, b = m, ax + by = c 对于x <= min(Li,Lj)是否有解, 有解就挂掉  
    FOR(b, lim, 1e6)
    {
        bool valid = true;
        LL x, y, mod;
        FOR(i, 1, n)
        {
            FOR(j, i + 1, n)
            {
                a = P[i] - P[j], c = C[j] - C[i];
                d = exgcd(a, b, x, y);
                if(c % d) continue;
                x = 1ll * x * (c / d), mod = abs(b / d);
                x = (x % mod + mod) % mod;
                if(x <= min(L[i], L[j])) 
                {
                    valid = false;
                    break;
                }   
            }
            if(!valid) break;
        }
        if(valid) 
        {
            printf("%d\n", b);
            break;
        }
    }
    return 0;
}

二、中国剩余定理（孙子定理CRT）

答案为：x=\sum_{i=1}^n{a_iM_it_i}，其中M_i为M=\prod_{i=1}^nm_i，M_i=M/m_{i}，M_it_i\equiv 1(mod\ m_i)

证明挺显然的。看一看注意事项：

中国剩余定理使用必要条件为：m_i两两互素；
最后的x应对M取模。
注意long\ long，以及求逆元时注意简化范围（mod\ m_i）

例题：[TJOI2009]猜数字

几乎就是crt了。只不过使用“龟速乘”。

const int N = 15;

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

int n, a[N], b[N];

int mul(int x, int y, int mod)
{
    int res = 0;
    while(y)
    {
        if(y & 1) res = (res + x) % mod;
        y >>= 1;
        x = (x + x) % mod;
    }
    return res;
}

main()
{
    scanf("%lld", &n);
    int m = 1;
    FOR(i, 1, n) scanf("%lld", &a[i]);
    FOR(i, 1, n) scanf("%lld", &b[i]);
    FOR(i, 1, n) m = m * b[i];
    FOR(i, 1, n) a[i] = (a[i] % b[i] + b[i]) % b[i];

    int ans = 0;
    FOR(i, 1, n)
    {
        int x, y; 
        exgcd(m / b[i], b[i], x, y);
        x = (x % b[i] + b[i]) % b[i];
        ans = (ans + mul(mul(m / b[i], a[i], m), x, m)) % m;
    }
    printf("%lld\n", ans);
    return 0;
 }

中国剩余定理能解决的问题不仅仅局限于上述。具体来说，我们可以解决一些诸如以下的问题：

x\equiv a(mod\ p)

对于这个等式，如果p不是质数，且a非常的大（允许动态求解）那么我们可以考虑从将模数p分解质因数：

p=p_1^{\alpha_1}*p_2^{\alpha_2}*……*p_k^{\alpha_k}

然后对于每一个p_i^{\alpha_i}当作一个模数，对a取模，构成k个方程。

考虑到模数两两互素，我们可以使用CRT求解。即x=\sum {a_i*M_i*t_i}+k*m，这个x就是a关于p的同余。

证明：若x\equiv a(mod\ p_i)，那么x\equiv a(mod\ \prod{p_i})

反证法。

考虑x=\sum {a_i*M_i*t_i}+k*m，假设存在x\not\equiv a(mod\ \prod{p_i})，那么x-a\not\equiv 0(mod\ \prod{p_i})

由于满足方程组，所以有：

x-a\equiv 0(mod\ p_i)

故p_i|x-a

由于p_i具有任意性，所以\prod{p_i}|x-a，矛盾。

故结论必定成立。

证毕。

该证明就是上述成立的基础。

例题：[SDOI2010]古代猪文

这道题首先将底数拎出来，然后就转化为对指数取模的问题。

又由欧拉定理可得：

a^{b}\equiv a^{b\ mod\ \ phi(n)}(mod\ n)|(a,n)=1

套用刚才的结论：

求： g^{\sum_{k | n}C_n^k}\ mod\ p

<=>$ 求:$ \sum_{k | n}C_n^k \mod \phi(p)

相当于求解:

x \equiv C_n^k (mod\ 2)

x \equiv C_n^k (mod\ 3)

x \equiv C_n^k (mod\ 4679)

x \equiv C_n^k (mod\ 35617)

考虑到模数m_i互素, 所以我们只需要求解CRT即可。

枚举n的因子，按照上式跑CRT；

至于出来的组合数，我们需要用到Lucas定理：

当p为质数时，有C_n^k\equiv C_{n/p}^{k/p}*C_{n/mod\ p}^{k/mod\ p}(mod\ p)

特别地，当g为模数的因数时，需要格外小心注意此前情况。

……
int m[5] = {0, 2, 3, 4679, 35617}, t[5] = {0, 1, 1, 1353, 31254};

int n, g, jc[5][N] = {}, jc_inv[5][N] = {};
int power(int x, int y, int k)
{
    int res = 1 % k;
    while(y)
    {
        if(y & 1) res = 1ll * res * x % k;
        y >>= 1;
        x = 1ll * x * x % k;
    }   
    return res;
}
int CRT(int *a)
{
    int x = 0;
    FOR(i, 1, 4)
    {
        x += 1ll * M / m[i] * a[i] % M * t[i] % M;
        x %= M;
    }
    return x;
}
int C(int x, int y, int l)
{
    if(y < x) return 0;
    return 1ll * jc[l][y] * jc_inv[l][x] % m[l] * jc_inv[l][y - x] % m[l];
}
int Lucas(int x, int y, int l)
{
    if(y < x) return 0; 
    if(!x || x == y) return 1;
    return 1ll * Lucas(x / m[l], y / m[l], l) * C(x % m[l], y % m[l], l) % m[l];
}
int solve(int p, int q)
{
    static int a[5];
    CLR(a, 0);
    FOR(i, 1, 4)
    {
        a[i] = Lucas(p, q, i);
    }
    return CRT(a);
}
void init()
{
    FOR(i, 1, 4)
    {
        jc[i][0] = 1;
        FOR(j, 1, m[i]) jc[i][j] = 1ll * jc[i][j - 1] * j % m[i];   
        jc_inv[i][0] = 1, jc_inv[i][m[i]] = 0;
        jc_inv[i][m[i] - 1] = power(jc[i][m[i] - 1], m[i] - 2, m[i]);
        ROF(j, m[i] - 2, 1) jc_inv[i][j] = 1ll * jc_inv[i][j + 1] * (j + 1) % m[i];
    }
    return;
}

int calc()
{
    int res = 0;
    for(Re int i = 1; 1ll * i * i <= n; ++ i)
    {
        if(n % i == 0)
        {
            res = (res + solve(i, n)) % M;
            if(i * i < n) res = (res + solve(n / i, n)) % M;
        }
    }
    return res;
}

int main()
{
    init();
    scanf("%d %d", &n, &g);
    if(g == mod) 
    {
        putchar('0');
        return 0; 
    } 
    g %= mod;
    printf("%d\n", power(g, calc(), mod) % mod);
    return 0;
}

三、高次同余方程

长这模样：

a^x\equiv b(mod\ p)

思路类似于中途相遇法（meet in the middle）

令t=\sqrt p，然后x=kt-r，接着移项，变成(a^t)^k\equiv b\times a^r，然后右式扔进hash(map)表里，枚举k，找是否有相等的值。

这就是有名的Baby Step Giant Step。

须注意的要点：

使用条件为：(a,p)=1；
判断一下a^t\equiv 0\ (mod\ p)的情况（即若b=0返回1否则返回-1）。

例题：[SDOI2011]计算器

这道题类型1时利用欧拉定理去做，询问类型2是用线性同余方程求解，类型三就是Baby Step Giant Step算法。

…
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 power(LL x, LL y, LL mod)
{
    LL res = 1 % mod;
    while(y)
    {
        if(y & 1) res = 1ll * res * x % mod; 
        y >>= 1;
        x = 1ll * x * x % mod;
    }
    return res;
}

LL Baby_Step_Giant_Step(LL a, LL b, LL mod)
{
    map <LL, LL> hash; hash.clear();
    b %= mod;
    LL t = (int)sqrt(mod) + 1, tmp = 1;
    for(int i = 0; i < t; ++ i)
    {
        LL val = 1ll * b * tmp % mod;
        hash[val] = i;
        tmp = 1ll * tmp * a % mod;
    }
    a = tmp % mod;
    if(!a) return b == 0 ? 1 : -1;
    tmp = 1;
    for(int i = 0; i <= t; ++ i)
    {
        LL val = tmp, res;
        int j = hash.find(val) == hash.end() ? -1 : hash[val];
        res = 1ll * i * t - j;
        if(j >= 0 && res >= 0) return res;
        tmp = 1ll * tmp * a % mod;
    }
    return -1;
}

int main()
{
    int T, k;
    LL a, b, c;
    scanf("%d %d", &T, &k);
    while(T --)
    {
        scanf("%lld %lld %lld", &a, &c, &b);
        if(k == 1)
        {
            a %= b, c %= b - 1;
            printf("%lld\n", power(a, c, b));
        }
        else if(k == 2)
        {
            LL x, y, d = exgcd(a, b, x, y);
            if(c % d) puts("Orz, I cannot find x!");
            else
            {
                x *= c / d;
                LL t = abs(b / d);
                x = (x % t + t) % t;
                printf("%lld\n", x);
            }
        }
        else
        {
            LL ans = Baby_Step_Giant_Step(a, c, b);
            if(ans == -1) puts("Orz, I cannot find x!");
            else printf("%lld\n", ans);
        }
    }
    return 0;
}

例题：随机数生成器

这道题看到是数列递推我就想到了”不动点“法求解一阶线性递推数列的公式。具体来说，”不动点“法求解的核心是待定系数法，也就是将左右两式配凑成一个等比数列的形式：

x_{i+1}+\frac b{a-1}≡a×(x_i+\frac b{a-1})(mod\ p)

换元，将这个式子推成高阶线性递推的式子，然后就不用说了吧！

交上去：Wrong Answer!

这道题阴险的地方在于：a可能为0或1，所以，对于这两种情况要特别判断。

另外，如果x_1=t则直接输出1，不用计算。

……
typedef long long LL;
int power(int x, int y, int mod)
{
    int res = 1 % mod;
    if((x %= mod) == 0) return y == 0 ? 1 : 0;
    while(y)
    {
        if(y & 1) res = 1ll * res * x % mod;
        y >>= 1;
        x = 1ll * x * x % mod;
    }
    return res; 
}
LL Baby_Step_Giant_Step(LL a, LL b, LL mod)
{
    map<LL, LL> hash; hash.clear();
    LL t = sqrt(mod) + 1, tmp = 1, val;
    for(LL i = 0; i < t; ++ i)  
    {
        val = 1ll * b * tmp % mod;
        hash[val] = i;
        tmp = 1ll * tmp * a % mod;
    }
    a = tmp, tmp = 1;
    if(!a) return b == 0 ? 1 : -1;
    for(LL i = 0; i <= t; ++ i)
    {
        int j = hash.find(tmp) == hash.end() ? -1 : hash[tmp];
        if(j >= 0 && 1ll * i * t - j >= 0) return 1ll * i * t - j + 1;
        tmp = 1ll * tmp * a % mod;
    }
    return -1;
}
LL p, a, b, x_1, t;
int main()
{
    int T;
    scanf("%d", &T);
    while(T --)
    {
        scanf("%lld %lld %lld %lld %lld", &p, &a, &b, &x_1, &t);
        if(x_1 == t) puts("1");
        else if(!a)
        {
            if(b == t) puts("2");
            else puts("-1");
        }
        else if(a == 1)
        {
            if(!b) puts("-1");
            else 
            {
                LL ans = 1ll * power(b, p - 2, p) * (((t - x_1) % p + p) % p) % p + 1;
                printf("%lld\n", ans);
            }
        }
        else 
        {
            LL tmp = 1ll * b * power(a - 1, p - 2, p) % p;
            printf("%lld\n", Baby_Step_Giant_Step(a, 1ll * (t + tmp) * power((x_1 + tmp) % p, p - 2, p) % p, p));
        }
    }
    return 0;
}

扩展BSGS

【模板】扩展BSGS

正如我们所知道的那样，BSGS算法使用条件为：(a,p)=1；

若不满足此条件该怎么办？链接：https://oi-wiki.org/math/bsgs/

接下来一些注意事项：

写扩展BSGS，建议写一个单独的函数；
最开始a和b都要对p取模。
在判断0->k处是否符合题意（即判断是否已有答案时），其模数可以为原始p，也可以是当前的p’，不过注意一点：p和b的对应关系不要混为一谈。
不要忘记最后给BSGS计算的b/D时乘以关于p/D的逆元！！

LL ex_Baby_Step_Giant_Step(LL a, LL b, LL p)
{
    a %= p, b %= p;
    LL x, y, k = 0, d, P = p, B = b, mul = 1;
    while((d = gcd(a, P)) != 1) 
    {
        if(B % d) return -1;
        ++ k, P /= d, B /= d, mul = mul * (a / d) % P;
        if(mul == B) return k;
    }
    exgcd(mul, P, x, y);
    x = (x % P + P) % P;
    LL res = Baby_Step_Giant_Step(a, B * x % P, P);
    return res == -1 ? -1 : res + k;
}