基于Dixon结式和逐次差分代换的多项式秩函数探测方法

袁月 李轶

摘 要：秩函數探测是循环程序终止性分析的重要方法，目前，已有很多研究者致力于为线性循环程序探测对应的线性秩函数，然而，针对具有多项式循环条件和多项式赋值的多项式型的循环，现有的秩函数探测方法还有所不足，解决方案大多是不完备的、或者具有较高的时间复杂度。针对现有工作对于多项式秩函数探测方法不足的问题，基于扩展Dixon结式（KSY方法）和逐次差分代换（SDS）方法，提出一种为多项式循环程序探测多项式型秩函数的方法。首先，将待探测的秩函数模板看作带参数系数的多项式，将秩函数的探测转换为寻找满足条件的参数系数的问题;然后，进一步将问题转换为判定相应的方程组是否有解的问题，至此，利用KSY方法中的扩展的Dixon结式，将问题更进一步简化为带参系数多项式（即结式）严格为正的判定问题;最后，利用SDS方法，找到一个充分条件，使得得到的结式严格为正，此时，可以获取满足条件的参数系数的取值，从而找到一个满足条件的秩函数，通过实验验证该秩函数探测方法的有效性。实验结果表明，利用该方法，可以有效地为多项式循环程序找到多项式秩函数，包括深度为d的多阶段多项式秩函数，与已有方法相比，该方法能够更高效地找到多项式秩函数，对于基于柱形代数分解（CAD）方法的探测方法因时间复杂度问题无法而应对的一些循环，利用所提方法能够在几秒内为这些循环找到秩函数。

关键词：循环程序终止性;多项式循环程序;多项式秩函数;多阶段秩函数;Dixon结式;逐次差分代换

Abstract： Ranking function detection is one of the most important methods to analyze the termination of loop program. Some tools have been developed to detect linear ranking functions corresponding to linear loop programs. However， for polynomial loops with polynomial loop conditions and polynomial assignments， existing methods for detecting their ranking functions are mostly incomplete or with high time complexity. To deal with these weaknesses of existing work， a method was proposed for detecting polynomial ranking functions for polynomial loop programs， which was based on extended Dixon resultants （the KSY （Kapur-Saxena-Yang） method） and Successive Difference Substitution （SDS） method. Firstly， the ranking functions to be detected were seen as polynomials with parametric coefficients. Then the detection of ranking functions was transformed to the problem of finding parametric coefficients satisfying the conditions. Secondly， this problem was further transformed to the problem of determining whether the corresponding equations have solutions or not. Based on extended Dixon resultants in KSY method， the problem was reduced to the decision problem whether the polynomials with symbolic coefficients （resultants） were strictly positive or not. Thirdly， a sufficient condition making the resultants obtained strictly positive were  found by SDS method. In this way， the coefficients satisfying the conditions were able to be obtained and thus a ranking function satisfying the conditions was found. The effectiveness of the method was proved by experiments. The experimental results show that polynomial ranking functions including d-depth multi-stage polynomial ranking functions are able to be detected for polynomial loop programs. This method is more efficient to find polynomial ranking functions compared with the existing methods. For loops whose ranking functions cannot be detected by the method based on Cylindrical Algebraic Decomposition （CAD） due to high time complexity， their ranking functions are able be found within a few seconds with the proposed method.

Key words： termination of loop program; polynomial loop program; polynomial ranking function; multi-stage ranking function; Dixon resultant; Successive Difference Substitution （SDS）

0 引言

循环程序的终止性分析是程序验证中一个重要且具有挑战性的研究领域。对于大多数软件来说，保证循环程序是可终止的是十分重要的。目前，已有很多出色的研究工作[1-4]基于各种计算模型来自动地进行程序的终止性分析。秩函数法是证明循环程序终止性证明的主要方法，秩函数的存在是循环程序终止的一个充分条件。秩函数的探测在程序的正确性证明方面起着核心的作用，并且也被认为是证明中最具挑战性的一部分之一[5]。

