Julia: Typen, Multimethoden und Arithmetik über Polynome

In dieser Veröffentlichung konzentrieren wir uns meiner Meinung nach auf das Hauptmerkmal der Julia-Sprache - die Darstellung von Funktionen in Form von Methoden mit Mehrfachversand. Auf diese Weise können Sie die Leistung von Berechnungen steigern, ohne die Lesbarkeit des Codes zu beeinträchtigen und einerseits die Abstraktheit zu beeinträchtigen, und andererseits können Sie mit mathematischen Konzepten in einer bekannteren Notation arbeiten. Als Beispiel wird die Frage der Einheitlichkeit (unter dem Gesichtspunkt linearer Operationen) bei der Darstellung der Koeffizientenliste mit Polynomen und bei Interpolationspolynomen betrachtet.

Grundlegende Syntax

Eine kurze Einführung für Unbekannte. Julia ist eine skriptähnliche Sprache mit REPL (Read-Evaluate-Print-Schleife, d. H. Eine interaktive Shell). Auf den ersten Blick sieht es zum Beispiel Python oder MATLAB ziemlich ähnlich.

Arithmetische Operationen

Die Arithmetik ist ungefähr gleich wie überall: +, -, *, /, ^ für die Potenzierung usw.
Vergleich:>, <,> =, <=, ==,! = Etc.
Zuordnung: =.
Features: Division durch / immer eine Bruchzahl zurück; Wenn Sie den ganzzahligen Teil der Division zweier Ganzzahlen benötigen, müssen Sie die Operation div(m, n) oder das Infix-Äquivalent m ÷ n .

Typen

Numerische Typen:

• Ganzzahlen ( Int ) - 2 , 3 , -42
• Ganzzahlen ohne UInt ( UInt ) - 0x12345
• Gleitkomma ( Float32 , Float64 ) - 1.0 , 3.1415 , -Inf , NaN
• Rational ( Rational ) - 3//3 , 7//2
• Real ( Real ) - alles oben Genannte
• Komplex ( Complex ) - 3+4*im , 2//3+2//3*im , 3.0+0.0*im ( im ist eine imaginäre Einheit, nur eine Zahl mit einem explizit geschriebenen Imaginärteil wird als komplex angesehen)
• Number - alle oben genannten

Zeichenfolgen und Zeichen:

• 'a' - Zeichen ( Char )
• "a" ist ein String ( String )

NB: Zeichenfolgen sind, wie jetzt in vielen Sprachen, unveränderlich.
NB: Zeichenfolgen (sowie Variablennamen) unterstützen Unicode, einschließlich Emoji.

Arrays:

• x = [1, 2, 3] - Angabe eines Arrays durch direkte Aufzählung von Elementen
• spezielle Konstruktoren: zeros(length) für ein Array von Nullen, ones(length) für ein Array von Einsen, rand(length) für ein Array von Zufallszahlen usw.
• Unterstützung für mehrdimensionale Arrays
• Unterstützung für lineare Algebraoperationen (Hinzufügen von Arrays, Skalarmultiplikation, Matrixvektormultiplikation und vieles mehr) in der Standardbibliothek

NB: Alle Sammlungen werden ab einer indiziert.
NB: weil Die Sprache ist für Rechenaufgaben gedacht. Arrays sind einer der wichtigsten Typen. Sie müssen mehr als einmal zu den Prinzipien ihrer Arbeit zurückkehren.

Tupel (geordneter Satz von Elementen, unveränderlich):

• (2, 5.3, "k") ist ein reguläres Tupel
• (a = 3, b = 4) - benanntes Tupel

NB: Auf die Felder eines benannten Tupels kann sowohl nach Namen durch einen Punkt als auch nach Index über [] zugegriffen werden.

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

Wörterbücher:

 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 

Grundlegende Kontrollsprachenkonstrukte

1. Variablen werden bei der Zuweisung automatisch erstellt. Typ ist optional.

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

