倭算数理研究所

科学・数学・学習関連の記事を、「倭マン日記」とは別に書いていくのだ!

spire で多項式 Polynomial を使ってみる (2) : 数学関数としてのメソッド

spire を使ってみるシリーズ(目次)。 今回は Polynomial に定義されている、数学関数に関連するメソッドを見ていきます。 記事中のサンプルコードでは、以下の import 文が書かれているものとします:

import spire.math._
import spire.implicits._

関数の評価

まずは apply() メソッド。 引数に変数の値を与えて多項式の値を計算します:

  val f = Polynomial("3x^3 - 2x^2 + 6")

  // apply() メソッド
  assert( f(0) == 6 )  // = 3*0^3 - 2*0^2 + 6
  assert( f(1) == 7 )  // = 3*1^3 - 2*1^2 + 6
  assert( f(2) == 22 ) // = 3*2^3 - 2*2^2 + 6

まぁ、これ以外の動作をしたら驚きますねぇ。 少し注意が必要なのは、上記の Polynomial オブジェクト f は Polynomial[Rational] 型なので(Rational は spire の有理数型)、apply() メソッドの返り値は Rational 型になります(Rational 型は Int 値と等値評価できるので、上記の assert 文は通ります)。

式の評価を Double 値として行いたい場合などは、evalWith() メソッドによって係数に変換を施して評価を行います:

  val f: Polynomial[Rational] = Polynomial("3x^3 - 2x^2 + 6")

  // evalWith() メソッド
  assert( f.evalWith(3.0)(_.toDouble) == 69.0 )

2つ目の引数リストに変換関数を渡し、1つ目の引数リストの値を使って式の評価を行います。 今回は扱いませんが、Polynomial 型には係数に変換を施す map() メソッドが定義されていて、evalWith() メソッドはこの map() メソッドを施した後に上記の apply() メソッドを呼んでいるだけです。

当然ですが、evalWith() メソッドは係数の型変換以外の変換でも使えます。

多項式の性質

次は多項式に関連するメソッド。 次数や係数を取得するメソッドを見ていきます。 係数を取得するメソッドはここで見るもの以外にたくさんありますが、それは別記事で見ていくことにします。

degree多項式の次数を返します:

  // degree
  assert( Polynomial("3x^3 - 2x^2 + 6").degree == 3 )

nth() メソッドは n 次の項の係数を返します:

  val p = Polynomial("3x^3 - 2x^2 + 6")

  // nth()  メソッド
  assert( p.nth(2) == -2 )
  assert( p.nth(0) == 6 )

maxOrderTermCoeff は最高次の係数を返します:

  val p = Polynomial("3x^3 - 2x^2 + 6")

  // maxOrderTermCoeff
  assert( p.maxOrderTermCoeff == 3 )

isZero, isConstant によって多項式が 0 や定数であることをテストできます:

  val p0 = Polynomial.constant(0)
  val p1 = Polynomial("3x^3 - 2x^2 + 6")

  // isZero
  assert( p0.isZero )
  assert( !p1.isZero )

  // isConstant
  assert( p0.isConstant )
  assert( !p1.isConstant )

signVariations は、高次の項から順に(もしくは低次の項から順に)係数の符号を見ていって、符号が変わった回数を返します:

  // signVariations
  assert( Polynomial("3x^3 - 2x^2 + 6").signVariations == 2 )
  assert( Polynomial("4x^2 + 3").signVariations == 0)

係数が0になった場合はカウントしません(その前後での符号の変化を見ます)。

他の多項式を返すメソッド

ここでは、既にある多項式から別の多項式を取得するメソッドを見ていきます。 ただし、和差積のような代数に関連するものは別記事で。

monic は最高次の係数が1になるように多項式を定数倍します:

  // monic 
  assert( Polynomial("3x^3 - 2x^2 + 6").monic == Polynomial("x^3 - 2/3x^2 + 2"))

reciprocal は係数を逆順にします。 言い換えれば n 次式では i 次の項の係数を n-i 次の項の係数にします:

  // reciprocal
  assert( Polynomial("3x^3 - 2x^2 + 6").reciprocal == Polynomial("6x^3 - 2x + 3"))

flip は最高次の係数のみ符号を逆転します:

  // flip
  assert( Polynomial("3x^3 - 2x^2 + 6").flip == Polynomial("-3x^3 - 2x^2 + 6"))

reductum は最高次の項のみを削除します:

  // reductum
  assert( Polynomial("3x^3 - 2x^2 + 6").reductum == Polynomial("-2x^2 + 6") )

