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;
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, использующей поле разложения многочлена, функции, предлагаемые нами, используют минимальный набор теоретических сведений из стандартного курса алгебры и теории чисел, и могут являться полезным учебным примером, показывающим, что нужно на самом деле сделать для того, чтобы запрограммировать алгоритм, базирующийся всего лишь на на нескольких простых утверждениях.