该题题义是给定如下一个方程组:
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; }