Polynomials.jl

Polynomials.jl — это пакет Julia, который предоставляет базовые арифметические операции, возможности интегрирования, дифференцирования, вычисления, нахождения корней и аппроксимации данных для одномерных многочленов.

Пакет Polynomials размещен на GitHub и устанавливается так же, как и другие пакеты Julia. Начиная с версии v3.0.0 требуется версия Julia не ниже 1.6.

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

using Polynomials

Создание и вычисление

Создание многочлена на основе коэффициентов, начиная с самого низкого порядка:

julia> Polynomial([1,0,3,4])
Polynomial(1 + 3*x^2 + 4*x^3)

Можно добавить необязательный параметр переменной.

julia> Polynomial([1,2,3], :s)
Polynomial(1 + 2*s + 3*s^2)

Создание многочлена на основе его корней:

julia> fromroots([1,2,3]) # (x-1)*(x-2)*(x-3)
Polynomial(-6 + 11*x - 6*x^2 + x^3)

Вычисление многочлена p в точке 1 с использованием нотации вызова:

julia> p = Polynomial([1, 0, -1])
Polynomial(1 - x^2)

julia> p(1)
0

Конструктор Polynomial сохраняет все коэффициенты с использованием стандартного базиса с вектором. Другие типы (например, ImmutablePolynomial, SparsePolynomial или FactoredPolynomial) используют иные внутренние контейнеры, которые могут иметь преимущества в ряде ситуаций.

Арифметические операции

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

julia> p = Polynomial([1,2])
Polynomial(1 + 2*x)

julia> q = Polynomial([1, 0, -1])
Polynomial(1 - x^2)

julia> 2p
Polynomial(2 + 4*x)

julia> 2 + p
Polynomial(3 + 2*x)

julia> p - q
Polynomial(2*x + x^2)

julia> p * q
Polynomial(1 + 2*x - x^2 - 2*x^3)

julia> q / 2
Polynomial(0.5 - 0.5*x^2)

julia> q ÷ p  # `div`, а также `rem` и `divrem`
Polynomial(0.25 - 0.5*x)

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

julia> p = Polynomial([1, 2, 3], :x)
Polynomial(1 + 2*x + 3*x^2)

julia> q = Polynomial([1, 2, 3], :s)
Polynomial(1 + 2*s + 3*s^2)

julia> p + q
ERROR: ArgumentError: Polynomials have different indeterminates
[...]

Однако это не касается операций с константными многочленами.

julia> p = Polynomial([1, 2, 3], :x)
Polynomial(1 + 2*x + 3*x^2)

julia> q = Polynomial(1, :y)
Polynomial(1)

julia> p + q
Polynomial(2 + 2*x + 3*x^2)

Смешение типов многочленов

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

julia> p, q = ImmutablePolynomial([1,2,3]), Polynomial([3,2,1])
(ImmutablePolynomial(1 + 2*x + 3*x^2), Polynomial(3 + 2*x + x^2))

julia> p + q
Polynomial(4 + 4*x + 4*x^2)

julia> p, q = ImmutablePolynomial([1,2,3]), SparsePolynomial(Dict(0=>1, 2=>3, 10=>1))
(ImmutablePolynomial(1 + 2*x + 3*x^2), SparsePolynomial(1 + 3*x^2 + x^10))

julia> p + q
LaurentPolynomial(2 + 2*x + 6*x² + x¹⁰)

Интегралы и производные

Многочлен p интегрируется почленно с добавлением при необходимости свободного члена C. Для ненулевых многочленов степень итогового многочлена на единицу больше степени p.

julia> integrate(Polynomial([1, 0, -1]))
Polynomial(1.0*x - 0.3333333333333333*x^3)

julia> integrate(Polynomial([1, 0, -1]), 2)
Polynomial(2.0 + 1.0*x - 0.3333333333333333*x^3)

Многочлен p дифференцируется почленно. Для ненулевых многочленов степень итогового многочлена на единицу меньше степени p.

julia> derivative(Polynomial([1, 3, -1]))
Polynomial(3 - 2*x)

Нахождение корней

Возвращает корни d (или нули) многочлена p степени d.

