#include <iostream>
#include <string.h>
#include <stdio.h>
#include <math.h>

using namespace std;
typedef long long LL;
const double eps = 1e-8;
const int MOD = 999983;
const int N = 55;

struct Poly
{
    int n;
    LL a[N];
};

Poly p[25];

LL gcd(LL a,LL b)
{
    return b? gcd(b,a%b):a;
}

Poly delete_zero(Poly x)
{
    int i,j;
    Poly tmp;
    for(i=0;i<x.n && x.a[i] == 0;i++);
    for(j=0;i<x.n;i++,j++) tmp.a[j] = x.a[i];
    tmp.n = j;
    return tmp;
}

Poly poly_gcd(Poly x,Poly y)
{
    x = delete_zero(x);
    y = delete_zero(y);
    Poly yy = y,tmp;
    tmp.n = 0;
    int i=0;
    if(x.n == y.n)
    {
        double k = x.a[0] / y.a[0];
        for(i=1;i<x.n;i++)
            if(fabs(x.a[i]*1.0 - k*y.a[i]) > eps) break;
        if(i == x.n) return y;
    }
    LL g = gcd(x.a[0],y.a[0]);
    LL tx = y.a[0] / g;
    LL ty = x.a[0] / g;
    for(i=0;i<x.n;i++)
    {
        x.a[i] *= tx;
        x.a[i] %= MOD;
    }
    for(i=0;i<y.n;i++)
    {
        y.a[i] *= ty;
        y.a[i] %= MOD;
    }
    if(x.n < y.n) swap(x,y);
    for(i=1;i<y.n;i++)
        tmp.a[i-1] = ((x.a[i] - y.a[i])%MOD + MOD)%MOD;
    for(;i<x.n;i++)
        tmp.a[i-1] = x.a[i];
    tmp.n = x.n - 1;
    tmp = delete_zero(tmp);
    if(tmp.n == 0) return yy;
    return poly_gcd(y,tmp);
}

信息论作业中有一道题目要求判断两个生成函数g1(x)和g2(x)构成的(2,1,m)卷积码是否为恶性卷积码。解题的关键是判断两多项式g1(x)和g2(x)是否有非常数的公因式。一种办法是：在GF(2)域上求两多项式的最大公因式。若最大公因式不是常数，则两多项式构成的卷积码是恶性卷积码。
在matlab的help文件里查了下，没有找到现成的GF(2)域求两多项式最大公因式的函数，于是自己动手写了一个。代码附在下面。
被注释掉的代码是第一个版本。当时不知道有个现成的gfdeconv函数可以直接在GF(2)计算多项式除法，因此写得笨笨的。知道以后就改写现在的样子。不过值得注意的是，与一般的翻卷积函数不同，gfdeconv函数把向量看做升幂排列的多项式。




clear;
%在GF(2)域中求两多项式的最大公因式
%注意多项式是降幂排列的还是升幂排列的
% b(x) = x^2, a(x) = x^3 + x^2 + 1
%降幂排列表示
% b = [1 0 0];
% a = [1 0 1 1];
%升幂排列表示
b = 
[ 0 
0 1];
a = 
[1 
1 0 
1];
while(sum(b)
> 0)
    t = b;
%注意！conv和deconb函数认为多项式按降幂排列
%     [q,r] =deconv(a,b);
%     b =mod(r,2);
%     ind =find(b,1,'first');
%     b =b(ind:end);
 
%而gfconv和gfdeconv等函数认为多项式按升幂排列
%例如[0 0 1]代表x^2
    [quot,remd]
= gfdeconv(a,b);
    b = remd;
    a = t;
end
    

大家一般用熟知的欧几里德算法来算最大公约数和公因式，下面介绍一种利用矩阵算最大公因式的方法：
一，在开始前先来说几个定理，你可以先看下面的部分，等充满疑惑后再来看这一部分：
定理1：（f（x），g（x））=（f（x）±g（x），g（x）±f（x））=（f（x）±g（x），g（x））=（f（x），g（x）±f（x））
定理2：如果f（x）=ax^n+bx^n-1+……+cx，那f（x）=x*（ax^n-1+bx^n-2+……+c）。而如果同时（x，g（x））=1，那易得（f（x），g（x））=（f（x）/x，g（x））
二，
用矩阵算多项式的最大公因式中最重要的2点：

