POJ 3526 高斯消元

POJ 3526 高斯消元

One of the tasks students routinely carry out in their mathematics classes is to solve a polynomial equation. It is, given a polynomial, say X2 − 4X + 1, to find its roots (2 ± √3).
If the students’ task is to find the roots of a given polynomial, the teacher’s task is then to find a polynomial that has a given root. Ms. Galsone is an enthusiastic mathematics teacher who is bored with finding solutions of quadratic equations that are as simple as a + b√c. She wanted to make higher-degree equations whose solutions are a little more complicated. As usual in problems in mathematics classes, she wants to maintain all coefficients to be integers and keep the degree of the polynomial as small as possible (provided it has the specified root).
Please help her by writing a program that carries out the task of the teacher’s side.You are given a number t of the form:
t = m√a + n√b
where a and b are distinct prime numbers, and m and n are integers greater than

1.In this problem, you are asked to find t’s minimal polynomial on integers, which is the polynomial F(X) = adXd + ad−1Xd−1 + ⋯ + a1X + a0 satisfying the following conditions.Coefficients a0, …, ad are integers and ad > 0.
2.F(t) = 0.
3.The degree d is minimum among polynomials satisfying the above two conditions.
4.F(X) is primitive. That is, coefficients a0, …, ad have no common divisors greater than one.

For example, the minimal polynomial of √3 + √2 on integers is F(X) = X4 − 10X2 + 1. Verifying F(t) = 0 is as simple as the following (α = √3, β = √2).

F(t)= (α + β)4 − 10(α + β)2 + 1
= (α4 + 4α3β + 6α2β2 + 4αβ3 + β4) − 10(α2 + 2αβ + β2) + 1
= 9 + 12αβ + 36 + 8αβ + 1 − 10(3 + 2αβ + 2) + 1
= (9 + 36 + 4 − 50 + 1) + (12 + 8 − 20)
αβ= 0Verifying that the degree of F(t) is in fact minimum is a bit more difficult. Fortunately, under the condition given in this problem, which is that a and b are distinct prime numbers and m and n greater than one, the degree of the minimal polynomial is always mn. Moreover, it is always monic. That is, the coefficient of its highest-order term (ad) is one.

Input
The input consists of multiple datasets, each in the following format.a m b nThis line represents m√a + n√b. The last dataset is followed by a single line consisting of four zeros. Numbers in a single line are separated by a single space.Every dataset satisfies the following conditions.m√a + n√b ≤ 4.mn ≤ 20.The coefficients of the answer a0, …, ad are between (−231 + 1) and (231 − 1), inclusive.
Output
For each dataset, output the coefficients of its minimal polynomial on integers F(X) = adXd + ad−1Xd−1 + ⋯ + a1X + a0, in the following format.ad ad−1 … a1 a0Non-negative integers must be printed without a sign (+ or −). Numbers in a single line must be separated by a single space and no other characters or extra spaces may appear in the output.
Sample Input
3 2 2 2
3 2 2 3
2 2 3 4
31 4 2 3
3 2 2 7
0 0 0 0
Sample Output
1 0 -10 0 1
1 0 -9 -4 27 -36 -23
1 0 -8 0 18 0 -104 0 1
1 0 0 -8 -93 0 24 -2976 2883 -32 -3720 -23064 -29775
1 0 -21 0 189 0 -945 -4 2835 -252 -5103 -1260 5103 -756 -2183

题意:

给定一个x=a(1/m)+b(1/n)。求原方程组

题解:

F(X) = adXd + ad−1Xd−1 + ⋯ + a1X + a0

防止系数a 与 无理数a 弄混,先写成

F(X) = kdXd + kd−1Xd−1 + ⋯ + k1X + k0

而 F(t) = 0, t = a m + b n t = \sqrt[m]{a} +\sqrt[n]{b} t=ma +nb

每 一 项 , k d t d = ∑ i = 0 d C d i ( a ) i / m ( b ) ( d − i ) / n 每一项, k _dt_d= \sum_{i=0}^d {C{^i_d} } {\left(a\right)}^{i/m } {\left(b\right)}^{(d - i)/n } {} kdtd=i=0dCdi(a)i/m(b)(di)/n

令 A = a ( 1 / m ) , B = b ( 1 / n ) 令A =a{^{(1/m)}} , B = b{^{(1/n)}} A=a(1/m),B=b(1/n)

则 , k d t d = ∑ i = 0 d C d i ( A ) i ( B ) ( d − i ) 则, k _dt_d= \sum_{i=0}^d {C{^i_d} } {\left(A\right)}^{i} {\left(B\right)}^{(d - i)} {} kdtd=i=0dCdi(A)i(B)(di)

所 以 , F ( t ) = 0 = ∑ j = 0 d k j t j 所以,F(t) = 0 = \sum_{j = 0}^d{ k _jt_j} F(t)=0=j=0dkjtj

//接下去,我也不是很理解

AB不同指数项的系数为0,则列出对应的mn个方程;变量为多项式的系数ai;方程系数按照(A+B)i展开一项项算。

参考博客
\quad

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <vector>
#include <iostream>
#include <ctime>
#define ll long long
using namespace std;
const int MAX_N = 25;
const double EPS = 1e-8;
int a, m, b, n;
ll C[MAX_N][MAX_N];

typedef vector<double> vec;
typedef vector<vec> mat;

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 <= MAX_N; i++) {
        C[i][0] = C[i][i] = 1;
        for(int j = 1; j < i; j++) {
            C[i][j] = C[i - 1][j] + C[i - 1][j - 1];
        }
    }
}

ll qpow(ll A, ll B) {
    ll res = 1;
    while(B) {
        if(B & 1) res *= A;
        A *= A;
        B >>= 1;
    }
    return res;
}

//四舍五入
double round(double x) {
    return x >= 0 ? floor(x + 0.5) : ceil(x - 0.5);
}

int main() {
    init();
    while(scanf("%d%d%d%d", &a, &m, &b, &n) != EOF && a, m, b, m) {
        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;
                ll y = i;
                ll 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 ", int(round(res[i])));
        }
        printf("%d\n", int(round(res[0])));
    }
}
  • 1
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值