秩函数将该循环的每个状态映射为一个有界递减集合中的元素，使得每当循环完成一次迭代时，对应的集合中的值减小。若能为一个给定的循环找到一个这种对应的秩函数，则意味着该循环是可终止的。近来，秩函数的探测已成为一个热门的研究方向。已有很多研究者致力于为线性循环程序探测线性秩函数，例如Colon等在2002年前后的一些早期工作[6-7]。在2004年，Podelski等[8]首次提出了一种完备有效地解决方案，基于Farkas引理[9]，为有理数域上的单分支线性循环探测线性秩函数。Ben-Amram等[10]最终在2014年，解决了针对整数域上的单分支线性循环程序的完备的线性秩函数探测方案的时间复杂度的问题。至此，从某种意义上说，如果对于一个单分支线性循环程序，存在一个对应的线性秩函数，则一定能找到该秩函数。然而，对于一个（线性）循环程序，可能存在对应的秩函数，却不存在线性秩函数，并且，一般的循环程序终止性是不可判定的[11]，即使对于线性循环程序也是如此[12]。一些研究者[11，13]探测了终止性可判定的一些线性循环程序子类。Yang等[14]旨在寻找一些线性循环程序子类可终止的充分条件，使用计算机代数工具DISCOVERER来进行辅助计算。此外，为了应对更多的循环程序，最近一些研究者[10，15-20]构造多个线性秩函数的组合（作为多阶段秩函数）来应对更加复杂的线性循环程序。

然而，具有线性秩函数的线性循环程序只是循环程序中的一小部分，为了捕捉更多的循环行为，一个研究方向是为具有多项式循环条件和多项式赋值语句的多项式型循环程序寻找多项式型的秩函数。Cousot[21]提出了一种针对半代数循环程序的多项式秩函数探测方法，他的方法基于的是拉格朗日松弛定理和半正定规划（Semi-Definite Programming， SDP），但由于采用浮点计算，其得到秩函数并不一定是真正的秩函数，同时，拉格朗日乘子的引入也导致了该方法是不完备的;而Shen等[22]则在此方法的基础上，利用有理恢复技术试图找到更加精确的候选秩函数。Chen等[5]则是将探测多项式秩函数的问题转化为了判断半代数系统是否有解的问题，并基于Collins早年提出的量词消去技术——柱形代数分解（Cylindrical Algebraic Decomposition， CAD），提出了一种完备的多项式秩函数探测方案，但由于CAD算法在糟糕情形下的复杂度是双指数的，故限制了他们的方法只适用于较小规模的系统。

为了解决时间复杂度以及精确性的问题，本文提出一种为多项式循环程序探测多项式秩函数的方法，利用杨路等的研究成果，即扩展的Dixon结式（KSY（Kapur-Saxena-Yang）方法[23]）和逐次差分代换（Successive Difference Substitution， SDS）方法[24]，来进行秩函数的探测。具体而言，本文将循环程序的终止性证明的问题转换为不等式的证明问题，首先，将带探测的秩函数模板看作带参数系数的多项式，将秩函数的探测转换为寻找满足条件的参数系数的问题，此时，秩函数应当满足的条件被转化为一系列带参系数多项式不等式之间应当满足一定的蕴含关系;然后，将蕴含关系的判定问题转换为判定相应的方程组是否有解的问题，至此，利用KSY方法中的扩展的Dixon结式，将问题进一步简化为带参系数多项式（即结式）严格为正的判定问题;最后，利用逐次差分代换（SDS）方法，找到一个充分条件，使得得到的结式严格为正。通过这种方式，可以获取满足条件的参数系数的取值，从而找到一个满足条件的秩函数。

实验结果表明，利用本文提出的方法可以有效地为多项式循环程序找到多项式秩函数，与现有基于CAD技术的完备的探测方案相比更加高效。此外，利用该方法，还能为多项式循环程序探测深度为d的多阶段秩函数，更具体而言，即嵌套秩函数（2017年由Ben-Amram等[16]提出的多阶段秩函数的一个更加普适性的定义）。与基于SDP方法的多项式秩函数探测方案不同的是，该方法可以找到精确的秩函数，不需要对找到的秩函数进行再一次的验证。