2. Der bedingte Sprungblock beginnt mit dem Ausdruck if <condition> und endet mit dem Worte end . Sie können auch ein else Licht oder elseif Lichter haben:

 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. Es gibt zwei Schleifenkonstrukte: while und for . Das zweite funktioniert wie in Python, d.h. Iteriert über die Sammlung. Eine häufige Verwendung ist das Durchlaufen eines Wertebereichs, dessen Syntax start[:increment]:end lautet. Im Gegensatz zu Python umfasst ein Bereich sowohl Start- als auch Endwerte, d. H. Der leere Bereich ist nicht 1:1 (dies ist ein Bereich von 1), sondern 1:0 . Das Ende des Schleifenkörpers ist mit dem Worte end markiert.

 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. Funktionen werden durch die Schlüsselwortfunktion definiert, die Definition der Funktion endet ebenfalls mit dem Worte end . Argumente mit Standardwerten und benannten Argumenten werden unterstützt.

 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 

Im Allgemeinen ist dies alles Python ziemlich ähnlich, mit Ausnahme geringfügiger Unterschiede in der Syntax und der Tatsache, dass die Codeblöcke nicht mit Leerzeichen, sondern mit Schlüsselwörtern belegt sind. In einfachen Fällen übersetzen Python-Programme sogar fast eins zu eins in Julia.
Es gibt jedoch einen signifikanten Unterschied in der Tatsache, dass Sie in Julia explizit Typen für Variablen angeben können, wodurch Sie Programme kompilieren und schnellen Code erhalten können.
Der zweite wesentliche Unterschied besteht darin, dass Python ein „klassisches“ OOP-Modell mit Klassen und Methoden implementiert, während Julia ein Multi-Dispatch-Modell implementiert.

Geben Sie Anmerkungen und Mehrfachversand ein

Mal sehen, was eine eingebaute Funktion ist:

 julia> sqrt sqrt (generic function with 19 methods) 

Wie REPL zeigt, ist sqrt eine generische Funktion mit 19 Methoden. Welche Art von verallgemeinerter Funktion und welche Art von Methoden?

Dies bedeutet, dass es mehrere sqrt Funktionen gibt, die für verschiedene Arten von Argumenten gelten und dementsprechend die Quadratwurzel mit verschiedenen Algorithmen berechnen. Sie können sehen, welche Optionen verfügbar sind, indem Sie eingeben

 julia> methods(sqrt) 

Es ist ersichtlich, dass die Funktion sowohl für verschiedene Arten von Zahlen als auch für Matrizen definiert ist.

Im Gegensatz zum „klassischen“ OOP, bei dem die konkrete Implementierung der Methode nur von der aufrufenden Klasse bestimmt wird (Dispatching durch das erste Argument), wird bei Julia die Auswahl einer Funktion durch die Typen (und die Anzahl) aller ihrer Argumente bestimmt.
Wenn Sie eine Funktion mit bestimmten Argumenten aus all ihren Methoden aufrufen, wird eine ausgewählt, die den spezifischen Satz von Typen, mit denen die Funktion aufgerufen wird, am genauesten beschreibt und der verwendet wird.

Eine Besonderheit besteht darin, dass ein Ansatz angewendet wird, der von den Autoren der Sprache als "nur vorzeitig" zusammengestellt wird. Das heißt, Funktionen werden beim ersten Aufruf für die angegebenen Datentypen kompiliert, wonach die folgenden Aufrufe viel schneller ausgeführt werden. Der Unterschied zwischen dem ersten und den nachfolgenden Anrufen kann sehr bedeutend sein:

 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 

Im schlechten Fall ist jeder Funktionsaufruf eine Überprüfung des Typs der empfangenen Argumente und eine Suche nach der gewünschten Methode in der Liste. Wenn Sie dem Compiler jedoch Hinweise geben, können Sie die Überprüfungen eliminieren, was zu schnellerem Code führt.

Betrachten Sie beispielsweise die Berechnung der Summe

$$

$$\ 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 

Der Benchmark zeigt, dass die Funktion S_typed() Gegensatz zu S() nicht nur schneller ausgeführt wird, sondern auch nicht für jeden Aufruf eine Speicherzuweisung erfordert. Das Problem hierbei ist, dass der Typ des von mysqrt() zurückgegebenen mysqrt() nicht definiert ist, genau wie der Typ der rechten Seite des Ausdrucks

 sum = sum + mysqrt(sgn) 

