多项式除法的C语言实现

问题的引出

在该博客中我将着手解决这样一个问题，在求解模素数$p$的同余式：
$$f(x)=a_n{x}^n+\cdots +a_{1}x+a_0\equiv 0\left(modp\right)$$
时，其中$a_n$恒不为零. 通常需要借助费马小定理和多项式的欧几里得除法将上述一个复杂的问题转化为一个次数更低的模素数同余式。

具体来说，由费马小定理可知，多项式$x^p-x\left(modp\right)$对任何整数取值为零.进一步，利用多项式欧几里得除法可知存在整系数多项式$q(x)$,$r(x)$使得：
$$f(x)=q(x)(x^p-x)+r(x),deg(r(x))<deg(g(x)).$$
从而容易将上述同余式：
$$f(x)=a_n{x}^n+\cdots +a_{1}x+a_0\equiv 0\left(modp\right)$$
降次为同余式：
$$r(x)\equiv 0\left(modp\right)$$
不难证明上述两个模$p$同余式等价.这时我便将初始的问题转化为一个次数更低的模$p$同余式，进而方便原复杂问题的进一步求解。在这个化简过程中，起到核心作用的便是多项式的欧几里得除法，所以为了最终成功编程实现模素数$p$任意次同余式的求解，我需要首先编程实现整系数的多项式欧几里得除法.

多项式欧几里得除法

引理（多项式欧几里得除法）：
设
$$f(x)=a_n{x}^n+\cdots +a_{1}x+a_0$$
为$n$次整系数多项式，
$$g(x)=x^n+\cdots +b_{1}x+b_0$$
为$m\ge 1$次首一整系数多项式，则存在整系数多项式$q(x)$,$r(x)$使得：
$$f(x)=q(x)(x^p-x)+r(x)，deg(r(x))<deg(g(x)).$$
证 :分以下两种情形讨论：
$(i)$ $n<m$.取$q(x)=0$,$r(x)=f(x)$,结论成立.
$(i)$ $n\ge m$.对$f(x)$的次数$n$做数学归纳法.
对于 $n=m$,有
$$f(x)-a_{n}g(x)=(a_{n-1}-a_n{b}_{m-1})x^{n-1}+\cdots +(a_1-a_n{b}_1)x+a_0-a_{n}b$$
因此，$q(x)=a_n,r(x)=f(x)-a_{n}g(x)$为所求.
假设$n-1\ge m$时，结论成立.
对于$n>m$，有
$$f(x)-a_n{x}^{n-m}g(x)=(a_{n-1}-a_n{b}_{m-1})x^{n-1}+\cdots +(a_{n-m}-a_n{b}_1)x^{n-m}+a_{n-m-1}{x}^{n-m-1}+\cdots +a_0-a_{n}b$$
这说明$f(x)-a_n{x}^{n-m}g(x)$是次数$\le n-1$的多项式。对其运用归纳假设，存在整系数多项式$q_1(x)$,$r_1(x)$使得：
$$f(x)-a_n{x}^{n-m}g(x)=q_1(x)g(x)+r_1(x)，deg(r_1(x))<deg(g(x)).$$
因此，$q(x)=a_n{x}^{n-m}g(x)+q_1(x),r(x)=r_1(x)$为所求.

多项式欧几里得除法的C语言实现

我利用整形数组来代替多项式本身进行运算，进行整型数组的加减法和乘法运算，本质上是从右至左进行运算.而用多项式的系数进行运算时为了方便，我们总是默认将多项式的系数从高次至低次依次输入数组，这样直接调用数组元素会加重代码复杂度，降低代码的可读性. 因此，我想到一个解决这个问题简单的方法，就是在进行加减法和乘法之前，现将初始数组倒置，这样就可以直观的进行运算，最终将数组再次倒置，便于阅读.

倒置函数

int* Reverse(int *a, int Len)
{
	int tmp;
	for (int i = 0; i < Len / 2; i++)
	{
		tmp = a[i];
		a[i] = a[Len - 1 - i];
		a[Len - 1 - i] = tmp;
	}
	return a;
}

多项式加法