1 预备知识

本章将介绍与KSY方法和SDS方法相关的背景知识。

1.1 Dixon结式和KSY方法

给定由k+1个含k个变元的多项式方程组S：

该多项式δ则被称为Dixon多项式。将δ看作关于α的多项式，令s′为系数个数、关于x的多项式oi（x）为每个系数，进一步得到以下方程组：

显然，无论α取值如何，S的任意一个解xi都是δ的一个零点，因此，xi也是多项式oi（x）的所有系数的零点，也就是说xi是方程组S′的一个解。进一步，S′可以寫成如下形式：

其中vi（i=1，2，…，s）是δ中关于变元x1，x2，…，xk的所有幂积，D是系数矩阵。D则被称为Dixon矩阵，而D的行列式则被称为Dixon结式。当s′=s时，Dixon结式有零点是原方程组S有解的一个必要条件[25]。

然而，Dixon矩阵常常可能是奇异的，即其行列式总等于0，或者s′≠s无法计算Dixon结式，此时，不能提供与原方程组有解相关的有用信息，因此，为了解决这个问题，Yang等[23]提出了KSY方法，来计算扩展的Dixon结式。

令F={p1（x），p2（x），…，pk+1（x）}，给定对x的约束C，有如下形式x≠0。设Dixon矩阵的第i列为ui，它所对应的单项式为monom（ui），令nvcol（C）为所有的满足C monom（ui）≠0的列，而D1为所有通过删除Dixon矩阵中属于nvcol（C）的一列而得到s′×（s-1）的子矩阵集合。设a1，a2，…，am为m个参数，Q为有理数域的代数闭域，Φ：{a1，a2，…，am} → Q为一个映射，也就是把参数的值代入，获取Q上的一个值，令Φ（F），Φ（D）和Φ（D1）为把参数值分别代入到F、D和D1而得到的结果，再令R={Y|Y是D的一个r×r的非奇异子矩阵}。

定理1 如果X∈D1.rank（X）此处的X，是矢量、向量或矩阵吗？另外，rank前面（D1后面）的符号，是点还是逗号，请明确。回复：rank前面（D1后面）的符号，是点，也就是s.t.;

根据定理1，可以得到如下的算法：先检查是否有X∈D1.rank（X）

1.2 SDS方法

差分代换（Difference Substitution， DS）方法常被用于对称形式的不等式的证明中，Yang[24]定义了一个更加一般的DS的概念，用于证明对称以及非对称形式的不等式，并且提出逐次差分替换法（即SDS方法），此后，很多研究者[26-29]对该方法进行了进一步的扩展研究。利用SDS方法，可以证明多元多项式在正实数域上是半正定的。

一般地，给定一个n元多项式P，其变量x1，x2，…，xn有n！种排序结果，对诸如x1≥x2≥…≥xn的每个排序，对应着一个替代变换：

其中，ti≥0。通过代换之后，P可以转化为一个关于t1，t2，…，tn的多项式，而所有n！个代换得到的多项式组成的集合称为DS。进一步，如果该集合中的每个元素的所有系数均为非负的，则无论x1，x2，…，xn取何种非负值都有P≥0，也就是说P在正实数域上是半正定的，而如果该集合中存在一个元素有系数为负，则继续对该多项式的DS用同样的方式进行检查。以上这个过程被称为逐次差分代换法（SDS方法）。

更进一步，利用SDS方法，可以得到多项式P严格为正的一个充分条件，即如果DS中的每个元素的所有系数非负且常数项严格为正，则P在正实数域上是严格正定的。

2 不等式蕴含关系的证明问题及其解决方法

本章将解决的是带参数系数的多元多项式不等式蕴含关系的证明问题，该问题定义如下。