julia> roots(Polynomial([1, 0, -1]))
2-element Vector{Float64}:
 -1.0
  1.0

julia> roots(Polynomial([1, 0, 1]))
2-element Vector{ComplexF64}:
 0.0 - 1.0im
 0.0 + 1.0im

julia> roots(Polynomial([0, 0, 1]))
2-element Vector{Float64}:
 0.0
 0.0

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

Функция roots по умолчанию использует собственные значения сопровождающей матрицы для многочлена. Эта операция имеет сложность 𝑶(n^3).

Для многочленов с коэффициентами типа BigFloat можно легко использовать пакет GenericLinearAlgebra:

julia> p = fromroots(Polynomial{BigFloat}, [1,2,3])
Polynomial(-6.0 + 11.0*x - 6.0*x^2 + 1.0*x^3)

julia> roots(p)
ERROR: MethodError: no method matching eigvals!(::Matrix{BigFloat})
[...]

julia> using GenericLinearAlgebra

julia> roots(p)
3-element Vector{Complex{BigFloat}}:
 0.9999999999999999999999999999999999999999999999999999999999999999999999999999655 + 0.0im
  1.999999999999999999999999999999999999999999999999999999999999999999999999999931 - 0.0im
  2.999999999999999999999999999999999999999999999999999999999999999999999999999793 + 0.0im

Замечания о нахождении корней

Пакет PolynomialRoots.jl обеспечивает альтернативный подход к нахождению комплексных корней одномерных многочленов, более эффективный, чем roots. Он основан на алгоритме Скоурона (Skowron) и Гоулда (Gould).

julia> import PolynomialRoots # импорт как `roots` приводит к конфликту

julia> p = fromroots(Polynomial, [1,2,3])
Polynomial(-6 + 11*x - 6*x^2 + x^3)

julia> PolynomialRoots.roots(coeffs(p))
3-element Vector{ComplexF64}:
  3.000000000000001 - 0.0im
 1.9999999999999993 + 0.0im
 1.0000000000000002 + 0.0im

Корни всегда возвращаются как комплексные числа.

Пакет FastPolynomialRoots предоставляет интерфейс для кода на Фортране, реализующего алгоритм Ауренца (Aurentz), Маха (Mach), Робола (Robol), Вандребриля (Vandrebril) и Ваткинса (Watkins), который может обрабатывать очень большие многочлены (он имеет сложность 𝑶(n^2) и является обратно устойчивым). Пакет AMRVW.jl реализует этот алгоритм на Julia с возможностью использовать другие числовые типы.

julia> using AMRVW

julia> AMRVW.roots(float.(coeffs(p)))
3-element Vector{ComplexF64}:
 0.9999999999999997 + 0.0im
 2.0000000000000036 + 0.0im
 2.9999999999999964 + 0.0im

Корни возвращаются как комплексные числа.

Как PolynomialRoots, так и AMRVW являются универсальными пакетами и работают, например, с коэффициентами типа BigFloat.

Пакет AMRVW может работать с гораздо более широким спектром многочленов, чем roots или PolynomialRoots.roots. Например, он позволяет быстро и точно определить корни следующего случайного многочлена 1000-й степени:

julia> filter(isreal, AMRVW.roots(rand(1001) .- 1/2))
2-element Vector{ComplexF64}:
  0.993739974989572 + 0.0im
 1.0014677846996498 + 0.0im

В пакете Hecke есть функция roots. Пакет Hecke использует библиотеку Arb для эффективных, точных вычислений:

julia> import Hecke # импорт как `roots` приводит к конфликту

julia> Qx, x = Hecke.PolynomialRing(Hecke.QQ)
(Univariate Polynomial Ring in x over Rational Field, x)

julia> q = (x-1)*(x-2)*(x-3)
x^3 - 6*x^2 + 11*x - 6

julia> Hecke.roots(q)
3-element Vector{Nemo.fmpq}:
 2
 1
 3

У следующего многочлена три вещественных корня, два из которых входят в кластер, и Hecke быстро определяет их:

julia> p = -1 + 254*x - 16129*x^2 + 1*x^17
x^17 - 16129*x^2 + 254*x - 1

