多元多项式最大公因子及因子分解

隐藏目录

多元多项式插值方法

稠密插值
稀疏插值

Euclid算法和一般模算法

概述
$\mathbb{Z}_p[x_1,x_2,\ldots,x_n]$ 上最大公因子
多元多项式的"Mignotte"界
$\mathbb{Z}[x_1,\ldots,x_n]$ 上最大公因子

Zippel稀疏插值模算法

一个具体的例子
算法描述

求GCD的其它方法

启发式算法(Heuristic GCD)
EZ-GCD

多元多项式因子分解的Kronecker算法

利用Hensel提升的因子分解算法

概述
Extended Zassenhaus 算法
因子还原
预先确定因子的首项系数
EEZ算法

在前面的章节中,我们已经具体地讨论了各种一元多项式环上的GCD问题以及因子分解问题,本章我们讨论多元多项式的相关问题.由于我们研究的是“精确”符号代数,因此这里只考虑整系数环上的多元多项式,亦即在 $\mathbb{Z}[x_1,\ldots,x_n]$ 中.

对于多元多项式的最大公因子和因子分解问题,我们也主要采用各种同态象的方法,以求将多元问题转化为一元问题求解.

多元多项式插值方法

回忆我们在 $\mathbb{Z}[x]$ 中处理问题时所用的方法.为了有效避免系数膨胀问题,我们无论求最大公因子还是做因子分解,都采用了模素数的做法,将其化为有限域上多项式问题,并且最后用中国剩余定理以及若干特殊的技术（如因子组合,格中短矢量等）将其恢复.对于多元多项式环,不仅要模素数在 $\mathbb{Z}_p$ 中讨论问题,还要将多元问题化为一元问题,这里我们主要都是采用赋值同态的方法.即同态映射 $\Phi_{x_i-a}(f)=f\bmod (x_i-a)=f(x_1,\ldots,x_{i-1},a,x_{i+1},\ldots,x_n).$ 这样相当于给某些未定元取赋值点,最后我们需要用中国剩余定理,或者说是插值的方法还原多项式.

前面我们曾经提到过多点的Lagrange插值快速算法,鉴于我们后面的算法可能是逐点插值,甚至还有所谓的“稀疏插值”,因此本节将介绍稠密插值和稀疏插值算法[1].它们对于多元多项式最大公因子的模算法是很有用的.

稠密插值

插值问题要解决的问题是已知若干个点 $u_1,u_2,\ldots,u_k$ 和 $v_1,v_2,\ldots,v_k$ ,求多项式 $f$ 使得 $\forall i=1,2,\ldots,k$ 有 $f(u_i)=v_i$ .下面直接给出Newton插值算法,很容易验证算法的正确性.

算法1(Newton插值) 输入和输出同前描述.

赋初值 $f(x)=v_1$ , $q(x)=(x-u_1)$ ,
对于 $i=2,3,\ldots,k$ ,顺次计算 $f(x)=f(x)+\frac{q(x)(v_i-f(u_i))}{q(u_i)},\quad q(x)=(x-u_i)q(x),$
输出 $f(x)$ .

多元多项式的稠密插值基于的思想十分简单,例如我们有一个函数 $F(x_1,\ldots,x_n)$ ,可以求得它在某些点上的函数值,我们的任务是找到一个各个变元次数均不超过 $d$ 的多元多项式 $P(x_1,\ldots,x_n)$ ,使得它们在各整点上的值相同.注意到我们在这里的问题的提法是很具有一般性的,不仅可以将GCD问题化为这种形式,甚至也可以将诸如多项式的乘法等问题这样提出来.

稠密插值的基本思想就是递归地依次将各个变元插值回来,假设我们有一初值点 $(x_{10},x_{20},\ldots,x_{n0})$ ,我们首先可以固定 $x_{20},\ldots,x_{n0}$ ,而再取 $d$ 个 $x_1$ 的值,由一元插值方法(Newton法或Lagrange法)求得多项式 $P(x_1,x_{20},\ldots,x_{n0})$ ,依次再确定各变元 $x_2,\ldots,x_n$ 即可.

若我们设 $d$ 次一元多项式插值问题的复杂度为 $M(d)$ ,则易知本算法的复杂度为 $n(d+1)^{n-1}M(d)$ [2].

稀疏插值

问题引入

如果我们先用一个例子来介绍稀疏插值算法,将会对它有一个更好的了解.

从上节稠密插值算法过程我们可以看出,对 $n$ 个变元次数不超过 $d$ 的多项式的插值我们大约要对函数 $F$ 求值 $(d+1)^n$ 次.为了便于说明,我们用[2]中所用的例子:对于多项式 $P(x,y,z)=x^5z^2+x^5z+xy^4+xyz^5+y^5z,$ 其需要计算函数值 $(5+1)^3=216$ 次！事实上,这么多插值次数对于这样稀疏的多项式显然是极其不划算的.我们重复一下插值的过程以说明稀疏插值是如何减少插值次数的.

首先可以由赋值点 $(x_i,y_0,z_0),i=0,\ldots,5$ 通过稠密插值得到 $P(x,y_0,z_0)=ax^5+bx+c$ .下一步即要再选择5个 $y$ 的赋值点 $y_i$ ,计算 $P(x,y_i,z_0)$ 才能由此插值得到 $P(x,y,z_0)$ .如果是稠密插值法,此时对于每个 $y_i$ ,我们都需取6个 $x_i$ 来插值,但是如果 $y_i$ 选的值恰当(所谓恰当的意义,后文定义2和定理1将会给出说明),我们完全可以假设对于 $y_i(1\le i\le 5)$ ,也有 $P(x,y_i,z_0)=a_ix^5+b_ix+c,$ 于是,对每个 $y_i$ 只需取3个 $x$ 的赋值点来插值,通过解方程的方法计算系数 $a_i,b_i,c_i$ .

这样,我们总共计算函数值的次数为 $6+3\times 5=21$ ,从而得到了 $P(x,y,z_0)$ ,注意从这一步开始,我们已经节省了插值次数,从原来的需要求值 $6\times 6=36$ 次降到了21次.并且可设求得的 $P(x,y,z_0)$ 具有形式: $P(x,y,z_0)=ax^5+bxy^4+cxy+dy^5.$

接下来的步骤就比较顺理成章了,我们同样要再计算 $P(x,y,z_i),1\le i\le 5$ ,此时可假设各多项式具有形式 $P(x,y,z_i)=a_ix^5+b_ixy^4+c_ixy+d_iy^5,$ 则需再求值 $4\times 5=20$ 次,总共求值次数 $41$ 次.

稀疏插值算法

为了叙述方便,给出如下记号.设多项式 $P(X)=P(X_1,X_2,\ldots,X_n)\in F[X]$ , $F$ 为某个域,每个变元 $X_i$ 的次数都不超过 $d$ ,并且 $P$ 的的单项式项式为 $T$ (对于稀疏多项式有 $T\ll (d+1)^n$ ).设 $v=(v_1,\ldots,v_n)$ 为一 $n$ 元有序对,定义单项式记号 $X^v=X_1^{v_1}X_2^{v_2}\cdots X_n^{v_n},$ 则多项式 $P$ 可记为 $P=c_1X^{e_1}+c_2X^{e_2}+\cdots+c_TX^{e_T},$ 其中 $e_1,\ldots,e_T$ 均为 $n$ 元序对.

从某种意义上说,集合 $\{e_1,e_2,\ldots,e_T\}$ 反映了多项式 $P$ 的“结构”,因而我们定义下面的:

定义1(模板集) 记多项式指数的集合为 $\mathrm{skel}P=\{e_1,e_2,\ldots,e_T\},$ 简称其模板集([2]中称为skeleton).并记其在前 $k$ 维上的投影为 $\mathrm{skel}_kP=\{(e_1,\ldots,e_k)|\exists e=(e_1,\ldots,e_T)\in\mathrm{skel}P\}.$

定义2(精确求值点) 设 $(x_1,\ldots,x_n)\in F^n$ ,若 $\forall k(1<k<n)$ ,有 $\mathrm{skel}P(X_1,\ldots,X_k,x_{k+1},\ldots,x_n)=\mathrm{skel}_kP(X),$ 则称其为精确求值点(precise evaluation point).

注意到一般情况下只有 $\mathrm{skel}P(X_1,\ldots,X_k,x_{k+1},\ldots,x_n)\subset\mathrm{skel}_kP(X),$ 欲使两者相等,则必须 $x_{k+1},\ldots,x_n$ 的取值使得以 $X_1,\ldots,X_k$ 为主变元的多项式 $P$ 的各项系数( $\in F[X_{k+1},\ldots,X_n]$ )均在其上不为零.据此,我们可以给出精确求值点的概率估计.

定理1(精确求值点概率估计) 设 $P(X)$ 是一整环上的 $n$ 元多项式,其每个变元次数不超过 $d$ ,总共有T项,设 $S$ 为一有限赋值点集合, $\#S=s$ ,取求值点 $x=(x_1,\ldots,x_n)\in S^n$ ,则其不是精确求值点的概率不超过 $\frac{n(n-1)dT}{s}.$