Infolgedessen kann der Compiler nicht einmal herausfinden, welche sum bei jeder Iteration sein wird. Boxen (Typ Label Hooking) ist also eine Variable und Speicher wird zugewiesen.
Für die Funktion S_typed() weiß der Compiler im Voraus, dass sum ein komplexer Wert ist, sodass der Code optimierter ist (insbesondere kann der Aufruf von mysqrt() effektiv inline erfolgen und den Rückgabewert immer an Complex ).

Noch wichtiger ist, dass S_typed() Compiler für S_typed() weiß, dass der Rückgabewert vom Typ Complex , aber für S() Typ des Ausgabewerts nicht erneut definiert, wodurch alle Funktionen verlangsamt werden, bei denen S() aufgerufen wird.
Mit dem Makro @code_warntype können Sie überprüfen, ob der Compiler über die vom Ausdruck zurückgegebenen Typen @code_warntype :

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

Wenn eine Funktion irgendwo in der Schleife aufgerufen wird, für die @code_warntype den Rückgabetyp nicht ausgeben kann oder für die irgendwo im Body der Empfang eines Werts vom Typ Any @code_warntype , führt die Optimierung dieser Aufrufe höchstwahrscheinlich zu einer spürbaren Leistungssteigerung.

Verbindungstypen

Ein Programmierer kann zusammengesetzte Datentypen für seine Anforderungen mithilfe des 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]) 

Strukturen in Julia sind unveränderlich, d. H. Durch Erstellen einer Instanz der Struktur ist es nicht mehr möglich, die Feldwerte zu ändern (genauer gesagt, Sie können die Adresse von Feldern im Speicher nicht ändern - Elemente von veränderlichen Feldern, wie z. B. sv im obigen Beispiel, können geändert werden). Mutable Strukturen werden durch das mutable Strukturkonstrukt erstellt, dessen Syntax dieselbe ist wie für reguläre Strukturen.

Die Vererbung von Strukturen im „klassischen“ Sinne wird nicht unterstützt, es besteht jedoch die Möglichkeit, Verhalten zu „erben“, indem zusammengesetzte Typen zu Supertypen oder, wie sie in Julia genannt werden, abstrakten Typen kombiniert werden. Typbeziehungen werden ausgedrückt als A<:B (A ist ein Subtyp von B) und A>:B (A ist ein Subtyp von B). Es sieht ungefähr so ​​aus:

 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 

Fallstudie: Polynome

Ein Typensystem, das mit Mehrfachversand gekoppelt ist, ist praktisch, um mathematische Konzepte auszudrücken. Schauen wir uns ein Beispiel einer einfachen Bibliothek für die Arbeit mit Polynomen an.
Wir führen zwei Arten von Polynomen ein: "kanonisch", definiert durch Koeffizienten bei Potenzen, und "Interpolation", definiert durch eine Menge von Paaren (x, f (x)). Der Einfachheit halber werden nur gültige Argumente berücksichtigt.

Zum Speichern eines Polynoms in einer üblichen Notation ist eine Struktur mit einem Array oder einem Tupel von Koeffizienten als Feld geeignet. Um völlig unveränderlich zu sein, soll es eine Wagenkolonne geben. Daher ist der Code zum Definieren des abstrakten Typs, der Struktur des Polynoms und zum Berechnen des Werts des Polynoms an einem bestimmten Punkt recht einfach:

 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 

Interpolationspolynome benötigen eine andere Darstellungsstruktur und Berechnungsmethode. Insbesondere wenn der Satz von Interpolationspunkten im Voraus bekannt ist und geplant ist, dasselbe Polynom an verschiedenen Punkten zu berechnen, ist die Newtonsche Interpolationsformel zweckmäßig:

$$

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

wobei n _k ( x ) Grundpolynome sind, n ₀ ( x ) und für k > 0

$$

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

wobei x _i die Interpolationsknoten sind.