julia> filter(isreal, Hecke._roots(p, 200)) # `_roots`, а не `roots`
3-element Vector{Nemo.acb}:
 [0.007874015748031496052667730054749907629383970426203662570129818116411192289734968717460531379762086419 +/- 3.10e-103]
 [0.0078740157480314960733165219137540296086246589982151627453855179522742093785877068332663198273096875302 +/- 9.31e-104]
 [1.9066348541790688341521872066398429982632947292434604847312536201982593209326201234353468172497707769372732739429697289 +/- 7.14e-119]

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

Пакет IntervalRootFinding служит для определения вещественных нулей одномерных функций и может использоваться для нахождения изолирующих интервалов с вещественными корнями. Пример:

julia> using Polynomials, IntervalArithmetic

julia> import IntervalRootFinding # его метод `roots` конфликтует с `roots`

julia> p = fromroots(Polynomial, [1,2,3])
Polynomial(-6 + 11*x - 6*x^2 + x^3)

julia> IntervalRootFinding.roots(x -> p(x), 0..10)
3-element Vector{IntervalRootFinding.Root{Interval{Float64}}}:
 Root([0.999999, 1.00001], :unique)
 Root([1.99999, 2.00001], :unique)
 Root([2.99999, 3.00001], :unique)

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

Пакет RealPolynomialRoots предоставляет функцию ANewDsc для нахождения изолирующих интервалов с корнями свободного от квадратов многочлена, заданного через коэффициенты:

julia> using RealPolynomialRoots

julia> st = ANewDsc(coeffs(p))
There were 3 isolating intervals found:
[2.62…, 3.62…]₂₅₆
[1.5…, 2.62…]₂₅₆
[-0.50…, 1.5…]₂₅₆

Эти изолирующие интервалы можно уточнить для нахождения численной оценки корней в виде значений BigFloat.

julia> refine_roots(st)
3-element Vector{BigFloat}:
 2.99999999999999999999...
 2.00000000000000000000...
 1.00000000000000000000...

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

julia> p = SparsePolynomial(Dict(0=>-1, 1=>254, 2=>-16129, 17=>1))
SparsePolynomial(-1 + 254*x - 16129*x^2 + x^17)

julia> ANewDsc(coeffs(p))
There were 3 isolating intervals found:
[1.5…, 3.5…]₅₃
[0.0078740157480314960682066…, 0.0078740157480314960873178…]₁₃₉
[0.0078740157480314960492543…, 0.0078740157480314960682066…]₁₃₉

В этом примере у пакета IntervalRootFinding возникают проблемы с различением кластеризованных корней:

julia> IntervalRootFinding.roots(x -> p(x), 0..3.5)
7-element Vector{IntervalRootFinding.Root{Interval{Float64}}}:
 Root([1.90663, 1.90664], :unique)
 Root([0.00787464, 0.00787468], :unknown)
 Root([0.00787377, 0.00787387], :unknown)
 Root([0.00787405, 0.00787412], :unknown)
 Root([0.00787396, 0.00787406], :unknown)
 Root([0.00787425, 0.00787431], :unknown)
 Root([0.00787394, 0.00787397], :unknown)

В данном примере функция filter(isreal, Hecke._roots(p)) также выделяет три вещественных корня, но не так быстро.

Большинство алгоритмов нахождения корней сталкиваются с проблемами, если корни имеют кратность. Например, и ANewDsc, и Hecke.roots предполагают, что многочлен свободен от квадратов. Для многочленов, не свободных от квадратов:

доступна функция Polynomials.Multroot.multroot для нахождения корней многочлена, в том числе кратных. Она основана на работе Цзэна (Zeng).

В данном случае мы видим, что у IntervalRootFinding.roots возникают проблемы с выделением корней из-за наличия кратности:

julia> p = fromroots(Polynomial, [1,2,2,3,3])
Polynomial(-36 + 96*x - 97*x^2 + 47*x^3 - 11*x^4 + x^5)

julia> IntervalRootFinding.roots(x -> p(x), 0..10)
335-element Vector{IntervalRootFinding.Root{Interval{Float64}}}:
 Root([1.99999, 2], :unknown)
 Root([1.99999, 2], :unknown)
 Root([3, 3.00001], :unknown)
 Root([2.99999, 3], :unknown)
 ⋮
 Root([2.99999, 3], :unknown)
 Root([2, 2.00001], :unknown)

