Надежный алгоритм поиска корней полиномиальных уравнений

Предположим, требуется найти все вещественные корни полинома:

Тогда решатель полиномиальных уравнений общего назначения это алгоритм поиска корней этого полинома, удовлетворяющий следующим критериям:

• КРИТЕРИЙ 1: Алгоритм должен быть абсолютно надежен, тоесть для любых коеффициентов, являющихся вещественными числами, все корни полинома гарантированно должны быть найдены. Исключение могут составлять случаи, когда корни настолько велики, что значение полинома в них не может быть вычислено по тем или иным причинам.

• КРИТЕРИЙ 2: Количество операций (арифметических и прочих), требуемое для нахождения всех корней полинома с заданной точностью, как функция от коеффициентов, ограничена сверху. Другими словами, ситуация, в которой с приближением вектора коеффициентов к некоторой точке в , время работы алгоритма неограничено растет, недопустима. В зависимости от численного метода, такими точками может быть, например, множество всех точек, в которых полином имеет кратные корни (такое подмножество называется дискриминантной поверхностью), или точки, в которых один или несколько старших коеффициентов становятся равны нулю (с приближением к таким точкам некоторые корни устремляются в бесконечность).

• КРИТЕРИЙ 3: Алгоритм должен быть вычислительно устойчив. Переход от вещественных чисел к числам с ограниченной точностью не должен нарушать принципы работы алгоритма. Если в результаты работы алгоритма и вносятся некоторые погрешности (при таком переходе), то они должны уменьшаться с увеличением точности представления чисел.

• КРИТЕРИЙ 4: Асимптотическая сложность алгоритма, как функция от степени полинома, количества итераций, и размера чисел, в которых производятся вычисления, должна быть полиномиальной (или лучше) по любой из переменных.

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

Индекс точки

Основой метода, который предлагается в этой статье, является некоторая целочисленная ступенчатая функция от вещественного , соответствующая рассматриваемому полиному. Будем в дальнейшем называть ее «индекс точки».

Приведем принцип ее вычисления. Для этого введем последовательность:

Тогда функция в точке равна числу перемен знака в ряду значений элементов , взятых в точке . Если какие либо значения равны при этом нулю, они удаляются из списка перед подсчетом числа перемен знака.

Индекс точки обладает совокупностью свойств, ключевых в поиске корней полинома:

• В точке корня полинома кратности индекс точки убывает на единиц.

• В точке корня кратности производной порядка индекс точки убывает на четное число единиц, близжайшее кратности .

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

• Количество единиц, на которое (в сумме) убывает индекс точки на всей числовой прямой в точности равняется степени полинома. Другими словами, индекс точки, монотонно убывая, изменяется от до нуля.

• В точках корней полинома и его производных индекс точки доопределен так, что интервалы постоянства его значения открыты вправо.

Рассмотрим поведение индекса точки на примере полинома

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

Доказательство (корень полинома)

Докажем, что в точке кратного корня полинома индекс точки убывает на величину кратности этого корня. В этом случае справедливо представление:

Где полином, очевидно не равный нулю в точке , иначе кратность была бы выше по крайней мере на единицу. Найдем производную:

В текущем контексте полином можно заменить на полином . Справедливо это по той причине, что в точке , а значит и в некоторой ее окрестности эти полиномы не равны нулю и при этом совпадают по знаку. Значит, и перемена знака между и при такой замене останется прежней в некоторой окрестности корня .

После такой замены между и имеется общий множитель, равный:

Очевидно, что в перемене знака между двумя числами, имеющими общий множитель, его можно сократить без изменения такой перемены знака. Сократив и на указанный множитель, имеем перемену знака между и . Ясно, что такая перемена всегда убывает (пропадает) в точке .

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

Доказательство (корень производной)

Теперь докажем, что в точке корня кратности производной порядка индекс точки убывает на четное число, близшайшее величине кратности корня.

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

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

