Julia: types, multimethods and arithmetic over polynomials

In this publication, we will discuss the main, in my opinion, distinctive feature of the Julia language - the representation of functions in the form of methods with multiple dispatching. This allows you to improve the performance of calculations without reducing the readability of the code and not spoiling the abstraction, on the one hand, and allows you to work with mathematical concepts in a more familiar notation, on the other. For example, we consider the question of uniform (from the point of view of linear operations) work with polynomials in the representation of a list of coefficients and with interpolation polynomials.

Basic syntax

A brief introduction for those who do not know. Julia is a scripting language, has a REPL (read-evaluate-print loop, i.e. interactive shell). At first glance it looks quite similar to, for example, Python or MATLAB.

Arithmetic operations

The arithmetic is about the same as elsewhere: +, -, *, /, ^ for exponentiation, etc.
Comparison:>, <,> =, <=, ==,! = Etc.
Assignment: =.
Features: division by / always returns a fractional number; if you need an integer part from dividing two integers, you need to use the operation div(m, n) or the infix equivalent m ÷ n .

Types

Numeric types:

• Integers ( Int ) - 2 , 3 , -42
• Unsigned UInt ( UInt ) - 0x12345
• Floating point ( Float32 , Float64 ) - 1.0 , 3.1415 , -Inf , NaN
• Rational ( Rational ) - 3//3 , 7//2
• Valid ( Real ) - All of the above
• Complex ( Complex ) - 3+4*im , 2//3+2//3*im , 3.0+0.0*im ( im is an imaginary unit, only a number with an explicitly written imaginary part is considered complex)
• Number ( Number ) - all of the above

Lines and characters:

• 'a' - character ( Char )
• "a" - string ( String )

NB: strings, as they are now in many languages, are immutable.
NB: strings (as well as variable names) support Unicode, including emoji.
')
Arrays:

• x = [1, 2, 3] - set the array by direct enumeration of elements
• special constructors: zeros(length) for an array of zeros, ones(length) for an array of ones, rand(length) for an array of random numbers, etc.
• support for multidimensional arrays
• support of linear algebra operations (addition of arrays, multiplication by a scalar, multiplication of a matrix by a vector, and much more) in the standard library

NB: all collections are indexed starting from one.
NB: because The language is intended for computational tasks, arrays are one of the most important types, and we have to return to the principles of their work more than once.

Tuples (ordered set of elements, immutable):

• (2, 5.3, "k") - regular tuple
• (a = 3, b = 4) - the named tuple

NB: the fields of the named tuple can be accessed both by name through a dot, and by index through []

 julia> x = (a = 5, b = 12) (a = 5, b = 12) julia> x[1] 5 julia> sqrt(xa^2 + x[2]^2) 13.0 

Dictionaries:

 julia> x = Dict('a' => 5, 'b' => 12) Dict{Char,Int64} with 2 entries: 'a' => 5 'b' => 12 julia> x['c'] = 13 13 julia> x Dict{Char,Int64} with 3 entries: 'a' => 5 'c' => 13 'b' => 12 

Basic language control structures

1. Variables are automatically created upon assignment. Type is optional.

 julia> x = 7; x + 2 9 julia> x = 42.0; x * 4 168.0 

2. The conditional transition block begins with the expression if <condition> and ends with the word end . You can also have an else flower or else if a flower:

 if x > y println("X is more than Y") elseif x == y println("X and Y are equal") else println("X is less than Y") end 

3. There are two loop constructs: while and for . The second works as in Python, i.e. iterates over the collection. Frequent use is iteration over a range of values, the syntax of which is start[:increment]:end . Unlike Python, the range includes both initial and final values, i.e. the empty range will not be 1:1 (this is a range of one value 1), but 1:0 . The end of the loop body is marked with the word end .

 julia> for i in 1:3; print(i, " "); end #   1  3   1 ( ) 1 2 3 julia> for i in 1:2:3; print(i, " "); end #   1  3   2 1 3 

4. Functions are defined by the function keyword, the function definition is also completed by the word end . Arguments with default values ​​and named arguments are supported.

 function square(x) return x * x end function cube(x) x * square(x) #       ; return   end function root(x, degree = 2) #  degree     return x^(1.0/degree) end function greeting(name; times = 42, greet = "hello") #       println(times, " times ", greet, " to ", name) end julia> greeting("John") 42 times hello to John julia> greeting("Mike", greet = "wassup", times = 100500) #           100500 times wassup to Mike 

In general, this is all quite similar to Python, except for small differences in the syntax and the fact that blocks of code are not allocated with spaces, but still keywords. In simple cases, Python programs are even translated to Julia almost one to one.
But there is a significant difference in that in Julia for variables you can explicitly specify types, which allows you to compile programs by getting a quick code.
The second significant difference is that Python implements the “classical” OOP model with classes and methods, and Julia implements the multiple dispatching model.