Функция roots определяет корни, но не их кратность:

julia> roots(p)
5-element Vector{Float64}:
 1.000000000000011
 1.9999995886034314
 2.0000004113969276
 2.9999995304339646
 3.0000004695656672

В свою очередь, функция multroot правильно определяет корни вместе со структурой кратности:

julia> Polynomials.Multroot.multroot(p)
(values = [1.0000000000000004, 1.9999999999999984, 3.0000000000000018], multiplicities = [1, 2, 2], κ = 35.11176306900731, ϵ = 0.0)

Помочь может функция square_free:

julia> q = Polynomials.square_free(p)
ANewDsc(q)
Polynomial(-0.20751433915978448 + 0.38044295512633425*x - 0.20751433915986722*x^2 + 0.03458572319332053*x^3)

julia> IntervalRootFinding.roots(x -> q(x), 0..10)
3-element Vector{IntervalRootFinding.Root{Interval{Float64}}}:
 Root([0.999999, 1.00001], :unique)
 Root([1.99999, 2.00001], :unique)
 Root([2.99999, 3.00001], :unique)

Аналогичным образом:

julia> ANewDsc(coeffs(q))
There were 3 isolating intervals found:
[2.62…, 3.62…]₂₅₆
[1.5…, 2.62…]₂₅₆
[-0.50…, 1.5…]₂₅₆

Подбор многочлена под произвольные данные

Функция fit подбирает многочлен (степени deg) под данные x и y, используя многочленную интерполяцию или (взвешенную) аппроксимацию по методу наименьших квадратов.

Подбирает многочлен (степени deg или меньше) под x и y, используя аппроксимацию по методу наименьших квадратов.

julia> xs = 0:4; ys = @. exp(-xs) + sin(xs);

julia> p =  fit(xs, ys); map(x -> round(x, digits=4), p)
Polynomial(1.0 + 0.0593*x + 0.3959*x^2 - 0.2846*x^3 + 0.0387*x^4)

julia> p = fit(ChebyshevT, xs, ys, 2); map(x -> round(x, digits=4), p)
ChebyshevT(0.5413⋅T_0(x) - 0.8991⋅T_1(x) - 0.4238⋅T_2(x))

Вот наглядный пример:

using Plots, Polynomials

xs = range(0, 10, length=10)
ys = @. exp(-xs)

f = fit(xs, ys)     # степень = length(xs) - 1
f2 = fit(xs, ys, 2) # степень = 2

scatter(xs, ys, markerstrokewidth=0, label="Data")
plot!(f, extrema(xs)..., label="Interpolation")
plot!(f2, extrema(xs)..., label="Quadratic Fit")

Другие базисы

Многочлен, например a_0 + a_1 x + a_2 x^2 + ... + a_n x^n, можно представить как набор коэффициентов [a_0, a_1, ..., a_n] относительно некоторого полиномиального базиса. Самым привычным базисом является стандартный: 1, x, x^2, … Возможны и другие базисы. Например, реализованы многочлены ChebyshevT. Конструктором является ChebyshevT — псевдоним для MutableDensePolynomial{ChebyshevTBasis}.

julia> p1 = ChebyshevT([1.0, 2.0, 3.0])
ChebyshevT(1.0⋅T_0(x) + 2.0⋅T_1(x) + 3.0⋅T_2(x))

julia> p2 = ChebyshevT{Float64}([0, 1, 2])
ChebyshevT(1.0⋅T_1(x) + 2.0⋅T_2(x))

julia> p1 + p2
ChebyshevT(1.0⋅T_0(x) + 3.0⋅T_1(x) + 5.0⋅T_2(x))

julia> p1 * p2
ChebyshevT(4.0⋅T_0(x) + 4.5⋅T_1(x) + 3.0⋅T_2(x) + 3.5⋅T_3(x) + 3.0⋅T_4(x))

julia> derivative(p1)
ChebyshevT(2.0⋅T_0(x) + 12.0⋅T_1(x))