证明根据引理5我们知道对于 $P$ 以 $X_1,\ldots,X_k$ 为主变元的某个系数多项式来说, $x_{k+1},\ldots,x_k$ 为其零点的概率不超过 $(n-k)d/s$ .因为系数最多有 $T$ 项,则概率不超过 $(n-k)dT/s$ .再对 $k=1,2,\ldots,n-1$ 求和即有概率不超过 $\frac{(n-1)dT}{s}+\frac{(n-2)dT}{s}+\cdots+\frac{dT}{s}=\frac{n(n-1)dT}{s}.$ 证毕. □

取足够大的 $s$ 可以减小此概率,这正是我们每一步稀疏插值时假设目标多项式具有同样的模板的概率依据.

下面，我们把整个插值还原的过程描述如下:

算法2(稀疏插值算法) 输入:一个可以求值的函数 $f(X_1,\ldots,X_n)$ ，精确求值点 $(x_{10},\ldots,x_{n0})$ ,各变元次数均不超过的 $d$ ,

输出:多项式 $P(X_1,\ldots,X_n)\in F[X]$ ,使得 $P=f$ .

随机任取 $d$ 个值 $x_{11},x_{12},\ldots,x_{1d}$ ,利用求得的值 $f(x_{1i},x_{20},\ldots,x_{n0})(0\le i\le d)$ ,用一元插值算法(如算法1)求出多项式 $P(X_1,x_{20},\ldots,x_{n0})$ ,
$S=\mathrm{skel}P(X_1,x_{20},\ldots,x_{n0})$ ,设 $S=\{s_1,s_2,\ldots,s_q\}$ ,
将 $i$ 顺次由 $2$ 循环到 $n$ ,做如下4-6步,
随机任取 $d$ 个值 $x_{i1},x_{i2},\ldots,x_{id}$ ,对于每个 $x_{ij}(1\le j\le d)$ ,设 $P(X_1,\ldots,X_{i-1},x_{ij},x_{i+1,0},\ldots,x_{n0})$ 具有如下形式 $p_1X^{s_1}+p_2X^{s_2}+\cdots+p_qX^{s_q},$ 取 $q$ 组赋值点 $(x_{1k},x_{2k},\ldots,x_{i-1,k})(1\le k\le q)$ 构造独立的线性方程组来求出 $p_1,\ldots,p_q\in F$ ,于是得到 $d$ 个多项式 $P(X_1,\ldots,X_{i-1},x_{ij},x_{i+1,0},\ldots,x_{n0}),(1\le j\le d)$
利用第4步求出的 $d$ 个多项式对每个系数 $p_k(1\le k\le q)$ 进行一元插值算法求出 $p_k(X_i)$ ,则 $P(X_1,\ldots,X_i,x_{i+1,0},\ldots,x_{n0})=p_1(X_i)X^{s_1}+p_2(X_i)X^{s_2}+\cdots+p_q(X_i)X^{s_q},$
$S=\mathrm{skel}P(X_1,\ldots,X_i,x_{i+1,0},\ldots,x_{n0})$ ,设 $S=\{s_1,\ldots,s_q\}$ ,
输出 $P$ .

注1 在算法第4步取 $q$ 组赋值点时,要使得它们对于要求解的变元 $p_1,\ldots,p_q$ 构成非奇异的独立线性方程组,否则需要重新选取某些赋值点.

注2 我们可以选取一组赋值点 $(x_{1},x_{2},\ldots,x_{i-1})$ ,然后令 $q$ 组赋值点为 $(1,\ldots,1),(x_{1}^1,\ldots,x_{i-1}^1),(x_{1}^2,\ldots,x_{i-1}^2),\ldots,(x_{1}^{q-1},\ldots,x_{i-1}^{q-1}),$ 这样对需求解的 $q$ 个变元来说，可以得到如下Vandermonde线性方程 $\begin{pmatrix} 1 &1 &\cdots &1\\ X^{s_1} &X^{s_2} &\cdots &X^{s_q}\\ \vdots &\vdots & &\vdots\\ (X^{s_1})^{q-1} &(X^{s_2})^{q-1} &\cdots &(X_{s_q})^{q-1} \end{pmatrix} \begin{pmatrix} p_1\\p_2\\ \vdots\\p_q \end{pmatrix} = \begin{pmatrix} f(1,\ldots,1,x_{ij},x_{i+1,0},\ldots,x_{n0})\\ f(x_1,\ldots,x_{i-1},x_{ij},x_{i+1,0},\ldots,x_{n0})\\ \vdots\\ f(x_1^{q-1},\ldots,x_{i-1}^{q-1},x_{ij},x_{i+1,0},\ldots,x_{n0}) \end{pmatrix}.$

利用此方法的好处是求解线性方程时,Vandermonde方程其系数矩阵具有特殊性,因而有特别的线性代数的处理方法,求解复杂度低于一般的线性方程.另一方面,Vandermonde行列式的非奇异性很容易得到保证,只需 $X^{s_1},\ldots,X^{s_q}$ 各不相同即可.

对于特征为零的域,这一点是很容易做到的,只需将 $x_1,\ldots,x_{i-1}$ 取为前 $i-1$ 个素数即可,或者当有限域的特征足够大(特征 $p>(2\times 3\times\cdots\times p_n)^d$ )时也可如此取赋值点,此时算法仅可能当初始点非精确赋值点时失败,因此根据定理1,我们知道失败的概率为 $\varepsilon<\frac{n(n-1)dT}{s}.$

另一方面,当域的特征 $p$ 有限时,此时我们随机选取这样 $i-1$ 个赋值点,当其非精确赋值点或Vandermonde行列式奇异时均会失败,其概率([2]P240-242) $\varepsilon<\frac{n(n-1)dT}{p-1}+\frac{dT(T-1)}{2(p-1)}<\frac{dT(2n^2+T)}{p-1}.$

Euclid算法和一般模算法

概述

我们已经解决了一元多项式的GCD问题,对于多元情形,我们仍然只考虑整系数环上的多项式,即 $\mathbb{Z}[X]=\mathbb{Z}[x_1,\ldots,x_n]$ 中的多项式.首先根据Guass引理,我们有: $\begin{align*} \gcd(f,g)&=\mathrm{cont}_{x_n}(\gcd\nolimits_{x_n}(f,g))\mathrm{pp}_{x_n}(\gcd\nolimits_{x_n}(f,g))\\ &=\gcd(\mathrm{cont}_{x_n}(f),\mathrm{cont}_{x_n}(g))\gcd\nolimits_{x_n}(\mathrm{pp}_{x_n}(f),\mathrm{pp}_{x_n}(g)). \end{align*}$

于是,我们顺利地将 $n$ 元问题化为了 $n-1$ 元子问题.将此过程递归地进行,最终化为一元问题可求解.显而易见,这种算法系数的增长是十分迅速的,不宜采用.

回忆前面处理一元问题采用模算法的思想,我们希望利用 $\mathbb{Z}[x_1,\ldots,x_n]$ 上的模算术来简化问题的计算.若我们取一个一次多项式 $p=y-a(a\in\mathbb{Z})$ , $D$ 为一UFD,考虑 $R=D[y][x]$ 上的多项式,并记 $R$ 到 $R/\idea{p}$ 的同态像为 $\Phi_p(f)$ 或 $\overline{f}$ ,则有下面的定理:

定理2 $f,g\in R$ 均为本原多项式, $h=\gcd(f,g)$ ,设 $\overline{h}\neq 0$ ,则 $\overline{h}=h_p(=\gcd(\overline{f},\overline{g})$ 当且仅当 $p=y-a\not|\mathrm{res}_x(f/h,g/h)$ .

证明本定理基本上是定理15的推广,证明也与其类似.

首先根据模同态的性质,我们显然有 $\overline{h}|h_p$ ,若 $\deg_xh_p=0$ ,则 $\overline{h}=h_p$ 显然.设 $\deg_xh_p\ge 1$ ,此是有 $s,t$ 使得 $sf/h+tg/h=\mathrm{res}_x(f/h,g/h)$ ,因此 $\overline{s}\overline{f}+\overline{t}\overline{g}=\overline{\mathrm{res}_x(f/h,g/h)}\overline{h}$ ,由于 $\overline{\mathrm{res}_x(f/h,g/h)}\neq 0$ ,所有 $h_p|\overline{h}$ ,即 $\overline{h}=h_p$ ,充分性得证.(必要性证明略去) □

于是我们可以得到类似于一元情形的一种赋值同态模算法[3],这里不再详细将算法列出,然而我们需要注意的是这里会有一个领项系数的问题,例如两多项式的最大公因子为 $h(x,y)=(y-1)x$ ,对 $y$ 进行赋值同态时我们可能会取如下值: $y=5,\quad h(x,5)=4x,$ $y=7,\quad h(x,7)=6x,$ 这样,通过两次赋值进行插值即得 $h(x,y)=(y-1)x$ ,然而若是我们“不幸”取了一个负值: $y=-5,\quad h(x,-5)=6x,$ 注意在 $\mathbb{Z}[x]$ 中 $\pm 1$ 为可逆元,因而求其中的GCD问题时可相差正负号,一般情况我们取首项系数为正,此时将 $h(x,-5)=6x$ 和 $h(x,5)=4x$ 插值则得不到正确的结果,只有取 $h(x,-5)=-6x$ 时才是正确的.这个问题仅在最大公因子对主变元 $x$ 的领项系数是一非平凡多项式才会出现,对于这种情况,若要比较好地处理,则须在赋值之前对领项系数进行处理.

我们同时注意到 $\mathbb{Z}_p$ 作为一个域的简单性和其能有效抑制 $\mathbb{Z}$ 上系数膨胀问题的性质,我们可以取模素数和赋值同态的综合算法求解.

$\mathbb{Z}_p[x_1,x_2,\ldots,x_n]$ 上最大公因子

根据前面的分析,我们先要解决 $\mathbb{Z}_p[x_1,x_2,\ldots,x_n]$ 上求最大公因子的问题,此要先通过赋值同态将其转化为 $\mathbb{Z}_p[x_1]$ 上的问题.

为了说明算法的方便,本算法中提到的容度 $\mathrm{cont}$ ,本原部分 $\mathrm{pp}$ ,首项系数 $\lc$ 等都是就多项式环 $\mathbb{Z}_p[x_n][x_1,\ldots,x_{n-1}]$ 而言的,亦即 $\mathrm{cont}$ 及 $\lc$ 都应是 $\mathbb{Z}_p[x_n]$ 中的元素.

算法3(有限域上多元多项式最大公因子算法) 输入: $f$ , $g\in\mathbb{Z}_p[x_1,\ldots,x_n]$ ,

输出: $h=\gcd(f,g)$ .

若 $n=1$ 则是一元问题,直接调用相关算法求得首一的 $h=\gcd(f,g)$ ,输出 $h$ ,
初始化 $a=\gcd(\mathrm{cont}(f),\mathrm{cont}(g))\in\mathbb{Z}_p[x_n]$ (注意这是一元问题), $f=\mathrm{pp}(f)$ , $g=\mathrm{pp}(g)$ , $b=\gcd(\plc{f},\plc{g})\in\mathbb{Z}_p[x_n]$ ,
赋初值 $q=1$ , $h=1$ , $m=\min(\deg_1f,\deg_1g)+1$ , $vlist=\{\}$ ,
循环做下面5-10步,
随机任取 $v\in\mathbb{Z}_p$ 使得 $b(v)\neq 0$ 且 $v\not\in vlist$ ,将 $v$ 添入 $vlist$ 中,
$f_v=f\bmod (x_n-v)$ , $g_v=g\bmod (x_n-v)$ ,递归调用本算法求解 $n-1$ 元子问题 $h=\gcd(f_v,g_v)$ ,令 $m_1=\deg_1h_v$ , $b_v=b(v)$ ,
将 $h_v$ 正则化使得 $\plc{h_v}=b_v$ ,即令 $h_v=b_vh_v/\plc{h_v}$ ,
若 $m_1<m$ 则令 $q=x_n-v$ , $h=h_v$ , $m=m_1$ ,
若上一步不成立且 $m_1=m$ ,则利用中国剩余定理或插值算法其出 $h$ ,使得 $h\bmod q$ 不变且 $h\bmod (x_n-v)=h_v$ ,即令 $h=h+\frac{(h_v-h\bmod(x_n-v))q(x_n)}{q(v)},\quad q(x_n)=(x_n-v)q(x_n),$
若 $\plc{h}=b$ 则:令 $pph=\mathrm{pp}(h)$ ,若 $pph|f$ 且 $pph|g$ 则输出 $a pph$ ,否则若 $m=0$ 则输出 $a$ .

多元多项式的"Mignotte"界

在处理整系数一元多项式最大公因子时我们曾经引进所谓的Mignotte界,这是Mignotte于1974年提出的对一元多项式因子系数界的估计.这一节我们来讨论对于整系数多元多项式同样的问题.为了后文叙述的方便,我们先回忆一下对于一元情形的Mignotte界,并表述为如下定理(见定理13):

定理3(Mignotte) 设非零多项式 $f,g\in\mathbb{Z}[x]$ 且 $g|f$ , $\deg g\le d$ ,则有 $\|f\|_2\ge 2^{-d}\|g\|_{\infty}.$

为了讨论多元情形,我们给出:

定义3(多元多项式p-范数) 设有多元多项式 $f=\sum_{i_1,i_2,\ldots,i_n}a_{i_1,i_2,\ldots,i_n}x_1^{i_1}\cdots x_n^{i_n},$ 则定义其p-范数为: $\|f\|_p=\left(\sum_{i_1,\ldots,i_n}a_{i_1,\ldots,i_n}^p\right)^{1/p}.$

Coron于2004年提出了二元情形下对Mignotte界的推广,有如下定理(见[4]和[5]3.2节):

定理4(Coron) 设有非零多项式 $f,g\in\mathbb{Z}[x,y]$ ,且 $g|f$ , $d=\max(\deg_xf,\deg_yf)$ ,则有 $\|f\|_2\ge 2^{-(d+1)^2}\|g\|_{\infty}.$

证明令 $f^*(x)=f(x,x^{d+1})$ ,则有 $\deg f^*\le(d+1)^2$ 且 $f^*(x)$ 和 $f(x,y)$ 有相同的整系数,因而 $\|f^*\|_2=\|f\|_2$ ,同理令 $g^*(x)=g(x,x^{d+1})$ ,则有 $\|g^*\|_{\infty}=\|g\|_{\infty}$ .并且有 $g(x,y)|f(x,y)$ 可知 $g^*(x)|f^*(x)$ ,因此根据Mignotte界有 $\|f\|_2\ge 2^{-(d+1)^2}\|g\|_{\infty},$ 证毕. □

我们很容易再从二元情况推广到多元情况,有如下定理:

定理5(多元多项式因子系数界) 设有非零多项式 $f,g\in\mathbb{Z}[x_1,\ldots,x_n]$ ,且 $g|f$ , $d=\max\{\deg_1f,\ldots,\deg_nf\}$ ,则有 $\|f\|_2\ge 2^{-(d+1)^n+1}\|g\|_{\infty}.$

证明过程是类似的,只需作相应的替换 $f^*(x)=f(x,x^{d+1},\ldots,x^{(d+1)^{n-1}})$ ,并且注意到 $f^*$ 的次数实际上是不超过 $d+d(d+1)+\cdots+d(d+1)^{n-1}=d\frac{(d+1)^n-1}{(d+1)-1}=(d+1)^n-1$ 即可.

有了多元多项式的"Mignotte"界,我们在用模素数方法处理多元多项式时,至少多了一个用对系数界的估计的方法,来帮助判断多项式是否能还原到整系数中.当然,我们知道Mignotte界是对一元情形一个相当好的估计,本节定理对多元情形的估计虽从其推广而来,却不一定是最好的估计,而且其随多项式规模和未定元的个数增长有可能会变得很大,我们在实际算法中仅参考.

$\mathbb{Z}[x_1,\ldots,x_n]$ 上最大公因子

在我们计算最大公因子的步骤 $\mathbb{Z}[x_1,\ldots,x_n]\rightarrow\mathbb{Z}_p[x_1,\ldots,x_n]\rightarrow\mathbb{Z}_p[x_1]\rightarrow\mathbb{Z}_p[x_1,\ldots,x_n]\rightarrow\mathbb{Z}[x_1,\ldots,x_n]$ 中,第二环节和第三环节已由算法3完成,本节将要讨论首尾两个环节,如无特别说明,本节中 $\mathrm{cont}$ , $\mathrm{pp}$ 等均是就整系数而言的.

算法4(整系数多元多项式最大公因子算法) 输入: $f,g\in\mathbb{Z}[x_1,\ldots,x_n]$ ,

输出: $h=\gcd(f,g)$ .

初始化 $a=\gcd(\mathrm{cont}(f),\mathrm{cont}(g))$ , $f=\mathrm{pp}(f)$ , $g=\mathrm{pp}(g)$ , $b=\gcd(\plc{f},\plc{g})$ ,
赋初值 $h=0$ , $q=1$ , $m=\min(\deg_nf,\deg_ng)$ , $limit=2^m|b|\min(\|f\|_{\infty},\|g\|_{\infty})$ , $plist=\{\}$ ,
循环做下面4—9步,
任取比较大的素数 $p$ 直至 $p\not|b$ 且 $p\not\in plist$ ,
令 $f_p=f\bmod p$ , $g_p=g\bmod p$ , $b_p=b\bmod p$ ,调用算法3计算 $h_p=\gcd(f_p,g_p)\in\mathbb{Z}_p[x_1,\ldots,x_n]$ ,并令 $m_1=\deg_nh_p$ ,
$h_p=b_p h_p/\plc{h_p}$ ,
若 $m_1<m$ 则令 $q=p$ , $h=h_p$ , $m=m_1$ ,
若上步判断不成立且 $m_1=m$ ,则用中国剩余定理计算 $h$ 使得 $h\bmod q$ 不变且 $h\bmod p=h_p$ ,再令 $q=pq$ ,
若 $q>limit$ 则:令 $pph=\mathrm{pp}(h)$ ,若 $pph|f$ 且 $pph|g$ 则输出 $a pph$ ,否则若 $m=0$ 则输出 $a$ .

关于本算法需要做一点说明.首先在 $m$ 和 $m_1$ 的计算中都是取了对变元 $x_n$ 的次数[6],实际上对所有的变元都取次数来比较也是可行的,并且可能更准确.其次, $limit$ 的计算本身并不是多元多项式因子系数的界,我们可以用前一节的"Mignotte"界来代替,也可以就用此限制,因为当 $q>limit$ 后还有判断,这只是预判断.我们还可以加上一则预判断,就是 $h$ 在中国剩余定理计算前后是否变化,当不变时再继续后面的判断.

Zippel稀疏插值模算法

Zippel稀疏模算法是求多元多项式最大公因子一个相当有效的算法,有关这方面的文献可以参见[7]P80-82,[1],[6]P312-313等.下面我们以一个文献中的具体的例子将这个问题简要地阐述一下.

一个具体的例子

例1 求多项式 $f,g$ 的最大公因子 $h$ .其中 $\begin{align*} f=&x^5+2yzx^4+(13yz^2-21y^3z+3)x^3\\ &+(26y^2z^3-42y^4z^2+2)x^2+(39yz^2-63y^3z+4yz)x+6,\\ g=&x^6+(13yz^2-21y^3z+z+y)x^4+3x^3\\ &+(13yz^3+13y^2z^2-21y^3z^2-21y^4z)x^2\\ &+(13yz^2-21y^3z+2z+2y)x+2. \end{align*}$ ◇

解首先我们取素数 $p_1=11$ ,并作如下计算: $h_{11}(x,1,2)\equiv x^3-x+2\pmod{11},$ $h_{11}(x,3,2)\equiv x^3+x+2\pmod{11},$ $h_{11}(x,5,2)\equiv x^3+4x+2\pmod{11},$ $h_{11}(x,-4,2)\equiv x^3+5x+2\pmod{11},$ $h_{11}(x,4,2)\equiv x^3-5x+2\pmod{11},$

用普通的稠密插值算法可以求得 $h_{11}(x,y,2)=x^3+(-3y+2y^3)x+2\pmod{11}.$ 然后我们再对 $z$ 取其它赋值点来计算,此时我们假定 $h_{11}(x,y,z_i)$ 具有形式 $x^3+(ay+by^3)x+2$ ,不必对每个 $z$ 的赋值点 $z_i$ 都取5个 $y$ 的赋值点来算,只需取2个点来计算即可.例如当 $z=-5$ 时,计算下面两个式子: $h_{11}(x,-3,-5)\equiv x^3-4x+2\pmod{11},$ $h_{11}(x,2,-5)\equiv x^3+5x+2\pmod{11}.$

于是我们得到下面的方程组: $-3a-5b\equiv -4\pmod{11},$ $2a-3b\equiv 5\pmod{11},$ 解得 $a=b=-5$ ,因此得到 $h_{11}(x,y,-5)=x^3+(-5y-5y^3)x+2\pmod{11}$ .同理我们还可以得到 $h_{11}(x,y,-3)\equiv x^3+(-4y-3y^3)x+2\pmod{11},$ $h_{11}(x,y,5)\equiv x^3+(-5y+5y^3)x+2\pmod{11}.$

此时利用 $z$ 在4个赋值点计算的结果,利用稠密插值可以得到 $h_{11}(x,y,z)\equiv x^3+(2yz^2+y^3z)x+2\pmod{11}.$

其次,我们再取素数 $p_2=17$ ,此时也不必取前面那么多赋值点,利用稀疏性假设,我们可以认为 $h_{17}(x,y,z)$ 具有形式 $x^3+(cyz^2+dy^3z)x+2$ ,只需取两组 $(y,z)$ 进行插值,例如: $h_{17}(x,7,-4)\equiv x^3+8x+2\pmod{17},$ $h_{17}(x,-2,4)\equiv x^3+x+2\pmod{17},$ 于是解方程可得 $h_{17}(x,y,z)=x^3+(-4yz^2-4y^3z)x+2\pmod{17}$ .

注3 利用中国剩余定理将 $\mathbb{Z}_{11}$ 和 $\mathbb{Z}_{17}$ 中的结果合起来,即为 $h=x^3+(13yz^2-21y^3z)x+2$ ,经验算,其恰为所欲求的最大公因子.这里只用了13次一元多项式GCD问题求解,而若用通常的模算法,则需将近40次.由此可见,稀疏插值模算法确能有效地减少计算次数.

算法描述

前面我们曾经介绍过了Zippel稀疏插值算法的思想以及具体操作的方法,现在我们只需稍加改动,具体给出求值函数的形式,就可以用来计算最大公因子.本节中将要给出的算法描述可以参考[2]15.3节.

先对记号做一些说明,在下面将要介绍的插值算法中,我们保留变元 $X_n$ ,并记之为 $Z$ .对集合 $\{(u_1,v_1),\ldots,(u_n,v_n)\}$ 的插值即求多项式 $f$ 使得 $\forall i(1\le i\le n),f(u_i)=v_i$ .

稀疏插值算法1和稀疏插值算法2两者互相递归调用.

算法5(稀疏插值算法1) 输入: $l(X_1,\ldots,X_k)$ , $f(X_1,\ldots,X_k,Z)$ , $g(X_1,\ldots,X_k,Z)$ ,其中 $0\le k<n$ ,

输出:多项式 $f,g$ 的最大公因子的某个倍数,使得其关于 $Z$ 的首项系数为 $l$ .

若 $k=0$ 则输出 $l\gcd(f,g)$ ,
令 $d=\min(\deg_kf,\deg_kg)+\deg_kl$ ,并任取一赋值点 $x_k$ ,
递归调用本算法计算 $f(X_1,\ldots,X_{k-1},x_k,Z),g(X_1,\ldots,X_{k-1},x_k,Z)$ 的最大公因子 $I(X_1,\ldots,X_{k-1},Z)$ ,其中第一个参数输入 $l(X_1,\ldots,X_{k-1},x_k)$ ,
令 $T$ 为多项式 $I$ 关于 $Z$ 的系数( $\in F[X_1,\ldots,X_{k-1}]$ )中含最多单项式的系数的单项式个数,同时对 $j$ 从 $0$ 循环到 $D$ 顺次命 $\mathcal{H}_j=\{(x_k,I\text{关于}Z^j\text{的系数})\}$ ,
对 $i$ 从 $1$ 循环到 $d$ 顺次做下面6,7步,
随机取不重复的赋值点 $y_i$ ,调用算法6,输入 $\mathrm{skel}I$ , $l(X_1,\ldots,X_{k-1},y_i)$ , $F(X_1,\ldots,X_{k-1},y_i,Z)$ , $G(X_1,\ldots,X_{k-1},y_i,Z)$ , $T$ ,得到稀疏插值求得的多项式 $W$ ,
对 $j$ 从 $0$ 循环到 $D$ 顺次命 $\mathcal{H}_j=\mathcal{H}_j\cup\{(y_i,W\text{关于}Z^j\text{的系数})\}$ ,
对 $j$ 从 $0$ 循环到 $D$ 顺次利用稠密插值算法由 $\mathcal{H}_j$ 计算到到多项式 $h_j\in F[X_1,\ldots,X_k]$ ,令 $h=h_DZ^D+\cdots+h_0$ ,
若 $h$ 不能同时整除 $f$ 和 $g$ 则退出整个算法重新选取赋值点（整个算法失败）,否则输出 $h$ .

算法6(稀疏插值算法2) 输入:模板集 $S$ ,多项式 $l(X_1,\ldots,X_k)$ , $f(X_1,\ldots,X_k,Z)$ , $g(X_1,\ldots,X_k,Z)$ , $T$ ,

输出:由它们计算得到稀疏插值的结果 $h$ ,也即 $f$ , $g$ 的最大公因子的某个倍数,其关于 $Z$ 的首项系数为 $l$ .

若 $k=0$ 则输出 $l\gcd(f,g)$ ,
任取 $k$ 个赋值点 $y_1,\ldots,y_k$ ,
对 $i$ 从 $1$ 循环到 $T$ 顺次做下面4,5步,
令 $W=l(y_1^i,\ldots,y_k^i)\times\gcd(f(y_1^i,\ldots,y_k^i,Z),g(y_1^i,\ldots,y_k^i,Z))$ ,
对 $j$ 从 $0$ 循环到 $D$ 顺次命 $\mathcal{H}_j=\mathcal{H}_j\cup\{W\text{关于}Z^j\text{的系数}\}$ ,
对 $j$ 从 $0$ 循环到 $D$ ,由 $\mathcal{H}_j$ ,以及 $S$ 中含 $Z^j$ 的项和 $y_1,\ldots,y_k$ 计算出Vandermonde矩阵各元素,解线性方程组可得 $h_j$ (具体方法见注5),
输出 $h_DZ^D+\cdots+h_0$ .

注4 算法5第1步,算法6第1步,第4步等计算一元GCD时均是指求得其首一化的GCD.

注5 算法6第6步计算的具体方法是:令 $S_j$ 为输入模板集 $S$ 中最后一项指标为 $j$ ,即 $Z$ 的系数为 $j$ 的指标在前 $k$ 维上投影的集合.例如对于指标集 $\{(1,1),(1,0),(0,0)\}$ 有 $S_0=\{(1),(0)\}$ , $S_1=\{(1)\}$ .设 $S_j$ 的元素个数为 $s_j$ ,由赋值点 $y_1,\ldots,y_k$ 和指标集 $S_j$ 可以计算得到一个 $s_j$ 阶的Vandermonde矩阵(若此矩阵奇异,则返回算法第2步重新选取赋值点).再取 $\mathcal{H}_j$ 的前 $s_j$ （肯定小于 $T$ ）项,以此可得一线性方程组.其解配上相应的模板集 $S_j$ 即到多项式 $h_j$ .

如果 $f$ 和 $g$ 都是关于 $Z$ 的首一多项式,则直接调用算法5,输入参数 $1,f,g$ 即可,但是对于关于 $Z$ 的首项系数不为1的多项式,我们可如下给出它们的最大公因子.其中为了使我们每一步最后的试除法能够成功,调用5时我们用 $lf$ , $lg$ .此时,为了使稠密插值的结果其首项系数为 $l$ ,我们在算法5第2步中令 $d=\min(\deg_kf,\deg_kg)+\deg_kl$ .

算法7(稀疏插值算法) 输入:多项式 $f,g\in F[X_1,\ldots,X_n]$ ,

输出: $f,g$ 的最大公因子.

递归调用本算法计算 $f_1=\mathrm{cont}_Z(f)$ , $g_1=\mathrm{cont}_Z(g)$ , $a=\gcd(f_1,g_1)$ ,并令 $f=f/f1$ , $g=g/g1$ ,
递归调用本算法计算 $l=\gcd(\mathrm{lc}_Z(f),\mathrm{lc}_Z(g))$ ,调用算法5,输入参数 $l,lf,lg$ ,得到 $h$ ,
输出 $a\mathrm{pp}(h)$ .

求GCD的其它方法

关于多元多项式的GCD问题,前述的一般算法以及Zippel稀疏插值模算法已经能够相当有效地处理,因此这一节所说的方法仅是作一些补充说明,算法的细节不再详细说明,有兴趣可以参考列出的一些文献.

启发式算法(Heuristic GCD)

启发式算法是一种将多元多项式GCD问题转化为大整数问题的算法.此算法的优势在于它将问题化为整数计算问题,而一般说来,对于一个计算机代数系统,其整数计算的部分都是在系统的Kernel中,因而由此提高计算效率,并且其前提为我们有一个非常高效处理整数的内核.其次,对于一般的小问题(即不是显得特别巨型冗长的多项式)来说,用启发式算法还是比插值算法可能来得快.

本节描述参考[3].至于算法的具体实现,可参考[6]7.7节.

首先我们给出如下一个有用的引理:

引理1 令 $f\in\mathbb{Z}[x_1,\ldots,x_n]$ , $a\in\mathbb{Z}$ , $x_k-a|f$ ,则 $|a|\le\|f\|_{\infty}$ .

证明设 $f=(x_k-a)g$ ,因为 $f\neq 0$ ,显然 $g\neq 0$ .再设 $g=g_0+\cdots+g_dx_k^d$ ,其中 $g_i\in\mathbb{Z}[x_1,\ldots,x_{k-1},x_{k+1},\ldots,x_n]$ ,则 $f=-ag_0+(g_0-ag_1)x_k+\cdots+(g_{d-1}-ag_d)x_k^d+g_dx_k^{d+1},$ 于是我们有 $\|-ag_0\|_{\infty}\le\|f\|_{\infty}$ , $\|g_{i-1}-ag_i\|_{\infty}\le\|f\|_{\infty}$ , $\|g_d\|_{\infty}\le\|f\|_{\infty}$ .若 $|a|>\|f\|_{\infty}$ ,则可推出 $g=0$ ,矛盾. □

利用上面的引理,下面的定理给出了一种多项式GCD模算法.以下本节中的范数均指无穷范数.

定理6 设 $f,g\in\mathbb{Z}[X]$ 非零, $r>2$ 且大于 $f,g$ 在 $\mathbb{Z}[X]$ 内任何因子的范数的2倍,令 $\overline{f}=\Phi_{x_n-r}(f)$ , $\overline{g}=\Phi_{x_n-r}(g)$ , $\overline{h}=\gcd(\overline{f},\overline{g})=h_0+h_1r+\cdots+h_dr^d$ ,其中 $h_i\in\mathbb{Z}[x_1,\ldots,x_{n-1}]$ ,令 $H=h_0+h_1x_n+\cdots+h_dx_n^d$ ,则 $H=\gcd(f,g)\Leftrightarrow H|f\wedge H|g.$

证明我们证明 $\Leftarrow$ .设 $h=\gcd(f,g)$ ,则由 $H|f\wedge H|g$ 知 $H|h$ ,设 $h=HG$ ,于是 $\Phi_{x_n-r}(h)=\Phi_{x_n-r}(H)\Phi_{x_n-r}(G)=\overline{h}\Phi_{x_n-r}(G),$ 而 $h=\gcd(f,g)\Rightarrow\Phi_{x_n-r}(h)|\gcd(\overline{f},\overline{g})=\overline{h}$ ,从而 $\Phi_{x_n-r}(G)=\pm 1\Rightarrow G\pm 1=0\vee (x_n-r)|(G\pm 1)$ .若 $G\pm 1\neq 0$ ,则由引理1可知 $\|G\pm 1\|\ge r$ .另一方面 $G$ 是 $f,g$ 的因子,则 $\|G\|<r/2\Rightarrow\|G\pm 1\|\le r/2$ ,矛盾，故 $G\pm 1=0$ ,命题成立. □

因为在该定理中要估计 $f,g$ 各因子的范数的上界,因此对于一元情形,可以利用Mignotte界来给出,但是对于多元情形则无法给出.下面给出一例说明用此方法做一元多项式GCD问题.

例2 设 $f=x^4+25x^3+145x^2-171x-360$ , $g=x^5+14x^4+15x^3-x^2-14x-15$ ,由Mignotte界给出的上界 $M=\sqrt{6}\cdot 2^5\times 15=1175.76$ ,因此我们取 $r=2400$ ,则 $f(2400)=33524034789240$ , $g(2400)=80090933754206385$ ,于是 $\gcd(f(r),g(r))=5793615=r^2+14r+15,$ 经检验, $H=x^2+14x+15$ 确是 $\gcd(f,g)$ . ◇

由于上面的方法一则依赖Mignotte界,其系数绝对值仍然可能很大,二则对于多元情形并不适用,因此仍需对其进行改进.

下面的Cauchy不等式给出了 $\mathbb{C}[x]$ 上多项式根的模的一个估计.

定理7 设 $f=\sum_{0\le i\le n}f_ix^i\in\mathbb{C}[x]$ , $r_1,\ldots,r_n$ 是它的根,则对任何一个根 $r$ 均有 $r_<1+\displaystyle\frac{\max(|f_0|,|f_1|,\ldots,|f_{n-1}|)}{|f_n|}.$

证明若 $|r|\le 1$ ,则无需证明.下面假设 $|r|>1$ ,由 $\sum_{0\le i\le n}f_ir^i=0$ 可得 $\begin{align*} |f_nr^n|&=\left|\sum_{0\le i\le n-1}f_ir^i\right|\le\max_{0\le i\le n-1}|f_i|\cdot (1+|r|+\cdots+|r|^{n-1})\\ &=\max_{0\le i\le n-1}|f_i|\cdot \frac{|r|^n-1}{|r|-1}<\max_{0\le i\le n-1}|f_i|\frac{|r|^n}{|r|-1}, \end{align*}$ 因此得到估计 $|r|<1+\displaystyle\frac{\max_{0\le i\le n-1}|f_i|}{|f_n|}$ . □

引理2 设 $p\in\mathbb{Z}[x]$ 非常数, $s\ge 1$ 为给定正实数,整数 $r$ 满足 $|r|\ge\|p\|_{\infty}+s+1$ ,若 $q\in\mathbb{Z}[x]$ 非常数且 $q|p$ ,则 $\|\Phi_{x-r}(q)\|_{\infty}>s$ .

证明在 $\mathbb{C}$ 中,设 $p=c\prod_{1\le i\le d}(x-r_i)$ ,则由Cauchy不等式估计有: $|r_i|<1+\frac{\|p\|_{\infty}}{|c|},$ 而 $\|\Phi_{x-r}(p)\|_{\infty}\ge|c|\left|\prod_{1\le i\le d}(|r|-|r_i|)\right|$ ,由于 $|r|\ge\|p\|_{\infty}+s+1$ ,则 $\big||r|-|r_i|\big|>1+\|p\|_{\infty}+s-(1+\frac{\|p\|_{\infty}}{|c|})\ge 0$ .故 $\|\Phi_{x-r}(p)\|_{\infty}>|c|s^d\ge s$ .同理,也可得到 $\|\Phi_{x-r}(q)\|_{\infty}\ge s$ . □

引理2可以用来引入下面的一元多项式模GCD算法.

定理8 设 $f,g\in\mathbb{Z}[x]$ 非零, $r\in\mathbb{Z}$ 满足 $r>1+\min(\|f\|_{\infty},\|g\|_{\infty})$ ,设 $\overline{h}=\gcd(\overline{f},\overline{g})=h_0+\cdots+h_dr^d$ ,令 $H=h_0+\cdots+h_dx^d$ ,则 $H=\gcd(f,g)\Leftrightarrow H|f\wedge H|g$ .

证明我们仍然只需要证明 $\Leftarrow$ .同样的,我们可以设 $h=HG$ ,容易得到 $\Phi_{x-r}(G)=\pm 1$ .设 $p$ 为 $f,g$ 中范数较小者,则 $r>1+\|p\|_{\infty}\Rightarrow r\ge 1+1+\|p\|_{\infty}$ .若 $H$ 非常数,则由引理2知 $\|\Phi_{x-r}(H)\|_{\infty}>1$ ,此与 $\Phi_{x-r}(H)=\pm 1$ 矛盾.故 $H$ 为常数, $H=\pm 1$ . □

当我们用此算法对例2再次讨论时,发现赋值点的选取要小得多.由于 $\|g\|_{\infty}=15$ ,因此我们取 $r=20$ ,比前面取的 $2400$ 要小得多.此时有 $f(20)=414220$ , $g(20)=5559305$ , $\gcd(f(20),g(20))=695=20^2+14*2=+15$ ,因此 $H=x^2+14x+15$ .

上述算法并不保证一定能够得到最大公因子.例如取 $f=(x-1)x$ , $g=(x+1)(x+2)$ ,则对二者取赋值点后求得的GCD必为偶数,但事实上 $\gcd(f,g)=1$ ,因此进一步我们必须考虑多项式的本原部分.下面的定理可以解决这一问题.

定理9 设 $f,g\in\mathbb{Z}[x]$ 是本原多项式, $r\in\mathbb{Z}$ 满足 $r>1+2\min(\|f\|_{\infty},\|g\|_{\infty})$ ,令 $\overline{h}=\gcd(\overline{f},\overline{g})$ ,定义多项式 $H$ 同前,则同样地有 $\mathrm{pp}(H)=\gcd(f,g)\Leftrightarrow\mathrm{pp}(H)|f\wedge\mathrm{pp}(H)|g$ .

证明我们仍然只证 $\Leftarrow$ .令 $a=\gcd(f,g)$ ,则有 $G$ 使得 $a=\mathrm{pp}(H)G$ .且 $\Phi_{x-r}(\mathrm{pp}(H))\Phi_{x-r}(G)=\overline{a}|\overline{h},$ 记 $c=\mathrm{cont}(H)$ ,则 $\mathrm{pp}(H)=H/c$ ,故 $\Phi_{x-r}(\mathrm{pp}(H))=\Phi_{x-r}(H/c)=\Phi_{x-r}(H)/c=\overline{h}/c,$ 于是 $\Phi_{x-r}(G)|c$ ,即 $\Phi_{x-r}(G)\in\mathbb{Z}$ 且 $\|\Phi_{x-r}(G)\|_{\infty}\le c\le r/2$ .

再设 $p$ 为 $f,g$ 中无穷范数较小者,于是 $r/2\ge 1+\|p\|_{\infty}$ .故 $\Phi_{x-r}(G)$ 非常数,于是由引理得 $\|\Phi_{x-r}(G)\|_{\infty}>s=r/2$ .因此, $G$ 为常数,由 $a=\mathrm{pp}(H)G$ 知 $G=\pm 1$ . □

对于多元情况,我们先引入一些记号.本节后面的命题不再重述.设 $\alpha_1,\ldots,\alpha_{n-1}$ 是随机选取的一些整数, $p\in\mathbb{Z}[x_1,\ldots,x_n]$ 是非常数多项式,且使得 $\deg_{x_n}(p)=\deg_{x_n}(\Phi_I(p))$ ,其中理想 $I=\idea{x_1-\alpha_1,\ldots,x_{n-1}-\alpha_{n-1}}$ .记 $\|p\|_I=\max(\|p\|_{\infty},\|\Phi_I(p)\|_{\infty})$ .

引理3 设有整数 $r$ 满足 $|r|\ge\|p\|_I+s+1$ , $s$ 是某一正实数且 $s\ge 1$ ,若 $q\in\mathbb{Z}[x_1,\ldots,x_n]$ 非常数且是 $p$ 的因子, $\Phi_{x_n-r}(q)\in\mathbb{Z}$ ,则 $\|\Phi_{x_n-r}(q)\|_{\infty}>s$ .

证明 $\Phi_I(q)|\Phi_I(p)$ .若 $\Phi_I(q)$ 是常数,则设 $p=qt$ ,有 $\deg_{x_n}\Phi_I(p)=\deg_{x_n}\Phi_I(t)\Phi_I(q)=\deg_{x_n}\Phi_I(t)\le\deg_{x_n}(t)\le\deg_{x_n}(q)+\deg_{x_n}(t)=\deg_{x_n}(p)$ ,由此 $\deg_{x_n}(q)$ ,则 $\Phi_{x_n-r}(q)=q\in\mathbb{Z}$ ,矛盾.故 $\Phi_I(q)$ 非常数.此时由引理2知 $\|\Phi_{x_n-r}(\Phi_I(q))\|_{\infty}>s$ .

因为 $\Phi_{x_n-r}(q)\in\mathbb{Z}$ ,则有 $\Phi_{x_n-r}(\Phi_I(q))=\Phi_I(\Phi_{x_n-r}(q))=\Phi_{x_n-r}(q),$ 因此 $\|\Phi_{x_n-r}(q)\|_{\infty}>s$ . □

下面两个定理的证明和前面一元情形两个定理的证明类似,这两个定理可以用来求多元多项式的最大公因子.

定理10 $f,g\in\mathbb{Z}[x_1,\ldots,x_n]$ 是非零多项式, $\alpha_1,\ldots,\alpha_{n-1}$ 的选取同前所述,使得 $\deg_{x_n}(f)=\deg_{x_n}(\Phi_I(f)),\quad \deg_{x_n}g=\deg_{x_n}(\Phi_I(g)),$ $r\in\mathbb{Z}$ 满足 $r>1+\min(\|f\|_I,\|g\|_I)$ ,令 $\overline{h}=\gcd(\overline{f},\overline{g})$ , $H$ 的引入同前,则 $H=\gcd(f,g)\Leftrightarrow H|f\wedge H|g.$

定理11 $f,g\in\mathbb{Z}[x_1,\ldots,x_n]$ 是非零本原多项式, $\alpha_1,\ldots,\alpha_{n-1}$ 的选取同定理10所述, $r\in\mathbb{Z}$ 满足 $r>1+2\min(\|f\|_I,\|g\|_I)$ . $\overline{h},H$ 的引入同前,则 $\mathrm{pp}(H)=\gcd(f,g)\Leftrightarrow\mathrm{pp}(H)|f\wedge\mathrm{pp}(H)|g.$

EZ-GCD

EZ-GCD算法即Extended Zassenhaus算法([8]),它利用Hensel提升来计算多元多项式的GCD.此算法的具体实现可以参考[6]7.6节.

多元多项式因子分解的Kronecker算法

对于给定的多项式 $f\in\mathbb{Z}[x_1,x_2,\ldots,x_n]$ ,我们取主变元为 $x_1$ ,并记之为 $x$ .设 $d=\max_{1\le i\le n}\deg_{x_i}f,$ 如果我们取多项式 $\tilde{f}(x)=f(x,x^{d+1},x^{(d+1)^2},\ldots,x^{(d+1)^{n-1}}),$ 易知如果有不可约因子分解 $f=f_1f_2\cdots f_r$ ,这里我们假设 $f$ 关于 $x$ 是无平方因子本原多项式,那么有 $\tilde{f}=\tilde{f_1}\cdots\tilde{f_r}.$ 然而当我们对 $\tilde{f}$ 进行因子分解时,得到的不一定是上式,因为 $\tilde{f_i}$ 不一定均是不可约的.设 $\tilde{f}$ 的不可约因子分解为 $\tilde{f}=g_1g_2\cdots g_t,$ 则由因子组合算法可以还原在 $\mathbb{Z}[x_1,\ldots,x_n]$ 中的分解.

由此可见,Kronecker算法的思想本身是简单的,下面给出具体的算法.

算法8(Kronecker因子分解算法) 输入:关于 $x=x_1$ 本原且无平方因子的多项式 $f\in\mathbb{Z}[x_1,\ldots,x_n]$ ,

输出:其各不可约因子 $\{f_1,\ldots,f_r\}$ .

令 $d$ 为各变元次数的上界,求得多项式 $\tilde{f}(x)=f(x,x^{d+1},x^{(d+1)^2},\ldots,x^{(d+1)^{n-1}})$ 的因子分解 $\tilde{f}=g_1g_2\cdots g_t.$
令 $T=\{1,2,\ldots,t\}$ , $s=1$ , $result=\{\}$ , $h=f$ ,
若 $2s\le\# T$ ,则循环做下面4-6步,否则转第7步,
枚举 $T$ 的所有 $s$ 元子集 $S$ ,并做下面第5步,
由多项式 $\prod_{i\in S}g_i\in\mathbb{Z}[x]$ 我们可以还原得到多项式 $g\in\mathbb{Z}[x_1,\ldots,x_n]$ (见注6),若 $g|h$ 则令 $result=result\cup\{g\}$ , $h=h/g$ , $T=T\setminus S$ ,并转第3步,
$s=s+1$ ,转第3步,
输出 $result\cup\{h\}$ .

注6 由某一多项式 $\tilde{f}(x)$ 还原得到多项式 $f(x_1,\ldots,x_n)$ 的方法是:对 $\tilde{f}$ 中的任何一个单项式 $ax^b$ ,由连续对 $d+1$ 的除法可以得到 $n$ 元序对 $(b_1,\ldots,b_n)$ ,其中 $0\le b_i\le d$ 使得 $b=b_1+b_2(d+1)+b_3(d+1)^2+\cdots+b_n(d+1)^{n-1},$ 将此单项用 $ax_1^{b_1}x_2^{b_2}\cdots x_n^{b_n}$ 代替即可还原多项式.

利用Hensel提升的因子分解算法

概述

正如我们前面所说的以及多元GCD求解的方法,利用赋值同态的方法,我们也可以将多元因子分解问题转化为一元问题.我们很容易会产生如下一般性的想法,这里假设我们考虑 $f\in \mathbb{Z}[x_1,x_2,\ldots,x_n]$ 的分解,并将变元 $x_1$ 视为主变元 $x$ .

$\begin{itemize} \item 将$f$化为关于主变元$x$本原以及无平方因子的多项式.这一点是很容易做到的,由多元GCD算法可以求出$f$关于$x$各系数多项式的最大公因子,而由$\gcd(f,\partial f/\partial x)$可以将其化为无平方因子多项式的分解问题. \item 利用赋值同态$I=\idea{x_2-a_2,\ldots,x_n-a_n}$将$f$化为$\tilde{f}=f\bmod I$,即$\tilde{f}(x)=f(x,a_2,\ldots,a_n)$,使得$\tilde{f}$无平方因子并且$\deg_x\tilde{f}=\deg_xf$. \item 处理一元分解问题,得到不可约分解$\tilde{f}=g_1g_2\cdots g_r$,这里可以将$\mathrm{cont}\tilde{f}$任意分配到不可约因子$g_i$上. \item 利用类似于Hensel提升的方法,将模$I$下的分解提升到足够高的模$I^k$下的分解,得到多元因子,最后采取诸如因子组合的方法将其还原到整系数多项式中的因子. \end{itemize}$

关于赋值点和主变元的选取,这里做一些补充.赋值点应优先选择 $\pm 1,0$ ,以保证 $\tilde{f}$ 无平方因子且关于 $x$ 次数不变.选择 $\pm 1,0$ 的好处是使得系数较小,由后面的算法我们看出对于非零赋值点,我们将对变量做一平移,以使提升算法的模运算也便于进行.

其次要考虑使得 $\tilde{f}$ 的不可约因子数尽可能少, $\tilde{f}$ 关于 $x$ 的首项系数尽可能小,如果是 $1$ 更好,那么所有的不可约因子的首项系数都将为 $1$ .

Extended Zassenhaus 算法

Zassenhaus算法基本上就是我们曾经提及的Hensel提升算法,具体也可参考[2]以及[9]第7节.

现在我们要介绍的扩展Zassenhaus算法,是指如下类似的问题.设 $\tilde{f}\in Z[x]$ 有分解 $\tilde{f}=gh$ ,其中 $\gcd(g,h)=1$ ,亦即 $f\equiv gh\pmod{I},$ 现在需找到 $g_k$ , $h_k$ 使得 $f\equiv g_kh_k\pmod{I^k},\quad f_k\equiv f\pmod{I},\quad g_k\equiv g\pmod{I}.$ 对于 $i=2,3,\ldots,n$ ,记 $y_i=x_i-a_i$ ,则此时 $I=\idea{y_2,\ldots,y_n}$ ,并记 $f^*(x_1,y_2,\ldots,y_n)=f(x_1,y_2+a_2,\ldots,y_n+a_n)$ .基于这些符号上的说明,我们可以提出一种归纳的算法如下:

算法9(EZ算法) 输入:多项式 $f^*$ (即 $f$ ), $g_{k-1}$ , $h_{k-1}$ ,并且满足前述相应条件,

输出:提升后的 $g_k$ , $h_k$ .

计算 $r_{k-1}$ , $e_{k-1}$ 使得 $r_{k-1}=g_{k-1}h_{k-1}-f^*,\quad e_{k-1}\equiv r_{k-1}\pmod{I^k},$
利用扩展Euclid算法计算唯一的 $\alpha_i(x)$ , $\beta_i(x)$ 使得 $\alpha_i(x)g(x)+\beta_i(x)h(x)=x^i,$ 并且 $\deg\alpha_i<\deg h$ .(命 $\alpha_i=\alpha_0x^i\bmod h$ ,即可.此时若 $i<\deg g+\deg h$ 则有 $\deg\beta_i<\deg g$ .并且此步可略去,因为可以在整个算法开始提升的第一步计算足够多的 $\alpha_i$ , $\beta_i$ ,保存起来并供后面各步提升时使用.)
将 $e_{k-1}$ 中的 $x^i$ 用 $\alpha_i(x)$ 代替得到多项式 $A(e_{k-1})$ ,同理将 $e_{k-1}$ 中的 $x^i$ 用 $\beta_i(x)$ 代替得到多项式 $B(e_{k-1})$ .
令 $g_k=g_{k-1}-B(e_{k-1})$ , $h_k=h_{k-1}-A(e_{k-1})$ 并输出 $g_k$ , $h_k$ .

证明(算法有效性) 首先由 $A(e_{k-1})$ 和 $B(e_{k-1})$ 的定义可知有 $A(e_{k-1})g(x)+B(e_{k-1})h(x)=e_{k-1}$ 成立.由 $e_{k-1}\bmod I^{k-1}=r_{k-1}\bmod I^{k-1}=0$ 知 $e_{k-1}$ 是 $y_2,y_3,\ldots,y_n$ 的 $k-1$ 次齐次多项式.因而 $A(e_{k-1})$ 和 $B(e_{k-1})$ 也是它们的 $k-1$ 次齐次式.由 $g_{k-1}=g-B(e_1)-B(e_2)-\cdots-B(e_{k-2})$ 可知 $A(e_{k-1})g_{k-1}(x)\equiv A(e_{k-1})g(x)\pmod{I^k},$ 同理有 $B(e_{k-1})h_{k-1}(x)\equiv B(e_{k-1})h(x)\pmod{I^k},$ 因而 $r_{k-1}\equiv e_{k-1}\equiv A(e_{k-1})g(x)+B(e_{k-1})h(x)\pmod{I^k}.$ 并且我们有 $A(e_{k-1})B(e_{k-1})\equiv 0\pmod{I^k}$ .

命 $r_k=g_kh_k-f^*$ ,则 $r_k\equiv r_{k-1}-A(e_{k-1})g_{k-1}-B(e_{k-1})h_{k-1}+A(e_{k-1})B(e_{k-1})\equiv 0\pmod{I^k}.$ 亦即 $f^*\equiv g_kh_k\pmod{I^k}$ . □

在向 $I^k$ 提升的过程中,倘若有某一步 $r_l=0$ ,则之后无须提升.

本算法是在整系数多项式下实现的.如果在 $\mathbb{Z}/m\mathbb{Z}$ 中,其中 $m$ 为一素数幂 $p^l$ ,此算法也是可行的.我们有下面几个命题来保证.

引理4 $a\in\mathbb{Z}_{p^l}$ 可逆当且仅当 $p\not|a$ .

证明充分性:若 $p\not|a$ ,则 $\gcd(a,p^l)=1$ ,由Bezout等式可知 $a$ 在 $\mathbb{Z}_{p^l}$ 中可逆.

必要性:设 $kp(0\le k<p^{l-1})$ 在 $\mathbb{Z}_{p^l}$ 中可逆,则 $p|\gcd(kp,p^l)|1$ ,矛盾. □

下面给出的定理揭示了 $\mathbb{Z}_{p^l}[x_1]$ 类似于Euclid整环的性质.

定理12 设 $g,h\in\mathbb{Z}_{p^l}[x_1]$ 满足 $p\not|\plc{g,x_1}$ , $p\not|\plc{h,x_1}$ ,且 $\gcd(\Phi_p(g),\Phi_p(h))=1$ ,则 $\forall f\in\mathbb{Z}_{p^l}[x_1]$ 都存在唯一的多项式 $s,t\in\mathbb{Z}_{p^l}[x_1]$ 使得 $sg+th\equiv f\pmod{p^l}$ ,且 $\deg_{x_1}s<\deg_{x_1}h$ .若 $\deg_{x_1}f<\deg_{x_1}g+\deg_{x_1}h$ ,还有 $\deg_{x_1}t<\deg_{x_1}g$ .

证明存在性:首先由 $\mathbb{Z}_p[x_1]$ 是Euclid整环我们知道 $\exists s^{(1)},t^{(1)}\in\mathbb{Z}_p[x_1]$ 使得 $s^{(1)}g+t^{(1)}h\equiv 1\pmod{p}.$

设 $s^{(k)},t^{(k)}$ 已求得并使 $s^{(k)}g+t^{(k)}h\equiv 1\pmod{p^k}$ ,则定义迭代算法如下:

设 $s_k,t_k\in\mathbb{Z}_p[x_1]$ 是 $s_kg+t_kh\equiv\frac{1-s^{(k)}g-t^{(k)}h}{p^k}\pmod{p}$ 的解,再令 $s^{(k+1)}=s^{(k)}+s_kp^k,\quad t^{(k+1)}=t^{(k)}+t_kp^k.$

显然 $s^{(k+1)}g+t^{(k+1)}h=s^{(k)}g+t^{(k)}h+p^k(s_kg+t_kh)\equiv 1\pmod{p^{k+1}}.$

因此 $fs^{(l)}g+ft^{(l)}h\equiv f\pmod{p^l}$ .由题设条件和引理4知 $\plc{h,x_1}$ 是 $\mathbb{Z}_{p^l}$ 中可逆元,我们可以作除法: $fs^{(l)}=uh+s\pmod{p^l},$ 其中 $\deg_{x_1}s<\deg_{x_1}h$ .再定义 $t=ft^{(l)}+ug$ ,即知 $s,t$ 满足定理.

唯一性:由 $s^{(l)},t^{(l)}$ 的存在知道在 $\mathbb{Z}_{p^l}$ 中 $g$ 和 $h$ 也互素,设满足定理的还有 $s',t'$ ,则易得 $(s-s')g=(t'-t)h$ ,由 $\deg_{x_1}(s-s')<\deg_{x_1}$ 且 $h|(s-s')$ 知 $s-s'=t-t'=0$ . □

我们可以将多元Hensel提升问题写为如下定理.

定理13 设 $f\in\mathbb{Z}[x_1,\ldots,x_n]$ ,理想 $I=\idea{x_2-a_2,\ldots,x_n-a_n}$ ,且素数 $p$ 满足 $p\not|\plc{f(x_1,a_2,\ldots,a_n),x_1}=\plc{\Phi_I(f),x_1}.$ 现有多项式 $g^{(1)},h^{(1)}$ 满足 $f\equiv g^{(1)}h^{(1)}\pmod{I,p^l},$ 且 $\Phi_p(g^{(1)})$ 和 $\Phi_p(h^{(1)})$ 在 $\mathbb{Z}_p[x_1]$ 中互素.

则 $\forall k\in\mathbb{Z}(k\ge 1)$ , $\exists g^{(k)},h^{(k)}\in\mathbb{Z}_{p^l}[x_1,\ldots,x_n]/I^k$ ,使得 $f\equiv g^{(k)}h^{(k)}\pmod{I^k,p^l},$ 且 $g^{(k)}\equiv g^{(1)}\pmod{I,p^l},\quad h^{(k)}\equiv h^{(1)}\pmod{I,p^l}.$

证明命题当 $k=1$ 时显然成立,现在我们归纳假设命题对于 $k$ 是成立的,此时设 $e^{(k)}=f-g^{(k)}h^{(k)}\in I^k$ ,可以将其在 $I^k$ 的基下进行表示,设表示为 $e^{(k)}=\sum_{2\le i_1\le n}\cdots\sum_{i_{k-1}\le i_k\le n}C_{i_1\cdots i_k}\prod_{j=1}^k(x_{i_j}-a_{i_j}),$ 其中 $C_{i_1\cdots i_k}\in\mathbb{Z}_{p^l}[x_1]$ .于是由定理12我们知道 $\forall i=i_1\cdots i_k$ ,都可以找到满足上面定理的唯一的 $s_i,t_i\in\mathbb{Z}_{p^l}[x_1]$ 满足 $s_ig^{(1)}+t_ih^{(1)}\equiv C_i\pmod{p^l}.$

令 $g^{(k+1)}=g^{(k)}+\sum_{2\le i_1\le n}\cdots\sum_{i_{k-1}\le i_k\le n}t_i\prod_{j=1}^k(x_{i_j}-a_{i_j}),$ $h^{(k+1)}=h^{(k)}+\sum_{2\le i_1\le n}\cdots\sum_{i_{k-1}\le i_k\le n}s_i\prod_{j=1}^k(x_{i_j}-a_{i_j}),$ 由算法9的证明易知 $g^{(k+1)},h^{(k+1)}$ 满足条件. □

[7]中给出了如下一个多元Hensel提升求因子分解的简单的例子.

例3 设 $f=x^2-3xz^2+2x-yx+3yz^2-2y+zx-3z^3+2z$ .取 $p=7$ , $l=1$ , $I=\idea{y,z}$ ,有 $f\equiv x(x+2)\pmod{I,p^l},$ 设 $g^{(1)}=x$ , $h^{(1)}=x+2$ ,因为在 $\mathbb{Z}_{p^l}[x_1]$ 中有 $3g^{(1)}-3h^{(1)}\equiv 1\pmod{p^l}$ ,即 $s=3,t=-3$ .

计算知 $e^{(1)}=f-g^{(1)}h^{(1)}=(-x-2)y+(x+2)z\pmod{I^2,p^l}$ ,即 $C_{10}^{(1)}=-x-2$ , $C_{01}^{(1)}=x+2$ ,易得 $C_{10}^{(1)}\equiv C_{10}^{(1)}sg^{(1)}+C_{10}^{(1)}th^{(1)}\equiv 0g^{(1)}+(-1)h^{(1)}\pmod{p^l}$ 即 $s_{10}^{(1)}=0$ , $t_{10}^{(1)}=-1$ ,同理可得 $s_{01}^{(1)}=0$ , $t_{01}^{(1)}=1$ ,经计算可得 $g^{(2)}=x-y+z$ , $h^{(2)}=x+2$ .

经过再一次提升,可得 $g^{(3)}=x-y+z$ , $h^{(3)}=x+2-3z^2$ ,此时已有 $f=g^{(3)}h^{(3)}$ . ◇

因子还原

不妨设我们已得到因子分解 $f\equiv g_1g_2\cdots g_t\pmod{m,I^k},$ 其中模去一个整数 $m$ ,是指在计算提升之前可以选择模掉一足够大的 $m$ 进行运算,以减小系数膨胀.当然也可以直接在整数环中计算,此种情况我们统一用 $m=0$ 表示.

因子还原的过程也是一个因子组合算法,这一算法已经多次出现,因此这里直接给出还原的算法.

算法10(EZ算法结果的因子还原) 输入:已知分解 $f\equiv g_1\cdots g_t\pmod{m,I^k}$ ,

输出: $f$ 的不可约因子集合 $\{f_1,f_2,\ldots,f_r\}$ .

令 $T=\{1,2,\ldots,t\}$ , $s=1$ , $result=\{\}$ , $h=f$ , $h^*=\plc{h}h$ ,
若 $2s\le\#T$ ,则循环做下面3-5步,否则转第6步,
枚举 $T$ 的所有 $s$ 元子集 $S$ ,并做下面第4步,
$g=\prod_{i\in S}g_i\bmod (m,I^k)$ , $g^*=\plc{h}\plc{g}^{-1}g\bmod (m,I^k)$ ,若 $g^*|h^*$ 则令 $result=result\cup\{\mathrm{pp}(g^*)\}$ , $h=h/\mathrm{pp}(g^*)$ , $h^*=\plc{h}h$ , $T=T\setminus S$ ,并转第2步,
$s=s+1$ ,转第2步,
输出 $result\cup\{h\}$ .

预先确定因子的首项系数

前几节组合起来可以给出一个完整的因子分解算法,然而其仍有一些效率上的不足.其中,因子还原时实际上涉及到Taylor展开,这正是我们前面尽量选择较小的赋值点的原因,以减轻此处计算负担.另外,首项系数的不确定使得EZ算法提升时的中间多项式将会比较复杂,最后因子还原时因而也有对其的处理.

我们先给出一个因子分解的例子.

例4 考虑多项式 $f=(y+2)^2x^2-1$ 的因子分解. ◇

解这里取 $x$ 为主变元, $y$ 取赋值点 $0$ ,则有分解 $f\equiv (2x+1)(2x-1)\pmod{y}$ .设 $g=2x+1$ , $h=2x-1$ ,提升算法给出 $g_2=\frac{1}{2}(2+y+4x(1+y)),\quad h_2=\frac{1}{2}(-2+4x+y),$ $g_3=\frac{1}{2}(2+y)(1+x(2+y)),\quad h_3=\frac{1}{4}(-4+8x+2y-y^2).$ 这时再经过一次首项系数的处理并取本原部分方可得到真正的因子 $(1+x(2+y))$ 和 $(-1+x(2+y))$ .

对于上例而言,若我们能预先确定各因子的首项系数,在提升之前就将 $g$ , $h$ 中的首项系数代替为实际首项系数,即 $g=(2+y)x+1$ , $h=(2+y)x-1$ ,则提升算法大大简化(对于此特例恰好无需任何提升).

设多项式 $f$ 关于主变元的首项系数为 $J=\plc{f,x}\in\mathbb{Z}[x_2,\ldots,x_n]$ ,设其不可约因子分解为 $J=\Omega\prod_{1\le i\le k}J_i^{e_i},$ 其中 $\Omega=\mathrm{cont}_{\mathbb{Z}}(J)$ , $J_i(1\le i\le k)$ 为其各不可约因子. 此时对于 $\tilde{f}=f\bmod I$ 有分解 $\tilde{f}=\delta g_1g_2\cdots g_r,$ 其中 $\delta=\mathrm{cont}\tilde{f}$ .我们的目的即是要确定各 $g_i$ 的首项系数,并将 $\delta$ 合理地分配到各个 $g_i$ 上.我们在选择赋值点时,需由满足如下几个条件:

定理14 选取赋值点 $a_2,a_3,\ldots,a_n$ , $\tilde{f}=f(x,a_2,\ldots,a_n)$ 使得:

$\deg\tilde{f}=\deg f$ ,即 $J(a_2,\ldots,a_n)\neq 0$ ,
$\tilde{f}$ 无平方因子,即 $\mathrm{res}_x(f,\partial f/\partial x)(a_2,\ldots,a_n)\neq 0$ ,
对任意 $J_i$ , $\tilde{J_i}=J_i(a_2,\ldots,a_n)$ 至少一有个素因子 $p_i$ ,其不能整除 $J_j(\forall j<i)$ 和 $\Omega$ , $\delta$ .