immut/dense チュートリアル

このチュートリアルでは、一変数多項式をイミュータブルな値として扱う方法を学びます。構築、算術、評価、合成、微分を行い、長い入力に適した乗算を選ぶ方法を説明します。

クイックスタート

モジュールを追加し、DensePolynomial を再エクスポートしているイミュータブルなファサードをインポートします。

moon add Luna-Flow/luna-poly@0.2.0
import {
  "Luna-Flow/luna-poly/immut",
}
test "dense quick start" {
  let p = @immut.DensePolynomial::from_coefficients([1, 2, 3])
  inspect(p, content="1 + 2x^1 + 3x^2")
  inspect(p.eval(2), content="17")
}

係数は定数項から順に並べるので、[1, 2, 3] は 1+2x+3x21 + 2x + 3x^2 であり、p(2)=1+4+12=17p(2) = 1 + 4 + 12 = 17 です。Luna-Flow/luna-poly/immut/dense を直接インポートすることもでき、型は同じです。

日常的なタスク

多項式を構築する

多項式全体には from_coefficients を、演算子で組み合わせる部品には constant、variable、monomial を使います。

test "building polynomials" {
  let x : @immut.DensePolynomial[Int] = @immut.DensePolynomial::variable()
  let two = @immut.DensePolynomial::constant(2)
  let p = x * x - two * x + @immut.DensePolynomial::monomial(0, 1)
  inspect(p, content="1 + -2x^1 + 1x^2")
  debug_inspect(p.to_coefficients(), content="[1, -2, 1]")
  inspect(@immut.DensePolynomial::from_coefficients([0, 0, 0]).is_zero(), content="true")
}

末尾のゼロはどこでも取り除かれるので、[0, 0, 0] は零多項式であり、to_coefficients は常に最短の形を返します。

係数と次数を読む

test "reading" {
  let p = @immut.DensePolynomial::from_coefficients([5, 0, 7])
  debug_inspect(p.degree(), content="Some(2)")
  inspect(p.coefficient(0), content="5")
  inspect(p.coefficient(10), content="0")
  debug_inspect(p.leading_coefficient(), content="Some(7)")
  let zero : @immut.DensePolynomial[Int] = @immut.DensePolynomial::zero()
  debug_inspect(zero.degree(), content="None")
}

零多項式には次数がないので、degree は 0 ではなく None を返します。

古い値を失わずに算術を行う

すべての演算は新しい多項式を返し、オペランドは元のままです。

test "value semantics" {
  let p = @immut.DensePolynomial::from_coefficients([1, 1])
  let q = p * p
  let r = q.scale(1, 2)
  inspect(p, content="1 + 1x^1")
  inspect(q, content="1 + 2x^1 + 1x^2")
  inspect(r, content="2x^1 + 4x^2 + 2x^3")
  inspect(p.pow(4), content="1 + 4x^1 + 6x^2 + 4x^3 + 1x^4")
}

scale(k, c) は c xkc\,x^k を掛け、pow(e) は繰り返し二乗法による反復乗算です。

評価、合成、微分

test "calculus" {
  let p = @immut.DensePolynomial::from_coefficients([1.0, -3.0, 0.0, 1.0])
  inspect(p.eval(2.0), content="3")
  let shift = @immut.DensePolynomial::from_coefficients([1.0, 1.0])
  let moved = p.substitute(shift)
  inspect(moved.eval(1.0), content="3")
  let dp = p.derivative()
  debug_inspect(dp.to_coefficients(), content="[-3, 0, 3]")
  inspect(dp.eval(1.0), content="0")
}

p.substitute(q) は合成 p(q(x))p(q(x)) なので、moved は p(x+1)p(x + 1) であり、moved.eval(1.0) は p.eval(2.0) に等しくなります。derivative は形式微分で、ここでは p′=3x2−3p' = 3x^2 - 3 が 11 で消えます。

長い多項式を高速に掛ける

* は筆算式の積で O(mn)O(mn) です。長いオペランドには karatsuba を呼び出してください。同じ結果を約 O(n1.585)O(n^{1.585}) で返します。