Рассмотрим случай, когда четно. Нетрудно показать, что в этом случае -ная производная не изменит знак своего значения в точке корня, тоесть идентична по знаку константе в некоторой выколотой окрестности корня. Следовательно, и соответствующая перемена знака между двумя знаковыми константами не изменится в точке .

Рассмотрим случай, когда нечетно. Нетрудно показать, что в этом случае -ная производная всегда меняет знак своего значения в точке корня. Следовательно, соответствующая перемена знака всегда либо пропадает, либо появляется.

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

Доказательство (значение в корнях)

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

В точке корня кратности первые элементов равны нулю, значит, первые перемен знака попросту не будут рассматриваться, так как соответствующие нули будут извлечены. В то же время, в правом (относительно точки корня) интервале постоянства индекса те же перемен знака пропадают (убывают). Следовательно, значение индекса в точке корня полинома доопределяется его правым значением (его значением в правом интервале постоянства).

В точке корня кратности производной порядка соответствующие перемен знака также выпадают из рассмотрения. Зато перемена знака между производными и порядков заменяется на перемену знака между производными и порядков. Остается показать, что правое значение производной -ого порядка и значение (в точке корня) -той производной совпадают по знаку. Действительно, дифференцируя раз представление

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

Доказательство (суммарное изменение)

Покажем, что на всей вещественной прямой индекс точки, убывая, изменяется от до нуля, где это степень полинома.

Пусть произвольный полином и его производная имеют вид:

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

По итогу имеем равнозначную перемену знака между следующими полиномами:

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

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

Так как для полинома степени индекс точки состоит из перемен знака, имеем, что индекс точки полинома действительно убывает от до нуля.

Прочие теоретические выводы

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

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

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

Скорее наоборот, отсюда следует, что если какая-то производная имеет корень, например, второй кратности, то полином не может иметь вещественных корней более, чем степень полинома, уменьшенная на два. Если же какая-то производная вдовесок имеет корень пятой кратности, то допустимое количество корней уменьшается на «два плюс по крайней мере четыре» единицы, в сумме на шесть единиц.

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

Алгоритм поиска корней №1

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

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

Далее по тексту величина обозначает в каком-то смысле середину интервала . Для простоты вначале ее можно считать равной среднему арифметическому и . Величина обозначает в каком-то смысле ширину интервала. Ее можно считать равной .

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

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

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

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

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

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

По окончании процесса рекурсивного деления считается, что все корни найдены.

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

Алгоритм поиска корней №2

Для нахождения всех корней полинома, очевидно, достаточно найти все точки, в которых индекс точки изменяется. Такие точки мы будем называть «ступенями». Из них затем достаточно отобрать те, которые соответствуют вещественным корням полинома.

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

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

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

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

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

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

Частичный индекс точки

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

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

Оптимизация такого вида может быть применена к обоим ранее изложенным алгоритмам, она однозначно уменьшает вычислительные затраты, однако я не выяснял, изменяет ли такая оптимизация асимптотическую сложность алгоритмов.

Вырождение полинома

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

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

Считается, что при любых значениях внешних переменных «реальная» степень полинома не может превышать фиктивной.

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

Ограничение точности вычислений

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

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

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

На практике повсеместно используется формат представления вещественных чисел с ограниченной точностью, называемый «floating point», или ФПЗ (формат чисел с плавающей запятой). В таком формате (со значительными упрощениями) двоичное число определенной ширины взаимно-однозначно и монотонно отображается в вещественное число , так, что максимальная относительная ошибка (называемая ошибкой округления), возникающая в результате такого преобразования, не превышает .

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

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

Рассмотрим реальный пример: дан полином с корнями вида:

Требуется вычислить значение (прежде всего знак значения) полинома в очень малой окрестности третьего корня (), используя ФПЗ единичной точности (32 бита). Замечу, что коеффициенты этого полинома целые числа допустимой величины, могущие быть точно представлены в указанном формате. Имеем картину:

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

