Programming: theory and practice by alexBlack

Содержание

1 Представление и операции с полиномами
2 Порядок следования членов в полиноме
3 S-полином
4 Редукция полинома по базису Гребнера
5 Алгоритм Бухбергера
6 В заключение

В статье описывается алгоритм Бухбергера построения редуцированного базиса Гребнера. В частности, базис Гребнера может быть применен для решения системы полиномиальных уравнений от нескольких переменных.

Я пытался применить этот алгоритм для получения точек пересечения кривых Безье, но позже использовал более простой способ, основанный на результантах. Эта статья, возможно, будет интересна только для первого знакомства с алгоритмом. Я не математик, и не претендую на полноту изложения. В статье не будет сложных математических выкладок. Все термины, которые необходимы для понимания статьи, будут объясняться по мере их появления в тексте и их не будет очень много. Более обстоятельно с алгоритмом можно познакомиться по книге Кокса^[1] или почитать самого Бухбергера^[2]. Кроме того, существуют еще как минимум два алгоритма для нахождения базиса Гребнера (они называются F4 и F5), но в них я не вникал.

Здесь же мы рассмотрим последовательность действий сразу с реализацией алгоритма Бухбергера и начнем с того, что договоримся как мы будем представлять основные структуры.

Представление и операции с полиномами

pascal

type
   TPolynom = record
      // Матрица коэффициентов
      list : array [0..mxst,0..mxst] of Double;   
   end;

   TMonom = record
      k : Double;            // Коэффициент
      x, y : integer;        // Степени переменных
   end;

procedure TPolynom.Add(m:TMonom);
begin
   list[m.x, m.y] := list[m.x, m.y]+m.k;
end;

	0	1	2	3
0	7	6	5	4
1	3	0	0	0
2	2	0	0	0
3	1	0	0	0

В общем случае алгоритм Бухбергера может быть применен для полиномов от произвольного числа переменных. Меня же интересовала система полиномов вида

  A₁x³ + A₂x² + A₃x + A₄y³ + A₅y² + A₆y + A₇ = 0 
  B₁x³ + B₂x² + B₃x + B₄y³ + B₅y² + B₆y + B₇ = 0

Такая система получается при поиске точек пересечения кубических кривых Безье. В этом случае (всего две переменных) отдельный полином удобно представлять в виде двумерного массива его коэффициентов TPolynom. Коэффициенты в массиве хранятся так, что индесы i и j массива - это степени переменных x и y в членах полинома. Например, полином

  x³ + 2x² + 3x + 4y³ + 5y² + 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;

Важно с самого начала определиться с порядком следования членов в полиноме. Дело в том, что далее нам понадобится сравнивать старшие члены двух полиномов, а это сравнение возможно только при заданном порядке следования. Кстати, от выбранного порядка следования зависит получаемый базис Гребнера. Выберем лесикографический порядок следования членов полинома (lex). При таком порядке сначала сравниваются имена переменных, затем их степени. (Для примера все полиномы выше написаны именно в lex-порядке). Теперь мы можем написать процедуру, которая будет возвращать старший член полинома - она ищет в нашей матрице первый ненулевой коэффициент перебирая матрицу снизу вверх и в каждой строке справа налево.

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;

S-полином от двух полиномов f и g (обозначается как S(f, g)) специально "сконструирован" для сокращения старших членов двух полиномов. Не будем здесь приводить формулу для S-полинома, а просто покажем на примере, тем более что для объяснения формулы потребуется несколько страниц, а смысл очень простой. Пусть, например, есть два полинома:

  f = x³y² - x²y³ + x 
  g = 3x⁴y + y²

Умножим первый из них на 3x , а второй на y , получим

  f₁ = 3x⁴y² - 3x³y³ + 3x² 
  g₁ = 3x⁴y² + y³

Теперь вычтем из первого полинома второй и получим

  f₁ - g₁ = - 3x³y³ + 3x² - y³

Это и есть 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;

Настала пора определиться, что же такое базис Гребнера. Грубо говоря (и нам достаточно такого определения) это множество многочленов. Будем представлять его как простой массив полиномов list и их количество count в списке ( TGBasis ).

Еще одно определение - отношение делимости двух членов полинома. Член m полинома делит член k полинома, когда степени всех переменных члена m меньше или равны соответствующих степеней члена k . Например, 2xy делит 3x²y (получаем 3x/2 ), а 3x²y не делит 2xy .

Слева приведена процедура divide деления члена полинома на член m . Мы получаем частное и true, если член m делит текущий член.

Перейдем собственно к редукции. Редукция полинома по базису - это последовательное вычитание старших членов полинома, которые делятся на старшие члены полиномов, входящих в базис. Процедура Reduce выполняет редукцию полинома F по текущему базису. Алгортим следующий:

В начале процедуры результат R - полином, не содержащий ни одного члена.

Для каждого полинома из базиса F проверяем делит ли его старший член ( F.LT ) старший член ( P.LT ) исходного полинома P.