julia> integrate(p2)
ChebyshevT(- 1.0⋅T_1(x) + 0.25⋅T_2(x) + 0.3333333333333333⋅T_3(x))

julia> convert(Polynomial, p1)
Polynomial(-2.0 + 2.0*x + 6.0*x^2)

julia> convert(ChebyshevT, Polynomial([1.0, 2,  3]))
ChebyshevT(2.5⋅T_0(x) + 2.0⋅T_1(x) + 1.5⋅T_2(x))

Итерация

Если базис является неявным, то многочлен можно рассматривать просто как вектор коэффициентов. Векторы отсчитываются от 1, но для удобства при индексировании большинства типов многочленов (например, с помощью функций getindex, setindex!, eachindex) отсчет начинается с 0. При итерировании по многочлену выполняется проход по базовым коэффициентам.

julia> as = [1,2,3,4,5]; p = Polynomial(as);

julia> as[3], p[2], collect(p)[3]
(3, 3, 3)

Итератор pairs перебирает индексы и коэффициенты, пытаясь определить то, как pairs применяется к базовой модели хранения:

julia> v = [1,2,0,4]
4-element Vector{Int64}:
 1
 2
 0
 4

julia> p,ip,sp,lp = Polynomial(v), ImmutablePolynomial(v), SparsePolynomial(v), LaurentPolynomial(v, -1);

julia> collect(pairs(p))
4-element Vector{Pair{Int64, Int64}}:
 0 => 1
 1 => 2
 2 => 0
 3 => 4

julia> collect(pairs(ip)) == collect(pairs(p))
true

julia> collect(pairs(sp)) # неупорядоченный словарь, содержащий только ненулевые члены
3-element Vector{Pair{Int64, Int64}}:
 0 => 1
 3 => 4
 1 => 2

julia> collect(pairs(lp))
4-element Vector{Pair{Int64, Int64}}:
 -1 => 1
  0 => 2
  1 => 0
  2 => 4

Неэкспортируемый итератор monomials перебирает члены (p[i]*Polynomials.basis(p,i)) многочлена:

julia> p = Polynomial([1,2,0,4], :u)
Polynomial(1 + 2*u + 4*u^3)

julia> collect(Polynomials.monomials(p))
4-element Vector{Any}:
 Polynomial(1)
 Polynomial(2*u)
 Polynomial(0)
 Polynomial(4*u^3)

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

julia> p = Polynomial([24, -50, 35, -10, 1])
Polynomial(24 - 50*x + 35*x^2 - 10*x^3 + x^4)

julia> q = convert(FactoredPolynomial, p) # неточная форма `factor`:
FactoredPolynomial((x - 4.0000000000000036) * (x - 2.9999999999999942) * (x - 1.0000000000000002) * (x - 2.0000000000000018))

julia> map(x -> round(x, digits=10), q)
FactoredPolynomial((x - 4.0) * (x - 2.0) * (x - 3.0) * (x - 1.0))

Тип элементов

Связь между `T` и `P{T,X}`

Сложение многочлена и скалярного значения, например

julia> using Polynomials

julia> p = Polynomial([1,2,3], :x)
Polynomial(1 + 2*x + 3*x^2)

julia> p + 3
Polynomial(4 + 2*x + 3*x^2)

кажется естественным, однако поскольку в Julia значение 3 имеет тип Int, а p — тип Polynomial{Int,:x}, некоторые особенности сложения необходимо уточнить. Основная идея в том, что значение 3 продвигается до константного многочлена 3 с неопределенной величиной :x в виде 3*one(p), а затем производится сложение p + 3*one(p).

Такое определение скалярного значения как константного многочлена работает и в другую сторону. Если q — это константный многочлен типа Polynomial{Int, :y}, то следует ожидать, что определена операция p+q, означающая p плюс свободный член q. И это действительно так:

julia> q = Polynomial(3, :y)
Polynomial(3)

julia> p + q
Polynomial(4 + 2*x + 3*x^2)

Если многочлен q не является константным, например variable(Polynomial, :y), то произойдет ошибка из-за несоответствия символов. (Математический результат должен быть многомерным многочленом, а не одномерным, как в случае с данным пакетом.)

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

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