test "karatsuba" {
  let a = @immut.DensePolynomial::from_coefficients(Array::makei(200, i => i % 5 - 2))
  let b = @immut.DensePolynomial::from_coefficients(Array::makei(150, i => i % 3 - 1))
  let fast = a.karatsuba(b)
  assert_true(fast == a * b)
  inspect(fast.length(), content="349")
}

短いほうのオペランドの係数が 33 個未満なら karatsuba は単に * を呼ぶので、実用上 * より遅くなることはありません。

さらに進む

係数型に関するジェネリックなアルゴリズム

ファサードが再エクスポートしている luna-generic の制約でアルゴリズムを書けば、必要な能力を持つすべての係数型で動作します。

fn[A : @immut.AddMonoid + Mul + Eq] square_eval(p : @immut.DensePolynomial[A], a : A) -> A {
  (p * p).eval(a)
}

test "generic coefficients" {
  inspect(square_eval(@immut.DensePolynomial::from_coefficients([1, 1]), 2), content="9")
  inspect(square_eval(@immut.DensePolynomial::from_coefficients([1U, 1U]), 2U), content="9")
  inspect(square_eval(@immut.DensePolynomial::from_coefficients([0.5, 0.5]), 1.0), content="1")
}

UInt には符号反転がありませんが、* と eval はそれを要求しないので、UInt の係数でも動作します。

表現に関するジェネリックなアルゴリズム

DensePolynomial::ops() は演算をレコードにまとめるので、同じコードをミュータブルな表現でも実行できます。

fn[P, A] horner_check(ops : @immut.UnivariateOps[P, A], p : P, a : A) -> A {
  ops.eval(ops.pow(p, 3), a)
}

test "ops records" {
  let p = @immut.DensePolynomial::from_coefficients([1, 1])
  let m = @mutable.DensePolynomial::from_coefficients([1, 1])
  inspect(horner_check(@immut.DensePolynomial::ops(), p, 1), content="8")
  inspect(horner_check(@mutable.DensePolynomial::ops(), m, 1), content="8")
}

中断しないエラー処理

monomial、coefficient、scale は負の累乗に対して中断します。それらの _checked 版は代わりに None を返します。

fn safe_shift(p : @immut.DensePolynomial[Int], k : Int) -> @immut.DensePolynomial[Int] {
  p.scale_checked(k, 1).unwrap_or(p)
}

test "checked" {
  let p = @immut.DensePolynomial::from_coefficients([3])
  inspect(safe_shift(p, 2), content="3x^2")
  inspect(safe_shift(p, -2), content="3")
}

arithmetic による累乗

DensePolynomial は @arithmetic.PowNatChecked を実装しているので、Luna-Flow/arithmetic に対して書かれたコードで多項式を自然数乗できます。

test "pow nat checked" {
  let p = @immut.DensePolynomial::from_coefficients([1, 1])
  let ctx = @arithmetic.ArithmeticContext::new(0)
  match @arithmetic.PowNatChecked::pow_nat_checked(p, 2, ctx) {
    Ok(q) => inspect(q, content="1 + 2x^1 + 1x^2")
    Err(_) => fail("pow_nat_checked never fails for polynomials")
  }
}

よくある落とし穴

  • 係数の順序。 配列は昇順です。[a, b, c] は a+bx+cx2a + bx + cx^2 であり、ax2+bx+cax^2 + bx + c ではありません。

  • ゼロの degree。 None です。零多項式に対して 0 が欲しい場合は length() を使ってください。

  • 固定幅のオーバーフロー。 Int の係数はラップアラウンドし、先頭係数がゼロに回り込んで次数が下がることがあります。

    test "wrapping" {
      let p = @immut.DensePolynomial::from_coefficients([1, 65536])
      inspect((p * p).degree().unwrap(), content="1")
    }
  • 浮動小数点の切り詰め。 切り詰められるのは、ゼロと等しいと比較される係数だけです。0.1 + 0.2 - 0.3 はゼロではないので、そのような係数は残ります。

  • 微分には NatHomomorphism が必要。 derivative は Float、Double、BigInt の係数で動作します。Int の係数は微分できません。

  • 0^0。 pow(0) はゼロを含むすべての多項式に対して 1 を返します。

  • 表示。 to_string は x^1 を明示的に書き、ゼロの項を省きます。厳密な出力には to_coefficients を使ってください。

次のステップ