定义1 P1。设p（x）为一个带参数系数η1，η2，…，ηm的多项式，其中，x=（x1，x2，…，xn）T为一个含n元变量的向量。P1问题为：找到一组η1，η2，…，ηm的实数赋值，使得对任意的x∈Rn有：

证明 当k=n时，SY为一个关于n个变元x的由n+1个多项式等式组成的方程组，又有如果H（η，ξ）≠0则SY无解，进一步，也就是说SY无解，更进一步，也就意味着P1问题中的关于变量x∈Rn蕴含式（1）成立。证毕。

注1 本文将多项式H（η，ξ）称为式（1）在k=n的情况下的扩展Dixon结式。在这种情况下，为解决P1问题，只需找到一个η使得H（η，ξ）在ξ∈Rk+1+上严格正定即可。考虑1.2节中介绍的SDS方法，其中，需要确保DS中每一个替换后的多项式的系数非负，而对于H（η，ξ）而言，如果每个替换后的多项式的所有系数非负且常数项严格为正，那么则可确保H（η，ξ）在ξ∈Rk+1+上严格正定，也就是说，此时可以找到一组P1问题中所要寻找的η1，η2，…，ηm的赋值。

引理2 k

那么，关于变量x∈Rn蕴含式（1）成立。

证明 有如果H（η，，ξ）≠0则SY无解，进一步，也就是说SY无解，更进一步，也就意味着P1问题中的关于变量x∈Rn蕴含式（1）成立。证毕。

至此，只需找到一个η使得HS中的每个多项式Hi在β∈Rn+1+上严格正定即可。利用与注1类似的方法，可以找到一个η使得HS中的每个多项式Hi在β∈Rn+1+上严格正定，进一步，也就是说，此时，H（η，，ξ）在ξ∈Rk+1+，∈Rn-k上严格正定，更进一步，至此可以找到一组P1问题中所要寻找的η1，η2，…，ηm的赋值。

3 多项式秩函数探测问题

本章将介绍如何利用第2章的结果来为单路径多项式循环程序探测多项式秩函数。首先，将给出所针对的循环程序的定义;然后，将介绍为一类循环程序探测秩函数（即深度为1的嵌套秩函数）的方法;最后将介绍为更加一般的循环程序探测深度为d的嵌套秩函数的方法。

引理3 给定一个形如式（5）的循环程序P，P在实数域上不可终止的充分必要条件是其所对应的形如式（6）的循环程序P在实数域上不可终止。

证明 分为以下两种情况考虑。

1）k≤b。在这种情况下，有n=b，也就是说P即为P，此时，显然，P不可终止的充要条件为P不可终止。

2）k>b。在这种情况下，有n=k>b。

本节将针对形如式（6）的循环P的一种简单的情况，即k=n的情况，为这类循环探测秩函数。r（x）为一个Rq → R的函数，r（x）被称为P的秩函数，如果以下条件被满足：

设r（x）为一个关于x的多项式，显然，（x）=r（x）-r（M（x））-1也为一个关于x的多项式，令η=（η1，η2，…，ηq+1）T∈Rq+1为r（x）的参数系数组成的向量，η′为（x）的参数系数向量，显然有η′=Bη，其中B是一个确定的矩阵，由循环中的赋值M和r（x）确定。

其中，η′=Bη，H1和H2分别是蕴含式（7）和（8）的扩展Dixon结式，则P有对应的秩函数。

证明 已知在蕴含式（7）和（8）中，G（x）≥0是由k个多项式不等式组成的合取式，并且x为一个含k个变量的向量，因此，根据引理1，如果存在η使得：对任意的ξ1∈Rk+1+有H1（η，ξ1）≠0并且对任意的ξ2∈Rk+1+有H2（η′，ξ2）≠0，那么式（7）对任意的x∈Rk恒成立并且对式（8）对任意的x∈Rk恒成立，此时，根据秩函数的有界和递减的条件可知，P有对应的秩函数r（x），该秩函数的系数由η确定。证毕。

