zoukankan      html  css  js  c++  java
  • 数论——线性同余方程

    下文讲的解都是整数解

    裴蜀定理和扩展欧几里得算法

    裴蜀恒等式:方程$ax+by=gcd(a,b)$存在一组特解$x=x0,y=y0$。

    证明:

    其实我们可以通过欧几里得算法来证明。

    设$gcd(a,b)=g$。

    欧几里得算法的最后一步会得到两个数$g, 0$。

    然后令$x0=1,y0=0$,就得到了一组解$g imes 1 + 0 imes 0 = g$。

    假设当前我们要求一组解满足$ax+by=g$,且我们已经得到了一组解满足$bx0+(a-[a/b]*b)y0=g$。

    拆开来得到:$ay0+b(x0-[a/b]*y0)=g$,所以我们可以令$x1=y0, y1=x0-[a/b]*b*y0$而得到一组解$x1,y1$。

    一直这么推下去就可以得到最后的答案。

    证毕。

    扩展欧几里得算法:求$ax+by=gcd(a,b)$的一组通解。

    要想求通解求出一组特解(x0,y0)即可。

    通解可表示为$x=x0+frac{b}{g}, y=y0-frac{a}{g}$,即两边同时减去最小公倍数。

    要求特解我们发现就是刚才裴蜀恒等式的证明。

    代码:

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

    求解一般不定方程:

    给定一个不定方程$ax+by=c$,求其通解。

    容易发现如果$gcd(a,b)$不整除c则无解。

    如果整除可以先算出一组方程$ax+by=gcd(a,b)$的特解$x0,y0$,然后再乘上$frac{c}{d}$。

    例题:CF510D Fox And Jumping

    思路:

    根据裴蜀定理,ax+by能表示出所有数当且仅当$(a,b)=1$。不难知道扩展到n项也是成立的。

    所以这变成了一个背包问题,设$dp_i$代表选出的数的最大公约数为$i$的最小价值。最后的答案为$dp_1$。

    我们可以通过$dp_j$来更新$dp_{gcd(j,c[i])}$。

    但gcd可能很大,不能直接枚举。所以我们可以用平衡树把所有已经可能出现的gcd存下来,然后用其来更新。

    注意自己也是可以单独的一个gcd的,这是初始化。

    #include <bits/stdc++.h>
    using namespace std;
    int n, l[310], c[310];
    map<int, int> mp;
    int main() {
        scanf("%d", &n);
        for (int i = 1; i <= n; i++) scanf("%d", &l[i]);
        for (int i = 1; i <= n; i++) scanf("%d", &c[i]);
        mp[l[1]] = c[1];
        for (int i = 2; i <= n; i++) {
            if (mp.find(l[i]) == mp.end()) mp[l[i]] = c[i];
            else mp[l[i]] = min(mp[l[i]], c[i]);
            for (map<int, int>::iterator it = mp.begin(); it != mp.end(); it++) {
                int val = it -> first;
                val = __gcd(val, l[i]);
                if (mp.find(val) == mp.end()) mp[val] = mp[it -> first] + c[i];
                else mp[val] = min(mp[val], mp[it -> first] + c[i]);
            }
        }
        if (mp.find(1) != mp.end()) printf("%d", mp[1]);
        else puts("-1");
        return 0;
    }
    CF510D

    整除&同余

    整除的定义:若一个数x乘上另一个数y可以得到z,则称x整除z,记作$x|z$。

    一些性质(因为markdown不好打就直接放图了):

    总结下来就是:

    1,乘上一个自然数整除性不变(类似等式)。

    2,传递性。

    3,如果有互质的部分可以去除。

    4,整除不看多于的系数。

    5,整除也可以推到等量关系。

    模运算的性质:

    都挺好理解的。

    注意除法是不能直接模的。

    同余的性质:

    这和上面没啥区别啊。

    完全剩余系:

    定义模n余0的所有数在一个集合,模n余1的所有数再一个集合,以此类推。得到模n的完全剩余系。

    简化剩余系:所有与n互质的同余类组成的集合。

    简化剩余系关于模n乘法封闭。

    简化剩余系可以证明欧拉定理,这里就不再赘述。

    线性同余方程

    解形如$ax equiv b (mod{m})$的方程。

    根据上面同余的理论可以得到$ax-b equiv 0 (mod m)$。

    所以$m|ax-b$,设$ax-b=-mk$,可以得到一个不定方程$ax+mk=b$。

    解这个不定方程即可。

    我们从这可以得出如果$gcd(a,m) mid b$方程无解

    模板题:同余方程

    Code:

    #include <iostream>
    #include <cstdio>
    #include <cstring>
    using namespace std;
    typedef long long ll;
    ll a, b;
    ll x0, y0;
    ll exgcd(ll n, ll m, ll &x, ll &y) {
        if (m == 0) {
            x = 1;
            y = 0;
            return n;
        }
        ll g = exgcd(m, n % m, x, y);
        ll tmp = x;
        x = y;
        y = tmp - (n / m) * y;
        return g;
    }
    int main() {
        scanf("%lld%lld", &a, &b);
        ll g = exgcd(a, b, x0, y0);
        x0 = (x0 % b + b) % b;
        if (x0 == 0) x0 += b;
        printf("%lld", x0);
        return 0;
    }
    同余方程

    乘法逆元

    同过线性同余方程我们可以求乘法逆元。

    即解方程:$ax=1(mod b)$。

    次方程有解当且仅当$gcd(a,b)=1$。

    如果b是质数,则还可用费马小定理。

    例题:[LuoguP1593]因子和

    我们知道对于一个整数n,可以分解成$p1^{c1} imes p2^{c2} imes cdots imes pm^{cm}$。

    其因子和应为$Pi_{i=1}^{m}sum_{j=0}^{c_i}p^{j}$。

    然后这道题就很好做了。

    答案应为$Pi_{i=1}^{m}sum_{j=0}^{{c_i}^{b}}p^{j}$。我们发现$sum$里面即为等比数列求和。直接套公式为$s=frac{p^{{c_i}^b+1}-1}{p-1}$。

    因为有除法取模所以要套逆元。

    #include <iostream>
    #include <cstdio>
    #include <cmath>
    using namespace std;
    const long long mod = 9901;
    bool mark[1000010];
    long long prime[1000010], tot;
    long long p[1000010], c[1000010], cnt;
    long long a, b, ans = 1;
    long long ksm(long long x, long long y, long long m) {
        long long ret = 1;
        while (y) {
            if (y & 1) ret = (ret * x) % m;
            x = (x * x) % m;
            y >>= 1;
        }
        return ret;
    }
    long long ni(long long x) {
        return ksm(x, mod - 2, mod);
    }
    void init(long long n) {
        for (long long i = 2; i <= n; i++) {
            if (!mark[i]) prime[++tot] = i;
            for (long long j = 1; j <= tot; j++) {
                if (prime[j] * i > n) break;
                mark[prime[j] * i] = 1;
                if (i % prime[j] == 0) break;
            }
        }
    }
    void divide(long long n) {
        for (long long i = 1; i <= tot && prime[i] * prime[i] <= n; i++) {
            while (n % prime[i] == 0) {
                p[++cnt] = prime[i];
                while (n % prime[i] == 0) {
                    c[cnt]++;
                    n /= prime[i];
                }
            }
        }
        if (n > 1) {
            p[++cnt] = n;
            c[cnt] = 1;
        }
    }
    int main() {
        cin >> a >> b;
        if (a == 0) {
            puts("0");
            return 0;
        }
        if (b == 0) {
            puts("1");
            return 0;
        }
        init(sqrt(a));
        divide(a);
        for (long long i = 1; i <= cnt; i++) {
            if ((p[i] - 1) % mod == 0) ans = (ans * (b * c[i] + 1)) % mod;
            else ans = (ans * (ksm(p[i], c[i] * b + 1, mod) - 1 + mod) * ni(p[i] - 1)) % mod;
        }
        cout << ans;
        return 0;
    }
    因子和

    习题:排列计数

    扩展——逆元线性求法。

    求$1,2cdots,n$模$p$的逆元。

    设$inv_i$代表$i$模$p$的逆元。

    设$p=ki+j, 0 le j < i$。

    $ki+j equiv 0 (mod p)$

    $(ki+j)i^{-1}j^{-1} equiv 0 (mod p)$

    $kj^{-1}+i^{-1} equiv 0 (mod p)$

    $i^{-1} equiv -kj^{-1} (mod p)$

    $inv_i = i^{-1} equiv (p-[p/i]) imes inv_{(p mod i)} (mod p)$

    我们可以得到下面的代码:

    inv[1] = 1;
    for (int i = 2; i <= n; i++) {
        inv[i] = (p - p / i) * inv[p % i] % p;
    }

    线性同余方程组

    求解一个同余方程组:

    中国剩余定理(crt)

    物不知数问题:

    有物不知其数,三三数之剩二,五五数之剩三,七七数之剩二。问物几何?

    化解为线性方程组来表示即为:

    $egin{cases}
    x equiv 2 (mod 3)\
    x equiv 3 (mod 5)\
    x equiv 2 (mod 7)\
    end{cases}
    $

    对于这个方程组我们可以对于每个方程构造一个数使得其模其它模数为0,模自己的模数为答案。

    即对于上面的方程可以构造:

    $egin{cases}
    a_1=35\
    a_2=63\
    a_3=30\
    end{cases}
    $

    全部加起来即可得到一组特解$x0=128$。然后得到通解$x=x0+Mt, t in $。M是模数的最小公倍数。

    可以得到最小正整数解:$x=23$。

    通过这个例子我们可以得到下面的结论:

    如果模数两两互素,设$M=Pi_{i=1}^{n}m_i, M_i=frac{M}{m_i}$。根据古人的智慧我们可以构造出一个特解:$x0=sum_{i=1}^{n}b_iM_i{M_i}{-1}$。

    容易证明这个特解是对的,然后我们就可以得到通解:$x=x0+M$。(M即为lcm(m1, m2, ..., mn))。

    多数情况下模数不是质数,所以要用exgcd求逆元。

    习题:[TJOI2009]猜数字(模板)

    Code:

    #include <bits/stdc++.h>
    using namespace std;
    typedef long long ll;
    ll n;
    ll a[11], b[11];
    ll M = 1, m[11], t[11];
    ll ans;
    ll exgcd(ll u, ll v, ll &x, ll &y) {
        if (v == 0) {
            x = 1;
            y = 0;
            return u;
        }
        ll g = exgcd(v, u % v, x, y);
        ll tmp = x;
        x = y;
        y = tmp - (u / v) * y;
        return g;
    }
    ll gsc(ll x, ll y, ll mod) {
        ll ret = 0;
        while (y) {
            if (y & 1) ret = (ret + x) % mod;
            x = (x + x) % mod;
            y >>= 1;
        }
        return ret;
    }
    int main() {
        cin >> n;
        for (int i = 1; i <= n; i++) {
            scanf("%lld", &b[i]);
        }
        for (int i = 1; i <= n; i++) {
            scanf("%lld", &a[i]);
            M *= a[i];
        }
        for (int i = 1; i <= n; i++) {
            m[i] = M / a[i];
            ll x;
            exgcd(m[i], a[i], t[i], x);
            t[i] = (t[i] % a[i] + a[i]) % a[i];
            ans += gsc(gsc(b[i], m[i], M), t[i], M);
        }
        ans = (ans % M + M) % M;
        cout << ans;
        return 0;
    }

    扩展中国剩余定理(excrt)

    模数不互质的情况,刚刚的构造数错的(想一想,为什么?)。

    这个时候,我们就要每个方程轮流按顺序求解。

    假设我们已经求出一个解$x=x0+Mt, t in $。这个解满足前$i-1$个方程。现在我们将其和第$i$个方程合并。

    具体操作就是将其代入第$i$个方程来解出得到未知数$t$的一个通解。

    设第$i$个方程是形如$x equiv a (mod b)$的。

    解方程的过程:

    $x0+Mt equiv a (mod b)$

    $Mt equiv a - x0 (mod b)$

    如果这个方程无解则整个方程组无解。

    设$g=gcd(M,b)$

    解得:$t=t0+frac{b}{g}t', t' in $

    带回$x$的解,得:$x=x0+M(t0+frac{b}{g}t')=x0+Mt0+lcm(M,b)t',t' in $。

    将$x0$更新为$x0+Mt0$,$M$更新为$lcm(M,b)$即可。

    第一个方程显然可以得到一组解$x=a+bt, t in $

    Code(扩展中国剩余定理):

    #include <iostream>
    #include <cstdio>
    using namespace std;
    typedef long long ll;
    ll n;
    ll a[100010], b[100010];
    ll exgcd(ll u, ll v, ll &x, ll &y) {
        if (v == 0) {
            x = 1;
            y = 0;
            return u;
        }
        ll g = exgcd(v, u % v, x, y);
        ll tmp = x;
        x = y;
        y = tmp - u / v * y;
        return g;
    }
    //(u * v) % m;
    ll Mul(ll u, ll v, ll m) {
        ll ret = 0;
        while (v) {
            if (v & 1) ret = (ret + u) % m;
            u = (u + u) % m;
            v >>= 1;
        }
        return ret;
    }
    ll excrt() {
        ll ans = b[1], M = a[1];
        for (int i = 2; i <= n; i++) {
            ll x = M;
            ll y = a[i];
            ll z = ((b[i] - ans) % y + y) % y;
            ll u, v;
            ll g = exgcd(x, y, u, v);
            if (z % g != 0) return -1;
            u = Mul(u, z / g, y / g);
            u = (u % (y / g) + y / g) % (y / g);
            ans += u * M;
            M = M * (y / g);
        }
        return (ans % M + M) % M;
    }
    int main() {
        cin >> n;
        for (int i = 1; i <= n; i++) cin >> a[i] >> b[i];
        cout << excrt();
        return 0;
    }

    例题:屠龙勇士(有意思的模板,转化为标准形式 | NOI竟然也考模板!!!和2017一样)。

  • 相关阅读:
    VIJOS P1057盖房子 (动态规划)
    RQNOJ PID57 / 找啊找啊找GF
    RQNOJ PID302 / [NOIP2001]统计单词个数 (动态规划)
    hdu 3829 Cat VS Dog 最大独立集
    并查集 找k颗树使节点数最多
    在 Sublime Text 3 中配置编译和运行 Java 程序
    在 Sublime Text 3 中配置编译和运行 Java 程序
    StarUML license key
    StarUML license key
    测试对于list的sort与sorted的效率
  • 原文地址:https://www.cnblogs.com/zcr-blog/p/13136734.html
Copyright © 2011-2022 走看看