Многочлены с рациональными коэффициентами

по материалам курсовой работы О. Пиценко

Хорошо известно, что если рациональное число p/q является корнем уравнения с целыми коэффициентами a_nxⁿ+ a_n-1x^n-1+ ... a₁x+ a₀=0, то a₀ делится на p и a_n делится на q. Это свойство дает возможность по коэффициентам a₀ и a_n данного уравнения найти все рациональные числа, которые могут быть корнями этого уравнения. Подставляя затем все эти числа в уравнение, можно установить, какие из них действительно являются его корнями. При этом переборе обычно используют еще одно необходимое условие: для того, чтобы рациональное число p/q являлось корнем уравнения f(x)=0 с целыми коэффициентами, необходимо, чтобы значение функции f(1) делилось на p-q, а значение f(-1) делилось на p+q.

Имея уравнение с рациональными коэффициентами, легко получить эквивалентное ему уравнение с целыми коэффициентами. Следующая функция использует перечисленные выше факты для вычисления рациональных корней многочлена с рациональными коэффициентами, при этом выводя на печать необходимые промежуточные вычисления.

RationalRootsOfPolynomial:=function(f)
local x, a0, diva0, an, g, divan, roots, p, q, x0, val, a, i, c, k,
      fplus1, fminus1, t, d;
Print("f(x) = ", f,"\n");
x:=IndeterminateOfUnivariateRationalFunction(f);
# Создаем пустой список для хранения в нем корней многочлена
roots:=[]; 
# Oпределяем коэффициенты 
a:=PolynomialCoefficientsOfPolynomial(f,x);
# Переводим коэффициенты в числа 
a:=List(a, i -> Value(i,0) ); 
a:=List(a,
i -> DenominatorRat(i) );
# Находим общий делитель знаменателей всех коэффициентов
a:=Lcm(a); 
# Домножаем многочлен на общий знаменатель
f:=f*a; 
if a>1 then
  Print("f(x) = 1/", a, " * (", f, ")\n");
fi;
a0:=Value(f,0); 
# Находим свободный член a0 и проверяем его равенство нулю. 
# Если a0=0, то делим многочлен на x и добавляем к списку roots
# корень уравнения 0 столько раз, какова его кратность
if a0=0 then
  k:=PositionProperty( PolynomialCoefficientsOfPolynomial(f,x), 
                       c -> not IsZero(c) ) - 1;
  f:=Quotient( f, x^k );
  AddSet(roots, 0);
  a0:=Value(f,0);
  if a=1 then
    Print("f(x) = x^", k, " * (", f,")\n");
  else
    Print("f(x) = 1/", a, " * x^", k, " * (", f,")\n");
  fi;
fi;
Print("a_0 = ", a0, "\n");
# Находим делители a0
diva0:=DivisorsInt(a0); 
Print("diva0 = ", diva0, "\n");
an:=LeadingCoefficient(f);
Print("a_n = ", an, "\n");
# Находим делители an
divan:=DivisorsInt(an); 
Print("divan = ", divan, "\n");
# Находим значение многочлена в точкe 1
fplus1:=Value(f,1); 
# Находим значение многочлена в точкe -1
fminus1:=Value(f,-1); 
Print(" p      q    (p-q)|f(1)      (p+q)|f(-1)    f(p/q)", "\n");
# Проверяем делимость f(1) на (p-q) и f(-1) на (p+q). В случае положительного 
# результата подстановкой проверяем, является ли число p/q корнем данного многочлена.
# Если число p/q является корнем многочлена, то добавляем его в список roots
for t in diva0 do
  for q in divan do
    for p in [t, -t] do
      x0:=p/q;
      Print(" ",p,"      ",q, "\c");
      d:=p-q;
      if d<>0 then
        if (fplus1 mod (p-q)) = 0 then
          Print("       +       ", "\c");
        else
          Print("       -       ", "\n");
          continue;
        fi;
      else
        Print( " ", "\c");
      fi;
      d:=p+q;
      if p+q<>0 then
        if (fminus1 mod (p+q)) = 0 then
          Print("       +       ", "\c");
        else
          Print("       -       ", "\n");
          continue;
        fi;
      else
        Print( " ", "\c");
      fi;
      val:=Value(f,x0);
      Print(" f(", x0, ") = ", val, "\n");
      if val=0 then
        AddSet(roots, x0);
        # Найдя корень многочлена x0, делим исходный многочлен на (x-x0)
        f:=Quotient(f,x-x0); 
      fi;
    od;
  od;
