CC BY
37
УДК 512.93
РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ УРАВНЕНИЙ В КОММУТАТИВНЫХ ОБЛАСТЯХ
© r.H. MajiauionoK
Malashonok G.I. Solving systems of linear equations in commutative domains. The article looks at fast methods of linear systems solving in commutative domains and in their fields of fractions. It proposes realisation of these methods for the rings of principal ideals with the help of Hensel lifting techniques. The complexity of the deterministic algorithm to find all the rational solutioas for linear systems with coefficient matrix A of size n X m and rank n with p-adic lifting is 0(«-’/w(/M-«+l)(logn+Iog||/l l+log/?)2) bit operations. The randomised algorithm to find all the Diophantine solutions for this system has essentially the same complexity.__________________________
ВВЕДЕНИЕ
Пусть К - коммутативная область с единицей, К - поле частных И. Одной из важных задач вычислительной коммутативной алгебры является задача решения систем линейных уравнений Ах = Ь, коэффициенты которых берутся из коммутативной области Я: А Е ДпХш, Ь е Яп.
При этом могут ставиться разные задачи: (1°) найти одно решение в Ктп, (2°) найти все решения в Кт, (3°) найти одно решение в К"1, (4°) найти все решения Кт. Рассмотрим каждую из этих задач для разных областей.
Область К является полем. Основным вычислительным методом здесь является метод исключения Гаусса и его модификации со сложностью 0(п2т) операций в поле Я. Быстрые методы решения систем линейных уравнений над полем ведут свое начало от работы Штрассена [1] и имеют такую же сложность, как и матричное умножение -0(пР~1т). 0(г?/’) - это число операций,
необходимое для умножения квадратных матриц п-го порядка. Для обычного умножения матриц /3 = 3, для алгоритма Штрассена /3 = 1о§2 7, а для лучшего на сегодня алгоритма /3 < 2.376 [2].
Произвольная коммутативная область
К с единицей. Оценку алгоритмов будем
производить по количеству операций над
элементами 11. При этом предполагаем, что в Я существует алгоритм сокращения, то
есть если а = Ьс и известны а и с, с ф 0, то существует алгоритм вычисления 6. Основным вычислительным методом здесь является метод Доджсона, иногда называемый методом Гаусса с точным делением, со сложностью 0{п2т)
операций в R. Впервые такой метод был опубликован в работе Доджсона [3], а его улучшения получены в работах [4|, [5], [6]. Быстрый метод решения систем линейных уравнений над коммутативной областью опубликован в работе [7]. Сложность этого метода 0(nf3~1m) такая же, как и сложность матричного умножения. И теперь дальнейший прогресс в задачах (1°) и (2°) целиком связан с улучшениями алгоритма матричного умножения.
Вероятностный алгоритм для решения задач (3°) и (4°) при самых общих предположениях относительно коммутативной области обсуждается в настоящей работе. Предполагается, что в области R имеется алгоритм, который может установить, является ли конечнопорожденный идеал / = («!, а2,, as) единичным. Если / является единичным, то этот алгоритм вычисляет коэффициенты k{ € R в разложении единицы 1 = Yli=1 Например, достаточно иметь
алгоритм вычисления в области R базиса Г ребнера.
Коммутативная область главных идеалов R. Это основная прикладная область. Лучшим методом для задачи (1°) в случае решения определенных систем в поле рациональных чисел (R = Z, т = n, det А Ф 0) остается метод Диксона [9|, использующий линейный р-адический подъем. Его сложность O(n3(logn + log || Л|| -f- logp)2) бит-операций.
Если применять метод Гаусса с точным делений для R = Z, то для решения задачи (2°) в обычной арифметике потребуется O(mn4(logn + log ||у!||)2) бит-операций. || Л || обозначает модуль наибольшего из коэффициентов системы. Если применить китайскую теорему об остатках, то
сложность можно уменьшить до величины 0(nin3(\og7i + log ||А||)2) бит-операций [8|.
Мой подход к решению задачи (2°) в области главных идеалов состоит в применении р-адического подьема, как и в методе Диксона [9]. Предлагаемый в настоящей работе метод имеет сложность 0(п2т(т — п + l)(logm + log ||v4|| + logp)2) бит-операций.
В работе [10] предлагается вероятностный алгоритм для задачи (3°) - получения одного решения в кольце целых чисел Z и кольце F[x] многочленов над полем F, который опирается на алгоритм Диксона.
Я предлагаю вероятностный алгоритм для задачи (4°), который основан на р-адическом подьеме и сложность которого имеет существенно кубическую зависимость от порядка системы, как и [10], но количество вычислений обратных матриц становится в m — п + 1 раз меньше.
Общий план всей работы следующий. Во втором параграфе приводится обоснование метода нахождения всех решений в поле частных и соответствующий алгоритм для коммутативной области с сокращением. В третьем параграфе приводится алгоритм решения систем линейных уравнений в поле частных с помощью линейного р-адичес.кого подъема. В четвертом параграфе обсуждается вероятностный метод получения всех решений системы линейных уравнений в коммутативной области и соответствующий алгоритм. Пятый параграф посвящен оценкам сложности этих алгоритмов. В заключении сравниваются предлагаемые алгоритмы получения точных решений систем линейных уравнений со стандартными приближенными методами решения систем линейных уравнений, использующими представление чисел с плавающей точкой.
РЕШЕНИЕ СИСТЕМЫ В ПОЛЕ ЧАСТНЫХ
Пусть R - коммутативная область с единицей, К - поле частных области R, А в RnXm, с € Rn,
Ах = с (1)
- система линейных уравнений над R. Пусть гапкЛ = г, S и Т - матрицы перестановок, такие, что det Л о 0, где Ао верхний левый блок размера г х г матрицы
SAT = ( £ А, ) и пусть Sc - ( С° )' со £ R7. Ai = {х\х 6 Кт,Лд: = с} -
множество решений системы (1). Единичную матрицу порядка г будем обозначать /г, а также обозначим - квадратные ма-
трицы, у которых только один ненулевой элемент - (г, ^’), и этот элемент равен единице. Будем пользоваться этими обозначениями на протяжении всей статьи.
Нам потребуются некоторые факты из теории систем линейных уравнений.
1. Если гапк(Л, с) ф г, то М. =0. Если гапк(Л, с) = г, то М. - это гиперплоскость размерности 771 — 7’ в пространстве размерности ттг. Она определяется т — г + 1 точками, которые не лежат в одной гиперплоскости меньшей размерности.
2. Если система (1) однородная (с = 0) и х\, Х2>... >хт-г ее линейно независимые решения, ТО М = {Х)Г=ТГ+1 х1иг\щ € К}.
3. Если система (1) неоднородная (с ф 0) и
Х\, .г’2,... ,а:т_г+1 ее линейно независимые решения, то М. = {Ц1^1Г+1 %ОН\щ €
К, = 11-
Определение. Базисным множеством решений однородной системы линейных уравнений (1) называется множество, состоящее из т — г линейно независимых решений системы (1).
Базисным множеством решений неоднородной системы линейных уравнений (1) называется множество, состоящее ИЗ 771 — г + 1 линейно независимых решений системы (1).
В следующих двух теоремах задача вычисления базисного множества решений системы (1) сводится к нескольким задачам решения определенных систем. В первой теореме рассматриваются однородные системы, во второй теореме - неоднородные системы.
Теорема 1 Пусть (1) - однородная система линейных уравнений и А\ = (с^, а2,..., ат_г), а2 € R7. Тогда системы
AoXj = -aj> j = 1,... ,77г - г (2)
определены, и их решения X] € Кг задают базисное множество решений системы (1):
Т( ^ )’ =
где е^ Е Rm-r — это столбцы единичной матрицы 1т—г (б1, 62, • • • ) &т—т)'
Доказательство. Обозначим у = Т~1х. По условию теоремы (Ао,А\)у = 0. Если искать
решение в виде у = ^ ^, то придем к
системе (2).
Следствие 1 Пусть В = А^ А\ = (61,62,
6т_г), 6^- Є Кг, тогда е/
т( ~е22 )■ ■■■■ т{ ) 'Гаэис’“>е
множество региений однородной системы линейных уравнений (1).
Доказательство. В (3) подставим решения систем (2).
Теорема 2 Пусть (1) - неоднородная
система линейных уравнений, Р - матрица перестановок, такая, что последний
элемент вектора 6 = РАц’со не равен 0. Пусть В = РА$ХА\, 3 С {1,...,га — г} -номера столбцов матрицы В с нулевыми элементами в последней строке и V = 1т—г А-1 ~і~ т -^'1,7-4-! • Пусть } А$Р — (-<4 о>^о)> (а0, РА\)и = (ао,аьа2,... ,атп_г). Тогда системы
{А'0,аэ)хэ = Рсо, і =0,1,...,т- г, (4)
определены и их решения X] = ^ ^ )
Кг_1, є К задают базисное множество решений системы (1):
ТУ V 7 = 0,1,. ...т-г, (5)
где е'^- 6 Кт_Г+1 — э7до столбцы единичной матрицы /т_г+1 = (е'о, е'ь ..., е'т_г), <2 = сПаё(Р,1т_г), Ж = сИад(1г_ьи), V = С№.
Доказательство. Обозначим у = У~1Т~1х. По условию теоремы Р(Ао, А\)Уу = Рс0 и Р(Л0, А\)У = (^4'о, а0, «1,... ,ат-г).
Если искать решение в виде у =
(А >
то
придем к системам (4).
Покажем, что системы (4) определены и ^ ф
0. Для этого умножим их слева на РА$Х Р. Так как Р = Р-1, то получим F,AQlP(A,o,aj)xj = РАо'со, ^ = 0,1,... ,т - г.
Так как 6= РЛц'со и РА^1Р(А'о,ао) = /г, то первая из систем (4) принимает вид 1гхо =
6. Обозначим /г = (Р.е), е € Кг, тогда РАо1РА'0 = /', РАогРа0 = е.
Так как по определению РЛд 1 Р(«о,..., ат-г) = РА^1 Р(ао, РА\)и = (е, Р)ГУ, то обозначим (I] = РЛо 1 Р<4, 3 — 1, • • • , 771 — ?• столбцы матрицы (е, Р)1/. Так как С/ = /т_г+1 + 7 ^1,Я-1
и - номера столбцов матрицы В с нулевыми элементами в последней строке, то последние элементы всех столбцов dj матрицы (е, В) и не равны нулю. Системы (4) принимают вид
(/ , Ду^Ху 6, ] 0, 1, . . . ,771 7, (6)
Algorithm RationalBasis[A, с]
Input: A G RnXm,c 6 Rn Output: Решения системы Ax = с или NIL if с = 0 then
{B,S, Т, г} := DiagonalFormH(Л); Comment: В = (6Ь 62,..., 6m_r) Comment: /m_r = (ei, e2,..., em_r).
return (T( t1).........T( "et;))
else
{B) 6,5, T, P, r} := DiagonalForm(A, c); Comment: P = (6ls 62,...,6m_r),
Comment: 6 = ^ ^ ^, bj — y ^ j
Comment : Q = diag(P, Im-r), Comment: Im_r+i = (e'0, e'i,..., e'm_r). if r<rank(A,c) then return NIL; else J := 0 for j := 1 to m — r do
if /3j = 0 then J := JU {j}; for j := 1 to m — r do
if j*J
then := /3//3j, fj := e'j else £j := /3, fj := e'j + e'0
pe'o r
, j = l,...,m -r
return ( TQ
TQ
6' - ^b'j fj
Рис. 1. Алгоритм ’’Рациональный базис”
при этом det(//, ду) ф 0. Так как у вектора 6 последний элемент не равен нулю, то решения
систем Хп =
(S)
^ . . таковы, что ф 0.
V У? /
Поэтому векторы ^е'^, з = 1,...,ттг — г
линейно независимы, следовательно, и векторы (5) линейно независимы.
Следствие 2 Пусть В = (6Ь 62,..., 6т_г), ь = (^ )’ = ( & )• 0 = и = е'’ 3 ф. 3. ^ = /3, /у = е!у 4- е'о для 3 € J. Тогда
1=1,...,т-г,
- базисное множество решений неоднородной системы линейных уравнений (1).
Доказательство. В (5) подставим решения систем (6). При этом будем
А ( Ъ'о \
учитывать, что, по построению а^ = I ^ I,
6^ = 1 при ^ е 3 и = Рз при ^ ф 3. Затем умножим на матрицу ИЛ
Следствия теорем 1 и 2 позволяют записать алгоритм построения базисного множества решений системы линейных уравнений для произвольной коммутативной области R.
Алгоритм приведен на рис. 1. В этом алгоритме полагается следующее.
R - коммутативная область с единицей, в которой существует алгоритм сокращения (вычисляющий сомножитель 6, по известным а и с, а = Ьс,с ф 0), (1) - система линейных уравнений над R.
Diagonal Form(A, с) - алгоритм
диагонализации матрицы (Л, с) для неоднородной системы (1). Он вычисляет
{£?, 6, Т, S, Р, ?•}, при этом г = гапкЛ, а матрицы В, 6, Т, .S', Р описаны выше.
Diagonal FormH (А) - алгоритм
диагонализации матрицы Л, который
вычисляет {В,Т, S, г} для однородной системы (1).
Для этих алгоритмов можно взять любой из методов [3], [4], [5], [6] со сложность 0(п2т) или метод [7| со сложностью 0(п/3~1т).
РЕШЕНИЕ СИСТЕМЫ С ПОМОЩЬЮ ПОДЪЕМА ПО ГЕНЗЕЛЮ
Теоремы 1 и 2 сводят задачу построения базисного множества решений системы (1) к задаче решения нескольких определенных систем. Для вычисления базисного множества решений неоднородной системы необходимо решить т — г + 1 систем с матрицей коэффициентов размера г х г (для однородных систем нужно решить т — г систем). Для решения определенных систем можно воспользоваться 1>адическим подъемом.
Напомним общую схему р-адических методов. Выбирается подходящий простой элемент р кольца R. Для решения системы необходимо, чтобы элемент р не делил определитель матрицы коэффициентов. Например, в кольце Z этот выбор может оказаться неудачным с вероятностью не более 1 /р. Если проверка покажет, что решение неверное, то выбирается другой простой элемент р. Кольцо вычетов по простому модулю р является полем, и решение в этом поле можно искать, например, методом Гаусса. С помощью неравенства Адамара находится верхняя оценка для числителей и знаменателей решения системы, а по ним находится оценка для величины рк - границы подъема. Затем решение поднимается по модулю р до рк, и затем реконструируется рациональное решение. Для решения одной определенной системы используем алгоритм, приведенный
:= SAT;
Algorithm RationalBasisH[A]
Input: Л е RnXm Output: Решения системы Ах = 0. р := FirstPrime\ repeat
Initialization:
{Лц T, r} := Diagonal FormH \modp Л]; if m = r then return {0};
Ло A\ A 2 As (А1 о, a0) := Л0;
(®1 > • • • j r) • A i, j := 1; Test := TRUE\ repeat Lifting:
I := 1; x := 0; c* := — ay, repeat
x := modpAQ 1 c*; c* := (c* - A0x)/p\ x := x + lx\ I := lp\ until ||/|| > Ь?[А\; x' 4
RationalRec.onstruction[x, , D];
Test:
Comment: (eb ..., em_r) = Im-r ■■= ''■( ); if Axj Ф 0 then Test := FALSE; j :=j + 1; until j > m — r or Test = FALSE; p := iV ext Prime [pj; until Test = TRUE] return {^i,... ,xm_r}.
Рис. ‘2. Алгоритм ” Вычисление рационального базиса однородной системы линейных уравнений с помощью 1>адического подъема”
Диксоном [9] и использующий линейный i> адический подъем.
Алгоритмы построения базисного множества решений системы (1), использующие линейный р-адический подъем, приведены на рис. 2 и 3. Алгоритм для решения однородных систем приведен на рис.
2, алгоритм для решения неоднородных систем
— на рис. 3.
В алгоритмах используются следующие функции.
Функция || || имеет следующее значение:
||а|| = |а| для а 6 Z, ||а|| = dega для а € Р[а:], для матрицы Л ||Л|| = гпахо€л ||а||, для столбца к матрицы Л || afc || = max» ||ai,fc||-
Из неравенства Адамара следуют оценки для числителя и знаменателя решения системы.
Для кольца Z получим:
Algorithm RationalBasis[A, с]
Input: А € RnXm, с е Rn
Output: Решения системы Ах = с или NIL
р := FirstPrime] i := 1; repeat
Initialization:
{Ao\B,b,S,T,P,r,rc} :=
DiagonalForm\modp(A, c)|; if r < rc then return NIL] Comment: В = (6i, 62,..., bm-r)
f y. \
Comment : bj = f JJ Q = diag(P, Im_r) Comment : Im_r+i = (e'o.e'i,... ,e'm_r).
Ao ^1 \ .. ОЛФ. / co
A2 A3 J с := Pc0; (A'o, a0) := PA0P;
(ai, • • •, am_r) := PAi; J := 0; for j := 1 to m-r do if /3j = 0 then J := J U {j}; j := 0; Test = TRUE; repeat
if j € J then Aj := (A'0,aj + a0) else Aj:=(A/0,aj)
Lifting:
:= SAT; ( ^ ) := Sc;
A-* := I -W
-1
I := 1; x := 0; c* := c; repeat
x := modpAj *c*; c* := (c* - Ajx)/p; x := x + lx; 1 := lp; until ||1|| > Lj (A, cl x' 4
RationalReconstruction[x, Nj, Dj]; Test:
p
x'o
£0
if j = 0 then xo := T else if j € J
( p( XJ
then Xj := T I у £j
\ ^3e3
f p( X3 else Xj := TI \ 0
\ 0ei if Axj ф с then Test := FALSE;
j := j + i;
until j > m — r or Test = FALSE
i := i 4- 1; p := NextPrime[p];
until Test = TRUE or i = NumberOfPrimes
if Test = FALSE then return NIL
else return {xo,xi,..., xm_r}.
Рис. 3. Алгоритм ” Вычисление рационального базиса неоднородной системы линейных уравнений с помощью 1>-адического подъема”
7 = IlLi11М1» Л = minfc=i,...,n-i Ilafc||> = min(A, ||aj||),
Dj = nn/2/y\\aj || — оценка для знаменателя решения, Nj = nn/2\\ajЦтЦсЦ^”1 — оценка для числителя, Lj[A,c] = 2NjDj - оценка для границы р-адического подъема.
Для F[x] получим: Dj = \\aj\\ + llafcll
— оценка для знаменателя, Nj = Dj + ||с|| + minfc=ii...in-i,j ||afc|| — оценка для числителя, Lj\A,c\ = Nj + Dj — оценка для границы р-адического подъема.
NextPrime\p\ — генератор простых элементов кольца R, который возвращает следующий простой элемент, после р. FirstPrime и MinPrime обозначают первый простой элемент и наименьший по норме простой элемент.
N umber О / Primes — количество простых элементов [( max Dj)\\MinPrime\\~]J.
j=n,n-fl
Пусть R = R/pR фактор-кольцо R по модулю />R.
DiagonalForm\modv{ А, с)] алгоритм
диагонализации в кольце R. Результатом вычислений будет {Дц \ В, b, S, Т, Р, г, гс}. Здесь все обозначения имеют прежний смысл. Знак тильды обозначает матрицы над факто]> кольцом R. г и гс — это ранги матриц modpA и modp(A) с). Так как R является полем, то для вычислений можно воспользоваться методом Гаусса или любым другим методом решения систем линейных уравнений над полем.
Diagonal F or 7nH[modp А\ алгоритм диагонализации в R, вычисляет {Aq , S', Т, 7-}.
RationalReconstraction\x, N{, Di\ восстановление вектора решения в поле частных кольца R по решению х в поле R/p^R когда числитель и знаменатель ограничены числами Ni и D*- Алгоритм описан в работе [111-
РЕШЕНИЕ СИСТЕМЫ В КОММУТАТИВНОЙ ОБЛАСТИ
Теперь рассмотрим задачу решения системы линейных уравнений в коммутативной области. Отметим, что для однородных систем эта задача не представляет трудности, так как всякое решение в поле частных, после умножения на подходящий множитель, дает решение в области. Поэтому дальше будем рассматривать только неоднородные системы.
Пусть а — непустое конечное подмножество в области R. Пересечение главных
идеалов, порожденных элементами множества
а, является главным идеалом U реа(р). Порождающий элемент этого идеала обозначим НОК(а), т. е. (НОК(а)) = ир€а(р).
Пусть К - поле частных области її, х
- вектор из пространства Ктп, х С И -множество знаменателей компонент вектора х.
Определение. Знаменателем вектора х будем называть х = НОК(х).
Вектор х будем записывать в виде
произведения х
= хх х Е Кт, х ^ К.
Знаменатель вектора х будем обозначать ОЕМ(а-).
Пусть М. = {х\х Є Кт, Ах = 6}
множество всех решений системы (1) в Кт. Обозначим = М и Ыт -
множество решений системы, лежащих в
модуле Кт. Будем называть Мо множеством диофантовых решений.
Для х = {хих2,...,ха) и у = (г/і,у2, • • • ,У«) будем обозначать (х, у) = ХіУі.
Теорема 3 Пусть {ж» = ХіХЇ~1|г = 1,2,..., /г.}
- базисное множество решений неоднородной системы линейных уравнений (1) и пусть х = (хі,х2,...,хь), Х = (Хі,Х2,---,Хл)- Тогда
Л< = {|^>І9 Є ИЛ, <*,«)/ 0}
Доказательство: Пусть и = (иі, ..., и^) Є КЛ такой, что = 1. Обозначим Эг =
гцхГ » 9 ~ ПОК знаменателей всех чисел |г =
1,2,..., /і}, (]і Эг, (]ї Є І 1,2,..., /ї, (]
(<7і,92,- •• , Чк)-
Пусть х = (хі,х2,.. .,х^ и г — (х,и), следовательно, г Є М. Тогда получим
н н
г = (х,и) = ^ХіхГ1^* = ^Хі3і = (*’Ч)9~\
h
г= 1 /і
г= 1 /і
(7)
1 = щ = 5Zln* = XiSi = ^ 1 •
г=1
(8)
(*,<7)
г=1 »=1
Поэтому
(х,<7)
= Z Є М.
Наоборот, если </ є R/', (х,<?) = <7 Ф ^ и
г = , то, обозначая и* = Хг<7г<7 \ получим
2 = (а;, и) и Х!г=1 иг = 1- Поэтому г € Л4.
Следствие 3 Множество всех знаменателей {[)Е№(х)\Ах = с,х 6 К”1} решений системы Ах = с вместе с нулевым элементом образует идеал 1д в кольце К.
Следствие 4 Для того чтобы система (1) имела диофантовы решения, необходимо и достаточно, чтобы 1д = К.
-1
Следствие 5 Пусть идеал, порожденный знаменателями базисного множества решений
Xi = хгХГ*> * = 1,2, ...,/г системы (1), является единичным, х = (XI, Х2, • • •, Хн)-Тогда найдется такой вектор <7 =
(91, ®, • • •, Чь) € К\ что (х, д) = 1 и г = (х, <7)
- диофантово решение системы (1). Если при
этом (}§ ф- 0, то X] ^.. •, — 1) *^5+1) • ■ *, *^/1
базисное множество решений этой системы.
Следствие 6 Система (1) имеет диофантовы решения, если идеал,
порожденный знаменателями базисного множества решений, является единичным.
Следствие 7 Пусть идеал, порожденный знаменателями базисного множества решений Х{ = х»хгГ 1, г = 1,2, ...,/г системы (1), является единичным, х = (XI, Хг, • • •, Хл)* Тогда найдется такой вектор д =
(<7ь <72, • • • ,дл) € что (х, 9> = 1 и г = (х, щ)
- диофантово решение системы (1). Если при этом ф 0, то х\,...,гсв_ 1, г, жв+ь ..., хь -базисное множество решений этой системы.
Следствие 8 Пусть XI = ХгХгГ 1, г — 1, 2,..., /г базисное множество решений системы (1) и XI = 1- Тогда хь Х»-Хх(х4 - 1), г = 2, 3,...,/г
- линейно независимые диофантовы решения этой системы.
Доказательство: Так как все г* = х» -
Х1 (хг - 1), * = 2,3,..., /г принадлежат Ят и ЯВЛЯЮТСЯ ЛИНеЙНО НезаВИСИМЫМИ ВМеСТе С Х\, остается показать, что - решения системы (1). Так как XI = 1, то
Zi — Xj Xi(x< - 1) =
Х1 (~Хг + 1) + Xj • 1
Хі(“Хг + О + Хг * 1
Следовательно, по теореме, г = 2,3,..., /г, -решения системы Ах = 6.
Определение. Будем называть
диофантовым базисом решений системы Ах = с - базисное множество решений этой системы, которое целиком лежит в области R.
Другими словами - это т + I — г линейно независимых решений неоднородной (т — г -для однородной) системы, принадлежащих Rm.
Следствия последней теоремы позволяют записать алгоритм вычисления диофантова базиса решений системы. Алгоритм приведен на рис. 4.
В этом алгоритме используются следующие функции.
Функция Rational В asis[A, с] вычисляет базисное множество {xix71>x2X2 ■ • • >XhXh '}
Algorithm DiophantineBasisjA, с]
Input: А е RnXm,c € Rn
Output: Диофантовы решения или NIL
{/г,х1хГ1,х2Х2 :=
Rational В asis\A, с]; if h = 0 then return NIL;
(,gcd}ui,u2, ...,uh):=
ExdendedGCD[(xuX2, • • ■ ,Xh)]; if gcd ф 1 then return NIL else
A; := 1; while Uk = 0 do к := к +
zk ■= XiUi’
for г = 1,2,..., A; — 1, к +1..., h do 2* := -xk(xi - 1); return {z1,z2y...,zh}.
Рис. 4. Алгоритм ”Диофантов базис”.
решений системы Ах = с и количество h векторов в нем.
Функция ExdendedGCD[(xi,X2, • • •, Хл.)) по заданным h числам ХъХ2> • • • > Xh вычисляет порождающий элемент gcd идеала (хъ Х2, —>Хь)> этот идеал
является главным, а также множители t*i, и2>... }ин G R, такие, что gcd = Yli=i uiXi-Если этот идеал не является главным, то gcd принимает значение NIL.
Если результатом является NIL, то необходимо выполнять следующую итерацию. Для следующей итерации нужно выбрать случайным образом новые матрицы перестановок S и Т и вычислить новый рациональный базис. Алгоритм ExdendedGC D нужно применять к знаменателям всех рациональных решений, полученным ранее. Как доказано в [10], ожидаемое число случайных рациональных решений, которое необходимо для получения одного диофантова решения, будет O(logn + log log ||(А, с)||). Если гапкЛ = п, то за одну итерацию, будет получено rn — п + 1 рациональных решений. Поэтому ожидаемое количество итераций будет 0((m — п+ 1)— 1 (log п -f" log log || (Л, с) ||)).
ОЦЕНКА СЛОЖНОСТИ АЛГОРИТМОВ
Сложность вычисления диофантового базиса определяется количеством операций в алгоритмах 1 и 4. Сложность первого из них составляет 0(ттг7г^-1), а второго 0{рх{т — г) + Со) операций в кольце R. Со — количество операций, необходимых для разложения единицы по ттг — г образующим единичного идеала в кольце R. Для такого разложения единицы можно воспользоваться, например, алгоритмом вычисления базиса Гребнера. Оценка сложности такого алгоритма
в общем случае выходит за рамки настоящей статьи. Отметим только, что для евклидовых колец С с — это количество операций в алгоритме Евклида, который вычисляет наибольший общий делитель для m — г чисел.
Более определенные оценки можно получить для алгоритмов, использующих линейный р-адический подъем.
Известны оценки сложности алгоритма Диксона [10]. Для целых чисел Z сложность алгоритма Диксона для решения системы (1) при тг = ттг и det А ф 0 ограничена числом Vz = 0(n3(log7i + log||A||+logp)2 + 7i(log||c||)2) бит-операций. Для кольца многочленов F[x] над полем F сложность алгоритма Диксона ограничена числом 2V[x] =
0(Ti3(degj4 + degp)2 + n(degb)2) операций в поле F.
Для вычисления базиса, состоящего из ттг — г решений потребуется не более, чем Cz = (ттг - r)T>z и Ср[х] = (т - г)Т>р[х\ операций, соответственно для каждого из этих случаев.
Перейдем к оценке сложности вычисления диофантова базиса при использовании р-адического подъема.
Количесто операций при вычислении диофантова базиса определяется главным образом количеством операций необходимых для вычисления рационального базиса и количеством итераций. Ожидаемое количество итераций равно О ((ттг — тг + l)-1(logn + log log || (Л, с)К)), а количество операций будет СООТВеТСТВеННО Cz И Ср |ж] для Z и F[x] соответственно. Следовательно, при m — тг << п количество операций для вычисления диофантова базиса от размеров системы зависит существенно по кубическому закону.
По сравнению с алгоритмом [10], среднее число вычислений обратной матрицы теперь в 771 — г раз меньше.
ЗАКЛЮЧЕНИЕ
В заключение сравним по количеству операций приведенные алгоритмы со стандартными приближенными методами решения систем линейных уравнений при использовании чисел с плавающей точкой.
С помощью масштабирования числа из формата с плавающей точкой переводятся в целые числа. А затем решение ищется в поле частных. Такой перенос в область точных решений позволяет избавиться от проблемы неустойчивости решения, накопления погрешности и т.п. Вопрос состоит в том, как при этом увеличивается вычислительная сложность.
Если использовать арифметику с s-битными словами, то метод Гаусса для решения системы с матрицей размера п х т ранга п потребует 0{n2ms2) бит-операций.
Приведенный алгоритм построения рационального базиса решений требует для R = Z примерно 0((m — п + 1 )n2mw2) бит-операций. Где w = logn -f log ||Л|| + logp. Их отношение равно О ((га — п + 1 )(w/s)2. В реальных задачах logTi + log||A|| меньше s, а р обычно выбирается так, что logp меньше s. Поэтому w/s - величина меньше 1. А искомое отношение - это число 0(т — п + 1).
Отсюда видно, что при т — п « п алгоритм вычисления рационального решения имеет такую же сложность, что и метод Гаусса для приближенной арифметики при использовании чисел с плавающей точкой. Однако с ростом разности т — п отношение растет линейно, как т — п.
ЛИТЕРАТУРА
1. Strassen. V. Gaussian Elimination is not optimal // Numerische Mathematik. 1969. V. 13. P. 354-356.
2. Coppersmith D. and Winograd S. Matrix multiplication via arithmetic progressions // Journal of Symbolic Computation. 1990. V.
9. P. 251-280.
3. Dodgson C.L. Condensation of determinants, being a new and brief method for computing their arithmetic values // Proc. Royal Soc. Lond. 1866. V. A. 15. P. 150-155.
4. Bareiss Е.И. Sylvester’s Identity and Multistep Integer-Preserving Gaussian Elimination // Math. Comp. 1968. V. 22 (103). P. 565-578.
5. Malaschonok G.I. Solution of a system of linear equations in an intergal domain // USSR J. of computational Mathematics and Mathematical Physics. 1983. V. 23. P. 1497-1500.
6. Malaschonok G.I. Algorithms for the solution of systems of linear equations in commutative rings // Effective Methods in Algebraic Geometry / Ed. by T. Mora and C. Traverso, Progress in Mathematics 94, Birkhauser, 1991. P. 289-298.
7. Malaschonok G.I. Recursive Method for the Solution of systems of Linear Equations // Computational Mathematics ( Sydow A. Ed, Proceedings of the 15th IMACS World Congress, Vol. I, Berlin, August 1997 ). Wissenschaft & Technik Verlag. Berlin, 1997. P. 475-480.
8. Bareiss E.H. Computational solutions of matrix problems over an integral domain // J. Inst. Maths Applies. 1972. V. 10. P. 68-104.
9. Dixon J. Exact solution of linear equations using p-adic expansions // Numer. Math. 1982. V. 40. P. 137-141.
10. Mxdders T. and Storjohann A. Diophantine Linear System Solving // Proceedings of IS-SAC’99: ACM International Symposium on Symbolic and Algebraic Computation. July 1999, Vancouver, Canada.
11. Collins G.E. and Encarnacion M.J. Efficient rational number reconstruction. Journal of Symbolic Computation. 1995. V. 20. P. 287-297.
Поступила в редакцию 1 февраля 2000 г.