julia> one(Polynomial, :x) + one(Polynomial, :y)
Polynomial(2.0)

julia> one(Polynomial, :y) + one(Polynomial, :x)
Polynomial(2.0)

Хотя оба многочлена являются константными типа Int, в первом неизвестная величина — :y, а во втором — :x.

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

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

julia> p = Polynomial([1,2],:x)
Polynomial(1 + 2*x)

julia> q = Polynomial([1,2],:y)  # неконстантные многочлены с разными неизвестными величинами
Polynomial(1 + 2*y)

julia> [1 p]
1×2 Matrix{Polynomial{Int64, :x}}:
 Polynomial(1)  Polynomial(1 + 2*x)

julia> [1 one(q)]
1×2 Matrix{Polynomial{Int64, :y}}:
 Polynomial(1)  Polynomial(1)

Однако поскольку в следующем случае результат будет неоднозначным,

julia> [one(p) one(q)]
ERROR: ArgumentError: Polynomials have different indeterminates
[...]

возникнет ошибка.

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

При использовании one(q) для константного многочлена с неизвестной величиной :y получается следующее.

julia> P = typeof(p)
Polynomial{Int64, :x} (alias for Polynomials.MutableDensePolynomial{Polynomials.StandardBasis, Int64, :x})

julia> P[one(p) one(q)]
1×2 Matrix{Polynomial{Int64, :x}}:
 Polynomial(1)  Polynomial(1)

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

julia> [one(p), one(q)+one(p)]
2-element Vector{Polynomial{Int64, :x}}:
 Polynomial(1)
 Polynomial(2)

но не такую матрицу:

julia> [one(p), one(p) + one(q)]
ERROR: ArgumentError: Polynomials have different indeterminates
[...]

Кроме того, смешение типов может приводить к неопределенным символам, как показано в этом примере:

julia> [1 p; p 1] + [1 2one(q); 3 4] # array{P{T,:x}} + array{P{T,:y}}
2×2 Matrix{Polynomial{Int64}}:
 Polynomial(2)        Polynomial(3 + 2*x)
 Polynomial(4 + 2*x)  Polynomial(5)

Хотя если в примере выше заменить 2one(q) неконстантным многочленом с неизвестной величиной y, сложение выдаст ошибку.

Нечисловые типы для `T`

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

Например, многочлен с матричными коэффициентами можно создать так:

julia> using Polynomials

julia> a,b,c = [1 0;2 1], [1 0; 3 1], [1 0; 4 1]
([1 0; 2 1], [1 0; 3 1], [1 0; 4 1])

julia> p = Polynomial([a,b,c])
Polynomial([1 0; 2 1] + [1 0; 3 1]*x + [1 0; 4 1]*x^2)

julia> q = derivative(p)
Polynomial([1 0; 3 1] + [2 0; 8 2]*x)

Доступны различные операции. Для примера выше была показана операция derivative, а вот операции в векторном пространстве:

julia> 2p
Polynomial([2 0; 4 2] + [2 0; 6 2]*x + [2 0; 8 2]*x^2)

julia> p + q
Polynomial([2 0; 5 2] + [3 0; 11 3]*x + [1 0; 4 1]*x^2)

Умножение многочленов:

julia> p * q
Polynomial([1 0; 5 1] + [3 0; 18 3]*x + [3 0; 21 3]*x^2 + [2 0; 16 2]*x^3)

Вычисление многочлена, в данном случае со скалярным значением или матрицей:

julia> p(2)
2×2 Matrix{Int64}:
  7  0
 24  7

julia> p(b)
2×2 Matrix{Int64}:
  3  0
 18  3

Однако если тип T не поддерживает некоторые универсальные функции, такие как zero(T) и one(T), возможны проблемы. Например, при T <: AbstractMatrix результатом вызова p[degree(p)+1] будет ошибка, так как реализация предполагает, что определено zero(T). Для статических массивов это не проблема, так как они поддерживают zero(T). Другие типы многочленов, такие как SparsePolynomial, имеют более ограниченную поддержку, так как некоторые специализированные методы предполагают более полную реализацию универсального интерфейса.