od;
return roots;
end;

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

Пример 1.

gap> Q:=Rationals;
Rationals
gap> x:=Indeterminate(Q,"x");
x
gap> f:=(x^2-4/9)*x^2;
x^4-4/9*x^2
gap> RationalRootsOfPolynomial(f);
f(x) = x^4-4/9*x^2
f(x) = 1/9 * (9*x^4-4*x^2)
f(x) = 1/9 * x^2 * (9*x^2-4)
a_0 = -4
diva0 = [ 1, 2, 4 ]
a_n = 9
divan = [ 1, 3, 9 ]
 p      q    (p-q)|f(1)      (p+q)|f(-1)    f(p/q)
 1      1        -       
 -1      1       -       
 1      3       -       

-1      3       -       
 1      9       -       
 -1      9       -       
 2      1       +              -       
 -2      1       -       
 2      3       +              +        f(2/3) = 0
 -2      3       +              +        f(-2/3) = 0
 2      9       -       
 -2      9       -       
 4      1       -       
 -4      1       +              -       
 4      3       +              -       
 -4      3       -       
 4      9       + 
            -       
 -4      9       -       
[ -2/3, 0, 2/3 ]
gap>

Пример 2.

gap> f:=(x^2-1/4)*(x+5/3);
x^3+5/3*x^2-1/4*x-5/12
gap> RationalRootsOfPolynomial(f);
f(x) = x^3+5/3*x^2-1/4*x-5/12
f(x) = 1/12 * (12*x^3+20*x^2-3*x-5)
a_0 = -5
diva0 = [ 1, 5 ]
a_n = 12
divan = [ 1, 2, 3, 4, 6, 12 ]
 p      q    (p-q)|f(1)      (p+q)|f(-1)    f(p/q)
 1      1        +        f(1) = 24
 -1      1       +         f(-1) = 6
 1      2       +              +        f(1/2) = 0
 -1      2       +              +        f(-1/2)
= 0
 1      3       +              -       
 -1      3       +              +        f(-1/3) = 16
 1      4       +              -       
 -1      4       -       
 1      6       -       
 -1      6       -       
 1      12       -       
 -1      12       -       
 5      1       +              +        f(5) = 80
 -5      1       +              -       
 5      2       +              -       
 -5      2       -       
 5      3       +     
        -       
 -5      3       +              +        f(-5/3) = 0
 5      4       +              -       
 -5      4       -       
 5      6       +              -       
 -5      6       -       
 5      12       -       
 -5      12       -       
[ -5/3, -1/2, 1/2 ]
gap>

Созданную нами функцию RationalRootsOfPolynomial удобно использовать для проверки заданий на нахождение рациональных корней многочленов в курсе алгебры и теории чисел, однако она не является универсальной. В частности, эта функция не учитывает кратность корней. Поэтому мы разработаем еще одну функцию MyRootsOfUPol, которая будет возвращать рациональные корни многочленов с учетом их кратности:

MyRootsOfUPol:=function(f) 
local x, a0, diva0, an, g, divan, roots, p, q, x0, val, a, i, c, k, 
      fplus1, fminus1, t, d, var; 
if IsZero(f) then  
  return [0]; 
elif DegreeOfUnivariateLaurentPolynomial(f)=0 then 
  return []; 