int* PA(int * a, int Len1, int * b, int Len2)
{
	a = Reverse(a, Len1);//倒置数组a
	b = Reverse(b, Len2);//倒置数组b
	int Len3 = Max(Len1, Len2);//获得用于表示和多项式的数组的长度
	int *arr = (int*)malloc(Len3 * sizeof(int));//为和多项式的数组开辟动态内存空间
	if (Len1 == Len3)//将多项式最高次的次数高的数组赋给新的数组
		for (int i = 0; i < Len1; i++)
			arr[i] = a[i];
	else
		for (int i = 0; i < Len3; i++)
			arr[i] = b[i];
	if (Len1 != Len3)//将多项式最高次的次数低的数组与新数组相加
		for (int i = 0; i < Len1; i++)
			arr[i] += a[i];
	else
		for (int i = 0; i < Len2; i++)
			arr[i] += b[i];
	arr = Reverse(arr, Len3);//倒置和多项式的数组，以便该数组直接表示和多项式从高至低的系数
	a = Reverse(a, Len1);//保持输入数组a不变
	b = Reverse(b, Len2);//保持输入数组b不变
	return arr;
}

多项式减法

int* PS(int * a, int Len1, int * b, int Len2)
{
	a = Reverse(a, Len1);//倒置数组a
	b = Reverse(b, Len2);//倒置数组b
	int Len3 = Max(Len1, Len2);//初步获得用于表示余式的数组的长度（可能比这个长度低，因为数组a和数组b的前多个髙次项系数和次数完全一样）
	int *arr = (int*)malloc(Len3 * sizeof(int));//为余式的数组开辟动态内存空间
	if (Len1 == Len3)//若数组a的长度大于等于数组b的长度，只需将数组a赋给新数组，再将数组a的前Len2个元素分别减掉数组b对应的元素
	{
		for (int i = 0; i < Len1; i++)//
			arr[i] = a[i];
		for (int i = 0; i < Len2; i++)
			arr[i] -= b[i];
	}
	else//若数组a的长度小于数组b的长度，先将数组a赋给新数组的前Len1个元素，再用数组a的元素分别减掉数组b对应的前len1个元素，最终将数组b的后Len3-Len1个元素置负后赋给新数组的后半部分
	{
		for (int i = 0; i < Len1; i++)
			arr[i] = a[i];
		for (int i = 0; i < Len1; i++)
			arr[i] -= b[i];
		for (int i = Len1; i < Len2; i++)
			arr[i] = -b[i];
	}
	arr = Reverse(arr, Len3);//倒置余式的数组，以便该数组直接表示余式从高至低的系数
	//由于数组a和数组b的前多个髙次项系数和次数完全一样，所以上述操作结束后余式数组arr可能前几项元素为0，即余式系数不再等于输入多项式的最高次项。
	//而输出的数组应该保留至首个非零的元素，即将系数为零的高次项消除。所以我们需要将余式从高次至低次检查系数是否为零，确保返回的余式首系数不为零。
	int count = 0;//用于记录余式数组arr需要出掉多少个零元素
	for (int i = 0; i < Len3; i++)//从左至右判断arr元素是否为零
	{
		if (arr[i] == 0)
			count++;
		else
			break;//一旦出现非零元素便跳出循环
	}
	if (count != 0)//判断数组arr所代表的余式是否需要保留首个系数不为零及之后的项
	{
		int *brr = (int*)malloc((Len3 - count) * sizeof(int));//为首个系数不为零及之后的项分配内存
		for (int i = 0; i < Len3 - count; i++)//从数组arr首个不为零的元素开始依次赋值给数组brr
			brr[i] = arr[i + count];
		free(arr);//释放最初数组arr
		arr = brr;//将数组brr的首地址赋给arr
		brr = NULL;//将数组brr置空
	}
	a = Reverse(a, Len1);//保持输入数组a不变
	b = Reverse(b, Len2);//保持输入数组b不变
	return arr;
}

多项式乘法

int* PM(int * a, int Len1, int * b, int Len2)
{
	a = Reverse(a, Len1);//倒置数组a
	b = Reverse(b, Len2);//倒置数组b
	int *arr = (int*)calloc(Len1 + Len2, sizeof(int));//为积多项式的数组开辟动态内存空间,并顺便将每个元素赋值为0
	for (int i = 0; i < Len2; i++)//用数组b的元素分别去乘数组a
		for (int j = 0; j < Len1 ; j++)
			arr[j+i] += b[i] * a[j];//并向左依次移一维
	arr = Reverse(arr, Len1 + Len2 - 1);//倒置积多项式的数组，以便该数组直接表示积多项式从高至低的系数
	a = Reverse(a, Len1);//保持输入数组a不变
	b = Reverse(b, Len2);//保持输入数组b不变
	return arr;
}

多项式除法