至此，再利用类似于注1中的方法，即基于SDS方法，可以找到一个η，使得H1（η，ξ1）和H2（Bη，ξ2）对任意的ξ1∈Rk+1+和任意的ξ2∈Rk+1+均严格为正，此时即可找到满足条件的秩函数。

3.3 深度为d的嵌套秩函数

多阶段秩函数的概念可以进一步扩展为一个更一般的概念（可以应对更多的循环程序），即嵌套秩函数[16]。接下来，将介绍为形如式（6）的循环探测深度为d的嵌套秩函数的方法。

显然，在3.2节中定义的秩函数为一个深度为1的嵌套秩函数。

4 算法和实验

深度为d的嵌套秩函数的探测算法如下算法1和算法2所示。在算法1中，SY指的是第3章中的式（2），Hi指的是定理3中的扩展的Dixon结式。

算法1 DETC。带参数系数的多元多项式不等式蕴含关系的探测。

输入：含d个关于n元变量的多项式的蕴含式的合取式I（在每个蕴含式中，前件为已知的k个多项式的合取式，后件为一个带参系数的多项式）。

输出：返回满足条件的参数系数的一组赋值。

过程：

1）为每个蕴含式中加入（k+1）个松弛变量，将I转化为形如SY的d个多项式方程组。

2）获取每个方程组的扩展Dixon结式Hi（i=1，2，…，d）。

3）如果存在一组参系数赋值使得H1，H2，…，Hd均严格为正，则返回该赋值;否则，返回空。

算法2 探测深度为d的嵌套秩函数。

输入：循环程序P（含n个变量，循环条件中有k个不等式，且k不大于n），秩函数模板r。

输出：如果存在深度为d的嵌套秩函数则返回True。

过程：

1）根据秩函数需满足的条件，将P和r转化为含d+1个蕴含式的合取式I;

讨论 两个算法的时间复杂度主要集中在算法1中的扩展Dixon结式的计算以及SDS方法的计算。首先是关于Dixon结式计算的复杂度，其中主要包括构造Dixon矩阵的复杂度（O（m21n！3∏ni=2m3i）[30]和计算Dixon结式的复杂度（即计算Dixon矩阵的行列式，O（n3n log n log log n））[31-32]。其次，关于SDS方法的计算复杂度，因为其过程是一个不可终止的过程，所以只能在给定差分代换次数上界的前提下来讨论其复杂度，给定一个n元多项式，对其进行1次差分代换，将产生n！个多项式，若对其进行m次差分代换，将产生（n！）m个多项式。这就需要对（n！）m个多项式进行分别处理（即从每个多项式中抽取出所有系数构成线性不等式系统，然后用线性规划工具进行求解，这个求解过程是多项式时间复杂度的），所以，复杂度较高。尽管如此，实验表明，对某些程序实例，与现有工具（如基于CAD的秩函数探测方法[5]）比较，本文方法能有效计算它们的秩函数。

实验环境为：一台处理器为Intel Core i5-6300HQ 2.3GHz、内存为16GB、操作系统为64位的Windows操作系统的计算机，利用的是Maple 18軟件。实验中将本文所提出的秩函数探测方法与基于CAD的秩函数探测方法[5]（利用Maple中的量消去工具RegularChains）进行了比较，表1是实验结果，其中循环程序10～13来自文献[16，20]。实验结果表明，利用本文提出的方法能够更加高效地探测到秩函数（仅需用几秒钟的时间）。特别地，对于表1中的循环程序11，利用本文的方法也无法计算其秩函数，这是因为本文方法找到的是秩函数存在的一个充分条件，而不是充要条件。下面将给出一些例子来具体说明秩函数探测的过程。

5 结语

