POJ 3526 组合数 + 高斯消元

题意

传送门 POJ 3526

题解

对于 $F(X)=a_d{X}^d+a_{d-1}{X}^{d-1}+\cdots +a_{1}X+a_0$，已知 $a_d=1$，则未知数为 $a_{d-1},a_{d-2},\cdots ,a_0$。考虑到 $a,b$ 是不同素数，则对于解 $\alpha +\beta (\alpha =\sqrt[m]{a},\beta =\sqrt[n]{b})$，只用考虑 $\alpha$ 指数模 $m$ 以及 $\beta$ 指数模 $n$ 的情况，共 $n\times m$ 种，方程右边系数为 $0$。

对于每一个指数项

$$X^i=\sum C_i^j{\alpha }^j{\beta }^{i-j}$$

组合数打表，联立关于 $a_{d-1},a_{d-2},\cdots ,a_0$ 的 $m\times n$ 元线性方程组

$$\{\begin{array}{l}{\sum }_{(i-j)\%n\times m+j\%m=0}{a}_i{C}_i^j{a}^{j/m}{b}^{(i-j)/n}=0\\ \dots \\ {\sum }_{(i-j)\%n\times m+j\%m=d-1}{a}_i{C}_i^j{a}^{j/m}{b}^{(i-j)/n}=0\end{array}$$

其中 $a_d=1$，方程中为常数项；高斯消元求解，四舍五入答案即可。

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <vector>
#define min(a, b) (((a) < (b)) ? (a) : (b))
#define max(a, b) (((a) > (b)) ? (a) : (b))
#define abs(x) ((x) < 0 ? -(x) : (x))
#define INF 0x3f3f3f3f3f3f3f3f
#define delta 0.85
using namespace std;

#define maxn 25
typedef long long ll;
const double EPS = 1e-8;
int a, m, b, n;
ll C[maxn][maxn];

typedef vector<double> vec;
typedef vector<vec> mat;
// Ax = b
vec gauss_jordan(const mat &A, const vec &b)
{
    int n = A.size();
    mat B(n, vec(n + 1));
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            B[i][j] = A[i][j];
        }
    }
    for (int i = 0; i < n; i++)
    {
        B[i][n] = b[i];
    }
    for (int i = 0; i < n; i++)
    {
        int pivot = i;
        for (int j = i; j < n; j++)
        {
            if (abs(B[j][i]) > abs(B[pivot][i]))
            {
                pivot = j;
            }
        }
        swap(B[i], B[pivot]);
        if (abs(B[i][i]) < EPS)
            return vec();
        for (int j = i + 1; j <= n; j++)
        {
            B[i][j] /= B[i][i];
        }
        for (int j = 0; j < n; j++)
        {
            if (i != j)
            {
                for (int k = i + 1; k <= n; k++)
                    B[j][k] -= B[j][i] * B[i][k];
            }
        }
    }
    vec x(n);
    for (int i = 0; i < n; i++)
        x[i] = B[i][n];
    return x;
}
// 组合数打表
void init()
{
    for (int i = 0; i < maxn; i++)
        C[i][0] = 1;
    for (int i = 1; i < maxn; i++)
    {
        for (int j = 1; j < maxn; j++)
        {
            C[i][j] = C[i - 1][j - 1] + C[i - 1][j];
        }
    }
}
// 快速幂
ll qpow(ll x, ll n)
{
    ll res = 1;
    while (n > 0)
    {
        if (n & 1)
            res *= x;
        x *= x;
        n >>= 1;
    }
    return res;
}
//四舍五入
double round(double x)
{
	return x >= 0 ? floor(x + 0.5) : ceil(x - 0.5);
}

void solve()
{
    int d = m * n;
    mat A(d, vec(d, 0));
    vec B(d, 0);
    for (int i = 0; i <= d; i++)
    {
        for (int j = 0; j <= i; j++)
        {
            ll fa = j, fb = i - j;
            ll x = fb % n * m + fa % m, y = i, z = C[i][j] * qpow(a, fa / m) * qpow(b, fb / n);
            if (i == d)
            {
                B[x] -= z;
            }
            else
            {
                A[x][y] += z;
            }
        }
    }
    vec res = gauss_jordan(A, B);
    printf("1 ");
    for (int i = d - 1; i >= 0; i--)
        printf("%d%c", int(round(res[i])), i == 0 ? '\n' : ' ');
}

int main()
{
    init();
    while (~scanf("%d%d%d%d", &a, &m, &b, &n) && (a | m | b | n))
    {
        solve();
    }
    return 0;
}