I would like to implement a Bernstein polynomial class in Julia for learning purposes. The goal in the end is to implement the quadratic clipping algorithm from this paper. I am completely new to Julia coming from Python + Numpy + Numba, so I want to make sure that I program in the canonical Julia way when using Julia.

The naming of variables follows the aforementioned paper.

using Polynomials

type Bernstein
Bernstein(n, i) = Bernstein(n, i, 0., 1.)

Base.string(B::Bernstein) = "Bernstein $(get_p(B))"
Base.show(io::IO, B::Bernstein) = print(io, string(B))
Base.display(io::IO, B::Bernstein) = print(io, string(B))

function get_p(n, i, α=0, β=1)
    @assert β >= α "Invalid interval"
    return (binomial(n, i)
            * poly([α for _ in 1:i])
            * poly([β for _ in 1:(n - i)])
            * (-1)^(n - i)
            / (β - α)^n)
get_p(B::Bernstein) = get_p(B.n, B.i, B.α, B.β)

function ∫(p::Poly, α=NaN, β=NaN)
    integrated = polyint(p)
    if isequal(α, NaN) || isequal(β, NaN)
        return integrated
        return integrated(β) - integrated(α)

function inner_product(p::Poly, q::Poly, α=0., β=1.)
    @assert β >= α "Invalid interval"
    ∫(p * q, α, β)

function inner_product(B::Bernstein, Q::Bernstein)
    @assert B.α >= Q.α "Invalid interval"
    @assert B.β >= Q.β "Invalid interval"
    m, i, n, j, α, β = B.n, B.i, Q.n, Q.i, B.α, B.β
    (β - α)binomial(m, i)binomial(n, j) / ((m + n + 1)binomial(m + n, i + j))

# QUESTION: Is it possible to automatically account for commutativity?
(inner_product(p::Poly, B::Bernstein)
 = inner_product(B::Bernstein, p::Poly)
 = inner_product(B.p, p::Poly))

function Base.LinAlg.norm(q::Poly, α=0., β=1.)
    @assert β >= α "Invalid interval"
    √inner_product(q, q) / (β - α)

function Base.LinAlg.norm(B::Bernstein)
    √inner_product(B, B) / (B.β - B.α)

B = Bernstein(3, 2)
Q = Bernstein(4, 3)

inner_product(B, Q)

1 Answer 1


The below code uses Julia 1.0. Full code as a Pkg3 project is available here.


There's enough code for it to get its own module. This also makes testing easier.

module Bernstein

using Polynomials

import Base: convert, promote_rule, show
import LinearAlgebra: dot, norm
import Polynomials: poly

Since we're in a module, we must define the exported things. Since most of the things is just overloaded methods, only the type remains.

export BernsteinPoly

You do this check so often, we can factor it out. Alternative: using a separate type for intervals.

macro checkinterval(α, β)
    :(@assert $(esc(α)) <= $(esc(β)) "Invalid interval")

Parametrization by a constrained type is better than fields of an abstract type. But in the case of n and i, concrete Ints make more sense, since they represent natural numbers always (this makes conversions below easier). And I'm calling it BernsteinPoly, since it is a variant of Poly.

struct BernsteinPoly{T<:Number}

BernsteinPoly(n, i) = BernsteinPoly(n, i, 0.0, 1.0)

Just implementing show using print is the recommended way for custom pretty-printing. Here, we can reuse Polynomials.printpoly for nicer and more consistent formatting:

function show(io::IO, b::BernsteinPoly)
    print(io, "BernsteinPoly(")
    Polynomials.printpoly(io, poly(b))
    print(io, ")")


get_p was essentially the conversion from BernsteinPoly to Poly. We can replace this by the implementation of proper convert methods. The default parameters are unnecessary, since already occuring defaulted in the BernsteinPoly constructor.

function convert(::Type{Poly{S}}, b::BernsteinPoly) where {S<:Number}
    n, i, α, β = b.n, b.i, convert(S, b.α), convert(S, b.β)
    @checkinterval α β

    return (binomial(n, i)
            * poly(fill(α, i))     # `fill` instead of list comprehension
            * poly(fill(β, n - i))
            * (-1)^(n - i)
            / (β - α)^n)

This allows you to just say convert(Poly, b), automatically reusing the inner type of b

convert(::Type{Poly}, b::BernsteinPoly{T}) where {T<:Number} =
    convert(Poly{T}, b)

While we're at it: conversion between different BernsteinPoly values

convert(::Type{BernsteinPoly{S}}, b::BernsteinPoly) where {S<:Number} =
    BernsteinPoly(b.n, b.i, convert(S, b.α), convert(S, b.β))

If we're handling different representations, we sometimes need to determine the "most general" form, which is called promotion:

promote_rule(a::Type{BernsteinPoly{S}}, b::Type{Poly{T}}) where {S<:Number, T<:Number} =
    Poly{promote_type(S, T)}

promote_rule(::Type{BernsteinPoly{S}}, ::Type{BernsteinPoly{T}}) where {S<:Number, T<:Number} = 
    BernsteinPoly{promote_type(S, T)}

Also add a method to the poly smart constuctor, which is now trivial:

poly(b::BernsteinPoly) = convert(Poly, b)

Linear Algebra

Now to the linear algebra part. We could add methods to dot and norm from LinearAlgebra:

function dot(p::Poly{T}, q::Poly{T}, α = zero(T), β = one(T)) where {T<:Number}
    @checkinterval α β
    polyint(p * q, α, β)

And the norm induced by that inner product:

function norm(q::Poly{T}, α = zero(T), β = one(T)) where {T<:Number}
    @checkinterval α β
    √dot(q, q, α, β) / (β - α)

But such "overloads from outside" are frowned upon. As of writing this, there's a norm method in Polynomials, but not dot. There's an issue about that; basically, dot is not unique, hence we need to specify the intervals each time here.

For Bernstein polynomials, on the other hand, the inner product is defined uniquely, if I understood correctly.

function dot(b::BernsteinPoly{T}, q::BernsteinPoly{T}) where {T<:Number}
    @checkinterval b.α q.α
    @checkinterval b.β q.β
    m, i, n, j, α, β = b.n, b.i, q.n, q.i, b.α, b.β
    (β - α) * binomial(m, i) * binomial(n, j) / ((m + n + 1) * binomial(m + n, i + j))

And the induced norm, as before.

function norm(b::BernsteinPoly)
    √dot(b, b) / (b.β - b.α)

As the last remaining question, how to do cross-type inner products. If dot is defined for Poly as written above, we can use dot(promote(p, q)..., α, β). But it's difficult to get that working as a method. I tried

dot(p::Union{Poly{T}, BernsteinPoly{T}}, q::Union{Poly{T}, BernsteinPoly{T}},
     α = zero(T), β = one(T)) where {T<:Number} =
     dot(promote(p, q)..., α, β)

but that doesn't work. It's probably not recommended to do that anyway. Use promote explicitely where necessary.

Example output

julia> B = BernsteinPoly(3, 2)
BernsteinPoly(3.0*x^2 - 3.0*x^3)

julia> Q = BernsteinPoly(4, 3)
BernsteinPoly(4.0*x^3 - 4.0*x^4)

julia> using LinearAlgebra

julia> dot(B, Q)

julia> norm(B)

end # module

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.