Если делит, полином F умножаем на частное от этого деления ( m ) и вычитаем из полинома P . Цикл в этом случае начинаем снова с первого полинома из базиса.

Если не делит, берем следующий полином из базиса и выполняем деление для него.

Когда мы прошли все полиномы базиса и деления не было, старший член исходного полинома P переносим в результат R . С новым полиномом P повторяем всю процедуру. Действия выполняются до тех пор, пока полином P содержит хотя бы один член.

Если после редукции получается полином, не содержащий ни одного члена, говорят, что полином редуцируется к нулю по базису.

Посмотрим выполнение алгоритма на примере. Пусть базис состоит из двух полиномов:

  y² - 1 
  xy - 1

Выполним редукцию полинома x²y + xy² + y² по этому базису. Берем старший член первого полинома y² . Он не делит старший член исходного полинома, поэтому переходим ко второму полиному. Старший член второго полинома xy делит старший член исходного. Выполняем умножение на частное от этого деления и вычитание:

  x²y + xy² + y² - x(xy - 1) = xy² + x + y²

Цикл начинаем сначала (с первого полинома). Теперь его старший член делит старший член xy² . Выполняем еще одну операцию:

  xy² + x + y² - x(y² - 1) = 2x + y²

Теперь старший член 2x не делится ни одним старшим членом из базиса. Переносим старший член 2x в результат, а с остатком y² повторяем деление. Он делится на старший член первого полинома из базиса. Выполняем операцию:

  y² - (y² - 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).

Для этих кривых получаем систему уравнений:

  - 6x² + 6x - 6y² + 6y - 2 
  4x³ - 6x² - 4y³ + 6y²

Получим для этих полиномов редуцированный базис Гребнера. Приведем здесь первые четыре полинома из базиса:

  24xy² - 24xy + 20x + 24y³ - 48y² + 12y - 4 
  128xy - 64x + 192y⁴ - 384y³ + 256y² - 128y + 26.67 
  -298.67x + 768y⁶ - 1536y⁵ + 768y⁴ - 256y³ + 704y² - 213.33y + 78.22 
  3072y⁸ - 6144y⁷ + 8192y⁶ - 15360y⁵ + 19968y⁴ - 11776y³ + 5006.22y² - 2588.44y + 400.59

Видно как из полиномов постепенно вычитаются члены, содержащие переменную x , в чем, собственно, и состоит смысл алгоритма Бухбергера.

Последний полином 8-ой степени содержит только одну переменную y . Можно приравнять его нулю и решить полученное уравнение. В результате получим два значения y

  y₁ = 0.788675 
  y₂ = 0.211325

Подставив эти значения в первое уравнение нашей системы, получим значения x .

  x₁ = 0.788675 
  x₂ = 0.211325

Вычисленные значения нужно поставить во второе уравнение, чтобы проверить и отбросить ненужные значения. Эти значения - параметры наших кривых Безье в точках пересечения.

В заключение

Поскольку мне нужно было получить точки пересечения кривых Безье, мне показалось более простым вывести один раз полином девятого порядка, чем каждый раз получать и исследовать базис. Поэтому я остановился на методе с результантами. Некоторое время я раздумавал стоит ли публиковать этот текст, но все-же решил это сделать. Для первого знакомства он вполне сгодится. А уже потом можно почитать Кокс'а.

1↑Д.Кокс, Дж.Литтл, Д.О'Ши Идеалы, многообразия и алгоритмы.

2↑Б.Бухбергер, Ж.Калме, Э.Калтофен. Компьютерная алгебра. Символьные и алгебраические вычисления.

3↑А.Ю.Рысцова Интеллектуализация программных средств решения систем нелинейных уравнений.

09 марта 2011

Downloads

Объявление и реализация классов (p9.zip) (4.7 Кб, просмотров: 1353 )

Comments

15.09.2012 03:07 Tink

 Доброе время суток, alexBlack. Большое спасибо за статью. Сам начал недавно интересоваться
 кривыми Безье. Правда пока мало опыта. Решил переписать твой код на С++ и поэкспереминтировать. 
 Но выяснилось, что нет функций умножения и вычитания (mul и sub), также функция sub
 перегруженна, так как принимает ввиди аргумента то объект m класса TMonom, то объект p2 класса TPolynom.

15.09.2012 03:08 Tink

 Реализация для объекта m я понимаю аналогична функции Add. Также не понятно , что содержат
 переменные x, y, k в функции TMonom.divide. Какое они должны содержать значение? Где их 
 объявление? И инициализация?
 
 P.S. Большое спасибо за статью и прошу ответь на мои вопросы, пожалуйста. Удачи в твоих делах!!!

01.10.2012 11:06 alexBlack

 Добавил исходники

22.12.2012 06:37 АВ

 очень хотелось посмтреть как будет работать,запускаю исходник,
 а в результате куча ошибок. почему?

22.12.2012 10:42 alexBlack

 :) Исходник только для того, чтобы уточнить объявления упоминаемых в статье классов. 
 Считаю, что этого будет достаточно чтобы разобраться.

Вы можете оставить комментарий или задать вопрос