Type annotations and multiple dispatching

Let's see what a built-in function is:

 julia> sqrt sqrt (generic function with 19 methods) 

As the REPL shows us, sqrt is a generic function with 19 methods. What kind of generic function and what methods?

And this means that there are several sqrt functions that are applied to different types of arguments and, accordingly, they calculate the square root using different algorithms. See what options there are, by typing

 julia> methods(sqrt) 

It can be seen that the function is defined for different types of numbers, as well as for matrices.

Unlike the “classical” OOP, where a specific implementation of a method is determined only by the calling class (dispatching by the first argument), in Julia the choice of a function is determined by the types (and number) of all its arguments.
When calling a function with specific arguments, from all its methods, the one that most accurately describes the specific set of types with which the function is called is selected, and it is used.

A distinctive feature is that the approach used, called by the authors of the language “just ahead-of-time” compilation. Those. the functions are compiled for the given data types on the first call, after which the next calls are made much faster. The difference between the first and subsequent calls can be quite significant:

 julia> @time sqrt(8) #  @time -      0.006811 seconds (3.15 k allocations: 168.516 KiB) #   ,        2.8284271247461903 julia> @time sqrt(15) 0.000002 seconds (5 allocations: 176 bytes) # 5   -     @time 3.872983346207417 

In the bad case, each function call is a check of the types of the arguments received and a search for the desired method in the list. However, if you give hints to the compiler, checks can be eliminated, which will lead to faster code.

For example, consider calculating the sum

$$

$$\ sum_ {k = 1} ^ N \ sqrt {(- 1) ^ k}$$

 function mysqrt(num) #    -     #   -           if num >= 0 return sqrt(num) else return sqrt(complex(num)) end end function S(n) #    sum = 0 sgn = -1 for k = 1:n sum += mysqrt(sgn) sgn = -sgn end return sum end function S_typed(n::Integer) # ..     ,      #     sum::Complex = 0.0 sgn::Int = -1 for k = 1:n sum += mysqrt(sgn) sgn = -sgn end return sum end 

The benchmark shows that the S_typed() function not only executes faster, but also does not require memory allocations on every call, unlike S() . The problem here is that the type of the mysqrt() returned from mysqrt() is not defined, as is the type of the right side of the expression

 sum = sum + mysqrt(sgn) 

As a result, the compiler cannot even understand what type the sum will be at each iteration. So, boxing (hitch type labels) variable and memory allocation.
For the S_typed() function, the compiler knows in advance that the sum is a complex value, so the code is more optimized (in particular, the mysqrt() call can be effectively zainlaynit, leading the return value to always Complex ).

More importantly, for S_typed() compiler knows that the return value is of type Complex , but for S() type of output value is not defined again, which will slow down all functions where S() will be called.
You can check what the compiler thinks about the types returned from the expression using the macro @code_warntype :

 julia> @code_warntype S(3) Body::Any #     ,      ... julia> @code_warntype S_typed(3) Body::Complex{Float64} #      ... 

If a function is called for you somewhere in the loop, for which @code_warntype cannot deduce the returned type, or for which it is in the body somewhere, it shows the receipt of the Any value - then optimizing these calls will most likely give a very noticeable performance increase.

Composite types

A programmer can define composite data types for his needs with a struct :

 julia> struct GenericStruct #   struct    name b::Int c::Char v::Vector end #       #       ,        julia> s = GenericStruct("Name", 1, 'z', [3., 0]) GenericStruct("Name", 1, 'z', [3.0, 0.0]) julia> s.name, sb, sc, sv ("Name", 1, 'z', [3.0, 0.0]) 

Structures in Julia are immutable, that is, by creating a copy of the structure, you cannot change field values ​​(more precisely, you cannot change the address of fields in the memory — elements of mutable fields, such as sv in the example above, can be changed). Mutable structures are created by the mutable struct construction, the syntax of which is the same as for ordinary structures.

Inheritance of structures in the “classical” sense is not supported, but there is a possibility of “inheriting” behavior by combining composite types into supertypes or, as they are called in Julia, abstract types. Type relationships are expressed as A<:B (A - subtype B) and A>:B (A - supertype B). It looks like this:

 abstract type NDimPoint end #   -     # ,    -     N  struct PointScalar<:NDimPoint x1::Real end struct Point2D<:NDimPoint x1::Real x2::Real end struct Point3D<:NDimPoint x1::Real x2::Real x3::Real end #     ;   Markdown """ mag(p::NDimPoint) Calculate the magnitude of the radius vector of an N-dimensional point `p` """ function mag(p::NDimPoint) sqrmag = 0.0 # ..   ,       #     T   fieldnames(T) for name in fieldnames(typeof(p)) sqrmag += getfield(p, name)^2 end return sqrt(sqrmag) end """ add(p1::T, p2::T) where T<:NDimPoint Calculate the sum of the radius vectors of two N-dimensional points `p1` and `p2` """ function add(p1::T, p2::T) where T<:NDimPoint #  -  , ..       #     list comprehension sumvector = [Float64(getfield(p1, name) + getfield(p1, name)) for name in fieldnames(T)] #     ,    #  ...      , .. # f([1, 2, 3]...) -   ,  f(1, 2, 3) return T(sumvector...) end 

