Julia に翻訳--222 多項式の値を求める

Julia に翻訳--222　多項式の値を求める

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

algorithm as241  appl. statist. (1988) vol. 37, no. 3, 477- 484.
https://jblevins.org/mirror/amiller/as181.f90

ファイル名: poly.jl　　関数名: poly

翻訳するときに書いたメモ

元は FORTRAN プログラム。書き方も古風なので（GO TO があったり），少し新しめに。

==========#

function poly(coefficients, x)

  #        algorithm as 181.2   appl. statist.  (1982) vol. 31, no. 2

  #        calculates the algebraic polynomial of order nored-1 with
  #        array of coefficients coefficients.  zero order coefficient is coefficients[1]

  nord = length(coefficients)
  result = coefficients[1]
  nord == 1 && return result
  p = x * coefficients[nord]
  if nord != 2
    n2 = nord - 2
    j = n2 + 1
    for i = 1:n2
      p = (p + coefficients[j]) * x
      j -= 1
    end
  end
  result + p
end

元のプログラムは古い FORTRAN で書かれているので，もっとわかりにくいが，書き換えてもなお，
プログラムが長く，一見何をやっているかわかりづらいが，nord が 1 であるかとか，2 であるとか，定数項を特別扱いしたりとかしなくてもよい。
たとえば，以下の5次式
0.1 + 0.221157∙x - 0.147981∙x^2 - 2.07119∙x^3 + 4.434685∙x^4 - 2.706056∙x^5
を ^ を使わないで
((((-2.706056∙x + 4.434685)∙x - 2.07119)∙x - 0.147981)∙x + 0.221157)∙x + 0.1
として計算する。

以下のように，短く簡潔でわかりやすいプログラムになる。

function poly2(coefficients, x)

  #        algorithm as 181.2   appl. statist.  (1982) vol. 31, no. 2

  #        calculates the algebraic polynomial of order nored-1 with
  #        array of coefficients coefficients.  zero order coefficient is coefficients[1]

  nord = length(coefficients)
  result = 0
  for i = nord:-1:2
    result = (result + coefficients[i]) * x
  end
  result + coefficients[1]
end

c6 = [0.1, 0.221157, -0.147981, -2.07119, 4.434685, -2.706056];
# 0.1 + 0.221157∙x - 0.147981∙x^2 - 2.07119∙x^3 + 4.434685∙x^4 - 2.706056∙x^5
# ((((-2.706056∙x + 4.434685)∙x - 2.07119)∙x - 0.147981)∙x + 0.221157)∙x + 0.1 これを計算する

x = 1.2
poly(c6, x)  # -0.9644910099199991
poly2(c6, x) # -0.9644910099199991

Julia では以下のようにすればよい。
using Polynomials
func = Polynomial(c6)
typeof(func) # Polynomial{Float64, :x}; func は関数
func(x)      # -0.9644910099199991

2次式，1次式，0次式の場合もちゃんと動くのを確認しておく。

c3 = [1.2, 3.4, 5.6]; # (5.6∙x + 3.4)∙x + 1.2
poly(c3, 3)       # 61.79999999999999
poly2(c3, 3)      # 61.79999999999999
Polynomial(c3)(3) # 61.79999999999999

c2 = [1.2, 3.4]; # 3.4∙x + 1.2
poly(c2, 3)       # 11.399999999999999
poly2(c2, 3)      # 11.399999999999999
Polynomial(c2)(3) # 11.399999999999999

c1 = [1.2]; # 1.2; x に何を指定しても無関係（定数が返される）
poly(c1, 5)       # 1.2
poly2(c1, 5)      # 1.2
Polynomial(c1)(5) # 1.2