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;
}