1，由定理1可得，2个多项式加减不改变最大公因式。
2，由定理2可得，移位不变最大公因式。（此时不理解没关系看完下面例子就知道了）
例1：求4x^5+4x^4+4x^3+4x^2+4x-20与5x^4+5x^3+5x^2+5x-20的最大公因式
系数矩阵
4 4 4 4 4 -20
0 5 5 5 5 -20
运用定理1，将第一行减去第二行
4 -1 -1 -1 -1   0
0  5  5   5  5 -20
再运用定理2将第一行向右移一位，在行首补上0，上下对齐

0  4 -1 -1  -1  -1
0  5  5   5  5 -20
再运用定理1使其中一行末尾出现0

0     4   -1   -1   -1  -1
0  -75  25   25  25  0   
移位

0   4    -1   -1    -1   -1
0   0   -75  25   25  25   （当然，可能已经有人发现了。其实一开始，第一行，和第二行有个公倍数4，5，，貌似把它们提出了也不要紧。对，对于一行的公倍数提出来对结果没影响，只是如果提出来，下面的运算就会出现分数，个人不太喜欢分数，所以没提）
就这样反复运用定理1，2，将一行的末尾化为0，再向右移位，最后上述矩阵化为：
0  0  0  0  1  -1
0  0  0  0  1  -1 这个就说明f（x），g（x）最后可以化成x-1，和x-1，，，公因式即为x-1
如果算到最后有一行全为0，但另一行却不全为0，，那就说明2个多项式互素。

实例四十八：最大公因式问题

问题描述：

编制利用辗转相除法求两个多项式的最大公因式。

算法思路：
设多项式的系数按照幂次由高到低的顺序存于一维数组中，多项式的最高幂次存于一变量中。

辗转相除法求两多项式去除另一个多项式，然后将所得余式变成除式，原除式变为被除式。如此反复相除，直到余式为零，最后的除式即为最大公因式。
#include<stdio.h>
#define N 10
int main()
{
	float a[N+1],b[N+1],c[N+1];
	int i,j,n,m,k,flag;
	float x,y;
	printf("Please input the highest power of two polynomials(n,m):");
	scanf("%d,%d",&n,&m);
	printf("Please enter the coefficients of the terms of the polynomial(x) in order of power from high to low:\n");
	for(i=ni>=0;i--)
		scanf("%f",&a[i]);
		printf("请按幂次由高到低输入多项式B(x)各项的系数:\n");
		for(j=m;j>=0;j--)
			scanf("%f",&b[j]);
		if(n<m)
		{
			for(i=n;i>=0;i--)
			{
				c[i] = a[i];
				a[i] = b[i];
				b[i] = c[i];
			}
			for(i=n+1;i<=m;i++)
				a[i] = b[i];
			flag = n;n = m;m = flag;
		}
		flag = 1;
		while(flag)
		{
			for(k=n-m;k>=0;k--)
			{
				x = a[k+m]/b[m];
				for(j=m-1;j>=0;j--)
				{
					y = a[k+j] - x * b[j];
					a[k+j] = (y<0.0005 && y>-0.0005)?0.0:y;
				}
			}
			for(i=m;i>=0;i--}
				c[i] = b[i];
			for(flag=0,i=m-1;i>=0;i--)
			{
				if(a[i]!=0)
					flag = 1;
				b[i] = a[i];
			}
			for(i=m;i>=0;i--)
				a[i] = c[i];
			n = m--;
			}
			printf("这两个多项式的公因式是:");
			for(i=n;i>=0;i--)
				if(i!=0)
					printf("%4.1fx^%d+",a[i],i);
				else
					printf("%4.1f\n",a[i]);
				return 0;
			}