int* PD(int * a, int Len1, int * b, int Len2)
{
	if (b[0] != 1)//除式首系数不为1，则该程序无法解决这样的多项式除法
		return 0;
	int count;//记录每次做差时，被减式与减式从左至右连续相同的位数
	int Len11 = Len1;//初始化被减式的位数
	int count2 = 0;//用于记录除法进行到哪一位
	int *Q = NULL;//初始化商式
	int *tmp1 = NULL;//初始化减式
	int* tmp2 = a;//初始化被减式（刚开始即为数组a）
	int *tmp3 = a;//初始化中间被差式(也是每一步的余式）
	while (Len11 >= Len2)//判断余式的最高次数是否小于除式的最高次数，若小于，那么多项式除法结束并返回最新的余式；若大于，那么多项式除法继续进行
	{
		printf("第 %d 次循环\n", count2);
		count = 0;//每次循环需要重新初始化被减式与减式从左至右连续相同的位数
		if (a[count2] == 0 && tmp2[1]==0)//判断这一循环位是否商零，若商零，用于记录除法进行到哪一位的count2自动加一，且直接进入下一个循环
		{
			count2++;
			continue;
		}
		else
		{
			printf("Len11 is %d\n", Len11);
			printf("count2 is %d\n", count2);
			Q = (int*)calloc(Len11 + 1 - Len2, sizeof(int));//分配本次循环的商式的内存。
			Q[0] = tmp2[0];//将被减式的首元素赋给商式的首元素
			printf("Q is \n");
			for (int j = 0; j < Len11 + 1 - Len2; j++)
				printf("%d", Q[j]);
			printf("\n");
			tmp1 = (int*)malloc((Len11) * sizeof(int));//为本次循环的减式分配内存
			printf("LenQ is %d\n", Len11 + 1 - Len2);
			tmp1 = PM(Q, Len11 + 1 - Len2, b, Len2);//计算本次循环的减式
			printf("tmp1 is \n");
			for (int j = 0; j < Len11; j++)
				printf("%d", tmp1[j]);
			printf("\n");
			for (int i = 0; i < Len11; i++)//记录本次循环被减式与减式从左至右连续相同的位数
				if (tmp1[i] == tmp2[i])
					count++;
				else
					break;

			printf("count is %d\n", count);
			printf("count2 is %d\n", count2);
			tmp3 = (int*)calloc(Len11 - count, sizeof(int));//为本次循环的余式分配内存
			printf("Len11 is %d\n", Len11);
			printf("Lentmp3 is %d\n", Len11 - count);
			tmp3 = PS(tmp2, Len11, tmp1, Len11);//计算本次循环的余式
			printf("tmp3 is \n");
			for (int i = 0; i < Len11 - count; i++)
				printf("%d", tmp3[i]);
			printf("\n");
			
			tmp2 = NULL;//每次循环到达这一步时，我们已经计算出了下一次循环的被除式，即tmp3，因此我们把旧的被除式置空
			tmp2 = tmp3;//并将新的被除式tmp3赋给tmp2
			printf("tmp2 is \n");
			for (int i = 0; i < Len11 - count; i++)
				printf("%d", tmp2[i]);
			printf("\n");
			Len11 = Len11 - count;//计算下次循环减式的长度
			
			count2++;
			printf("\n");
			printf("\n");
		}
	}
	return tmp2;//返回最新的余式，
}

主程序调用

我在这里会举出一个实际的例子：
求同余式：
$$f(x)=3x^{14}+4x^{13}+2x^{11}+x^9+x^6+x^3+12x^2+x\equiv 0(mod5),$$
等价的次数小于5的同余式：$r(x)$.

int main()
{
	int a[] = { 3,4,0,2,0,1,0,0,1,0,0,1,12,1,0 };
	int b[] = { 1,0,0,0,-1,0 };
	int Len1 = sizeof(a) / sizeof(a[0]);
	int Len2 = sizeof(b) / sizeof(b[0]);
	int *f = PD(a, Len1, b, Len2);
	for (int i = 0; i < 4; i++)
		printf("%d ", f[i]);
	free(f);
	_getch();
	return 0;
}

运行结果

也就是说：
$$f(x)=3x^{14}+4x^{13}+2x^{11}+x^9+x^6+x^3+12x^2+x\equiv 0(mod5),$$
在模5的意义下与:
$$3x^3+16x^2+6x$$
等价.在这个问题中我们的除式选择为：
$$x^5-x$$
原因只需要用费马小定理解释即可！

注：在解决这个问题的过程中，我使用的编译器为Visual Studio 2017,
不同版本的编译器在使用我的代码的过程中可能需要简单调试，这里因给大家造成不便而表示歉意.