本文基于杨路等的研究成果，即KSY方法和SDS方法，提出了一种为多项式循环程序探测多项式秩函数的方法，给定秩函数模板，寻找满足条件的秩函数模板中的参数，从而确定满足条件的秩函数。实验结果表明，对于一些多项式循环程序，可以有效地为其找到多项式秩函数（即深度为1的多阶段秩函数）和深度为d的多阶段秩函数。此外，与现有的基于SDP方法的解决方案不同，利用本文提出的方法，可以发现精确的秩函数，无需进行再次验证。本文方法能在较短时间内为多项式循环找到对应的多项式秩函数。在后续的工作中，本研究将针对多项式秩函数的探测问题进行进一步探索，以应对更多的多项式循环程序。

參考文献 （References）

[1] LEIKE J， HEIZMANN M. Geometric nontermination arguments [C]// TACAS 2018： Proceedings of the 2018 International Conference on Tools and Algorithms for the Construction and Analysis of Systems. Berlin： Springer， 2018： 172-186.

[2] FEDYUKOVICH G， ZHANG Y L， GUPTA A. Syntax-guided termination analysis[C]// CAV 2018： Proceedings of the 2018 International Conference on Computer Aided Verification. Berlin： Springer， 2018： 124-143.

[3] CHEN Y F， HEIZMANN M， LENGAL O， et al. Advanced automata-based algorithms for program termination checking[C]// PLDI 2018： Proceedings of the 2018 ACM SIGPLAN Conference on Programming Language Design and Implementation. New York： ACM， 2018： 135-150.

[4] LI Y. Witness to non-termination of linear programs[J]. Theoretical Computer Science， 2017， 681： 75-100.

[5] CHEN Y H， XIA B C， YANG L， et al. Discovering non-linear ranking functions by solving semi-algebraic systems[C]// ICTAC 2007： Proceedings of the 2007 International Colloquium on Theoretical Aspects of Computing. Berlin： Springer， 2007： 34-49.

[6] COLON M， SPIMA H. Synthesis of linear ranking functions[C]// TACAS 2001： Proceedings of the 2001 International Conference on Tools and Algorithms for the Construction and Analysis of Systems. Berlin： Springer， 2001： 67-81.

[7] COLON M， SPIMA H. Practical methods for proving program termination[C]// CAV 2002： Proceedings of the 2002 International Conference on Computer Aided Verification. Berlin： Springer， 2002： 442-454.

[8] PODELSKI A， RYBALCHENKO A. A complete method for the synthesis of linear ranking functions[C]// VMCAI 2004： Proceedings of the 2004 International Conference on Verification， Model Checking and Abstract Interpretation. Berlin： Springer， 2004： 239-251.

[9] SCHRIJVER A. Theory of Linear and Integer Programming [M]. New York： John Wiley and Sons， Inc.， 1986： 87-90.

[10] BEN-AMRAM A M， GENAIM S. Ranking functions for linear-constraint loops[J]. Journal of the ACM， 2014， 61（4）： Article No. 26.

[11] TURING A M. On computable numbers， with an application to the entscheidungs problem[J]. Proceedings of the London Mathematical Society， 1937， 42（2）： 230-265.

[12] BRAVEMAN M. Termination of integer linear programs[C]// CAV 2006： Proceedings of the 2006 International Conference on Computer Aided Verification. Berlin： Springer， 2006： 372-385.

[13] TIWARI A. Termination of linear programs[C]// CAV 2004： Proceedings of the 2004 International Conference on Computer Aided Verification. Berlin： Springer， 2004： 70-82.

[14] YANG L， ZHAN N J， XIA B C， et al. Program verification by using DISCOVERER[C]// VSTTE 2005： Proceedings of the 2005 Working Conference on Verified Software： Theories， Tools， and Experiments. Berlin： Springer， 2008： 528-538.

[15] BAGNARA R， MESNARD F. Eventual linear ranking functions[C]// PPDP 2013： Proceedings of the 2013 International Symposium on Principles and Practice of Declarative Programming. New York： ACM， 2013： 229-238.

[16] BEN-AMRAM A M， GENAIM S. On multiphase-linear ranking functions[C]// CAV 2017： Proceedings of the 2017 International Conference on Computer Aided Verification. Berlin： Springer， 2017： 601-620.