fi; 
x:=IndeterminateOfUnivariateRationalFunction(f); 
roots:=[]; 
a:=PolynomialCoefficientsOfPolynomial(f,x); 
a:=List(a, i -> Value(i,0) ); 
a:=List(a, i -> DenominatorRat(i) ); 
a:=Lcm(a);  
if a<>1
then 
  f:=f*a; 
fi;  
a0:=Value(f,0); 
if a0=0 then 
  k:=PositionProperty( PolynomialCoefficientsOfPolynomial(f,x), 
                       c -> not IsZero(c) ) - 1; 
  f:=Quotient( f, x^k ); 
  Append(roots, List([1..k], i -> 0) ); 
  a0:=Value(f,0); 
fi; 
diva0:=DivisorsInt(a0); 
an:=LeadingCoefficient(f); 
divan:=DivisorsInt(an); 
fplus1:=Value(f,1); 
fminus1:=Value(f,-1); 
var:=[]; 
for t in diva0 do 
  for
q in divan do 
    for p in [t, -t] do 
      AddSet(var, p/q); 
    od; 
  od; 
od; 
for x0 in var do 
  p:=NumeratorRat(x0); 
  q:=DenominatorRat(x0); 
  d:=p-q; 
  if d<>0 then  
    if (fplus1 mod (p-q)) <> 0 then 
      continue; 
    fi; 
  fi; 
  d:=p+q; 
  if p+q<>0 then  
    if (fminus1 mod (p+q)) <> 0 then 
      continue; 
    fi; 
  fi; 
  val:=Value(f,x0); 
  if
val=0 then  
    Add(roots, x0); 
    if DegreeOfUnivariateLaurentPolynomial(f)>0 then 
      f:=Quotient(f,x-x0); 
      Append(roots, MyRootsOfUPol(f)); 
      return roots; 
    fi;  
  fi;  
od; 
return roots; 
end; 

Алгоритм данной функции практически совпадает с предыдущей, однако она работает быстрее, так как не тратит время на вывод промежуточных вычислений. Интересно сравнить производительность функции MyRootsOfUPol со стандартной функцией системы GAP RootsOfUPol, а также проверить правильность расчетов сравнением результатов, возвращаемых обеими функциями. Зададим многочлен десятой степени, случайным образом сформировав его коэффициенты, и найдем его рациональные корни двумя способами:

gap> f:=Product(List([1..10],i->x-Random(Rationals)));
x^10-17/3*x^9+43/16*x^8+247/16*x^7+197/64*x^6-279/64*x^5-77/64*x^4+41/192*x^3+1/16*x^2
gap> RootsOfUPol(f);               
[ 4, 3, 1/2, 1/4, 0, 0, -1/4, -1/3, -1/2, -1 ]
gap> time;
283
gap> MyRootsOfUPol(f);
[ 0, 0, -1, -1/2, -1/3, -1/4, 1/4, 1/2, 3, 4 ]
gap> time;
9
gap> SortedList(MyRootsOfUPol(f))=SortedList(RootsOfUPol(f));
true
gap>

Из приведенного выше примера видно, что результаты расчетов совпадают. Как показывает этот и аналогичные примеры, при небольшой степени многочлена функция MyRootsOfUPol находит рациональные корни многочлена быстрее, чем стандартная функция RootsOfUPol. Возьмем теперь более сложный многочлен, например, многочлен 30-й степени и сравним результаты и время работы обеих функций:

gap> f:=Product(List([1..30],i->x-Random(Rationals)));
x^30+1051/60*x^29+141331/1200*x^28+3082613/8640*x^27+58522771/259200*x^26-436162163/259200*x^2\
5-35773349/7200*x^24-7211778329/2073600*x^23+33235831763/4147200*x^22+125514755/6912*x^21+7199\
839177/1036800*x^20-9780069169/518400*x^19-34659602443/1382400*x^18-3695473133/2073600*x^17+98\
3525023/51840*x^16+9485611669/691200*x^15-9502389071/4147200*x^14-7894762583/1036800*x^13-1052\
660579/345600*x^12+190257823/207360*x^11+4904751197/4147200*x^10+216549797/691200*x^9-1813583/\
28800*x^8-595997/9600*x^7-54551/3200*x^6-3513/1600*x^5-9/80*x^4
gap> RootsOfUPol(f);
[ 2, 1, 1, 1, 2/3, 2/3, 3/5, 0, 0, 0, 0, -1/5, -1/4, -1/2, -1/2, -1/2, -1/2, -3/4, -3/4, -1, 
  -1, -1, -1, -1, -4/3, -3/2, -5/3, -2, -3, -6 ]