Видим, что знак значения меняется отнюдь не раз в окрестности выбранного корня, более того, оказывается полином имеет в этой окрестности аж восемь «корней». Для справки, ширина окрестности колебания знака на изображении составляет .

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

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

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

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

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

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

Рассмотрим теперь корень кратности производной порядка . В этом случае из перемен знака перемен убывают в точке корня всегда. Оставшаяся перемена всегда меняет значение в любом направлении ( нечетно) или не изменяется ( четно).

В изображенном ряду значения производных, которые (в переходной зоне точки ) могут иметь произвольный знак (в том числе быть равными нулю). Значения и обозначают граничащие с указанными производными значения элементов . Таковые могут иметь произвольный знак (кроме нуля), но сохраняют его в переходной зоне. Обозначим за сумму перемен знака в этом ряду чисел.

Приведу несколько фактов, доказательство которых элементарно:

• Если и совпадают по знаку, а четно, то .

• Если и совпадают по знаку, а нечетно, то .

• Если и отличаются по знаку, а четно, то .

• Если и отличаются по знаку, а нечетно, то .

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

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

Рассмотрим на втором алгоритме. Пусть искомая ступень заведомо содержится в интервале. Если точка оказалась вне переходных зон, то решение о том, в каком из интервалов или находится искомая ступень, всегда принимается верно.

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

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

С первым алгоритмом размышления аналогичны.

Почему не теорема Штурма ?

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

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

Пусть полиномы и , расположенные последовательно, имеют вид:

Тогда, применяя алгоритм Евклида, имеем, что полином , в последовательности следующий за полиномами и , вычисляется по несложной формуле:

Очевидно, сложность вычисления коеффициентов полинома пропорциональна степени , (предполагая, что каждая арифметическая операция производится точно и занимает фиксированное время). Следовательно, вычисление всех полиномов последовательности займет время, пропорциональное .

Теорема Штурма является важным теоретическим результатом алгебры. Однако, несмотря на некоторые преимущества, она не сыскала широкого применения на практике (по крайней мере, применительно к рассматриваемой задаче) ни в 19 веке, ни сейчас.

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

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

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

Сепаратор полинома

Индекс точки в текущем определении несовершенен в том смысле, что не все точки, в которых он убывает, соответствуют вещественным корням полинома. Такой недостаток хорошо проявляется, например, в случае, когда требуется найти минимальный корень полинома. Придется найти несколько побочных корней, чтобы «добраться» до первого вещественного корня. В целом с увеличением степени полинома статистически все меньшее количество корней полинома являются вещественными.

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

Здесь для решения той же проблемы я предлагаю поступить немного иначе:

Cепаратором полинома будем называть полином, который, если и имеет корни, то только по одному в каждом промежутке между корнями . Тем самым, корни , чередуясь с корнями , как-бы отделяют корни полинома друг от друга.

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

Выясним, как себя ведет перемена знака между полиномом и его сепаратором:

Видим, что в каждой точке корня перемена убывает, а в точке корня , снова появляется (на изображении в верхних строках таблиц указан знак полинома в окрестности корня, в нижних указан знак сепаратора).

Допустим, индекс сепаратора известен. Покажем, что для того, чтобы получить искомый индекс полинома достаточно прибавить к перемену знака между полиномом и его сепаратором:

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

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

Таким образом, для вычисления для полинома требуется знать для сепаратора (обозначим его ) . В свою очередь, требуется знать

Пусть полином пятой степени имеет три корня. Тогда должен иметь два корня, являясь полиномом четвертой степени. Далее должен иметь один корень, являясь полиномом третьей степени. Продолжая, не должен иметь кореней, являясь полиномом второй степени. Но в случае , он также не должен иметь корней, вопреки тому, что имеет первую (нечетную!) степень.