[17] BRADLEY A R， MANNA Z， SIPMA H B. Linear ranking with reachability[C]// CAV 2005： Proceedings of the 2005 International Conference on Computer Aided Verification. Berlin： Springer， 2005： 491-504.

[18] BRADLEY A R， MANNA Z， SIPMA H B. The polyranking principle[C]// ICALP 2005： Proceedings of the 2005 International Conference on Automata， Languages and Programming. Berlin： Springer， 2005： 1349-1361.

[19] LEIKE J， HEIZMANN M. Ranking templates for linear loops[C]// TACAS 2014： Proceedings of the 2014 International Conference on Tools and Algorithms for the Construction and Analysis of Systems. Berlin： Springer， 2014： 172-186.

[20] LI Y， ZHU G， FENG Y. The L-depth eventual linear ranking functions for single-path linear constraint loops[C]// TASE 2016： Proceedings of the 2016 International Symposium on Theoretical Aspects of Software Engineering. Piscataway， NJ： IEEE， 2016： 30-37.

[21] COUSOT P. Proving program invariance and termination by parametric abstraction， Lagrangian relaxation and semidefinite programming[C]// VMCAI 2005： Proceedings of the 2005 International Conference on Verification， Model Checking and Abstract Interpretation. Berlin： Springer， 2005： 1-24.

[22] SHEN L Y， WU M， YANG Z F， et al. Generating exact nonlinear ranking functions by symbolic-numeric hybrid method[J]. Journal of Systems Science and Complexity， 2013， 26（2）： 291-301.

[23] KAPUR D， SAXENA T， YANG L. Algebraic and geometric reasoning using Dixon resultants[C]// ISSAC 1994： Proceedings of the 1994 International Symposium on Symbolic and Algebraic Computation. New York： ACM， 1994： 99-107.

[24] YANG L. Solving harder problems with lesser mathematics[C]// ATCM 2005： Proceedings of the 2005 Asian Technology Conference in Mathematics. Radford， VA： Mathematics and Technology， 2005： 37-46.

[25] DIXON A L. The eliminant of three quantics in two independent variables[J]. Proceedings of the London Mathematical Society， 1909， s2-7（49）： 49-69.

[26] 杨路，姚勇.差分代换矩阵与多项式的非负性判定[J].系统科学与数学，2009，29（9）：1169-1177.（YANG L， YAO Y. Difference substitution matrices and decision on nonnegativity of polynomials [J]. Journal of Systems Science and Mathematical Sciences， 2009， 29（9）：1169-1177.）

[27] HOU X R， SHAO J W. Bounds on the number of steps of WDS required for checking the positivity of integral forms[J]. Applied Mathematics and Computation， 2011， 217（24）： 9978-9984.

[28] 韓京俊.基于差分代换的正半定型判定完备方法[J].北京大学学报（自然科学版），2013，49（4）：545-551.（HAN J J. A complete method based on successive difference substitution method for deciding positive semi-definiteness of polynomials[J]. Acta Scientiarum Naturalium Universitatis Pekinensis， 2013， 49（4）： 545-551.）

[29] 杨路.差分代换与不等式机器证明[J].广州大学学报（自然科学版），2006，5（2）：1-7.（YANG L. Difference substitution and automated inequality proving[J]. Journal of Guangzhou University （Natural Science edition）， 2006， 5（2）： 1-7.）

[30] QIN X， WU D， TANG L， et al. Complexity of constructing Dixon resultant matrix[J]. International Journal of Computer Mathematics， 2016， 94（10）： 1-16.

[31] LI Y. An effective hybrid algorithm for computing symbolic determinants[J]. Applied Mathematics and Computation， 2009， 215（7）： 2495-2501.

[32] KALTOFEN E. On computing determinants of matrices without divisions [C]// ISSAC 1992： Proceedings of the 1992 International Symposium on Symbolic and Algebraic Computation. New York： ACM， 1992： 342-349.