gap> time;
102
gap> MyRootsOfUPol(f);
[ 0, 0, 0, 0, -6, -3, -2, -5/3, -3/2, -4/3, -1, -1, -1, -1, -1, -3/4, -3/4, -1/2, -1/2, -1/2, 
  -1/2, -1/4, -1/5, 3/5, 2/3, 2/3, 1, 1, 1, 2 ]
gap> time;
993
gap> SortedList(MyRootsOfUPol(f))=SortedList(RootsOfUPol(f));
true
gap> 

Данный пример показывает, что для сложных многочленов стандартная функция RootsOfUPol все-таки находит рациональные корни значительно быстрее. Таким образом, наша разработка не может претендовать на замещение стандартной функции системы. Однако, в отличие от RootsOfUPol, использующей поле разложения многочлена, функции, предлагаемые нами, используют минимальный набор теоретических сведений из стандартного курса алгебры и теории чисел, и могут являться полезным учебным примером, показывающим, что нужно на самом деле сделать для того, чтобы запрограммировать алгоритм, базирующийся всего лишь на на нескольких простых утверждениях.

---

Предлагаем Вам программу для системы GAP, которая возвращает ряд Штурма для заданного многочлена c рациональными коэффициентами, при этом выводя на печать необходимые промежуточные вычисления. Кроме того, программа также содержит функции, которые используют ряд Штурма для определения числа положительных и отрицательных корней многочлена, числа его корней в заданном интервале, и вычисления приближенного значения корней многочлена.

Текст программы находится здесь. Для начала работы с ней загрузите предлагаемый файл и прочитайте его с помощью команды

Read("sturm.g");

Сначала зададим некоторый многочлен от одной переменной:

gap> x:=Indeterminate(Rationals,"x");
x
gap> f:=x^3-6*x+1;
x^3-6*x+1

Первой функцией, которую мы продемонстрируем, будет SturmSystem. Данная функция находит ряд функций Штурма для заданного многочлена. Заметим, что после вычисления система Штурма сохраняется в качестве атрибута многочлена, что исключает затраты времени на ее повторное вычисление. Чтобы убедиться в этом, воспользуемся функцией KnownAttributesOfObject(f), которая позволяет узнать, какими атрибутами обладает объект (в нашем случае многочлен f), до и после вычисления его системы Штурма:

gap> KnownAttributesOfObject(f);
[ "CoefficientsOfLaurentPolynomial", "IndeterminateNumberOfUnivariateRationalFunction" ]
gap> s:=SturmSystem(f);
f(x) = x^3-6*x+1
f'(x) = 3*x^2-6
x^3-6*x+1 = 3*x^2-6 * 1/3*x - 4*x-1
3*x^2-6 = 4*x-1 * 3/4*x+3/16 -93/16
4*x-1 = 93/16 * 64/93*x-16/93 -0
[ x^3-6*x+1, 3*x^2-6, 4*x-1, 93/16 ]

gap> time;
440
gap> s:=SturmSystem(f);
[ x^3-6*x+1, 3*x^2-6, 4*x-1, 93/16 ]
gap>
time;
0

gap> KnownAttributesOfObject(f);
[ "CoefficientsOfLaurentPolynomial", "IndeterminateNumberOfUnivariateRationalFunction", 
  "IndeterminateOfUnivariateRationalFunction", "DegreeOfLaurentPolynomial",
  "Derivative", "SturmSystem" ]

