Рассмотрим структуру класса для работы с многочленами, алгоритмы для работы с объектами типа «полином» и примеры методов полиномиального класса, реализующие эти алгоритмы:
template <class YourOwnFloatType> class polynom: public vector<YourOwnFloatType>
{
Если мы запишем коэффициенты полинома в порядке убывания, включая нулевые, мы получим упорядоченный кортеж длиной n+1: (an, an-1, …, a2, a1, a0), то есть не что иное, как вектор размерности n+1, полностью характеризующий заданный полином. Поэтому вполне естественным является то, что наш полином будет базироваться на векторе.
Как и векторный класс, он будет параметризованным. При этом тип-параметр будет у нас относиться как к коэффициентам многочлена, так и к подставляемым в него значениям. Внутренний формат для хранения нашего полинома будет a0+a1x+a2x2+a3x3+a4x4+a5x5+...+an-1xn-1, что обусловлено требованием удобства его индексирования, а внешнее представление будет в канонической форме –ввод и вывод будет производиться, начиная с коэффициента при наивысшей степени.
void optimize();
В процессе работы с полиномом мы можем прийти к ситуации, когда его порядок необходимо уменьшить в связи с тем, что обнулился коэффициент при старшей степени (или несколько подряд идущих коэффициентов). В связи с этим после каждой операции, которая может привести к изменению степени полинома, производится проверка, записан ли полином в канонической форме. Если нет, мы с помощью этой внутренней функции понижаем его порядок до тех пор, пока полином не будет преобразован в каноническую форму.
polynom<YourOwnFloatType> reverse();
В связи с различием внутреннего и внешнего представления полинома иногда бывает необходимо записать полином в обратном порядке.
void In(istream &);
void Out(ostream &);
По аналогии с векторами переопределим две виртуальные функции ввода и вывода. При этом отпадает необходимость в перегрузке операций потокового ввода-вывода: однажды определённые в векторном классе, они используют именно эти виртуальные функции, а определение, из какого именно класса необходимо их вызывать, осуществляется уже на этапе выполнения, в зависимости от того, к какому типу преобразуется базовый указатель. Так, любой полиномиальный объект можно преобразовать к векторному.
public:
В этом разделе мы разместим общедоступные дружественные функции и методы полиномиального класса.
polynom(char *);
Полином из файла? Почему бы и нет! Параметром этого конструктора является имя файла, в котором находятся данные в виде [Степень Многочлена-1] [Свободный Член] [Коэффициент При Первой Степени] … [Коэффициент При Старшей Степени].
polynom(long,YourOwnFloatType *);
Конструктор, принимающий два параметра –уменьшенную на единицу степень многочлена и указатель на данные –коэффициенты многочлена, начиная со свободного члена. Почему бы не передавать в этот конструктор степень многочлена? Просто показалось удобнее пользоваться понятием не «степень многочлена», а «размерность полинома», то есть количество коэффициентов, его составляющих.
polynom(polynom<YourOwnFloatType> &);
Конструктор копирования делает слепок с заданного полинома.
polynom();
Конструктор по умолчанию создаёт полином особого вида –нуль-полином, то есть многочлен размерности 1 (соответственно степени 0), свободный член которого равен нулевому элементу типа-параметра полиномиального класса.
polynom(long);
В некоторых случаях бывает полезно задать сначала степень полинома, но временно оставить неопределёнными его коэффициенты. Этот конструктор создаёт полином заданной размерности и обнуляет коэффициенты, нарушая каноническую форму записи многочлена. Это бывает необходимо редко и только в тех случаях, когда коэффициенты многочлена становятся известными после того, как задана его степень. Разумеется, и в этом случае параметром конструктора является степень многочлена-1.
polynom<YourOwnFloatType> operator-();
Полином, противоположный полиному P(x)=anxn+ +an-1xn-1+…+a2x2+a1x+a0, определяется как –P(x)=–anxn– –an-1xn-1–…–a2x2–a1x–a0, то есть знак коэффициентов при степенях многочлена меняется на противоположный.
polynom<YourOwnFloatType> operator+();
Унарный плюс –операция, просто возвращающая копию текущего многочлена.
Дружественная функция сложения многочленов принимает два параметра – многочлены-слагаемые, и возвращает результирующий многочлен-сумму. Суммой двух многочленов, f(x)=anxn+an-1xn-1+…+a2x2+a1x+a0 и g(x)= =bmxm+bm-1xm-1+…+b2x2+b1x+b0 степеней n и m соответственно (n≥m) называют многочлен cnxn+cn-1xn-1+…+ +c2x2+c1x+c0, где ci=ai+bi (i=0, 1, …, m), ci=ai (i=m+1, …, n). Разумеется, после получения многочлена-суммы его необходимо записать в канонической форме.
Дружественная функция, реализующая операцию сокращённого сложения двух многочленов, складывает первый многочлен со вторым и записывает результат в первый, возвращая его копию для дальнейших преобразований.
Так как операция сложения многочленов обратима, то мы можем определить вычитание многочленов, используя операции сложения и получения многочлена, противоположного к данному. Таким образом, разностью двух полином f(x) и g(x) является полином P(x) такой, что
Умножение числа на полином, равно как и умножение полинома на число – это операция, воздействующая на коэффициенты полинома, и соответствует умножению числа на вектор (вектора на скаляр), то есть:
friend long operator<(polynom<YourOwnFloatType> &, polynom<YourOwnFloatType> &);
friend long operator>(polynom<YourOwnFloatType> &, polynom<YourOwnFloatType> &);
friend long operator<=(polynom<YourOwnFloatType> &, polynom<YourOwnFloatType> &);
friend long operator>=(polynom<YourOwnFloatType> &, polynom<YourOwnFloatType> &);
Набор операций для сравнения полиномов. Два полинома считаются равными, если они одинаковой степени и коэффициенты при соответствующих степенях х равны. В противном случае полномы не равны. При этом считается, что первый полином меньше второго, если в канонической форме его степень ниже (размерность меньше). Если же эти полиномы одинаковой размерности (а, соответственно, и степени), то мы последовательно сравниваем их коэффициенты, начиная со старшей степени. Первое различие в коэффициентах и определяет, какой знак необходимо поставить между этими многочленами.
YourOwnFloatType &operator[](long);
Индексация полинома –операция, принимающая в качестве параметра степень одночлена, коэффициент при котором необходимо вернуть.
YourOwnFloatType operator()(YourOwnFloatType);
Получить значение полинома в заданной точке мы можем, используя эту функцию. Её задача – либо просто просуммировать произведения коэффициентов многочлена на соответствующую степень аргумента данной функции, либо вычислить это значение по теореме Безу.
friend long div(polynom<YourOwnFloatType> &, polynom<YourOwnFloatType> &, polynom<YourOwnFloatType> &, polynom<YourOwnFloatType> &);
Из курса алгебры известна теорема, гласящая о том, что для любых двух многочленов f(x) и g(x) существует единственная пара многочленов l(x) и r(x), удовлетворяющих условию
f(x)=g(x)l(x)+r(x),
где степень r(x) меньше степени g(x) или r(x) – нуль-многочлен.
Для определения вида этих многочленов рассмотрим многочлены
2. Пусть n≥m. В этом случае составим вспомогательный многочлен степени n1. Если n1≤m, то ограничиваемся одним этим многочленом, если же n1≥m, то строим (где – старший коэффициент многочлена f1(x)). Если степень n2 многочлена f2(x) меньше m, то ограничиваемся построенными двумя многочленами, если же n2≥m, то продолжаем (аналогично предыдущему) строить вспомогательные многочлены f3(x), f4(x), …, fs(x), …. При этом степень каждого следующего многочлена меньше степени предыдущего, поэтому на определённом k-м шаге имеем многочлен , степень которого меньше m (или fk(x)=0).
Сложим почленно полученные равенства:
f1(x)+f2(x)+…+fk(x)=f(x)+f1(x)+…+fk-1(x)–
.
Выражение в скобках представляет собой сумму многочленов, а потому (по определению суммы) также является многочленом. Обозначим его через l(x). Многочлен fk(x) обозначим r(x). По построению, степень r(x) меньше степени f(x), или r(x)=0.
По аналогии, назовём f(x) делимым, g(x) – делителем, l(x) – неполным частным, а r(x) – остатком от деления f(x) на g(x).
Используя описанную выше функцию, определим две вспомогательные операторные функции –для нахождения целой части и остатка от деления.
long IsEqual(void *);
Виртуальная функция с таким же именем, определённая в векторном классе, предназначена для сравнения текущего вектора с вектором, адрес которого передаётся через обобщённый указатель. Та же операция проделывается и для полиномов с одной единственной целью –использовать единожды, в векторном классе, определённые операторные функции == и !=, что увеличивает гибкость и универсальность данного класса за счёт использования механизма виртуальных функций.
Полином, как функция аналитическая, удобен тем, что для него всегда можно найти аналитическую производную. Как известно, производной степенной функции является степенная же функция, поэтому производная от полинома будет тоже полином:
.
Параметром данной функции является порядок производной –количество раз, которое полином дифференцируется.
По аналогии с аналитической производной для полинома можно определить аналитический интеграл. Первообразная степенной функции –степенная же функция степени на единицу выше, поэтому интеграл от полинома также будет полиномом:
.
Параметром данной функции является кратность интегрирования –количество раз, которое полином интегрируется.
От параметризованных полиномов, тип которых определяется подстановкой в классовые уголки, перейдём к комплексным полиномам вида
f(z)=anzn+an-1zn-1+…+a2z2+a1z+a0,
подставив тип «комплексное число»:
typedef polynom<complex> cpolynom;
typedef vector<complex> cvector;
typedef matrix<complex> cmatrix;
Значения z, при подстановке которых многочлен обращается в нуль, называются корнями этого многочлена. Таким образом, корни f(z) – это решения уравнения:
f(z)=anzn+an-1zn-1+…+a2z2+a1z+a0=0,
называемого алгебраическим уравнением n-й степени.
По теореме Безу, остаток, получаемый при делении многочлена f(z) на (z–a), равен f(а). В частности, для того чтобы многочлен f(z) делился на (z–а) без остатка, необходимо и достаточно условие f(a)=0, т. е. для того чтобы многочлен делился на двучлена (z–а) без остатка, необходимо и достаточно, чтобы z=a было корнем этого многочлена.
Таким образом, зная корень z=a многочлена f(z), мы можем выделить из этого многочлена множитель (z–a): f(z)=(z–a)f1(z), где f1(z)=bn-1zn-1+bn-2zn-2+…+b2z2+b1z+b0 (bn-1=an-1); нахождение остальных корней приводит к решению уравнения
bn-1zn-1+bn-2zn-2+…+b2z2+b1z+b0=0
(n–1)-й степени.
Продолжая этот процесс, мы получим окончательно следующее разложение f(z) на множители: f(z)=a0(z–zl)(z– –z2)...(z–zn), т.е. всякий многочлен n-й степени разлагается на (n+1) множителей, один из которых равен старшему коэффициенту, а остальные – двучлены первой степени вида (z–а). При подстановке z=zs (s=1, 2, ..., n) по крайней мере один из множителей в разложении обратится в нуль, т. е. значения z=zs суть корни f(z).
Рассмотрим частные методы решения алгебраических уравнений в радикалах, а также метод Ньютона для нахождения всех корней комплексного полинома.
Пусть задано алгебраическое уравнение степени n коэффициентами определяющего его многочлена. Количество корней этого уравнения будет равно степени многочлена, и в общем случае эти корни будут комплексными, поэтому возвращаемым значением функции, предназначенной для нахождения корней полинома, будет вектор размерности n, компоненты которого суть корни полинома.
Первая из рассматриваемых функций предназначена для решения алгебраического уравнения второй степени, или, как его ещё называют, квадратного уравнения. Параметрами его являются коэффициенты при второй, первой и нулевой степенях соответственно. Алгоритм решения такого уравнения определяется следующими соотношениями:
Для решения уравнений третьей и четвёртой степеней с комплексными коэффициентами общих методов нет. Однако существуют два частных метода для решения уравнений этих степеней с действительными коэффициентами, метод Кардано-Тартальи для решения уравнений третей степени и метод Феррари для решения уравнений четвёртой степени. Так как эти методы не охватывают все уравнения данного класса, мы на них останавливаться не будем, однако приведём ниже пример их реализации.
cvector newton(cpolynom);
Пусть необходимо найти все корни уравнения P(z)=0. Для численного решения можно воспользоваться итерационным методом Ньютона. Замечательной особенностью этого метода является его гарантированная сходимость на всей комплексной плоскости для выпуклых функций (а полиномы степени выше первой – функции выпуклые), причём скорость сходимости не зависит от начального приближения. Для начала итерации зададимся точностью ε, с которым нам необходимо найти решение, и начальным приближением к некоторому корню z0=x0+iy0. Алгоритм вычислений может быть таким:
1. Находим производную полинома в точке zk(k – номер итерации). Если она равна нулю, дифференцируем P(z) дальше до тех пор, пока P(j)(zk)не станет отличной от нуля.
2. Следующее приближение к корню находим по формуле , где параметр t подбирается так, чтобы выполнялось условие: |P(zk+1)|<|P(zk)|.
3. Проверяем выполнение условия |P(zk+1)|<ε; если условие не выполняется, переходим к пункту 1. Если условие выполняется, то zk+1 считаем корнем уравнения, полином P(z) делится на двучлен (z–zk) и получаем полином (n–1)-ой степени, для которого повторяем все вычисления, начиная с п. 1.
Описанные алгоритмы, к примеру, можно реализовать так:
#ifndef __EQUATION_H
#define __EQUATION_H
/*Этот файл включения содержит объявление двух типов - комплексного полинома и комплексного вектора, а также заголовки четырёх функций для решения уравнений 2-ой, 3-ей, 4-ой и высших степеней */
/*Общий вид уравнения, корни которого мы ищем - a3*x^3+a2*x^2+a1*x+a0=0. Так как уравнение имеет третью степень, то и корней у него три, что и определяет размерность вектора-результата*/
cvector res=3;
if(a3==0)//если коэффициент при x^3 нулевой,
{
//решаем соответствующее квадратное уравнение
cvector res2=square(a2,a1,a0);
/*в этом случае принимаем, что два корня являются совпадающими */
res[0]=res[1]=res2[0];
res[2]=res2[1];
}
else
{
/*иначе - приводим кубическое уравнение к каноническому виду, когда коэффициент при третьей степени равен 1 */
double a=a2/a3,b=a1/a3,c=a0/a3;
//находим слагаемые кубического дискриминанта
double p=b-a*a/3,q=2*a*a*a/27-a*b/3+c;
/*проведя переобозначение x=y-a/3, решаем в дальнейшем уравнение y^3+p*y+q=0 */
//находим кубический дискриминант
double diskr=pow(q/2,2)+pow(p/3,3);
/*если дискриминант равен 0, то имеем 3 действительных корня, из них два совпадающих */
if(diskr==0)
{
res[0]=3*q/p;
res[1]=res[2]=-res[0]/2;
}
/*один действительный корень и два комплексно сопряжённых*/
/*В этом цикле находится только один корень z. Для нахождения остальных делим исходный полином на х-z, понижая тем самым степень на единицу, и находим ещё один корень, и т.д. до первой степени*/
result[0]=z;
cpolynom z1=2;
z1[0]=-z,z1[1]=1;
cvector more=newton(p/z1);
for(long i=1;i<result.getm();i++)
result[i]=more[i-1];
return result;//возвращаем вектор результата
}
По приведённым алгоритмам работы с полиномиальными объектами можно составить, к примеру, такой интерфейс для данного класса:
/*Перегружать операторы ввода из потока и вывода в поток необходимости нет - достаточно перегрузить две виртуальные функции ввод полинома из потока, начиная с коэффициентов при старших степенях*/
template <class YourOwnFloatType>
void polynom<YourOwnFloatType>::In(istream &is)
{
for(long i=getm()-1;i>=0;i--)
is>>(*this)[i];
optimize();
}
/*вывод полинома в поток, начиная с коэффициентов при старших степенях*/
Из двух многочленов одинаковой длины меньше тот, у которого коэффициент при старшей степени меньше. При равенстве рассматриваем более низкую степень и т.д.
*/
template <class YourOwnFloatType>
long operator<(polynom<YourOwnFloatType> &f, polynom <YourOwnFloatType> &s)
{
polynom<YourOwnFloatType> test1=f,test2=s;
test1.optimize();
test2.optimize();
if(test1.getm()<test2.getm())
return 1;
if(test1.getm()>test2.getm()||test1==test2)
return 0;
//точно не равны - выясняем, кто же из них меньше
for(long i=test1.getm()-1;i>=0;i--)
if(test1[i]<test2[i])
return 1;
else
if(test1[i]>test2[i])
return 0;
return 0;
}
/*все остальные операции определяем через уже известные*/
template <class YourOwnFloatType>
long operator>(polynom<YourOwnFloatType> &f, polynom <YourOwnFloatType> &s)
{
return (!(f<s))&&(f!=s);
}
template <class YourOwnFloatType>
long operator<=(polynom<YourOwnFloatType> &f, polynom <YourOwnFloatType> &s)
{
return (f<s)||(f==s);
}
template <class YourOwnFloatType>
long operator>=(polynom<YourOwnFloatType> &f, polynom <YourOwnFloatType> &s)
{
return (f>s)||(f==s);
}
/* деление полинома на полином по алгоритму Евклида */
.
степени n1. Если n1≤m, то ограничиваемся одним этим многочленом, если же n1≥m, то строим 
(где 
– старший коэффициент многочлена f1(x)). Если степень n2 многочлена f2(x) меньше m, то ограничиваемся построенными двумя многочленами, если же n2≥m, то продолжаем (аналогично предыдущему) строить вспомогательные многочлены f3(x), f4(x), …, fs(x), …. При этом степень каждого следующего многочлена меньше степени предыдущего, поэтому на определённом k-м шаге имеем многочлен 
, степень которого меньше m (или fk(x)=0).
.
.
.
, где параметр t подбирается так, чтобы выполнялось условие: |P(zk+1)|<|P(zk)|.