Julia：多项式的类型，多重方法和算术

在我看来，该出版物将重点介绍Julia语言的主要特色-以具有多个分派的方法的形式表示功能。 一方面，这可以提高计算性能，而又不会降低代码的可读性，并且不会破坏抽象性，另一方面，您可以使用更熟悉的表示法来处理数学概念。 作为一个例子，考虑了均匀性的问题（从线性运算的角度来看）与系数列表中的多项式和插值多项式一起工作。

基本语法

对于不认识的人的简要介绍。 Julia是一种类似于脚本的语言，具有REPL（读取-评估-打印循环，即交互式外壳）。 乍一看，它看起来与Python或MATLAB非常相似。

算术运算

算术与其他地方几乎相同：+，-，*，/，^用于求幂，等等。
比较：>，<，> =，<=，== 、! =等。
作业：=。
特点：除以/总是返回小数； 如果您需要将两个整数相除的整数部分，则需要使用div(m, n)运算或等价的in等于m ÷ n 。

种类

数值类型：

• 整数（ Int ） -42 ， -42
• 无符号整数（ UInt ） 0x12345
• 浮点数（ Float32 ， Float64 ） -Inf ， -Inf ， NaN
• 有理（ Rational ） 3//3 7//2
• 实数（ Real ）-以上所有
• 复杂（ Complex ） 3+4*im ， 2//3+2//3*im ， 3.0+0.0*im （ im是一个虚数单位，只有具有明确书写的虚部的数字才被视为复数）
• Number -以上所有

字符串和字符：

• 'a' Char字符（ Char ）
• "a"是一个字符串（ String ）

注意：就像许多语言一样，字符串是不可变的。
注意：字符串（以及变量名）支持Unicode，包括表情符号。

数组：

• x = [1, 2, 3] -通过直接枚举元素指定数组
• 特殊构造函数： zeros(length)表示零数组，一ones(length)表示一数组， rand(length)表示随机数数组，等等。
• 多维数组支持
• 支持标准库中的线性代数运算（数组加法，标量乘法，矩阵向量乘法等）

注意：所有集合均从一个索引开始。
注意：因为 该语言用于计算任务，数组是最重要的类型之一，您将不得不不止一次地返回其工作原理。

元组（元素的有序集合，不可变）：

• (2, 5.3, "k")是一个常规元组
• (a = 3, b = 4) -命名元组

注意：命名元组的字段可以通过句点按名称访问，也可以通过[]通过索引访问

 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灯：

 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 。 第二个工作就像在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定义，函数的定义也以单词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程序甚至几乎一对一地转换为Julia。
但是，在Julia中您可以显式指定变量的类型，这使您可以编译程序并获得快速的代码，这是一个很大的不同。
第二个显着差异是Python使用类和方法实现“经典” OOP模型，而Julia实现了多调度模型。

类型注释和多次分派

让我们看看一些内置函数是什么：

 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()不同， S_typed()函数不仅运行速度更快，而且不需要为每个调用分配内存。 这里的问题是，从mysqrt()返回的mysqrt()的类型未定义，就像表达式右侧的类型一样

 sum = sum + mysqrt(sgn) 

结果，编译器甚至无法弄清楚每次迭代的类型sum 。 因此，装箱（类型标签挂钩）是一个变量，并分配了内存。
对于S_typed()函数，编译器事先知道sum是一个复数值，因此代码得到了更优化（特别是，可以有效地内联调用mysqrt() ，始终将返回值返回到Complex ）。

更重要的是，对于S_typed()编译器知道返回值的类型为Complex ，但是对于S()则不再定义输出值S()类型，这将减慢将调用S()所有函数的速度。
您可以使用@code_warntype宏来验证编译器是否考虑了从表达式返回的类型：

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

如果在循环中某处调用了某个函数，而@code_warntype无法@code_warntype打印返回类型，或者该@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]) 

Julia中的结构是不可变的，即通过创建该结构的实例不再可以更改字段值（更确切地说，您无法更改内存中字段的地址-可以更改可变字段的元素，例如上例中的sv ）。 可变结构是由mutable struct构造创建的，其语法与常规结构的语法相同。

不支持“古典”意义上的结构继承，但是，有可能通过将复合类型组合成超类型或抽象类型（如Julia中所称）来“继承”行为。 类型关系表示为A<:B （A是A<:B的子类型）和A>:B （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 

案例研究：多项式

带有多个分派的类型系统便于表达数学概念。 让我们看一个用于处理多项式的简单库的示例。
我们介绍两种类型的多项式：通过幂系数定义的“规范”和由一对对（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的形式方便地组织存储，并且可以类似于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 

在这两种情况下，用于计算多项式值的函数称为evpoly() ，但接受不同类型的参数。

除了计算函数外，最好编写一个从已知数据创建多项式的函数。

在Julia中，有两种技术可以使用：外部构造函数和内部构造函数。 外部构造函数只是一个函数，它返回适当类型的对象。 内部构造函数是在结构描述中引入的函数，它取代了标准构造函数。 建议使用内部构造函数构造插值多项式，因为

• 不通过插值节点和系数而是通过插值函数的节点和值来获得多项式更加方便
• 插值节点必须不同
• 节点数和系数必须匹配

编写一个内部构造函数（确保可以遵循这些规则）可确保至少所有创建的InterpPolynomial类型的变量都可以由evpoly()函数正确处理。

我们编写了一个普通多项式的构造函数，该构造函数将一维数组或系数的元组作为输入。 插值多项式的构造函数在其中接收插值节点和所需值，并使用除差法计算系数。

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

要添加两个普通多项式，只需将系数相加就足够了，将插值多项式与另一个相加时，您可以在几个点上找到和值并从中建立一个新的插值。

 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. 朱莉娅语言文档： docs.julialang.org
2. 朱莉娅语言讨论平台： discourse.julialang.org
3. J. Stoer，W。Bulirsch。 数值分析导论
4. 朱莉娅· 哈伯（ Julia Hub）： habr.com/en/hub/julia
5. 想想朱莉娅： benlauwens.imtqy.com/ThinkJulia.jl/latest/book.html