Julia: jenis, multimetode, dan aritmatika di atas polinomial

Dalam publikasi ini, kami akan fokus pada fitur utama, dalam pendapat saya, fitur khas dari bahasa Julia - representasi fungsi dalam bentuk metode dengan pengiriman ganda. Ini memungkinkan Anda untuk meningkatkan kinerja perhitungan tanpa mengurangi keterbacaan kode dan tanpa merusak abstraknya, di satu sisi, dan memungkinkan Anda untuk bekerja dengan konsep matematika dalam notasi yang lebih akrab, di sisi lain. Sebagai contoh, pertanyaan seragam (dari sudut pandang operasi linier) bekerja dengan polinomial dalam representasi daftar koefisien dan dengan polinomial interpolasi dipertimbangkan.

Sintaks dasar

Pengantar singkat untuk mereka yang tidak tahu. Julia adalah bahasa mirip naskah, memiliki REPL (loop baca-evaluasi-cetak, mis. Shell interaktif). Sepintas tampilannya mirip, misalnya, dengan Python atau MATLAB.

Operasi aritmatika

Aritmatika hampir sama dengan di mana-mana: +, -, *, /, ^ untuk eksponensial, dll.
Perbandingan:>, <,> =, <=, == ,! = Dll
Tugas: =.
Fitur: pembagian melalui / selalu mengembalikan angka pecahan; jika Anda memerlukan bagian integer dari pembagian dua bilangan bulat, Anda perlu menggunakan div(m, n) operasi div(m, n) atau infix setara m ÷ n .

Jenis

Jenis numerik:

• Integer ( Int ) - 2 , 3 , -42
• Bilangan bulat tak UInt ( UInt ) - 0x12345
• Floating point ( Float32 , Float64 ) - 1.0 , 3.1415 , -Inf , NaN
• Rasional ( Rational ) - 3//3 , 7//2
• Nyata ( Real ) - semua hal di atas
• Kompleks ( Complex ) - 3+4*im , 2//3+2//3*im , 3.0+0.0*im ( im adalah unit imajiner, hanya angka dengan bagian imajiner yang ditulis secara eksplisit dianggap kompleks)
• Number - semua hal di atas

String dan karakter:

• 'a' - karakter ( Char )
• "a" adalah string ( String )

NB: string, seperti sekarang dalam banyak bahasa, tidak dapat diubah.
NB: string (dan juga nama variabel) mendukung Unicode, termasuk emoji.

Array:

• x = [1, 2, 3] - menentukan array dengan penghitungan langsung elemen
• konstruktor khusus: zeros(length) untuk array nol, ones(length) untuk array yang, rand(length) untuk array angka acak, dll.
• dukungan array multidimensi
• dukungan untuk operasi aljabar linier (penambahan array, perkalian skalar, perkalian vektor matriks, dan banyak lagi) di perpustakaan standar

NB: semua koleksi diindeks mulai dari satu.
NB: karena bahasa ini dimaksudkan untuk tugas komputasi, array adalah salah satu tipe yang paling penting, Anda harus kembali ke prinsip kerja mereka lebih dari sekali.

Tuples (kumpulan elemen yang dipesan, tidak dapat berubah):

• (2, 5.3, "k") adalah tupel biasa
• (a = 3, b = 4) - bernama tuple

NB: bidang tuple bernama dapat diakses baik dengan nama melalui periode dan indeks melalui []

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

Kamus:

 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 

Konstruksi bahasa kontrol dasar

1. Variabel secara otomatis dibuat atas penugasan. Jenisnya opsional.

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

2. Blok melompat bersyarat dimulai dengan ekspresi if <condition> dan berakhir dengan kata end . Anda juga dapat memiliki lampu else atau lampu lain:

 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. Ada dua konstruksi loop: while dan for . Yang kedua berfungsi seperti di Python, yaitu Iterate atas koleksi. Penggunaan umum adalah iterasi pada rentang nilai yang sintaksnya start[:increment]:end . Tidak seperti Python, rentang mencakup nilai awal dan akhir, mis. rentang kosong tidak akan menjadi 1:1 (ini adalah kisaran 1), tetapi 1:0 . Akhir dari badan loop ditandai dengan kata 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. Fungsi didefinisikan oleh function kata kunci, definisi fungsi juga berakhir dengan kata end . Argumen dengan nilai default dan argumen bernama didukung.

 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 

Secara umum, ini semua sangat mirip dengan Python, kecuali untuk perbedaan kecil dalam sintaks dan fakta bahwa blok kode tidak dialokasikan dengan spasi, tetapi masih dengan kata kunci. Dalam kasus-kasus sederhana, program-program Python bahkan menerjemahkan ke dalam Julia hampir satu lawan satu.
Tetapi ada perbedaan yang signifikan dalam fakta bahwa di Julia Anda dapat secara eksplisit menentukan jenis untuk variabel, yang memungkinkan Anda untuk mengkompilasi program, mendapatkan kode cepat.
Perbedaan signifikan kedua adalah bahwa Python mengimplementasikan model OOP "klasik" dengan kelas dan metode, sementara Julia mengimplementasikan model multi-pengiriman.