Aus den obigen Formeln ist ersichtlich, dass die Speicherung zweckmäßigerweise in Form eines Satzes von Interpolationsknoten x _i und Koeffizienten c _{i organisiert} ist und die Berechnung auf ähnliche Weise wie das Horner-Schema durchgeführt werden kann.

 """ 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 

Die Funktion zur Berechnung des Wertes des Polynoms heißt in beiden Fällen gleich - evpoly() - akzeptiert jedoch unterschiedliche Arten von Argumenten.

Zusätzlich zur Berechnungsfunktion wäre es schön, eine Funktion zu schreiben, die aus bekannten Daten ein Polynom erstellt.

In Julia gibt es dafür zwei Techniken: externe Konstruktoren und interne Konstruktoren. Ein externer Konstruktor ist einfach eine Funktion, die ein Objekt des entsprechenden Typs zurückgibt. Ein interner Konstruktor ist eine Funktion, die in die Strukturbeschreibung eingeführt wird und den Standardkonstruktor ersetzt. Es ist ratsam, den internen Konstruktor zu verwenden, um Interpolationspolynome zu konstruieren, da

• Es ist bequemer, ein Polynom nicht durch die Interpolationsknoten und -koeffizienten zu erhalten, sondern durch die Knoten und Werte der interpolierten Funktion
• Interpolationsknoten müssen unterschiedlich sein
• Die Anzahl der Knoten und Koeffizienten muss übereinstimmen

Durch das Schreiben eines internen Konstruktors, in dem diese Regeln garantiert eingehalten werden, wird sichergestellt, dass alle erstellten Variablen vom Typ InterpPolynomial zumindest von der Funktion evpoly() korrekt verarbeitet werden können.

Wir schreiben einen Konstruktor gewöhnlicher Polynome, der ein eindimensionales Array oder ein Tupel von Koeffizienten als Eingabe verwendet. Der Konstruktor des Interpolationspolynoms empfängt die Interpolationsknoten und die darin enthaltenen gewünschten Werte und berechnet die Koeffizienten nach der Methode der geteilten Differenzen .

 """ 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 

Neben der eigentlichen Erzeugung von Polynomen wäre es schön, mit ihnen arithmetische Operationen durchführen zu können.

Da arithmetische Operatoren in Julia gewöhnliche Funktionen sind, zu denen eine Infixnotation als syntaktischer Zucker hinzugefügt wird (die Ausdrücke a + b und +(a, b) sind beide gültig und absolut identisch), erfolgt ihre Überladung auf die gleiche Weise wie beim Schreiben zusätzliche Methoden zu ihren Funktionen.

Der einzige subtile Punkt ist, dass der Benutzercode vom Hauptmodul (Namespace) gestartet wird und sich die Funktionen der Standardbibliothek im Basismodul befinden. Wenn Sie also überladen, müssen Sie entweder das Base importieren oder den vollständigen Namen der Funktion schreiben.

Also fügen wir die Addition eines Polynoms mit einer Zahl hinzu:

 # -   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 

Um zwei gewöhnliche Polynome hinzuzufügen, reicht es aus, die Koeffizienten zu addieren. Wenn Sie das Interpolationspolynom zum anderen hinzufügen, können Sie die Summenwerte an mehreren Punkten finden und daraus eine neue Interpolation erstellen.

 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) #         #       xmid = 0.5 * xmax + 0.5 * xmin dx = 0.5 * (xmax - xmin) / cos(0.5 * π / (deg + 1)) chebgrid = [xmid + dx * cos((k - 0.5) * π / (deg + 1)) for k = 1:deg+1] 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) xmid = 0.5 * xmax + 0.5 * xmin dx = 0.5 * (xmax - xmin) / cos(0.5 * π / (deg + 1)) chebgrid = [xmid + dx * cos((k - 0.5) * π / (deg + 1)) for k = 1:deg+1] fsample = [evpoly(p1, x) + evpoly(p2, x) for x in chebgrid] InterpPolynomial(chebgrid, fsample) end function Base.:+(p1::Polynomial, p2::InterpPolynomial) p2 + p1 end 

Auf die gleiche Weise können Sie andere arithmetische Operationen für Polynome hinzufügen, was zu deren Darstellung im Code in einer natürlichen mathematischen Notation führt.

Das ist alles für jetzt. Ich werde versuchen, weiter über die Implementierung anderer numerischer Methoden zu schreiben.

Zur Herstellung wurden folgende Materialien verwendet:

1. Julia Sprachdokumentation: docs.julialang.org
2. Julia Sprachdiskussionsplattform: diskurs.julialang.org
3. J. Stoer, W. Bulirsch. Einführung in die numerische Analyse
4. Julia Hub: habr.com/en/hub/julia
5. Denken Sie an Julia: benlauwens.imtqy.com/ThinkJulia.jl/latest/book.html