POJ 3526 The Teacher’s Side of Math 高斯消元

The Teacher’s Side of Math

Time Limit: 5000MS Memory Limit: 65536KTotal Submissions: 680 Accepted: 190

Description

One of the tasks students routinely carry out in their mathematics classes is to solve a polynomial equation. It is, given a polynomial, say X² − 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 + ⁿ√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) = a_dX^d + a_d−1X^d−1 + ⋯ + a₁X + a₀ satisfying the following conditions.

1. Coefficients a₀, …, a_d are integers and a_d > 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 a₀, …, a_d have no common divisors greater than one.

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

F(t)= (α + β)⁴ − 10(α + β)² + 1 = (α⁴ + 4α³β + 6α²β² + 4αβ³ + β⁴) − 10(α² + 2αβ + β²) + 1 = 9 + 12αβ + 36 + 8αβ + 1 − 10(3 + 2αβ + 2) + 1 = (9 + 36 + 4 − 50 + 1) + (12 + 8 − 20)αβ = 0

Verifying 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 (a_d) is one.

Input

The input consists of multiple datasets, each in the following format.

a m b n

This line represents ^m√a + ⁿ√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.

1. ^m√a + ⁿ√b ≤ 4.

2. mn ≤ 20.

3. The coefficients of the answer a₀, …, a_d are between (−2³¹ + 1) and (2³¹ − 1), inclusive.

Output

For each dataset, output the coefficients of its minimal polynomial on integers F(X) = a_dX^d + a_d−1X^d−1 + ⋯ + a₁X + a₀, in the following format.

a_d a_d−1 … a₁ a₀

Non-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 23 2 2 32 2 3 431 4 2 33 2 2 70 0 0 0

Sample Output

1 0 -10 0 11 0 -9 -4 27 -36 -231 0 -8 0 18 0 -104 0 11 0 0 -8 -93 0 24 -2976 2883 -32 -3720 -23064 -297751 0 -21 0 189 0 -945 -4 2835 -252 -5103 -1260 5103 -756 -2183

Source

Japan 2007

设t=a^(1/m)+b^(1/n),找一个项数最少的多项式使A(d)*t^d+A(d-1)*t^(d-1)+...+A0=0.

高斯消元

观察发现一共有n*m项。

把每一个次方项二项式展开，每一项有关a^(x/m)*b^(y/n) (x<m,y<n)的因为不是整数，系数之和必须为0，才能使得这些项正好消掉。

把所有系数当成未知数，二项式展开之后正好得到一个线性方程组。接着，用高斯消元解就是了。

第一次写高斯消元，写的非常痛苦。

#include <cstdio>#include <iostream>#include <string.h>#include <string> #include <map>#include <queue>#include <vector>#include <set>#include <algorithm>#include <math.h>#include <cmath>#include <stack>#define mem0(a) memset(a,0,sizeof(a))#define meminf(a) memset(a,0x3f,sizeof(a))using namespace std;typedef long long ll;typedef double db;const int maxn=25,inf=0x3f3f3f3f;  const ll llinf=0x3f3f3f3f3f3f3f3f;   int c[maxn][maxn]; db ans[maxn],r[maxn][maxn];void gauss(int size) {int i,j,k;for (i=0;i<size;i++) {/*for (k=0;k<size;k++) {for (j=0;j<=size;j++) {cout << r[k][j] << ' ';}cout << "\n";}cout << '\n';*/k=i;for (j=i+1;j<size;j++) {     //找最大的列主元 if (fabs(r[j][i])>fabs(r[k][i])) k=j;}if (fabs(r[k][i])<1e-9) continue; for (j=i;j<=size;j++) {db z=r[k][j];r[k][j]=r[i][j];r[i][j]=z;}for (j=0;j<size;j++) {if (i==j) continue;if (fabs(r[j][i])<1e-9) continue;db l=r[j][i]/r[i][i];for (k=i;k<=size;k++) {r[j][k]-=l*r[i][k];}}}/*for (i=0;i<size;i++) {for (j=0;j<=size;j++) {cout << r[i][j] << ' ';}cout << "\n";}cout << '\n';*/ans[size-1]=1.0;for (i=size-2;i>=0;i--) {ans[i]=r[i][size]/r[i][i];}}int main() {int m,n,i,j,a,b;mem0(c);c[0][0]=1;for (i=1;i<=20;i++) {c[i][0]=1;for (j=1;j<=i;j++) {c[i][j]=c[i-1][j]+c[i-1][j-1];}}scanf("%d%d%d%d",&a,&m,&b,&n);while (a!=0||m!=0||b!=0||n!=0) {mem0(r);r[0][0]=1;for (i=1;i<=n*m;i++) {for (j=0;j<=i;j++) {db z=c[i][j]*pow(a*1.0,j/m)*pow(b*1.0,(i-j)/n);r[(j%m)*n+(i-j)%n][i]+=z;}}r[n*m][n*m]=1;r[n*m][n*m+1]=1;gauss(n*m+1);for (i=n*m;i>=0;i--) {if (ans[i]<1e-9) ans[i]-=0.5; else ans[i]+=0.5;printf("%d ",(int)(ans[i]));}printf("\n");scanf("%d%d%d%d",&a,&m,&b,&n);}return 0;}