整数矩阵和多项式矩阵求逆的复杂性

1×|| |||| ||||||||||||×|||||| ||∈|||||| |||| |||| |||| ||||||∈||||||关于整数矩阵和多项式矩阵求逆的复杂性阿恩·斯托约翰David R.加拿大安大略省滑铁卢大学Cheriton计算机科学学院N2L 3G1二○ ○八年七月二十四日摘要提出了一种利用O<$（n3（logA+ logκ（A）位运算概率地计算非奇异nn整数矩阵A的精确逆的算法.这里，A= maxijA ij表示绝对值中的最大项，κ（A）：=A−1A是输入矩阵的条件数，软O符号O表示某些缺失的logn和log logA因子。对多项式矩阵给出了该算法的一种变形。 任何非奇异n-n矩阵的逆矩阵，其元素是域上的d次多项式，可以使用期望数量的O（n-3d）f ield运算来计算。这两种方法都是在拉斯维加斯的运行时间范围内进行的：失败的返回概率最多为1/2，如果没有返回，则在相同的运行时间范围内证明输出是正确的。1介绍设AZn×n是非奇异的.我们用A：= max A ij表示A中元素的最大值，用κ（A）：= A A−1表示矩阵关于最大范数的条件数。我们给出了一个LasVegas概率算法，该算法的期望运行时间为O <$（n3（log A + log κ（A）位运算来计算A的精确逆. 因此，对于一个条件良好的A，κ（A）由一个n logA的多项式函数限制，这个成本估计变为O（n3 log A）。为了比较，在非奇异AZn×n的逆中的所有条目的位长度的总和可以大于n3 log 2A位。据我们所知，整数矩阵求逆的最佳已知复杂度估计是O <$（n ω+1log A）位运算，由任何经典方法支持，例如同态成像和中国式迭代[6，第5.5节]，通过牛顿迭代的二次提升，或无分数高斯消除的递归版本[23，第2节]。这里，ω是环上矩阵乘法的指数[2，第1章]，软O符号O表示缺少log n和log log||一||因素2∈|| |||| |||| |||| |||| ||i=1i=1现在考虑多项式矩阵的情况。设K是一个域，且AK[x]n×n是非奇异的，其元素的次数以d >0为界.A的逆可能需要n个3d域元素的数量级来表示.与整数情况类似，各种经典的应用程序的polyn的omial矩阵的版本ex是t，所有的成本估计O<$（nω+1d）的字段操作从K。我们参考[10，第1节]对以前的方法进行了简要的调查，并讨论了在获得多项式矩阵问题的快速算法方面取得的一些进展。Jeannerod和Villard [10]提出了一种新的方法，对于n为2的幂的泛型A，将在O（n3d）场运算中计算A−1。在本文中，我们适应我们的方法整数矩阵求逆多项式的情况下。通过将一些代数预处理技术，特别是多项式，我们得到了拉斯维加斯概率算法，预期运行时间O（n3d）的任何非奇异输入。[10]中一般多项式矩阵的本质上最优的求逆算法的发现，以及最近在降低整数矩阵上的许多基本线性代数问题的复杂性方面所取得的进展，促使我们开发一种适用于整数矩阵的求逆的替代算法。 假设ω = 3，我们在接下来的两段中回顾了与本文特别相关的两个问题的一些关键结果：非奇异有理线性系统求解和行列式/史密斯形式计算。 对于计算这些和其他整数矩阵不变量以及结合快速矩阵乘法技术所做的工作的调查，我们参考[12，24]。早在四分之一个世纪以前，人们就已经证明，有理系统解A−1b可以用线性p进提升在O（n3logA）位运算中计算出来[1 5，3]。取b为单位矩阵的一列表明，逆矩阵的任何一列都可以在O（n3 log A）位运算中计算。线性p-adic提升是本文算法的一个关键子程序. 我们表明，提升可以用来计算所有n列的逆一个良好的条件矩阵在同一时间作为一个单一的列。现在考虑行列式的计算，假设ω=3。经典的方法，如同态成像[6，第5.5节]需要O<$（n4logA）位操作。一个突破是Kaltofen [11]的Krylov方法，它给出了O（n3. 5log A）Las Vegas算法。最近，Eberly等人给出了一种计算A的Smith形式的[5]的第10段。 史密斯表格诊断（s1，s2，. . . ，sn）是在幺模预乘和后乘下的典范对角化，参见第2节。 不变因子s i满足|det A|= s1s2···sn，因此，一旦已知了A的行列式，就可以很容易地求出A的行列式. 文[5]中的应用是建立在计算k∈O（n）A的Smith形式的不同不变因子，以及它们的重数，通过使用线性p-adic提升计算O<$（k）个有理线性系统解，每个解的代价为O<$（n3logA因此，[5]中的算法的总成本对k敏感，k是数字不同的不变因子。一方面，我们在这里提出的算法避免了这一点通过能够利用以下事实来提高对k的灵敏度：独立于A的不变结构。 另一方面，我们的算法是敏感的条件数为κ（A）。我们在时间上连续计算所有n个不变因子与大约n2<$n（logsi + logκ（A））位操作成比例3阿yXaIn−1u1阿yXa×−- × −- -A=⟨ ⟩≡−[||♩ ||♩−76−746992接下来，我们通过回忆关于秩为n的域上的n矩阵A的伴随的事实来激励我们的整数矩阵求逆的方法1，然后观察这个给出了Z上一类矩阵的伴随的一个简单公式。因为A的秩为n1，所以A至少有一个非零次维数n1。在不失一般性的情况下（直到行和列置换），假设A的前导（n 1）（n 1）子矩阵A<$是非奇异的，并且部分A为A=100%。如果我们设u：=−xA<$−1，arowvectorr，且v：=−A<$−1y，acolumnvector，则Σ Σ ΣΣ Σ=- 是的（一）在（1）的右手边的矩阵的伴随将具有除最后一行和最后一行中的元素之外的所有元素零，其中h ich将等于detA'。将（1）的b或h侧adj（detA<$）vdetA<$拉吉乌1美元。（二）注意，当A的条目来自于一个主理想环，甚至一个零因子环，如剩余类环的此外，提供了detA=0和detA′是来自环的单元（即，Aisivertible环）。 特别考虑一个非奇异矩阵A∈Zn×n，（n−1）×（n−1）子矩阵A<$满足detA<$detA。则在剩余类环上Z/detA，我们有一个detA0（moddetA）和detA′auni t，所以方程（2）有模detA.换句话说，Z上A的伴随与（2）中的秩1矩阵逐元素全等。此外，如果detA足够大，则可以通过将（2）中的外积相乘并在对称范围[（detA 1）/2， detA / 2]中减少所有模detA的条目来获得A在Z例如，矩阵−93 −32 8 44A=无菌−72 −4 99 −31−2 27 29 67In−1v1阿04−∈−||≡具有detA=6 4 33 4 0 45，其相对于detA<$=5 3 36 66而言是时间。公式（2）给出−28408A形容词861667181262533666澳门银河-27 25 95 0611 40 47 41-895 16 51澳门银河（mod6433 40 45）−853217 368292 −179420 −28408409178 −744548 233455 861667−584392 295157 438250 18126262584 183281 −289125 533666(mod64334045）。（三）对于这个例子，detA足够大，可以捕获Aadj中的所有元素;（3）中的矩阵，通过乘以外积并在对称范围中以64334045为模减少元素而获得，是A在Z上的伴随。为了使刚刚描述的方法适应于计算任意输入矩阵的逆，需要处理当A中的维度为n_1的所有子矩阵与detA具有公因子时的情况（即，A的Smith形式是非平凡的）。我们的解决方案是将公式（2）推广到具有非平凡不变结构的A，通过给出Aadjmod detA的表达式作为标度外积之和。Aadjadj.形容词Av uSdetA+vuSdetA+···+v u−S(mod detA）。n nnn−1n−1n1 1 11每个vi是列向量，每个ui是行向量。我们称这种结构为外部结构A的乘积伴随公式实际上，由于Aadj中的所有元素都可以被（detA）/sn整除，所以我们在第4节中对公式的定义如上所示，但将detA替换为sn，左手边替换为snA−1。我们注意到，任何主理想整环R上的任何非奇异ARn×n的外积伴随公式的存在性直接由R上的Smith形式的存在性导出。我们在本文中的贡献是双重的。首先，我们确定并利用公式的一些有用的性质：（1）当R=Z或R=K[x]时，外积伴随公式可以使用与表示A大致相同的空间来表示;（2）外积伴随公式给出了计算Rem（Aadjv，detA）的本质上最优的方法，对于给定的向量v，其具有模detA约简的条目。其次，我们提出了快速的概率算法来计算外积伴随公式。当应用外积伴随公式来计算Aadj时，如果detA太小而不能捕获对称范围模detA中的Aadj的条目，则会出现问题。在这种情况下，我们将Aadj的计算分为两部分。考虑分解sn A−1= Rem（sn A−1，sn）+sn R其中Rem表示模sn的余数。首先，我们通过外积伴随公式计算Rem（sn A−1，sn）然后我们使用线性p-adic提升计算剩余部分R：=（1/sn）（sn A−1Rem（sn A−1，sn））不幸的是，计算这两部分的成本对logκ（A）很敏感。5∈∗×∗∗∈ ∗∗∈−→−→本文的其余部分组织如下。在第2节中，我们修正了一些符号，并回忆了一些定义，包括史密斯标准形的定义。第三节考虑对非奇异的ARn×n，计算Smith分解：A =snf（A），其中n是R上的 n个矩阵。第4节定义并给出计算外积伴随公式的算法，外积伴随公式本质上是史密斯分解s n A−1= snf（s n A−1）。 第3节和第4节一般性地发展了结果，因此它们都适用于R=Z和R=K[x]。第五节给出了整数矩阵求逆的算法，第六节给出了多项式矩阵的变分算法。偶然的是，我们可以应用一些简单的代数预处理技术，将原始d次输入矩阵A K [x] n×n变换为一个新的矩阵B，它的行列式为nd。一旦通过外积伴随公式确定了B的逆，则可以在分配的时间内反转预处理以恢复A的逆。成本函数设M：Z>0R>0使得可以使用t个最多M（t）位运算来乘以大小为2t的整数。S_h_nha ge_Str a ss en a lg o rithm [2 1 ]都有M（t）= O（t（log t）（log logt））. 我们假设M（a）+M（b）≤ M（a + b）和M（ab）≤ M（a）M（b），其中a，bN≥2。我们参考[6，8.3节]以获得更多关于整数乘法的参考和讨论。定义一个额外的函数B来限制与整数gcd相关的计算的成本将是有用的。我们可以取B（t）=M（t）logt。然后，扩展的gcd问题的两个整数的大小有界的2t，和有理数重建问题[6，5.10节]的模有界的2t，可以解决O（B（t））域操作[20]（与[19]相比）。我们将稍微重载符号，并使用M：Z≥0R>0作为多项式乘法的代价函数：两个次数为d的K[x]中的多项式最多可以使用M（d）个域运算相乘。类似地，B将被用作与gcd相关的问题的成本函数，如有理函数重构和扩展的gcd。与整数情况类似，我们可以取B（d）=M（d）log d。我们参考[6，11.1节]了解更多细节和参考。2定义、符号和术语设R是一个主理想环，一个有单位元的交换环，其中每个理想都是主-NR。本文主要研究整环R=Z和R=K[x].在[17]之后，我们规定了一个完整的非结合集A（R），并且对于每个非零s∈R，规定了一个完整的剩余集R（R，s），如下所示。（四）R=ZR=K[x]（ZR（Z，s）=−）={0Σ，1，2，、、、. .. }A（K[x]）={ 0|− 1|− 12、|S|2，R（K [x]，s）={f ∈ K [x]|度f<度s}}f ∈ K [x] |f是monic6∈A∈∈AR∈R∈∈⟨ ⟩阿格拉河AR|≤≤−∈ A≤≤注意，我们对R（Z，s）的选择对应于通常的对于非零N，函数Rem（a，N）返回R（R，N）中与a模N全等的元素。下一个引理是我们对R（R，N）的选择的结果。引理1. 设a，N ∈ R，其中N非零. 则a = Rem（a，N），如果R = Z：|N |≥ 2|一|+2个R=K[x]：度N≥度a+1对于a，sR，我们用gcd（a，s）表示由a和s生成的理想在（R）中的唯一主生成元。我们允许gcd取任意数量的参数，包括矩阵和向量以及R的单个元素。例如，如果B是R上的矩阵，s是R的元素，则gcd（B，s）表示s和B中所有元素的gcd。我们可以使用R上的gcd的定义来归纳和的定义。对于残留物类环R/ s的定义.这将在下面有用，当我们回忆如何通过对R/s进行处理来计算R上的史密斯形式时，对于一个良好选择的s。对于非零sR，我们用元素集（R，s）标识剩余类环R/ s，并定义A（R/s）={gcd（a，s）|a ∈ R（R，s）}和R（R/gcd，b）= R（R，gcd（b，s））.A和R的这些选择使我们能够容易地获得基本运算的算法，关于R上基本运算的算法，参见[23，第1节]。史密斯标准形对应于每个ARn×n，存在R上可逆的么模矩阵P，QRn×n，使得snf（A）= S = UAV = Diag（s1，s2，. . . ，Sn），其中S是Smith标准型[17，第二章]，即对于1in1，si（R）为si + 1. 当R是主理想整环时，Sn是（detA）/gcd（Aadj）的伴随.因此，最大不变因子sn是R的最小非零元素（K[x]上的最小度和Z上的最小幅度），使得sn A−1在R上。计算环R=Z或R=K[x]上Smith型的经典方法是应用一系列初等行和列运算。一个众所周知的问题是，工作矩阵中的条目可能会变得过大。为了避免这种现象，许多作者（例如，[4，9，8]）已经使用了以下引理的思想与引理1的思想相结合。该引理独立于和的选择而成立，并且从任何主理想环上的史密斯形式的存在性和唯一性得出[13]，特别是在R/S上。引理2. 设ARn×n是非奇异的，sR是A的最大不变因子Sn的非零倍. 如果S是A在R上的史密斯形式，并且S<$是A<$=Rem（A，s）在R/S <$s上的史密斯形式，则S<$=Rem（S，s）。7x1 x2 x3 ···xn11. ..11−u/s1In−1s1uv1−v/s1In−1S1阿Σ∈∈∈为了帮助阐明在随后的章节中所发展的算法，我们在这里简述了[7]中将非奇异A∈Zn×n变换为Smith形式的方法。假设我们已经预先计算了s，sn的非零倍数，如引理2所示。计算c2，c3，. . .，c n∈ Z使得s的gcd与列向量A1c2gcd（s，A）.现在计算x1，x2，. . .，x n，其中x1是，使得c3 ···cnT等于1一1个用户. ..=、吉 尔科夫1其中s1等于s和A中所有条目的gcd。接下来，应用一些更基本的行和列操作，将s1下面和右边的条目清零。ΣΣ ΣΣ ΣΣ=ΣΣ.通过应用引理2的思想，A的n次尝试可以被简化为模。在一个类似的fa shion中，通过对A的诅咒，重新计算了remaining inv a riatfactors。 实际上，由于s1可以整除Di ag（s1，A<$）中的所有n次尝试，因此我们可以用（1/s）A<$和模s/s1来 代替。3Smith分解设BRn×n是非奇异的，且snf（B）=Diag（g1，g2，. . .，g n）。在本节中，我们采用上一节中描述的史密斯形式计算方法，不仅计算B的史密斯形式，而且还构造列向量v1，v2，. . .，vnRn×1和行向量u1，u2，. . .，unR1 ×n，使得B可以表示为标度外积之和：B= T n+1= g1v n u n+ g2vn−1un−1+···+g n v1u1。（五）请注意，最初我们将直接在R上工作。之后，我们将应用引理2的思想来避免表达式膨胀的问题。为方便起见，定义g0= 1。定义ej=gj/gj−1使得gj=e1e2···ej，1≤j≤n。我们将构造矩阵B1，B2，. . . ，Bn+1∈ Rn×n使得TJ与1Bj= gj−1（B−n（g1vnun+g2v n−1n−1 +x`·· ·+gj−1vn−j+2un−j+2））（6）gj/gj−1 gj+1/gj−1gn/gj−110） A = A（你好，ejex`j+ ，的。 . . ，n∈j∈j+x1，0，0，. . . ，0）。 （七）.CnC21C3181xIn1y In1−un−i+1In埃伊乌v Bi∈∈输入：·nonsingularB∈Rn×n设B1=B.对于i从1到n，设x∈R1×n，y∈Rn×1使得xBi y=gcd（Bi）.v：=Bi y;u：=xBi;ei：=xv;设Bgi：=e1e2. . . ei; vn−i+1：=v/ei;1un−i+1：=u/ei;n一期+1=（B−e vei伊因u×n断言：对于j=i+1，等式（6）和（7）在R上成立。−i+1n−i+1）∈R.od断言：B = g1v n u n+ g2vn−1un−1+···+g nv1u1。考虑（7），snf（Bn+1）将是零矩阵，表明（5）将从（6）得出，j = n +1。现在我们解释如何通过对j，1 ≤ j ≤ n +1的归纳来构造B j，其中基情形j= 1。初始化B1=B。则B1满足（6）和（7）。假设满足（6）和（7）的Bj已经被计算出来，对于1≤j≤i，对于某些i，其中1≤j≤n。我们展示了如何计算满足（6）和（7）的Bi+1 因为R是主理想环，存在向量xR1×n和yRn×1使得ei：=xBiy，Bi的所有元素的gcd. (We推迟到第5节和第6节来描述如何计算x和y。 如果我们嵌入到一个较大的零矩阵中，矩阵的非零不变结构不会改变。设v：= B i y和u：= xB i，并考虑从B i通过增加初始行和列的零而获得的矩阵的幺模变换。Σ Σ ΣΣ Σ Σ=Σ中国（8）设vn−i+1：=v/ei和un−i+1：=u/ei，并将另一个幺正矩阵n应用于（8）右侧的矩阵，以将e i下方和右侧的元素归零。Σ Σ Σ Σ Σ=中文（简体）在 （ 9 ） 的 右 手 边 的 矩 阵 的 所 有 元 素 都 可 被 ei 整 除 。 设 Bi+1= （ 1/ei ）（Bi−eivn−i+1un−i+1）∈Rn×n，满足（6）对j=i+1.由于（9）右边的m a与B i具有相同的非零不变因子，我们得出结论：（7）对j = i +1也成立。下面的代码片段总结了刚才描述的过程，通过归纳得出其正确性。图1：计算Smith分解图2调整图1中的算法以计算（5）中的分解，该分解仅保持模s，其中s是gn的非零倍数。所有操作都是显式0B我埃伊乌v Bi1−vn−i+1IneiB我 −eivn−i+1un−i +19≤≤≤≤输入：·nonsingularB∈Rn×n设B1=B.• s∈R，B对于i从1到n，v：= Rem（Bi y，s/gi−1）;u：= Rem（xBi，s/gi−1）;ei：= Rem（xv，s/gi−1）;设x∈R1×n和y∈Rn×1使得Rem（xBi y，s/gi−1）= Rem（gcd（Bi，s/gi−1），s/gi−1）.如果ei= 0，则断开fi;设Bgi：=e1e2. . . ei; vn−i+1：=v/ei;un−i+1：=u/ei;1n一期+1=（B−e vei伊因u断言：等式（6）和（7）对于j=i+1保持模s/gi−i+1n−i +1 ）∈R.×nod作为序列：Bg1vnun+g2vn−1un−1+···+givn−i+1un−i+1（m ods）。图2：计算Smith分解模s在R上，但引理2在迭代i时应用以计算模s/gi−1。有两个微妙之处需要注意。首先，在迭代i时，矩阵Bi只有模s/gi−1才是正确的。通过引理2，我们需要计算ei作为R/s/gi−1上Bi的所有元素的gcd：在R上我们计算s/gi−1的gcd，其中所有元素在Bi中，然后约简模s/gi−1。其次，如果s是gn的n个类，则存在极小k≤n，使得gk=gk+1=···= g n，在这种情况下，B k与零矩阵模s/gk−1全等：这可以通过退出循环来处理，因为在（5）中的sumgkvn−k+1+gk+1vn−kun−k+···+g n v1u1的剩余部分已知在这种情况下与零模s全等。4外积伴随公式设R是主理想整环，A∈Rn×n是非奇异的，满足Smith形式S=UAV = Diag（ s1，s2，. . . ，s n）。设vi和ui分别是幺模矩阵V和U的第i列和第i行，1i n。将方程S =UAV的两边取反，乘以s n，然后求解s n A−1， 得到s A−1=V（s S−1）U=snv uSSn+vu+···+snv u.（十）−S无无无无无无无nn−1n−1n1 1 11请注意，snf（s n A−1）= Diag（s n/s n，s n/sn−1，.. . . ，s n/s1）。现在考虑取方程（10）模sn。由于每个外积vi ui都被sn/si缩放，如果vi和ui中的条目以si，1in为模进行约简，则等式仍将以sn为模保持。这一概念通过以下定义得到了明确。定义3.非奇异A ∈ Rn×n的一个外积伴随公式是元组集（sn−i +1，vn−i +1，un−i +1）1 ≤i≤k 使得：10R···.Snn−1sn−k+1我n我外积伴随（A，sn）输入：·nonsingularA∈Rn×n输出：A的外积伴随公式。设B1=snA−1且sn+1=sn。对于i到n做• Sn∈R，A的最大不变因子如果ei=0，则断开fi;u：=Rem（xBi，sn−i+2）;v：=Rem（Biy，sn−i+2）;ei：=Rem（xv，sn−i+2）;设x∈R1×n，y∈Rn×1满足Rem（xBiy，sn−i+2）=Rem（gcd（Bi，sn−i+2），sn−i+2）.sn−i+1：=sn−i+2/ei;设B=（B-e1un−i+1：=u/ei;vn−i+1：=v/ei;n一期+1ei伊因 vuod;−i+1n−i +1 ）∈R.×nreturn（sn−j+1，vn−j +1，un−j+1） 1≤j≤i− 1• A的Smith形式是Di ag（1，1，. . . ，1，sn−k+1，sn−k+2，. . . ，sn），其中sn−k+1=/1，• vn−i+1∈R（R，sn−i+1））n×1和un−i+1∈R（R，sn−i+1）1×n，其中1≤i≤k，• snA−1 <$snvnun + snvn−1un−1+···+snvn−k+1un−k+1（mod sn）。图4所示的算法OuterProductAdjoint将计算外积伴随公式。该算法与图2所示的算法相同，除了s=s n和gi= sn/sn−i+1。图3：算法OuterProductAdjoint当我们专门研究环R=Z和R=K[x]时，上面的大多数循环都将完全按照所写的那样实现。注意，需要在每次迭代的开始只需要构造v。我们将求助于已知的结果来构造v作为B i Y的R-线性组合，其中Y是随机选择的，使得gcd（BiY，sn−i+2）=gcd（Bi，sn−i+1），其中hp rb b i ty为高。 一旦v是已知的，则x被计算为扩展的欧几里德问题的解。该技术的主要目的是有效地计算Rem（BiY，sn−i+2），其中Y∈（R，sn−i+1）n×m。 （计算Rem（xBi，sn−i+2）的向量矩阵的问题是类似的，所以 我 们 在 这 里 关 注 矩 阵 向 量 的 情 况 。 与 第 3 节 完 全 相 同 ， 但 这 里 s=sn ，gi=sn/sn−i+1，如果我们可以的话T=snv u snSn+vn−1n−1n−1Sn++vn−i+2n−i−1un−i−1 ∈Rn×n，（11）则在循环迭代i= 1，2，. . .，k，我们使用矩阵B=sn−i+2sSnA−1−Ti∈Rn×n（12）nn11. Σ−s“我的天，i=1n−j+1n−j +1n−j+1n−j+1SSSSnO ute rP rod Mu l（s，（sn−j+1，vn−j+1，un−j+1）1≤j≤i−1，Y）输入：·nonzeros∈ A（R）• vn−j+1∈R（R，sn−j+1）n×1和un−j+1∈R（R，sn−j+1）1×n，1≤j≤i−1输出：Rem（（• Y∈ R（R，s）n×m• sn−i+2|sn−i+3|···|sn，其中sn=s且sn−i+2为非n单位，如果i>1Σ我j=1−1（秒/秒）n−j+1n）v−j+1nu −j+1）Y，s）。Yn+1：=Y;sn+1：=sn;如果i≤1，则返回n×m零矩阵fi;对于j从1到i- 1，od;Yn−j+1：= Rem（Yn−j +2，sn−j +1）; Xn−j+1：= Rem（vn−j+1（un−j+1Yn−j+1），sn−j+1）V：=Xn−i+2;返回V对于j，从mi−2downtoo1doV：=Rem（Xn−j+1+（sn−j+1/sn−j）V，sn−j+1）od;i=1和m odulussn−i+2，其中sn+1被定义为sn。对于给定的Y∈ R（R，Sn）n×m，利用嵌套乘法可以有效地计算R em（Ti Y，Sn）c.TiYI1svn−j+1 un−j+1Y（modsn）j=1i−1n−j+1Yn−j+1nRem（v（uRem（Yx，`s），s）j=1sn−j+1Xn−j+1无意义的;Xn+···+sn−i+5。X+sn−i+4。X+sn−i+3XΣΣ···Σ图4中所示的算法Outer Mul实现了上述方案。 在第一个l op中，i teratio njdoesO（nm）arithmeticoperi onsmodulosn−j+1，其中操作符来自R（R，sn−j+1）和R（R，sn−j+2）。第二个循环的迭代j是 O（nm）的运算 模 n− j +1 ， 运 算 模来 自R（R，sn−j+1）。 Lemma4从M.图4：算法外接Mul引理4. 算法4（外接Mul）是正确的。算法的代价是R=Z：O（nmM（l ogQksn−i+1））位操作。nn−i+4n−i+4n−i+3n−i+3n−i+2n−i+212R=K[x]：O（nmM（numk degsn−i+1））域运算。13..|||| ||||||||∈现在考虑Rem（BiY，sn−i+2）的计算。请注意，thatmatricesBi并没有明确计算。相反，根据（12）中Bi的定义，我们有：Rem（BiY，sn−i+2）= Remsn−i+2sSnA−1Y−Rem（Ti Y，sn），sn−i+2设N∈R与Sn互素.如果Rem（Bi Y，N）=Bi Y，则还Rem（Rem（BiY，N），sn−i+1）=Rem（BiY，sn−i+1），我们可以应用以下配方。W1：= Rem（A−1Y，N）;Rem（BiY ，sn−i+2）=W2：=O ute rP rod Mu l（s，（sn−j+1 ，vn−j+1 ，un−j+1）1≤j≤i−1，Y）;（十三）返回nRem（Rem（sn−i+2W1−（sn−i+2/sn）W2，N），sn−i+2）5. 如果Y∈R（R，sn−i+2）n×m，则recipe（13）返回Rem（BiY，sn−i+2），如果R=Z：N≥sn−i+2（nsn−i+2||+1）+ 2。||+1)+2.R=K[x]：degN≥max x（degsn−i+2，degsnA−1+degsn−i+2）。证据 矩阵s n A−1Y和Rem（T i Y，s n）在R上，它们的差可被sn/sn−i+2整除。如果我们表明N满足引理1的界，则结果将如下，其中B i Y扮演a的角色。在K[x]上，我们有degsn A−1Y≤ degsn A−1+ degY≤ degsn A−1+ degsn− 1和degRem（Ti Y，sn）≤degsn− 1。定理中给出的度N的界是通过在两个界的最大值上加1并减去度sn/sn−i+2得到的。O v er Zweh ave||≤nsn||A −1||Y||≤ ns n||s n − i + 2 / 2 ， ||Rem （ T i Y ， s n ）||A−1||sn−i+2/2and||Rem (TiY,sn)||≤sn/2，因此，由ysn−i+2/sn相乘的矩阵的差的大小由ynsn−i+2A−1有界Y+sn−i+2/2;定理中N的界是通过将这个界乘以2再加上2得到的。5计算整数矩阵的逆算法IntInverse如图5所示。阶段1使用随机化来概率地计算输入矩阵AZn×n的最大不 变 因 子 ， 以 及A-1 。 阶 段 2 通 过 引 入 随 机 化 来 调 整 图 3 所 示 的 算 法OuterProductAdjoint，以获得用于计算A的外积伴随的概率算法。第三阶段证明了第二阶段计算的外积伴随公式的正确性，并计算了A的显式逆。在下面的小节中，我们填写该实现详细描述并证明了每个阶段的概率结果和复杂度估计。首先，我们回顾一些已知的结果。通过一个有理矩阵或向量的分母，我们意味着最小的正整数倍，将清除的分解器，假设分数减少。n14∈我j=1ei∈≥/-∈IntInverse（A）输入：非奇异AZn×n。输出：A-1或失败。 失败将以1/ 2的概率返回。<1. [编辑]从第一个12[logn2（nn/2）]开始，选择唯一且唯一||n）n素数;||n)♩primes;p：=当p≥log时，p的整数幂中的最小psitve如果p<$/pdetAreturnnfailseC：=Rem（A−1，p）fi;n+ log||一||;α：= ||A−1X||其中X∈ {− 1，1}n×4是均匀随机选取的;M：= 6 +[2 n（log 2n + log 2||一||）|;L：={0，. . . ，M − 1};s：=A−1X的分母，其中X∈Ln×12是均匀随机选择的;m：= 2 [（2 + log 2 n）/log 2 3|;2. [计算外积伴随公式]设B1表示esA−1，且sn+1=s。对于i到n做均匀随机地选取Y∈Ln×mY：=Rem（Y，sn−i+2）;集合N：=pk，其中k∈Z>0是极小的. t. pk≥sn−i+2（nαsn−i+2+1）+2;计算V：=Rem（BiY，sn−i+2）使用recipe（13）;使用引理9来计算x和y使得gcd（xVy，sn−i+2）=gcd（V，sn−i+2）;v：=Rem（Vy，sn−i+2）;计算u：=Rem（BTxT，sn−i+2）T使用recipe（13）;ei：=Rem（xv，sn−i+2）;ifei=0thenbreakelsesn−i+1：=sn−i+2/eifi;if（i=1andsn/=s）OrQisn−j+1>nn/2||n或e i /||norei/|uthenreturnFailfi;un−i+1：=u/ei; vn−i+1：=v/ei;设Bi+1表示1（Bieivn−i+1un−i+1）Zn×n。od3. [计算反演和测定正确性]若O ute rP rod uc t Mul（s，（sn−j+1，vn−j+1，un−j+1）1≤j≤i−1，A） 0n×n，则n个环n个域fi;C0：=O ut e r P ro d u c tM ul（s，（sn−j+1，vn−j+1，un−j+1）1≤j≤i−1，In）;N1：=p k，其中kZ>0是最小s. t。 pks2αs+2;C1：= Rem（sRem（A−1，N1），N1），使用p-adic提升计算C：= Rem（N1Rem（1/N1，s）C0+sRem（1/s，N1）C1，N1s）;N2：= pk，其中k ∈ Z>0是极小s.t. p k≥（2 n/s）||一||||C||+2;如果Rem（A Rem（（1/s）C，N2），N2）=I n，则返回失败fi;返回（1/s）C，分数减少图5：IntInverse算法15{−}[|||||∈≤||||∈∈|| ||≤≤||||∈Qj=1|| ||× × ∗ × ∗ ∗×引理6（[5，定理2.1]）. 设A ∈ Zn×n是非奇异的. 若X ∈ Ln×2是均匀随机选取的，其中L：= 0，1，. . . ，M 1且M：= 6 + 2 n（log 2n +log 2A），则A−1X的分母是A的最大不变因子，概率至少为1/3。对[5，定理2.1]的证明的考察揭示了引理6中M的定义是由A的最大不变因子Sn的大 小 驱 动 的 ， 并 且 与 A 和 n 无 关 。 由 于 sn 是 一 个 因 素 的 detA ， 我 们 有sndetAnn/2An，第二个不等式是隐含的阿达玛下面的结果是[5，定理2.1]证明的推论推论7. 设sZ>0满足界sn n/2A n.对于任何BZn×n，如果X Ln×2是均匀随机选择的，则gcd（BX，s）= gcd（B，s）的概率至少为1/3。我们算法中的主要计算工具是线性p-adic提升[3，15]，用于计算给定X的Rem（A−1X，p k）Zn×m和k.作了介绍和分析在[3，15]中假设了标准的整数运算，但是如果p的位长选择得很好，那么通过呼吁用于算术运算（例如基数转换）的快速算法，可以没有太大困难地并入快速整数运算。我们参考[16，第5节]，以获得以下部分（a）和（b）的推导。引理8. 设A ∈ Zn×n是非奇异的. 假设我们有p ∈ Z>0，其中p ∈ A，logp ∈ Θ（log n + log||一||），加上C：= Rem（A−1，p）。(a) 若N：=pk且Y∈Zn×m满足log ||Y|| ∈O（log N），则Rem（A−1Y，N）的计算时间为O（n2km B（log n + log||一||）+nm B（k（log n + log||一||））位操作。(b) 若Y ∈ Zn×m满足log ||Y||时间复杂度为O（n（log n +log||一||）），则A−1Y可以计算为时间复杂度O（n）||一||））位操作。(c) 考虑recipe（13）。 如果i−1sn−i+1≤nn/2||n和log N ∈ O（l o g s n − i +2 + l o gn + log||），则Rem（B i Y，N）可以在下式中计算：||A−1||), then Rem(B i Y, N)can be computed inO（nm）B（n（10gn+10 g||））+ n 2 k n − i + 2 m B（l o g n + l o g||））||A||))位操作，其中k：=1+logsn−i+2+logκ（A）n−i+2其中k（A）：=||A−1||||一||.logn + log||一||logn + log||一||证据 关于部分（a）和（b）的详细推导，我们参考[16，第5节]。为了理解（c）部分，回忆一下（a）部分的成本分析会有所帮助。p-adic提升的每一步都涉及两个矩阵乘积：和C哪里是nm矩阵，其元素的位长度为O（log p）。这就产生了O（n2kmB（log n+log A））项。第二个任期限制了所有额外工作的成本，特别是在最后16∈∈∈R∈R|| |||| ||⊥≥||||∈||||≤||||||||||提升，基数转换从p-adic到标准表示的nm数界定的位长度由logN。现在考虑（c）部分。算法Outer Mul计算W2的时间复杂度为O（nmM（n（logn+log||一||））位操作（引理4）。这也限制了在recipe（13）的return语句中执行的算术运算的成本。计算W1：=Rem（A−1Y ，N）的成本是在（a）中陈述的，其中k=loggpN。 kn−i+2的定义是满足tkn−i+2≥1且kn−i+2∈Θ（logpN）.在p art（c）中估计的成本通过使用kn−i+2∈O（n）简化了从p a rt（a）估计的成本的第二项。通常的扩展欧几里德问题以一个列向量v Zn×1作为输入，并要求一个行向量xZ1 ×n，使得xv= gcd（v）。