ジュリア：型、マルチメソッド、多項式上の算術

この出版物は、私の意見では、ジュリア言語の主な特徴である、複数のディスパッチを伴うメソッドの形式での関数の表示に焦点を当てます。 これにより、コードの可読性を低下させることなく、抽象性を損なうことなく、計算のパフォーマンスを向上させることができます。また、数学的な概念をより馴染みのある表記で扱うことができます。 例として、（線形操作の観点から）均一の問題は、係数のリストの表現での多項式と補間多項式で動作します。

基本的な構文

知らない人のための簡単な紹介。 ジュリアは、いわば、スクリプト言語であり、REPL（読み取り-評価-印刷ループ、つまり対話型シェル）を備えています。 一見すると、たとえばPythonやMATLABと非常によく似ています。

算術演算

算術はどこでもほぼ同じです：累乗の場合は+、-、*、/、^など

比較：>、<、> =、<=、==、！=など

割り当て：=。

機能：除算/



常に小数を返します。 2つの整数の除算の整数部分が必要な場合、演算div(m, n)



または中置に相当するm ÷ n



を使用する必要があります。

種類

数値タイプ：

• 全体（ Int
  
  
  
  ） -42
  
  
  
  、 -42
  
  
  
  
• 符号なし整数（ UInt
  
  
  
  ） 0x12345
  
  
  
  
• 浮動小数点（ Float32
  
  
  
  、 Float64
  
  
  
  ） Float32
  
  
  
  、 -Inf
  
  
  
  、 NaN
  
  
  
  
• Rational（ Rational
  
  
  
  ） 3//3



7//2
  
  
  
  
• Real（ Real
  
  
  
  ）-上記のすべて
• Complex（ Complex
  
  
  
  ） 3+4*im
  
  
  
  、 2//3+2//3*im
  
  
  
  、 3.0+0.0*im
  
  
  
  （ im
  
  
  
  は虚数単位であり、虚数部が明示的に記述された数値のみが複素数と見なされます）
• Number
  
  
  
  -上記のすべて

文字列と文字：

• 'a'



Char
  
  
  
  キャラクター（ Char
  
  
  
  ）
• "a"
  
  
  
  は文字列です（ String
  
  
  
  ）

NB：現在、多くの言語で見られるように、文字列は不変です。

注意：文字列（および変数名）は、絵文字を含むUnicodeをサポートしています。



配列：

• x = [1, 2, 3]
  
  
  
  -要素の直接列挙による配列の指定
• 特別なコンストラクター： zeros(length)
  
  
  
  の配列の場合はゼロones(length)
  
  
  
  、1の配列の場合は1 rand(length)
  
  
  
  、 rand(length)
  
  
  
  配列の場合はrand(length)
  
  
  
  など。
• 多次元配列のサポート
• 標準ライブラリでの線形代数演算（配列の追加、スカラー乗算、行列ベクトル乗算など）のサポート

注意：すべてのコレクションは、1から始まるインデックスが付けられます。

NB：なぜなら この言語は計算タスクを対象としているため、配列は最も重要なタイプの1つであるため、作業の原則に何度も戻る必要があります。



タプル（要素の順序セット、不変）：

• (2, 5.3, "k")
  
  
  
  は通常のタプルです
• (a = 3, b = 4)
  
  
  
  -名前付きタプル

NB：名前付きタプルのフィールドには、ピリオドを介した名前と[]を介したインデックスの両方でアクセスできます

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

辞書：

 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
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

基本的な制御言語構造

1.変数は割り当て時に自動的に作成されます。 タイプはオプションです。

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

2.条件付きジャンプブロックは、式if <condition>



で始まり、単語end



ます。 else



ライトまたはelseif



ライトを使用することもできelse



。

 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. while



とfor



2つのループ構造があります。 2番目は、Pythonのように機能します。 コレクションを反復処理します。 一般的な使用法は、構文がstart[:increment]:end



である値の範囲を反復処理することです。 Pythonとは異なり、範囲には開始値と終了値の両方が含まれます。 空の範囲は1:1



（これは1の範囲）ではなく、 1:0



ます。 ループ本体の終わりには、単語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.関数はキーワードfunction



によって定義され、 function



の定義もend



という単語でend