Case study: Polynomials

The type system, coupled with multiple dispatching, is convenient for expressing mathematical concepts. Consider the example of a simple library for working with polynomials.
We introduce two types of polynomials: “canonical”, defined by coefficients with powers, and “interpolation”, defined by a set of pairs (x, f (x)). For simplicity, we will consider only valid arguments.

A structure that has an array or a tuple of coefficients as a field is suitable for storing a polynomial in a regular record. So that it is completely immutable, let it be a tuple. Thus, the code for defining an abstract type, the structure of a polynomial and calculating the value of a polynomial at a given point is quite simple:

 abstract type AbstractPolynomial end """ Polynomial <: AbstractPolynomial Polynomials written in the canonical form """ struct Polynomial<:AbstractPolynomial degree::Int coeff::NTuple{N, Float64} where N # NTuple{N, Type} -    N    end """ evpoly(p::Polynomial, z::Real) Evaluate polynomial `p` at `z` using the Horner's rule """ function evpoly(p::Polynomial, z::Real) ans = p.coeff[end] for idx = p.degree:-1:1 ans = p.coeff[idx] + z * ans end return ans end 

Interpolation polynomials require a different representation structure and method of calculation. In particular, if the set of interpolation points is known in advance, and the same polynomial is planned to be calculated at different points, Newton's interpolation formula is convenient:

$$

$$P (x) = \ sum_ {k = 0} ^ N {c_k n_k (x)},$$

where n _k ( x ) are basic polynomials, n ₀ ( x ) and for k > 0

$$

$$n_k (x) = \ prod_ {i = 0} ^ {k-1} {(x-x_i)},$$

where x _i - interpolation nodes.

From the above formulas it is clear that storage is conveniently organized as a set of interpolation nodes x _i and coefficients c _i , and the calculation can be done in a manner similar to Horner’s scheme.

 """ InterpPolynomial <: AbstractPolynomial Interpolation polynomials in Newton's form """ struct InterpPolynomial<:AbstractPolynomial degree::Int xval::NTuple{N, Float64} where N coeff::NTuple{N, Float64} where N end """ evpoly(p::Polynomial, z::Real) Evaluate polynomial `p` at `z` using the Horner's rule """ function evpoly(p::InterpPolynomial, z::Real) ans = p.coeff[p.degree+1] for idx = p.degree:-1:1 ans = ans * (z - p.xval[idx]) + p.coeff[idx] end return ans end 

The function for calculating the value of a polynomial in both cases is called the same - evpoly() - but it takes different types of arguments.

In addition to the calculation function, it would be nice to also write a function that creates a polynomial from known data.

For this, Julia has two methods: external constructors and internal constructors. An external constructor is simply a function that returns an object of the appropriate type. The inner constructor is a function that is entered inside the structure description and replaces the standard constructor. To construct interpolation polynomials, it is advisable to use the internal constructor, since

• it is more convenient to get a polynomial not through interpolation nodes and coefficients, but through nodes and values ​​of the interpolated function
• interpolation nodes must necessarily be different
• the number of nodes and coefficients must match

Writing an internal constructor, in which these rules are guaranteed to be observed, ensures that all created variables of the InterpPolynomial type InterpPolynomial at least be correctly processed by the evpoly() function.

