• POJ2065 SETI 高斯消元,扩展GCD


    该题题义是给定如下一个方程组:

    F(1) = C1 (mod) P
    F(2) = C2 (mod) P
    F(3) = C3 (mod) P ...

    其中F(1) = A(1,1)*x1 + A(1, 2)*x2 + A(1, 3)*x3...

    其中A(i, j) = i ^ (j-1).

    面对这样一个方程,我们的做法就是先不管方程右端的(mod)P,因为F(i) 和 Ci 是同余的,那么在方程两端乘以一个数是没有影响的,代入某一个方程同样是允许的。于是剩下的就是高斯消元。由于解是唯一的,因此解出来的系数矩阵一定是满秩。接下来就可一用扩展gcd求解未知元了。

    代码如下:

    #include <cstdlib>
    #include <cstring>
    #include <algorithm>
    #include <cstdio>
    using namespace std;
    
    typedef long long int Int64;
    
    char s[100];
    
    int length, MOD;
    
    Int64 x[100];
    
    int _pow(int a, int b)
    { // 快速幂
        int ret = 1;
        while (b) {
            if (b & 1) {
                ret *= a;
                ret %= MOD;
            }
            a *= a;
            a %= MOD;
            b >>= 1;
        }
        return ret;
    }
    
    struct Matrix
    {
        Int64 a[100][100];
        void init()
        {
            for (int i = 1; i <= length; ++i) {
                for (int j = 1; j <= length; ++j) {
                    a[i][j] = _pow(i, j-1);
                }
                a[i][length+1] = s[i];
            }
        }
        void rswap(int x, int y, int s)
        { // 用于将后面的行交换到第一行上面来
            Int64 t;
            for (int j = s; j <= length+1; ++j) {
                t = a[x][j];
                a[x][j] = a[y][j];
                a[y][j] = t;
            }
        }
        void multi(int x, Int64 M, int s)
        { // 同一行乘以某一数,当然这时为消元准备的
            for (int j = s; j <= length+1; ++j) {
                a[x][j] *= M;
            }
        }
        void relax(int x, int y, int s)
        { // 将某一方程整体迭代进另外一个表达式
            for (int j = s; j <= length + 1; ++j) {
                a[y][j] -= a[x][j];
                a[y][j] = (a[y][j] % MOD + MOD) % MOD;
                a[x][j] = (a[x][j] % MOD + MOD) % MOD;
            }
        }
        void print()
        {
            for (int i = 1; i <= length; ++i) {
                for (int j = 1; j <= length+1; ++j) {
                    printf("%5I64d ", a[i][j]);
                }
                puts("");
            }
            puts("");
        }
    }M;
    
    Int64 GCD(Int64 a, Int64 b)
    {
        Int64 t;
        while (b) {
            t = b;
            b = a % b;
            a = t;
        }
        return a;
    }
    
    Int64 LCM(Int64 a, Int64 b)
    {
        return a / GCD(a, b) * b;
    }
    
    Int64 ext_gcd(Int64 a, Int64 b, Int64 &x, Int64 &y)
    {
        Int64 temp, ret;
        if (!b) {
            x = 1, y = 0;
            return a;
        }
        ret = ext_gcd(b, a % b, x, y);
        temp = x, x = y, y = temp - a/b*y;
        return ret;
    }
    
    void Gauss()
    {
        int i = 1, k;
        Int64 a, b, sum;
        memset(x, 0, sizeof (x));
        for (int j = 1; j <= length; ++j) {
            for (k = i; k <= length; ++k) {
                if (M.a[k][j]) break; 
            }
            if (k > length) continue;
            if (k != i) {
                M.rswap(i, k, j);
            }
            for (k = i + 1; k <= length; ++k) {
                if (M.a[k][j]) {
                    Int64 lcm = LCM(M.a[i][j], M.a[k][j]);
                    a = lcm/M.a[i][j];
                    b = lcm/M.a[k][j];
                    if (a > 1) M.multi(i, a, j);
                    if (b > 1) M.multi(k, b, j);
                    M.relax(i, k, j);
                }
            }
            ++i;
        }
        // M.print();
        // 一定会有唯一解,所以没有进行无解判定
        for (k = i-1; k >= 1; --k) { // 当然这里可以暴力进行求解,毕竟给定的素数不是很大
            Int64 xx, yy, A, B, C = 0, L, G;
            for (int j = k+1; j <= length; ++j) {
                C += x[j] * M.a[k][j];
            }
            A = M.a[k][k], B = MOD;
            C = M.a[k][length+1] - C; 
            G = ext_gcd(A, B, xx, yy);
            L = B / G;
            xx *= C / G; 
            xx = (xx % L + L) % L;
            x[k] = xx;
        }
        for (int j = 1; j <= length; ++j) {
            printf(j == 1 ? "%I64d" : " %I64d", x[j]);
        }
        puts("");
    }
    
    int main()
    {
        int N;
        scanf("%d", &N);
        while (N--) {
            scanf("%d %s", &MOD, s+1);
            length = strlen(s+1);
            for (int i = 1; i <= length; ++i) {
                if (s[i] != '*') {
                    s[i] -= ('a' - 1);
                }
                else {
                    s[i] = 0;
                }
            }
            M.init();
            Gauss();
        }
        return 0;
    }
  • 相关阅读:
    学习 Message(12): 整合鼠标 Down 消息
    合并两个 Wav 文件流的函数 回复 "刘文强" 的问题
    “博客无双”第一期拍卖活动获奖名单公告
    [获奖公告]“博客无双”12月27日第一期获奖名单
    “博客无双”活动拍卖时间调整公告
    致歉
    祝大家新年快乐
    博客园电子期刊2010年12月刊发布啦
    “博客无双”拍卖活动将于14:00开始
    2011年4月微软最有价值专家(MVP)申请截止时间:2011年1月13日
  • 原文地址:https://www.cnblogs.com/Lyush/p/2611547.html
Copyright © 2020-2023  润新知