Jenis Anotasi dan Pengiriman Ganda

Mari kita lihat apa fungsi bawaan:

 julia> sqrt sqrt (generic function with 19 methods) 

Seperti yang ditunjukkan REPL kepada kami, sqrt adalah fungsi umum dengan 19 metode. Seperti apa fungsi umum dan metode apa?

Dan ini berarti ada beberapa fungsi sqrt yang berlaku untuk berbagai jenis argumen dan, karenanya, menghitung akar kuadrat menggunakan berbagai algoritma. Anda dapat melihat opsi apa yang tersedia dengan mengetik

 julia> methods(sqrt) 

Dapat dilihat bahwa fungsi didefinisikan untuk berbagai jenis angka, serta untuk matriks.

Berbeda dengan OOP "klasik", di mana implementasi konkret dari metode ini hanya ditentukan oleh kelas panggilan (dispatching oleh argumen pertama), di Julia pilihan fungsi ditentukan oleh jenis (dan nomor) dari semua argumennya.
Saat memanggil suatu fungsi dengan argumen spesifik dari semua metodenya, satu dipilih yang paling akurat menggambarkan set tipe spesifik yang dengannya fungsi dipanggil, dan inilah yang digunakan.

Fitur yang khas adalah bahwa suatu pendekatan yang disebut kompilasi "tepat waktu" oleh penulis bahasa diterapkan. Yaitu fungsi dikompilasi untuk tipe data yang diberikan pada panggilan pertama, setelah itu panggilan berikut dilakukan jauh lebih cepat. Perbedaan antara panggilan pertama dan selanjutnya bisa sangat signifikan:

 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 

Dalam kasus yang buruk, setiap pemanggilan fungsi adalah pemeriksaan jenis argumen yang diterima dan pencarian metode yang diinginkan dalam daftar. Namun, jika Anda memberikan petunjuk kepada kompiler, Anda dapat menghilangkan pemeriksaan, yang akan menyebabkan kode lebih cepat.

Misalnya, pertimbangkan untuk menghitung jumlahnya

$$

$$\ 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() menunjukkan bahwa fungsi S_typed() tidak hanya berjalan lebih cepat, tetapi juga tidak memerlukan alokasi memori untuk setiap panggilan, tidak seperti S() . Masalahnya di sini adalah bahwa jenis mysqrt() dikembalikan dari mysqrt() tidak didefinisikan, sama seperti jenis sisi kanan ekspresi

 sum = sum + mysqrt(sgn) 

Akibatnya, kompiler bahkan tidak bisa mengetahui sum tipe apa yang akan ada pada setiap iterasi. Jadi, tinju (tipe label hooking) adalah variabel dan memori dialokasikan.
Untuk fungsi S_typed() , kompiler tahu sebelumnya bahwa sum adalah nilai yang kompleks, sehingga kode lebih dioptimalkan (khususnya, memanggil mysqrt() dapat secara efektif sebaris, selalu mengembalikan nilai kembali ke Complex ).

Lebih penting lagi, untuk S_typed() kompiler tahu bahwa nilai kembali adalah tipe Complex , tetapi untuk S() tipe nilai output tidak didefinisikan lagi, yang akan memperlambat semua fungsi di mana S() akan dipanggil.
Anda dapat memverifikasi bahwa kompilator berpikir tentang tipe yang dikembalikan dari ekspresi menggunakan makro @code_warntype :

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

Jika suatu fungsi dipanggil di suatu tempat di loop dimana @code_warntype tidak dapat mencetak tipe kembali, atau yang di suatu tempat di dalam tubuh menunjukkan tanda terima dari nilai tipe Any , maka optimasi panggilan-panggilan ini kemungkinan besar akan memberikan peningkatan kinerja yang sangat nyata.

Jenis Senyawa

Seorang programmer dapat mendefinisikan tipe data komposit untuk kebutuhannya menggunakan 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]) 

Struktur di Julia tidak dapat diubah, yaitu, dengan membuat contoh struktur, tidak lagi mungkin untuk mengubah nilai bidang (lebih tepatnya, Anda tidak dapat mengubah alamat bidang dalam memori - elemen bidang yang bisa berubah, seperti sv dalam contoh di atas, dapat diubah). Struktur yang dapat berubah dibuat oleh konstruk mutable struct , sintaksisnya sama dengan struktur biasa.

Warisan struktur dalam arti "klasik" tidak didukung, namun ada kemungkinan perilaku "mewarisi" dengan menggabungkan tipe komposit menjadi supertipe atau, sebagaimana mereka disebut dalam Julia, tipe abstrak. Tipe hubungan dinyatakan sebagai A<:B (A adalah subtipe dari B) dan A>:B (A adalah subtipe dari B). Itu terlihat seperti ini:

 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 

Studi kasus: Polinomial