shift() メソッドは多項式の変数の値を指定した分だけずらします:

  val p = Polynomial("3x^3 - 2x^2 + 6")

  // shift() メソッド
  assert( p.shift(r"2") == Polynomial("3x^3 + 16x^2 + 28x + 22") )
    // = 3(x+2)^3 - 2(x+2)^2 + 6
    //   上記の r"2" は単に有理数 (Rational) の2のこと。 shift() メソッドとは直接関係なし。

  // compose() メソッドを使うと以下のように書ける
  assert( p.shift(r"2") == p.compose(Polynomial("x + 2")) )

shift() メソッドは、次に見る compose() メソッドを使って同じことができますが、shift() メソッドを使う方がパフォーマンスが良いそうです。 ちなみにグラフの平行移動とは異なるので注意。

compose() メソッドは関数の合成を行います:

  val f = Polynomial("2x^2 + 3x + 4")
  val g = Polynomial("x - 2")

  // compose() メソッド
  assert( (f compose g) == (2*:(g**2) + 3*:g + 4) )
    // = 2(x-2)^2 + 3(x-2) + 4
    //   *: は右結合の定数倍

  assert( (f compose g).degree == f.degree * g.degree )
    // 合成した多項式の次数は、それぞれの次数の積

  // もちろん、順序を変えれば結果は変わり得る
  assert( (g compose f) == (f - 2) )

Scala の関数オブジェクトに定義されている andThen() メソッド(compose と順序が逆な合成)は、Polynomial には定義されていないようです。

多項式微分積分は、それぞれ derivative, integral で計算できます:

  val p = Polynomial("3x^3 - 2x^2 + 6")

  // derivative
  assert( p.derivative == Polynomial("9x^2 - 4x") )

  // integral
  assert( p.integral == Polynomial("3/4x^4 - 2/3x^3 + 6x") )

積分定数は0になります。

解に関するメソッド

最後は多項式を = 0 とおいた代数方程式の解(根)に関連するメソッドを見ていきます。

解は、基本的には roots を呼び出すだけで計算してくれますが、係数の型が Rational や BigDecimal だと精度を丸めないといけないので、(必要なら)係数を実数型 Real や Number 型に変換しなければいけません*1。 係数の型の変換は map() メソッドを使います:

  val p: Polynomial[Real] = Polynomial("x^3 - 6x^2 + 11x - 6").map(_.toReal)
    // = (x-1)(x-2)(x-3)

  // roots
  assert( p.roots.toSet == Set(1, 2, 3) )
    // 返される Set は実際には Set[Real] 型。 Real 型と Int 値は等値評価できる。

返り値は Scala のコレクションではなく spire.math.poly.Roots 型ですが、toList や toSet などで Scala コレクションに変換できます。 Roots 型には解の個数を返す count がありますが、これは重解はダブって数えます:

  val p = Polynomial("x^3 - 4x^2 + 5x - 2").map(_.toReal) // = (x-1)^2(x-2)

  assert( p.roots.toSet == Set(1, 2) )  // 異なる解の個数は2個
  assert( p.roots.count == 3 )  // 解の個数は3個

Roots 型には他にもいろいろメソッドが定義されています。

方程式を Double 値で解きたい場合は Polynomial[Number] 型(この Number は spire.math.Number 型)にする必要があります*2

  val p: Polynomial[Number] = Polynomial("x^3 - 6x^2 + 11x - 6").map(_.toDouble)

  assert( p.roots.toSet == Set(1.0, 2.0, 3.0) )
    // Double 値の等値評価はよくない

実数型で解を求めた場合とパフォーマンスがどの程度違うのか気になるところですが、特に比べてません。

removeZeroRoots は、多項式 = 0 の代数方程式から x = 0 の解となるものを取り除いた多項式を返します。 要は最低次の項が定数項になりように x の冪で割るだけですが:

  assert(Polynomial("3x^3 + 2x^2").removeZeroRoots == Polynomial("3x + 2"))

まぁ、なんてことないですね。

今回は Polynomial に定義されている数学関数としてのメソッドを見てきました。 いくつか、直接自分では使わなさそうなメソッドもありましたが、どこかのアルゴリズムの中で使われているのでしょう。

次回は代数体としてのメソッドを見ていく予定。

Scalaスケーラブルプログラミング第3版

Scalaスケーラブルプログラミング第3版

*1:BigDecimal で精度を MathContext によって指定する方法もあるようですが、ここでは扱いません。

*2:これは spire の型クラスの宣言なり暗黙のオブジェクトの定義がおかしいせいだと思うのだけど、よくわかりません。 ドキュメントには Double 値では簡単に呼び出せると書いてるんですが・・・