We write the constructor of ordinary polynomials, which takes as input a one-dimensional array or a tuple of coefficients. The constructor of the interpolation polynomial takes as input the interpolation nodes and the desired values ​​in them and uses the method of divided differences to calculate the coefficients.

 """ Polynomial <: AbstractPolynomial Polynomials written in the canonical form --- Polynomial(v::T) where T<:Union{Vector{<:Real}, NTuple{<:Any, <:Real}}) Construct a `Polynomial` from the list of the coefficients. The coefficients are assumed to go from power 0 in the ascending order. If an empty collection is provided, the constructor returns a zero polynomial. """ struct Polynomial<:AbstractPolynomial degree::Int coeff::NTuple{N, Float64} where N function Polynomial(v::T where T<:Union{Vector{<:Real}, NTuple{<:Any, <:Real}}) #     /     P(x) ≡ 0 coeff = isempty(v) ? (0.0,) : tuple([Float64(x) for x in v]...) #   -   new #  -    return new(length(coeff)-1, coeff) end end """ InterpPolynomial <: AbstractPolynomial Interpolation polynomials in Newton's form --- InterpPolynomial(xsample::Vector{<:Real}, fsample::Vector{<:Real}) Construct an `InterpPolynomial` from a vector of points `xsample` and corresponding function values `fsample`. All values in `xsample` must be distinct. """ struct InterpPolynomial<:AbstractPolynomial degree::Int xval::NTuple{N, Float64} where N coeff::NTuple{N, Float64} where N function InterpPolynomial(xsample::X, fsample::F) where {X<:Union{Vector{<:Real}, NTuple{<:Any, <:Real}}, F<:Union{Vector{<:Real}, NTuple{<:Any, <:Real}}} #   ,    ,   f  ,   if !allunique(xsample) throw(DomainError("Cannot interpolate with duplicate X points")) end N = length(xsample) if length(fsample) != N throw(DomainError("Lengths of X and F are not the same")) end coeff = [Float64(f) for f in fsample] #     (Stoer, Bulirsch, Introduction to Numerical Analysis, . 2.1.3) for i = 2:N for j = 1:(i-1) coeff[i] = (coeff[j] - coeff[i]) / (xsample[j] - xsample[i]) end end new(N-1, tuple([Float64(x) for x in xsample]...), tuple(coeff...)) end end 

In addition to the actual generation of polynomials, it would be nice to be able to perform arithmetic operations with them.

Since in Julia, arithmetic operators are ordinary functions to which an infix notation is added as syntactic sugar (the expressions a + b and +(a, b) are both valid and absolutely identical), their overload is done in the same way as writing additional methods to its functions.

The only subtle point is that the custom code is run from the module (namespace) Main , and the functions of the standard library are in the Base module, so when overloaded, you must either import the Base module or write the full name of the function.

So, add the addition of a polynomial with a number:

 # -   Base.+  , #    Base.:+,   " :+   Base" function Base.:+(p::Polynomial, x::Real) Polynomial(tuple(p.coeff[1] + x, p.coeff[2:end]...)) end function Base.:+(p::InterpPolynomial, x::Real) # ..           - #          . #       - #        fval::Vector{Float64} = [evpoly(p, xval) + x for xval in p.xval] InterpPolynomial(p.xval, fval) end #       function Base.:+(x::Real, p::AbstractPolynomial) return p + x end 

To add two ordinary polynomials, it is enough to add the coefficients, and when adding an interpolation polynomial with another, you can find the sum values ​​at several points and build a new interpolation on them.

 function Base.:+(p1::Polynomial, p2::Polynomial) #    ,      deg = max(p1.degree, p2.degree) coeff = zeros(deg+1) coeff[1:p1.degree+1] .+= p1.coeff coeff[1:p2.degree+1] .+= p2.coeff Polynomial(coeff) end function Base.:+(p1::InterpPolynomial, p2::InterpPolynomial) xmax = max(p1.xval..., p2.xval...) xmin = min(p1.xval..., p2.xval...) deg = max(p1.degree, p2.degree) dx = xmax - xmin #         #       chebgrid = zeros(deg+1) for k = 1:deg-1 chebgrid[k+1] = xmax - 0.5 * dx * (1 + cos((k - 0.5) * π / (deg - 1))) end chebgrid[1] = xmin chebgrid[end] = xmax fsample = [evpoly(p1, x) + evpoly(p2, x) for x in chebgrid] InterpPolynomial(chebgrid, fsample) end function Base.:+(p1::InterpPolynomial, p2::Polynomial) xmax = max(p1.xval...) xmin = min(p1.xval...) deg = max(p1.degree, p2.degree) dx = xmax - xmin chebgrid = zeros(deg+1) for k = 1:deg-1 chebgrid[k+1] = xmax - 0.5 * dx * (1 + cos((k - 0.5) * π / (deg - 1))) end chebgrid[1] = xmin chebgrid[end] = xmax fsample = [evpoly(p1, x) + evpoly(p2, x) for x in chebgrid] InterpPolynomial(chebgrid, fsample) end function Base.:+(p1::Polynomial, p2::InterpPolynomial) p2 + p1 end 

In the same way, you can add other arithmetic operations on polynomials, as a result of getting a representation of them in the code in a natural mathematical notation.

That's all for now. I will try to write further about the implementation of other numerical methods.

In preparing the materials used:

1. Julia language documentation: docs.julialang.org
2. Julia language discussion platform: discourse.julialang.org
3. J.Stoer, W. Bulirsch. Introduction to Numerical Analysis
4. Julia Hub: habr.com/ru/hub/julia
5. Think Julia: benlauwens.imtqy.com/ThinkJulia.jl/latest/book.html