Suatu sistem tipe yang digabungkan dengan beberapa pengiriman sangat nyaman untuk mengekspresikan konsep matematika. Mari kita lihat contoh perpustakaan sederhana untuk bekerja dengan polinomial.
Kami memperkenalkan dua jenis polinomial: "kanonik", yang didefinisikan melalui koefisien pada kekuatan, dan "interpolasi", yang didefinisikan oleh seperangkat pasangan (x, f (x)). Untuk mempermudah, kami hanya akan mempertimbangkan argumen yang valid.

Untuk menyimpan polinomial dalam notasi biasa, struktur yang memiliki larik atau tupel koefisien sebagai bidang cocok. Agar benar-benar kekal, biarkan ada iring-iringan. Dengan demikian, kode untuk mendefinisikan tipe abstrak, struktur polinomial, dan menghitung nilai polinomial pada titik tertentu cukup sederhana:

 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 

Polinomial interpolasi membutuhkan struktur representasi dan metode perhitungan yang berbeda. Secara khusus, jika himpunan titik interpolasi diketahui sebelumnya, dan direncanakan untuk menghitung polinomial yang sama pada titik yang berbeda, rumus interpolasi Newton mudah digunakan:

$$

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

di mana n _k ( x ) adalah polinomial dasar, n ₀ ( x ) dan untuk k > 0

$$

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

di mana x _i adalah node interpolasi.

Dari rumus di atas dapat dilihat bahwa penyimpanan mudah diatur dalam bentuk seperangkat node interpolasi x _i dan koefisien c _i , dan perhitungan dapat dilakukan dengan cara yang mirip dengan skema Horner.

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

Fungsi untuk menghitung nilai polinomial dalam kedua kasus disebut sama - evpoly() - tetapi menerima berbagai jenis argumen.

Selain fungsi perhitungan, alangkah baiknya menulis fungsi yang membuat polinomial dari data yang diketahui.

Ada dua teknik untuk ini di Julia: konstruktor eksternal dan konstruktor internal. Konstruktor eksternal hanyalah sebuah fungsi yang mengembalikan objek dengan tipe yang sesuai. Konstruktor internal adalah fungsi yang diperkenalkan di dalam deskripsi struktur dan menggantikan konstruktor standar. Disarankan untuk menggunakan konstruktor internal untuk membangun polinomial interpolasi, karena

• lebih mudah mendapatkan polinomial bukan melalui interpolasi node dan koefisien, tetapi melalui node dan nilai fungsi interpolasi
• node interpolasi harus berbeda
• jumlah node dan koefisien harus cocok

Menulis konstruktor internal di mana aturan-aturan ini dijamin untuk diikuti memastikan bahwa semua variabel yang dibuat dari tipe InterpPolynomial , setidaknya, dapat diproses dengan benar oleh fungsi evpoly() .

Kami menulis konstruktor polinomial biasa yang menggunakan larik satu dimensi atau tupel koefisien sebagai input. Konstruktor polinomial interpolasi menerima node interpolasi dan nilai yang diinginkan di dalamnya dan menggunakan metode perbedaan terbagi untuk menghitung koefisien.

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

Selain generasi polinomial yang sebenarnya, alangkah baiknya untuk dapat melakukan operasi aritmatika dengan mereka.

Karena operator aritmatika di Julia adalah fungsi biasa, yang notasi infiksnya ditambahkan sebagai gula sintaksis (ekspresi a + b dan +(a, b) keduanya valid dan benar-benar identik), overloading mereka dilakukan dengan cara yang sama seperti menulis metode tambahan untuk fungsinya.

Satu-satunya titik halus adalah bahwa kode pengguna diluncurkan dari modul Main (namespace), dan fungsi-fungsi perpustakaan standar ada di modul Base , jadi ketika kelebihan, Anda harus mengimpor modul Base atau menulis nama lengkap fungsi.

Jadi, kami menambahkan penambahan polinomial dengan nomor:

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

Untuk menambahkan dua polinomial biasa, cukup menambahkan koefisien, dan ketika menambahkan polinomial interpolasi ke yang lain, Anda dapat menemukan nilai penjumlahan di beberapa titik dan membangun interpolasi baru dari mereka.

 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 

Dengan cara yang sama, Anda dapat menambahkan operasi aritmatika lainnya pada polinomial, menghasilkan representasi mereka dalam kode dalam notasi matematika alami.

Itu saja untuk saat ini. Saya akan mencoba untuk menulis lebih lanjut tentang penerapan metode numerik lainnya.

Dalam persiapan, bahan-bahan berikut digunakan:

1. Dokumentasi Bahasa Julia: docs.julialang.org
2. Platform Diskusi Bahasa Julia: discourse.julialang.org
3. J. Stoer, W. Bulirsch. Pengantar Analisis Numerik
4. Julia Hub: habr.com/en/hub/julia
5. Think Julia: benlauwens.imtqy.com/ThinkJulia.jl/latest/book.html