ます。 デフォルト値を持つ引数と名前付き引数がサポートされています。

 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
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

一般に、これはPythonと非常によく似ていますが、構文のわずかな違いと、コードブロックにはスペースが割り当てられず、キーワードが割り当てられるという事実があります。 単純な場合、Pythonプログラムはほぼ1対1でJuliaに変換されます。

しかし、ジュリアでは変数の型を明示的に指定できるため、プログラムをコンパイルして高速なコードを取得できるという点で大きな違いがあります。

2番目の重要な違いは、Pythonがクラスとメソッドを備えた「クラシックな」OOPモデルを実装するのに対し、ジュリアはマルチディスパッチモデルを実装することです。

型注釈と複数ディスパッチ

いくつかの組み込み関数を見てみましょう：

 julia> sqrt sqrt (generic function with 19 methods)
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

REPLが示すように、 sqrt



は19のメソッドを持つ汎用関数です。 どのような一般化された関数とどのようなメソッドですか？



そしてこれは、さまざまなタイプの引数に適用され、それに応じてさまざまなアルゴリズムを使用して平方根を計算するいくつかの sqrt



関数があることを意味します。 次のように入力して、使用可能なオプションを確認できます。

 julia> methods(sqrt)
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

この関数は、さまざまなタイプの数値と行列に対して定義されていることがわかります。



メソッドの具体的な実装が呼び出しクラスによってのみ決定される（最初の引数によるディスパッチ）「クラシック」OOPとは異なり、Juliaでは、関数の選択はそのすべての引数のタイプ（および数）によって決定されます。

すべてのメソッドから特定の引数を使用して関数を呼び出す場合、関数が呼び出される特定のタイプのセットを最も正確に記述するものが選択され、それが使用されます。



特徴的な機能は、言語の作者による「ジャスト・アヘッド・タイム」コンパイルと呼ばれるアプローチが適用されることです。 つまり 関数は、最初の呼び出しで指定されたデータ型用にコンパイルされ、その後の次の呼び出しははるかに高速です。 最初の呼び出しと後続の呼び出しの違いは非常に重要です。

 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
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

悪いケースでは、各関数呼び出しは受け取った引数のタイプのチェックであり、リスト内の目的のメソッドを検索します。 ただし、コンパイラにヒントを与えると、チェックを省略でき、コードの高速化につながります。



たとえば、合計を計算することを検討してください

$$

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

ベンチマークは、 S_typed()



関数がより高速に実行されるだけでなく、 S()



とは異なり、各呼び出しにメモリ割り当てを必要としないことを示しています。 ここでの問題は、式の右側の型のように、 mysqrt()



から返されるmysqrt()



の型が定義されていないことです。

 sum = sum + mysqrt(sgn)
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

その結果、コンパイラーは各反復でどのタイプのsum



なるかを把握することさえできません。 したがって、ボクシング（タイプラベルフッキング）は変数であり、メモリが割り当てられます。

S_typed()



関数の場合、コンパイラーはsum



が複素数値であることを事前に認識しているため、コードはより最適化されています（特にmysqrt()



呼び出しは事実上インラインで、常に戻り値をComplex



返します）。



さらに重要なことは、 S_typed()



コンパイラーは戻り値の型がComplex



であることを知っていますが、 S()



場合、出力値の型は再定義されないため、 S()



が呼び出されるすべての関数の速度が低下します。

@code_warntype



マクロを使用して、 @code_warntype



が式から返されたタイプを@code_warntype



することを確認できます。

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

@code_warntype



が戻り値の型を出力できないループのどこかで関数が呼び出された場合、または本体のどこかでAny



型の値の受信を示す関数が呼び出された場合、これらの呼び出しの最適化により、パフォーマンスが大幅に向上する可能性が高くなります。

複合タイプ

プログラマーは、 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])
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

ジュリアの構造体は不変です。つまり、構造体のインスタンスを作成することにより、フィールド値を変更できなくなります（より正確には、メモリ内のフィールドのアドレスを変更することはできません-上記の例のsv



などの可変フィールドの要素は変更できます）。 mutable struct



は、 mutable struct



によって作成されます。その構文は、通常の構造と同じです。



「古典的な」意味での構造の継承はサポートされていませんが、複合型をスーパータイプに結合したり、ジュリアでは「抽象型」と呼ばれるようにしたりすることで、動作を「継承」する可能性があります。 型の関係は、 A<:B



