Расширение многочленов

Тип AbstractUnivariatePolynomial предусматривает возможность расширения.

Коэффициенты многочлена соотносятся с каким-либо базисом. Тип Polynomial связывает коэффициенты [a0, a1, ..., an], допустим, с многочленом , посредством стандартного базиса . Новые типы многочленов обычно представляют многочлен через другой базис. Например, CheyshevT использует базис . Для этого типа коэффициенты [a0,a1,...,an] связаны с многочленом `a0\cdot T_0 + a_1 \cdot T_1 + \cdots + a_n\cdot T_n.

Тип-многочлен состоит из типа контейнера (с родительским типом AbstractUnivariatePolynomial) и типа базиса (с родительским типом AbstractBasis). Реализовано несколько разных типов для хранения данных.

Для реализации нового типа-многочлена P следует реализовать следующие методы:

Функция Обязательная Примечания

Функция	Обязательная	Примечания
Тип контейнера	x	Обычно выбирается из числа доступных.
Тип базиса	x
`variable`		Удобная функция для нахождения одночлена `x` в новом базисе.
`Base.evalpoly(x, p::P)`	x	Для вычисления многочлена в точке `x`
`*(::P, ::P)`		Умножение многочленов
`convert(::P, p::Polynomial)`	По умолчанию соответствует вычислению многочлена. Может использоваться для определения `*` путем кругового обхода типа `Polynomial`
`convert(::Polynomial, p)`		По умолчанию соответствует вычислению многочлена с использованием `evalpoly`, `variable`, `*`
`scalar_add(c::S, ::P)`		Скалярное сложение. По умолчанию требует определения `one`.
`one`	x	Удобная функция для нахождения константы в новом базисе.
`map(f, p)`	x	Используется для определения скалярного умножения
`divrem`		Требуется для `gcd`
`vander`		Требуется для `fit`
`companion`		Требуется для `roots`
`Polynomials.domain`		Должна возвращать экземпляр `Polynomials.Interval`

Тип контейнера

Обычно выбирается из числа доступных.

Тип базиса

variable

Удобная функция для нахождения одночлена x в новом базисе.

Base.evalpoly(x, p::P)

Для вычисления многочлена в точке x

*(::P, ::P)

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

convert(::P, p::Polynomial)

По умолчанию соответствует вычислению многочлена. Может использоваться для определения * путем кругового обхода типа Polynomial

convert(::Polynomial, p)

По умолчанию соответствует вычислению многочлена с использованием evalpoly, variable, *

scalar_add(c::S, ::P)

Скалярное сложение. По умолчанию требует определения one.

one

Удобная функция для нахождения константы в новом базисе.

map(f, p)

Используется для определения скалярного умножения

divrem

Требуется для gcd

vander

Требуется для fit

companion

Требуется для roots

Polynomials.domain

Должна возвращать экземпляр Polynomials.Interval

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

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

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

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

Большинство правил продвижения выполняют приведение к типу Polynomial, поэтому если предоставлена функция преобразования, реализовывать нужно не все методы.

Новый тип базиса

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

Доступны и другие определения. Ниже для вычисления применяется трехточечная рекурсия, описываемая посредством A, B и C.

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

Сначала мы загружаем пакет и импортируем несколько неэкспортированных функций:

julia> using Polynomials;

julia> import Polynomials: AbstractUnivariatePolynomial, AbstractBasis, MutableDensePolynomial;

Базис определяется следующим образом:

julia> struct LaguerreBasis{alpha} <: AbstractBasis end

julia> Polynomials.basis_symbol(::Type{<:AbstractUnivariatePolynomial{LaguerreBasis{α},T,X}}) where {α,T,X} =
           "L^$(α)"

Для вывода этого базиса мы добавили метод basis_symbol. Способ отображения символа базиса по умолчанию не самый подходящий. Приведенный выше метод требует полного типа, так как неопределенная величина X может входить в выводимые данные. Говоря шире, для Polynomials.printbasis могут быть добавлены методы для учета разных отображаемых типов.

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

julia> P = MutableDensePolynomial{LaguerreBasis{0}}
MutableDensePolynomial{LaguerreBasis{0}}

Теперь можно создавать экземпляры:

julia> p = P([1,2,3])
MutableDensePolynomial(1L^0_0 + 2*L^0_1 + 3*L^0_2)

Или использовать другие типы для хранения данных:

julia> Polynomials.ImmutableDensePolynomial{LaguerreBasis{1}}((1,2,3))
Polynomials.ImmutableDensePolynomial(1L^1_0 + 2*L^1_1 + 3*L^1_2)

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

julia> q = P([1,2])
MutableDensePolynomial(1L^0_0 + 2*L^0_1)

julia> p + q
MutableDensePolynomial(2L^0_0 + 4*L^0_1 + 3*L^0_2)

julia> 2p
MutableDensePolynomial(2L^0_0 + 4*L^0_1 + 6*L^0_2)

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

Для вычисления многочленов Лагерра можно использовать рекурсию методом Кленшо.

julia> function ABC(::Type{LaguerreBasis{α}}, n) where {α}
           o = one(α)
           d = n + o
           (A=-o/d, B=(2n + o + α)/d, C=(n+α)/d)
       end
ABC (generic function with 1 method)

julia> function clenshaw_eval(p::P, x::S) where {α, Bᵅ<: LaguerreBasis{α}, T, P<:AbstractUnivariatePolynomial{Bᵅ,T}, S}
           d = degree(p)
           R = typeof(((one(α) * one(T)) * one(S)) / 1)
           p₀ = one(R)
           d == -1 && return zero(R)
           d == 0 && return p[0] * one(R)
           Δ0 = p[d-1]
           Δ1 = p[d]
           @inbounds for i in (d - 1):-1:1
               A,B,C = ABC(Bᵅ, i)
               Δ0, Δ1 =
                   p[i] - Δ1 * C, Δ0 + Δ1 * muladd(x, A, B)
           end
           A,B,C = ABC(Bᵅ, 0)
           p₁ = muladd(x, A, B) * p₀
           return Δ0 * p₀ + Δ1 * p₁
       end
clenshaw_eval (generic function with 1 method)

На внутреннем уровне вызывается метод evalpoly, поэтому мы передаем его.

julia> Polynomials.evalpoly(x, p::P) where {P<:AbstractUnivariatePolynomial{<:LaguerreBasis}} =
               clenshaw_eval(p, x)

Для тестирования мы передаем переменную x в стандартном базисе:

julia> p = P([0,0,1])
MutableDensePolynomial(L^0_2)

julia> x = variable(Polynomial)
Polynomial(1.0*x)

julia> p(x)
Polynomial(1.0 - 2.0*x + 0.5*x^2)

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

julia> convert(ChebyshevT, p)
ChebyshevT(1.25⋅T_0(x) - 2.0⋅T_1(x) + 0.25⋅T_2(x))

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

julia> q = Polynomials.basis(MutableDensePolynomial{LaguerreBasis{0//1}, Int}, 5)
MutableDensePolynomial(L^0//1_5)

julia> x = variable(Polynomial{Int})
Polynomial(x)

julia> q(x)
Polynomial(1//1 - 5//1*x + 5//1*x^2 - 5//3*x^3 + 5//24*x^4 - 1//120*x^5)

Значения one и variable определяются очень просто, как и или :

julia> Polynomials.one(::Type{P}) where {B<:LaguerreBasis,T,X,P<:AbstractUnivariatePolynomial{B,T,X}} =
           P([one(T)])

julia> Polynomials.variable(::Type{P}) where {B<:LaguerreBasis,T,X,P<:AbstractUnivariatePolynomial{B,T,X}} =
           P([one(T), -one(T)])

Проверить правильность можно так:

julia> variable(P)(x) == x
true

Для скалярного сложения по умолчанию вызывается one(p), а такой вызов уже определен:

julia> 2 + p
MutableDensePolynomial(2L^0_0 + L^0_2)

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

julia> function Polynomials.scalar_add(c::S, p::P) where {B<:LaguerreBasis,T,X,
                                                          P<:MutableDensePolynomial{B,T,X},S}
           R = promote_type(T,S)
           iszero(p) && return MutableDensePolynomial{B,R,X}(c)
           cs = convert(Vector{R}, copy(p.coeffs))
           cs[1] += c
           MutableDensePolynomial{B,R,X}(cs)
       end

julia> p + 3
MutableDensePolynomial(3L^0_0 + L^0_2)

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

julia> function Base.:*(p::MutableDensePolynomial{B,T,X},
                        q::MutableDensePolynomial{B,S,X}) where {B<:LaguerreBasis, T,S,X}
           x = variable(Polynomial{T,X})
           p(x) * q(x)
       end

Если бы был определен метод convert для преобразования из Polynomial в LaguerreBasis, его можно было бы использовать для реализации умножения, так как мы определили метод variable.

Новый тип контейнера

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

У нас есть два метода конструктора. Первый представляет собой обычный путь выполнения кода. Он создает копию коэффициентов, а затем заключает их в тип контейнера многочлена. В целях повышения производительности, как правило, полезно передавать флаг, указывающий, что копирование или проверка входных данных не требуется (Val{false}). Он используется некоторыми унаследованными методами при специализации для типа StandardBasis. Обычно тип контейнера может принимать смещение, но в нашем случае это не так; неявно используется отсчитываемый от 0 вектор.

julia> using Polynomials

julia> struct AliasPolynomialType{B,T,X} <: Polynomials.AbstractDenseUnivariatePolynomial{B, T, X}
           coeffs::Vector{T}
           function AliasPolynomialType{B, T, X}(coeffs::AbstractVector{S}, o::Int=0) where {B, T, S, X}
               new{B,T,Symbol(X)}(convert(Vector{T}, copy(coeffs)))
           end
           function AliasPolynomialType{B, T, X}(::Val{false}, coeffs::AbstractVector{S}, o::Int=0) where {B, T, S, X}
               new{B,T,Symbol(X)}(convert(Vector{T}, coeffs))
           end
       end

julia> Polynomials.@poly_register AliasPolynomialType

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

Чтобы было возможным индексирование, необходимо определить несколько методов:

julia> Base.firstindex(p::AliasPolynomialType) = 0

julia> Base.lastindex(p::AliasPolynomialType) = length(p.coeffs) - 1

julia> Polynomials.constructorof(::Type{<:AliasPolynomialType{B}}) where {B} = AliasPolynomialType{B}

Необходимо также добавить операции в векторном пространстве:

julia> function Base.:+(p::AliasPolynomialType{B,T,X}, q::AliasPolynomialType{B,S,X}) where {B,S,T,X}
                       R = promote_type(T,S)
                        n = maximum(degree, (p,q))
                   cs = [p[i] + q[i] for i in 0:n]
                   AliasPolynomialType{B,R,X}(Val(false), cs)  # сохраняем копию
                       end

julia> function Base.:-(p::AliasPolynomialType{B,T,X}, q::AliasPolynomialType{B,S,X}) where {B,S,T,X}
                       R = promote_type(T,S)
                               n = maximum(degree, (p,q))
                               cs = [p[i] - q[i] for i in 0:n]
                               AliasPolynomialType{B,R,X}(Val(false), cs)
                       end

julia> function Base.map(fn, p::P) where {B,T,X,P<:AliasPolynomialType{B,T,X}}
                  cs = map(fn, p.coeffs)
                  R = eltype(cs)
                  AliasPolynomialType{B,R,X}(Val(false), cs)
              end

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

julia> AliasPolynomial = AliasPolynomialType{Polynomials.StandardBasis};

Посмотрим на этот новый тип многочлена в действии:

julia> xs = [1,2,3,4];

julia> p = AliasPolynomial(xs)
AliasPolynomialType(1 + 2*x + 3*x^2 + 4*x^3)

julia> q = AliasPolynomial(1.0, :y)
AliasPolynomialType(1.0)

julia> 2p - q
AliasPolynomialType(3.0 + 4.0*x + 6.0*x^2 + 8.0*x^3)

julia> (derivative ∘ integrate)(p) == p
true

julia> p(3)
142

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

julia> Base.:*(p::AliasPolynomialType{T,X}, q::AliasPolynomialType{S,X}) where {T,S,X} = Polynomials._standard_basis_multiplication(p,q)

julia> p * p
AliasPolynomialType(1 + 4*x + 10*x^2 + 20*x^3 + 25*x^4 + 24*x^5 + 16*x^6)

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

Для скалярного сложения p+c по умолчанию производится операция p + c*one(p), то есть сложение многочленов, которое не выполняется на месте, если не принять дополнительные меры. Поэтому мы создаем новый метод и инфиксный оператор:

julia> function scalar_add!(c::T, p::AliasPolynomial{T}) where {T}
           p.coeffs[1] += c
           p
       end;

julia> p::AliasPolynomial +ₛ c::Number = scalar_add!(c, p);

julia> c::Number +ₛ p::AliasPolynomial = scalar_add!(c, p);

Ввиду того факта, что многочлен представляет вектор коэффициентов, следует ожидать, что векторные операции должны быть возможны, где это имеет смысл. Скалярное умножение — это векторная операция, поэтому, по-видимому, есть смысл переопределить механизм трансляции для реализации выполняемой на месте операции (например, p .*= 2). По умолчанию типы многочленов не транслируются через свои коэффициенты. Потребуется внести изменение в функцию copyto!:

julia> Base.broadcastable(p::AliasPolynomial) = p.coeffs;

julia> Base.ndims(::Type{<:AliasPolynomial}) = 1

julia> Base.copyto!(p::AliasPolynomial, x) = (copyto!(p.coeffs, x); chop!(p));

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

julia> p .*= 2
AliasPolynomialType(2 + 4*x + 6*x^2 + 8*x^3)

julia> p ./= 2
AliasPolynomialType(1 + 2*x + 3*x^2 + 4*x^3)

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

Теперь при трансляции p рассматривается как вектор p.coeffs, поэтому возможны неожиданные последствия. Например, следующее выражение возвращает вектор, а не многочлен:

julia> p .+ 2
4-element Vector{Int64}:
 3
 4
 5
 6

Неэкспортируемый тип многочлена Polynomials.PnPolynomial реализует многие из описанных выше возможностей.