Исходя из такого противоречия можно сделать вывод о том, что в отличие от ряда Штурма, всегда имеющего полиномов, ряд

должен иметь полиномов ровно на один больше того, сколько корней у.

Можно заметить, что случай кратных корней у здесь не был рассмотрен. Определение сепаратора полинома легко расширить, например, так: если два корня совпадают в точке , то промежутком между ними будет считаться сама точка (соответственно, если совпадают три корня, возникает два одноточечных промежутка, и т.д). Тем самым, в точке корня кратности индекс убывает на величину кратности .

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

Поиск на всей числовой прямой

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

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

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

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

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

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

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

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

Полиномиальная оценка

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

В первом случае число хранимых коеффициентов растет квадратично со степенью полинома, значит уже степени на 10 никакого регистрового файла (даже с SIMD) для хранения этих коеффициентов не хватит. Значит, будут обращения к кешу первого уровня, каждое из которых занимает 5-7 циклов, в то время как всего за один такт процессоры способны выполнять до четырех двойных инструкций double FMA (умножение с накоплением). Во втором же случае усложнение с дополнительными умножениями способно замедлить алгоритм по крайней мере раза в два.

Существует более выгодный алгоритм вычисления всех производных полинома в одной точке, открытый, по-видимому, также Горнером (или Безу, источник мне теперь не известен) и найденный мной в книге «The art of computer programming» Дональда Кнута:

Для от до выполнить:
Для от до выполнить:
.
закончить.
закончить.

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

Связь между производными полинома , взятыми в точке и коеффициентами по полинома очень просто установить. Достаточно разложить по формуле Тейлора в соответствующую сумму:

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

Разберем алгоритм подробнее на примере полинома:

Тогда последовательность операций выглядит следующим образом:

На этой диаграмме хорошо заметны связи между узлами. Вычисление каждого из них, в соответствии с таким алгоритмом, происходит «слева-направо и сверху-вниз», тоесть по горизонталям, отделенным пунктирной линией:

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

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

Оригинальный алгоритм нетрудно доработать так, чтобы локальная параллельность была существенно выше. Достаточно переупорядочить последовательность вычисления значений узлов таким образом:

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

Переупорядочивая операции, тот же алгоритм можно записать таким образом:

Для от до выполнить:
Для от до выполнить:
.
закончить.
закончить.

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

Комплексные корни (формальное решение)

Еще в 1875 году Раусс нашел способ (называемый алгоритмом Раусса) по вещественным коеффициентам полинома определить количество корней, вещественные части которых положительны (читайте Гантмахера, теория матриц). Его выкладки могут быть основаны на индексе Коши рациональной функции, принципе агрумента и теореме Штурма.

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

Ползуясь алгоритмом №1 или №2, заменив индекс на индекс , можно найти все вещественные части корней полинома. Мнимые их части затем ищутся элементарно.

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

исключать переменную , пользуясь теорией исключения (результантами), или тем более применять базисы Гребнера.

Асимптотическая сложность

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

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

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

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

Заключение

Здесь я не привел сравнение этого метода с другими, достаточно известными методами, также теми, которые применяются в системах компьютерной алгебры (таких как Maxima, Mathematica, Wolfram, и др.). Не сравнил то, как общие и локальные оптимизации влияют на время вычислений. Не рассказал об истории проблемы отделения корней, тоесть о правиле знаков Декарта, о теореме Штурма, теореме Будана-Фурье, и о многом другом.

Если у меня будет время и, что немаловажно, желание читателей, я постараюсь написать вторую часть, в которой, надеюсь, изложу указанное. Надеюсь также, что многие из тех читателей, кто задавался вопросом «как численно и надежно, и за заранее ограниченное время решить такое-то уравнение такой-то степени», теперь имеют на него ответ.

Поправки, предложения и конструктивная критика приветствуются, с учетом того, что это моя первая статья на Хабре.