| Главная | Эксперименты | Утилиты | Компоненты | Что почитать | Контакты |
| ||||||||
|
Я пытался применить этот алгоритм для получения точек пересечения кривых Безье, но позже использовал более простой способ, основанный на результантах. Эта статья, возможно, будет интересна только для первого знакомства с алгоритмом. Я не математик, и не претендую на полноту изложения. В статье не будет сложных математических выкладок. Все термины, которые необходимы для понимания статьи, будут объясняться по мере их появления в тексте и их не будет очень много. Более обстоятельно с алгоритмом можно познакомиться по книге Кокса[1] или почитать самого Бухбергера[2]. Кроме того, существуют еще как минимум два алгоритма для нахождения базиса Гребнера (они называются F4 и F5), но в них я не вникал.
Здесь же мы рассмотрим последовательность действий сразу с реализацией алгоритма Бухбергера и начнем с того, что договоримся как мы будем представлять основные структуры.
Представление и операции с полиномами
|
|
A1x3 + A2x2 + A3x + A4y3 + A5y2 + A6y + A7 = 0 B1x3 + B2x2 + B3x + B4y3 + B5y2 + B6y + B7 = 0Такая система получается при поиске точек пересечения кубических кривых Безье. В этом случае (всего две переменных) отдельный полином удобно представлять в виде двумерного массива его коэффициентов TPolynom. Коэффициенты в массиве хранятся так, что индесы i и j массива - это степени переменных x и y в членах полинома. Например, полином
x3 + 2x2 + 3x + 4y3 + 5y2 + 6y + 7будет храниться как показано в таблице слева. Такое представление полинома удобно тем, что умножение полинома на x или y в этом случае равносильно сдвигу массива вниз или вправо.
С другой стороны, нам понадобится отдельный член полинома, который мы будем представлять записью TMonom, содержащей коэффициент и степени при переменных x и y .
При таком представлении сложение полинома и отдельного члена (также как и вычитание) сводится к сложению коэффициентов. Немного сложнее реализуется умножение полинома на отдельный член.
Порядок следования членов в полиноме
pascal
function TPolynom.LT : TMonom; var i, j:integer; begin // lex for i := mxst downto 0 do begin for j := mxst downto 0 do begin if (abs(list[i,j]) < eps) then continue; result.k := list[i, j]; result.x := i; result.y := j; exit; end; end;
S-полином
pascal
function SPolynom(f1, f2:TPolynom):TPolynom; var m1, m2, k1, k2 : TMonom; p2 : TPolynom; begin m1 := f1.LT; m2 := f2.LT; k1.k := m2.k; k1.x := max(m2.x-m1.x, 0); k1.y := max(m2.y-m1.y, 0); k2.k := m1.k; k2.x := max(m1.x-m2.x, 0); k2.y := max(m1.y-m2.y, 0); p2 := f2; p2.mul(k2); result := f1; result.mul(k1); result.sub(p2); end;
f = x3y2 - x2y3 + x g = 3x4y + y2Умножим первый из них на 3x , а второй на y , получим
f1 = 3x4y2 - 3x3y3 + 3x2 g1 = 3x4y2 + y3Теперь вычтем из первого полинома второй и получим
f1 - g1 = - 3x3y3 + 3x2 - y3Это и есть S-полином от полиномов f и g . Просто, не правда ли ? Похоже на то как мы сокращаем члены уравнений в методе Гаусса для систем линейный уравнений (кстати, этот метод является частным случаем алгоритма Бухбергера).
Еще раз. Мы подбираем коэффициенты для полиномов так, чтобы при вычитании сократились их старшие члены. (А они, как мы помним зависят от выбранного порядка следования). Отсюда совсем простая процедура нахождения S-полинома.
Редукция полинома по базису Гребнера
pascal
type TGBasis = record list : array [0..maxPolynoms] of TPolynom; count : integer; end; function TMonom.divide(m:TMonom):boolean; begin result := (x >= m.x) and (y >= m.y); k := k/m.k; x := x-m.x; y := y-m.y; end; function TGBasis.Reduce(F:TPolynom):TPolynom; var R, P:TPolynom; divided : boolean; i : integer; m : TMonom; begin P := F; R := TPolynom.Create; while not P.LT.zero do begin divided := false; i := 0; while (i <= Count - 1) and not divided do begin F := list[i]; m := P.LT; if not m.divide(F.LT) then begin inc(i); continue; end; F.mul(m); P.sub(F); divided := true; end; if not divided then begin m := P.LT; R.add(m); P.sub(m); end; end; result := R; end;
Еще одно определение - отношение делимости двух членов полинома. Член m полинома делит член k полинома, когда степени всех переменных члена m меньше или равны соответствующих степеней члена k . Например, 2xy делит 3x2y (получаем 3x/2 ), а 3x2y не делит 2xy .
Слева приведена процедура divide деления члена полинома на член m . Мы получаем частное и true, если член m делит текущий член.
Перейдем собственно к редукции. Редукция полинома по базису - это последовательное вычитание старших членов полинома, которые делятся на старшие члены полиномов, входящих в базис. Процедура Reduce выполняет редукцию полинома F по текущему базису. Алгортим следующий:
В начале процедуры результат R - полином, не содержащий ни одного члена.
Для каждого полинома из базиса F проверяем делит ли его старший член ( F.LT ) старший член ( P.LT ) исходного полинома P.
Если делит, полином F умножаем на частное от этого деления ( m ) и вычитаем из полинома P . Цикл в этом случае начинаем снова с первого полинома из базиса.
Если не делит, берем следующий полином из базиса и выполняем деление для него.
Когда мы прошли все полиномы базиса и деления не было, старший член исходного полинома P переносим в результат R . С новым полиномом P повторяем всю процедуру. Действия выполняются до тех пор, пока полином P содержит хотя бы один член.
Если после редукции получается полином, не содержащий ни одного члена, говорят, что полином редуцируется к нулю по базису.
Посмотрим выполнение алгоритма на примере. Пусть базис состоит из двух полиномов:
y2 - 1 xy - 1Выполним редукцию полинома x2y + xy2 + y2 по этому базису. Берем старший член первого полинома y2 . Он не делит старший член исходного полинома, поэтому переходим ко второму полиному. Старший член второго полинома xy делит старший член исходного. Выполняем умножение на частное от этого деления и вычитание:
x2y + xy2 + y2 - x(xy - 1) = xy2 + x + y2Цикл начинаем сначала (с первого полинома). Теперь его старший член делит старший член xy2 . Выполняем еще одну операцию:
xy2 + x + y2 - x(y2 - 1) = 2x + y2Теперь старший член 2x не делится ни одним старшим членом из базиса. Переносим старший член 2x в результат, а с остатком y2 повторяем деление. Он делится на старший член первого полинома из базиса. Выполняем операцию:
y2 - (y2 - 1) = 1Старший член 1 снова не делится ни одним старшим членом базиса. Переносим ее в результат и, т.к. остатка нет, завершаем выполнение. Итак, в результате редукции получили полином
2x + 1
Алгоритм Бухбергера
pascal
class function TGBasis.Create(f1, f2:TPolynom):TGBasis; var G:TGBasis; k, i : integer; s : TPolynom; begin G := TGBasis.Create; G.NewPolynom(f1); G.NewPolynom(f2); k := 1; while (k < G.count) and (G.count <= maxPolynoms) do begin f1 := G.list[k]; i := 0; while i <= G.count-1 do begin if i = k then begin inc(i); continue; end; f2 := G.list[i]; s := G.Reduce(SPolynom(f1, f2)); if not s.isConst then begin G.NewPolynom(s); if G.count > maxPolynoms then break; end; inc(i); end; inc(k); end; result := G; end;
Начальное состояние - базис из двух полиномов нашей системы уравнений ( f1, f2). Назовем его текущим базисом.
Для каждой пары из текущего базиса получаем S-полином. Если этот S-полином не редуцируется к нулю на текущем базисе, добавляем остаток от редукции в базис. Вычисления заканчиваются когда S-полином любой пары полиномов из текущего базиса редуцируется к нулю по этому базису.
Не обязательно каждый раз проверять S-полиномы всех возможных пар из базиса. Достаточно проверить S-полиномы каждого добавляемого в базис полинома со всеми остальными, т.к. остальные пары уже проверены и результат не изменится.
Возьмем, к примеру, две кривые с опорными точками (0,2), (2,2), (2,0), (0,0) и (2,2), (0,2), (0,0), (2,0).
|
- 6x2 + 6x - 6y2 + 6y - 2 4x3 - 6x2 - 4y3 + 6y2Получим для этих полиномов редуцированный базис Гребнера. Приведем здесь первые четыре полинома из базиса:
24xy2 - 24xy + 20x + 24y3 - 48y2 + 12y - 4 128xy - 64x + 192y4 - 384y3 + 256y2 - 128y + 26.67 -298.67x + 768y6 - 1536y5 + 768y4 - 256y3 + 704y2 - 213.33y + 78.22 3072y8 - 6144y7 + 8192y6 - 15360y5 + 19968y4 - 11776y3 + 5006.22y2 - 2588.44y + 400.59Видно как из полиномов постепенно вычитаются члены, содержащие переменную x , в чем, собственно, и состоит смысл алгоритма Бухбергера.
Последний полином 8-ой степени содержит только одну переменную y . Можно приравнять его нулю и решить полученное уравнение. В результате получим два значения y
y1 = 0.788675 y2 = 0.211325Подставив эти значения в первое уравнение нашей системы, получим значения x .
x1 = 0.788675 x2 = 0.211325Вычисленные значения нужно поставить во второе уравнение, чтобы проверить и отбросить ненужные значения. Эти значения - параметры наших кривых Безье в точках пересечения.
В заключение
Поскольку мне нужно было получить точки пересечения кривых Безье, мне показалось более простым вывести один раз полином девятого порядка, чем каждый раз получать и исследовать базис. Поэтому я остановился на методе с результантами. Некоторое время я раздумавал стоит ли публиковать этот текст, но все-же решил это сделать. Для первого знакомства он вполне сгодится. А уже потом можно почитать Кокс'а.1↑Д.Кокс, Дж.Литтл, Д.О'Ши Идеалы, многообразия и алгоритмы.
2↑Б.Бухбергер, Ж.Калме, Э.Калтофен. Компьютерная алгебра. Символьные и алгебраические вычисления.
3↑А.Ю.Рысцова Интеллектуализация программных средств решения систем нелинейных уравнений.
|