Как видно из примера, время повторного вычисления равняется нулю, а список известных атрибутов пополняется, среди прочих, атрибутом SturmSystem.

Далее рассмотрим функцию ChangeSigns(s,a), которая определяет число перемен знаков в ряде Штурма s в заданной точке a:

gap> ChangeSigns(s,0);
2
gap> ChangeSigns(s,-3);
3
gap> ChangeSigns(s,2);
1

Функция NrRoots(f,a,b) находит число действительных корней многочлена f, которые больше чем a и не превышают b:

gap> NrRoots(f,-3,2);
2

Для подсчета количества перемен знаков в ряде Штурма s на плюс и минус бесконечности применяются функции ChangeSignsMinusInfinity(s) и ChangeSignsPlusInfinity(s):

gap> ChangeSignsPlusInfinity(s);
0
gap> ChangeSignsMinusInfinity(s);
3

С помощью этих двух функций легко можно определить количество корней многочлена f в функции NrOfRealRoots(f):

gap> NrOfRealRoots(f);
3

Функции NrOfPositiveRealRoots(f) и NrOfNonPositiveRealRoots(f) определяют соответственно количество положительных и отрицательных корней многочлена f:

gap> NrOfNonPositiveRealRoots(f);
1
gap> NrOfPositiveRealRoots(f);
2

Для нахождения корней многочлена f, сначала требуется найти интервал, на котором содержатся все корни многочлена или, что то же самое, найти верхнее ограничение модуля корней. Для этого была разработана функция UpperRootsLimit(f):

gap> UpperRootsLimit(f);
7

Функция IntervalsOfRoots(f,eps) выводит интервал на котором содержится корень многочлена f, с заданной точностью eps:

gap> IntervalsOfRoots(f,10);
[ [ -7, 0 ], [ 0, 7 ] ]

gap> IntervalsOfRoots(f,6);
[ [ -7/2, 0 ], [ 0, 7/2 ] ]

gap> IntervalsOfRoots(f,1);
[ [ -21/8, -7/4 ], [ 0, 7/8 ], [ 7/4, 21/8 ] ]

gap> IntervalsOfRoots(f,10^-3);
[ [ -1295/512, -20713/8192 ], [ 1365/8192, 343/2048 ],
[ 19341/8192, 4837/2048 ] ]

gap> List(last, r -> (r[1]+r[2])/2 );
[ -41433/16384, 2737/16384, 38689/16384 ]

gap> IntervalsOfRoots(f,10^-50);
[ [ -946180540226875416749522184033111355909143980149321/
374144419156711147060143317175368453031918731001856,
-1892361080453750833499044368066222711818287960298635/
748288838313422294120286634350736906063837462003712 ],
[ 15662545086391001699568813492298469860337205786443/
93536104789177786765035829293842113257979682750464,
125300360691128013596550507938387758882697646291551/
748288838313422294120286634350736906063837462003712 ],
[ 441765179940655704975623465031958738233897578501771/
187072209578355573530071658587684226515959365500928,
1767060719762622819902493860127834952935590314007091/
748288838313422294120286634350736906063837462003712 ] ]

В качестве приближенных значений корней можно взять середину полученных интервалов:

gap> List(last, r -> (r[1]+r[2])/2 );
[ -3784722160907501666998088736132445423636575920597277/
1496577676626844588240573268701473812127674924007424,
250600721382256027193101015876775517765395292583095/
1496577676626844588240573268701473812127674924007424,
3534121439525245639804987720255669905871180628014175/
1496577676626844588240573268701473812127674924007424 ]

Как видно из примера, с помощью функции для вычисления системы Штурма и вспомогательных функций можно вычислить рациональные приближения корней многочлена с заданной точностью. Заметим, что рациональные корни многочлена можно найти с помощью функции RootsOfUPol (см. пример по теме "Многочлены с рациональными коэффициентами").

---