Аналогичным образом, в качестве T можно использовать многочлены:

julia> a,b,c = Polynomial([1],:y), Polynomial([0,1],:y), Polynomial([0,0,1],:y)
(Polynomial(1), Polynomial(y), Polynomial(y^2))

julia> p = Polynomial([a,b,c], :x)
Polynomial(Polynomial(1) + Polynomial(y)*x + Polynomial(y^2)*x^2)

julia> q = derivative(p)
Polynomial(Polynomial(y) + Polynomial(2*y^2)*x)

Опять-таки, многие возможности доступны:

julia> 2p
Polynomial(Polynomial(2) + Polynomial(2*y)*x + Polynomial(2*y^2)*x^2)

julia> p + q
Polynomial(Polynomial(1 + y) + Polynomial(y + 2*y^2)*x + Polynomial(y^2)*x^2)

julia> p(2)
Polynomial(1 + 2*y + 4*y^2)

julia> p(b)
Polynomial(1 + y^2 + y^4)

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

Рациональные функции

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

julia> P = FactoredPolynomial
FactoredPolynomial

julia> p,q = fromroots(P, [1,2,3,4]), fromroots(P, [2,2,3,5])
(FactoredPolynomial((x - 4) * (x - 2) * (x - 3) * (x - 1)), FactoredPolynomial((x - 5) * (x - 2)² * (x - 3)))

julia> pq = p // q
((x - 4) * (x - 2) * (x - 3) * (x - 1)) // ((x - 5) * (x - 2)² * (x - 3))

julia> lowest_terms(pq)
((x - 4.0) * (x - 1.0)) // ((x - 5.0) * (x - 2.0))

julia> d,r = residues(pq); r
Dict{Float64, Vector{Float64}} with 2 entries:
  5.0 => [1.33333]
  2.0 => [0.666667]

julia> x = variable(p);

julia> for (λ, rs) ∈ r # воссоздаем p/q на основе результата `residues`
           for (i,rᵢ) ∈ enumerate(rs)
               d += rᵢ//(x-λ)^i
           end
       end


julia> d
((x - 4.0) * (x - 1.0000000000000002)) // ((x - 5.0) * (x - 2.0))

Предоставляется базовый шаблон графика.

using Plots, Polynomials
P = FactoredPolynomial
p,q = fromroots(P, [1,2,3]), fromroots(P, [2,3,3,0])
plot(p//q)

Связанные пакеты

StaticUnivariatePolynomials.jl — одномерные многочлены фиксированного размера на основе кортежа.
MultiPoly.jl — разреженные многомерные многочлены.
DynamicPolynomals.jl — реализация многомерных многочленов с коммутативными и некоммутативными переменными.
MultivariatePolynomials.jl — многомерные многочлены и моменты коммутативных и некоммутативных переменных.
PolynomialRings.jl — библиотека для арифметических и алгебраических операций с многочленами нескольких переменных.
AbstractAlgebra.jl, Nemo.jl для универсальных колец многочленов, матричных пространств, полей частных, колец вычетов, степенных рядов; Hecke.jl для алгебраической теории чисел.
LaurentPolynomials.jl — пакет для работы с многочленами Лорана.
CommutativeAlgebra.jl — базовая система компьютерной алгебры, предназначенная для дискретных вычислений с поддержкой многочленов.
PolynomialRoots.jl — быстрое нахождение комплексных корней многочленов. Для задач более высокой степени также можно использовать пакеты FastPolynomialRoots и AMRVW.
SpecialPolynomials.jl — пакет, предоставляющий различные типы многочленов помимо многочленов в стандартном базисе, доступных в Polynomials.jl. Включает в себя интерполяционные многочлены, многочлены Бернштейна и классические ортогональные многочлены.
ClassicalOrthogonalPolynomials.jl — пакет Julia для работы с классическими ортогональными многочленами и расширениями. Включает в себя многочлены chebyshevt, chebyshevu, legendrep, jacobip, ultrasphericalc, hermiteh и laguerrel. В том же репозитории доступны пакеты FastGaussQuadrature.jl, FastTransforms.jl и ApproxFun.

Сотрудничество

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