（AはBのサブタイプ）およびA>:B



（AはA>:B



サブタイプ）として表されます。 次のようになります。

 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
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

事例研究：多項式

複数のディスパッチと組み合わせた型システムは、数学的な概念を表現するのに便利です。 多項式を扱うためのシンプルなライブラリの例を考えてみましょう。

2種類の多項式を導入します。累乗係数で定義される「正準」と、ペアのセット（x、f（x））で定義される「補間」です。 簡単にするために、有効な引数のみを考慮します。



通常の表記法で多項式を保存するには、フィールドとして係数の配列またはタプルを持つ構造が適しています。 完全に不変であるために、モーターカードがあります。 したがって、抽象型、多項式の構造を指定し、特定のポイントで多項式の値を計算するためのコードは非常に簡単です。

 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
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

補間多項式には、異なる表現構造と計算方法が必要です。 特に、補間点のセットが事前にわかっていて、異なる点で同じ多項式を計算することが計画されている場合、 ニュートンの補間式が便利です。

$$

$$P（x）= \ sum_ {k = 0} ^ N {c_k n_k（x）}、$$

ここで、 n _k （ x ）は基本多項式、 n ₀ （ x ）およびk > 0の場合

$$

$$n_k（x）= \ prod_ {i = 0} ^ {k-1} {（x-x_i）}、$$

ここで、 x _iは補間ノードです。



上記の式から、ストレージは補間ノードx _iおよび係数c _iのセットとして便利に編成でき、計算はホーナーのスキームと同様の方法で実行できることがわかります。

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

両方の場合に多項式の値を計算する関数は同じものと呼ばれますevpoly()



-異なる種類の引数を受け入れます。



計算関数に加えて、既知のデータから多項式を作成する関数を作成すると便利です。



Juliaには、外部コンストラクターと内部コンストラクターの2つのメソッドがあります。 外部コンストラクターは、適切なタイプのオブジェクトを返す単なる関数です。 内部コンストラクターは、構造記述内に導入され、標準コンストラクターを置き換える関数です。 内部コンストラクターを使用して補間多項式を構築することをお勧めします。

• 補間ノードと係数ではなく、補間関数のノードと値から多項式を取得する方が便利です
• 補間ノードは明確でなければなりません
• ノードと係数の数は一致する必要があります

これらのルールが遵守されることが保証されている内部コンストラクターを作成すると、作成されたInterpPolynomial



型のすべての変数が、少なくともevpoly()



関数によって正しく処理されることが保証されます。



入力として1次元配列または係数のタプルを受け取る通常の多項式のコンストラクターを記述します。 補間多項式のコンストラクターは、補間ノードとその中の目的の値を受け取り、 分割された差分の方法を使用して係数を計算します。

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

多項式の実際の生成に加えて、それらを使用して算術演算を実行できると便利です。



Juliaの算術演算子は通常の関数であり、それに中置記法が構文糖として追加されるため（式a + b



および+(a, b)



は両方とも有効であり、完全に同一です）、それらのオーバーロードは書き込みと同じ方法で行われます機能への追加メソッド。



唯一の微妙な点は、ユーザーコードがMain



モジュール（名前空間）から起動され、標準ライブラリの関数がBase



モジュールにあるため、オーバーロードする場合は、 Base



モジュールをインポートするか、関数のフルネームを記述する必要があります。



そのため、数のある多項式を追加します。

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

2つの通常の多項式を追加するには、係数を追加するだけで十分です。補間多項式を他の多項式に追加すると、いくつかのポイントで合計値を見つけて、それらから新しい補間を作成できます。

 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
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

同様に、多項式に他の算術演算を追加して、自然な数学表記でコードに表現することができます。



今のところすべてです。 他の数値的手法の実装についてさらに記述していきます。



準備では、次の材料が使用されました。

1. Julia言語のドキュメント： docs.julialang.org
2. Julia Language Discussion Platform： discourse.julialang.org
3. J.ストーア、W。ブリルシュ。 数値解析の概要
4. ジュリアハブ： habr.com/en/hub/julia
5. ジュリアを考える： benlauwens.github.io/ThinkJulia.jl/latest/book.html