第4章 · 多項式

式をデータにして、根の個数を調べる¶

√2 は有理数として作れませんが、 x²−2=0 を満たす数として指定できます。この章では式を係数列として保存し、代入とスツルム列（Sturm列）の実装から、根が1と2の間に一つだけあることを確かめます。

前章から使うもの

有理数の表現と演算、約分、作成時の不変条件、ログの表示範囲。

この章で学ぶこと

変数、係数列、ホーナー法（Horner法）、根、符号変化、スツルム列、高速化のための選択。

この章を終えると

x²−2 と係数列 (−2,0,1) を対応付けられる
Polynomial.__init__ と evaluate の処理を追える
根の存在と一意性を、別の根拠で説明できる

数学

有限個の有理数係数を持つ式へ値を代入し、0になる根を調べます。

データ表現

Polynomial.coefficients が、定数項から順に係数を保存します。

実装

コンストラクタが表現を正規化し、ホーナー法とスツルム列が根を調べます。

この章で横に置く実装

完全なクラスと根数え関数: peano/polynomial.py。本文では __init__、evaluate、 sturm_sequence、count_real_roots を抜粋します。

正規化・代入・根数の検査: tests/test_polynomial.py。ファイル索引は実装リファレンスにあります。

数学上の定義

`x²−2` の根として√2を指定する¶

x は、あとから数を代入する場所です。x² は x×x の短い表記です。数と x のべき乗を足して作る式を多項式と呼びます。

x²−2 = −2 + 0x + 1x²

−2、0、1 は各項の係数です。x²−2 へ x=1 を代入すると −1、x=2 なら 2 です。代入した結果が0になる x を根と呼びます。したがって、探している正の √2 は x²−2 の根です。

多項式の値は、入力を少し動かしたとき突然飛びません。この連続性と中間値の定理により、1で負、2で正なら、開区間 (1,2) に根が少なくとも一つあります。ただし、端点の符号だけでは根が一つとは限りません。開区間 (1,2) は、1と2そのものを含まず、その間にある実数全体です。ここでは中間値の定理を証明せず、連続な関数が負から正へ変わる途中では 0を通る、という保証として使います。

データの表し方

式を、定数項から始まる係数列として保存する¶

ライブラリは多項式の見た目の文字列を保存しません。係数を x⁰, x¹, x², … の順に並べた、変更できない列として保存します。

数学の式	保存する係数列
`x²−2`	`(−2,0,1)`
`1+x`	`(1,1)`
`3`	`(3)`

polynomial((-2,1),(0,1),(1,1)) の各組は、係数を (分子,分母) で書いたものです。内部では三つの Rational を作り、 Polynomial へ渡します。

このクラスは、これまでの数のクラスと違って init=False です。

@total_ordering
@dataclass(frozen=True, slots=True, init=False, eq=False, repr=False)
class Polynomial:
    _coefficients: tuple[Rational, ...]

    def __init__(self, *coefficients: Rational) -> None:
        if not coefficients:
            coefficients = (Q_ZERO,)
        if any(not isinstance(value, Rational) for value in coefficients):
            raise TypeError("Polynomial coefficients must be Rational")

        normalized = [value.reduction() for value in coefficients]
        while len(normalized) > 1 and normalized[-1] == Q_ZERO:
            normalized.pop()
        object.__setattr__(self, "_coefficients", tuple(normalized))

dataclassの生成を止め、正規化するコンストラクタを書く

係数を受け取ったまま保存すると、末尾の0や未約分の係数が残ります。そこで init=False で自動生成を止め、正規化してから一つのタプルへ保存する __init__ を手書きします。 dataclassはfrozenとslotsの制約には引き続き使われます。

Q_ZERO は前章で構成した有理数0を一度作って名前を付けた定数です。

すべての有理数係数を約分する
最高次側に続く0を取り除く
変更できない tuple として保存する

1+0x と 1 は同じ多項式なので、どちらも係数列 (1) へ正規化されます。整数と有理数では入力した代表元を残しましたが、多項式では表現を一意に寄せます。これは型ごとの設計上の違いです。

frozen=True では通常の self._coefficients = ... も拒否されるため、初期化中だけ基底の object.__setattr__ を明示的に呼びます。これは不変オブジェクトを手書きで初期化するためのPythonデータモデル上の処理で、数学的な操作ではありません。

Pythonでの実装

ホーナー法で係数列から値を計算する¶

係数が (a₀,a₁,a₂) の多項式は a₀+a₁x+a₂x² です。同じ式を次のように括ると、x の掛け算を一段ずつ進められます。

a₀ + x(a₁ + xa₂)

これをホーナー法と呼びます。Polynomial.evaluate は係数列を逆向きに読み、有理数の乗法と加法を繰り返します。

@log(log_level=31)
def evaluate(self, value: object) -> tuple[Rational, LogMessage]:
    point = cast2r(value)
    result = Q_ZERO
    for coefficient in reversed(self.coefficients):
        result = (result * point + coefficient).reduction()
    return result, lambda: f"{self!r}: x={point!r} -> {result!r}"

cast2r は、自然数・整数・有理数の入力を有理数へ揃える
reversed によって、最高次の係数から読む
毎回 reduction() して、途中の分数表現が膨らみすぎないようにする

x²−2 へ1を代入すると、内部状態は次のように変わります。

読んだ係数	更新式	`result`
1	`0×1+1`	1
0	`1×1+0`	1
−2	`1×1−2`	−1

この表は実装のループを手で実行したものです。

スツルム列の実装は、根を数えるために別の計算を組み立てます。次は入力検査と戻り値を含む完全な関数です。

def sturm_sequence(value: Polynomial) -> tuple[Polynomial, ...]:
    if not isinstance(value, Polynomial):
        raise TypeError("sturm_sequence expects a Polynomial")
    if value.degree <= 0:
        raise ValueError("a constant polynomial has no Sturm sequence")

    square_free = value.square_free()
    sequence = [square_free, square_free.derivative()]
    while sequence[-1]:
        remainder = sequence[-2] % sequence[-1]
        if not remainder:
            break
        sequence.append(-remainder)
    return tuple(sequence)

square_free() は、同じ場所に重なっている根を一つとして数えられる形へ直します。 x²−2 には重なった根がないため、ここでは元の式のままです。次の derivative() は導関数を作ります。係数と指数を掛けて次数を一つ下げるので、各項 aₖxᵏ は k×aₖxᵏ⁻¹ へ変わり、x²−2 の導関数は 2x です。

その後は前二つを多項式として割り、余りの符号を反転して追加します。 x²−2 なら、実装のループを次のように手で追えます。

段階	割る計算	余り	列へ追加するもの
初期化	—	—	`x²−2`, `2x`
1	`(x²−2) ÷ 2x`	`−2`	符号を反転した `2`
2	`2x ÷ 2`	`0`	追加せず終了

最初の余りは、x²−2 = (x/2)(2x)−2 と書けば −2 だと確認できます。したがってスツルム列は次の三つです。

x²−2,　2x,　2

実行して確かめる

代入ログと符号変化を、別の証拠として読む¶

evaluate には @log(log_level=31) が付いています。ログはループの各周ではなく、すべての係数を処理して戻り値ができたあとに一行だけ出ます。

実行前に、ホーナー法を一度手で追う

上に抜粋した evaluate のループを参照し、 x=2 の三つの result を書いてから実行してください。

⌘ / Ctrl + Enter

出力

二つの代入ログと、printの二つの結果を区別します。

<P(-2 + x^2)> は保存された多項式、<Q(1/1)> は代入した有理数、最後の <Q(-1/1)> はホーナー法で得た有理数です。代入ログの一行だけではループ途中は見えないため、実装を根拠に上の表を作ることに意味があります。

次にスツルム列の符号を端点で調べます。0を除いた符号が何回切り替わるかを sign_variations が数えます。

点	`x²−2, 2x, 2` の符号	変化回数
1	`−,+,+`	1
2	`+,+,+`	0

スツルムの定理により、左端の変化回数から右端の回数を引いた値が、開区間の相異なる実根の個数です。

実装も、その差をそのまま返します。次は二つの関数の完全な実装です。

def sign_variations(
    sequence: tuple[Polynomial, ...],
    point: Rational,
) -> int:
    signs: list[bool] = []
    for value in sequence:
        sign = value.sign_at(point)
        if sign == 0:
            continue
        signs.append(sign < 0)
    return sum(left != right for left, right in zip(signs, signs[1:]))


def count_real_roots(
    value: Polynomial,
    lower: Rational,
    upper: Rational,
) -> int:
    if not isinstance(value, Polynomial):
        raise TypeError("value must be a Polynomial")
    lower = cast2r(lower)
    upper = cast2r(upper)
    if lower >= upper:
        raise ValueError("lower must be less than upper")
    if value.sign_at(lower) == 0 or value.sign_at(upper) == 0:
        raise ValueError("interval endpoints must not be roots")
    sequence = sturm_sequence(value)
    return sign_variations(sequence, lower) - sign_variations(sequence, upper)

表示幅に合わせて関数のシグネチャと最後の式を改行していますが、処理は実ファイルと同じです。

存在と個数を分けて予想

上の count_real_roots のreturnを参照してください。中間値の定理が保証することと、符号変化の差 1−0 が保証することを、それぞれ一文で書いてください。

⌘ / Ctrl + Enter

出力

列、両端の符号変化、差として得る根の個数を順に読みます。

テストで確かめる

テストは表現・計算・根の個数を分離する¶

係数列の正規化は、末尾の0がある表現とない表現が等しいことで検査します。

self.assertEqual(
    Polynomial(rational(1, 1), rational(0, 1)),
    Polynomial(rational(1, 1)),
)

代入は、x²−2 へ 3/2 を入れた結果が 1/4 になることを確かめます。

self.assertEqual(
    self.sqrt_two_polynomial.evaluate(rational(3, 2)),
    rational(1, 4),
)

スツルム列については、(−2,2) に二つの実根 −√2 と √2 があること、列の長さが3になることを検査します。

self.assertEqual(
    count_real_roots(self.sqrt_two_polynomial, rational(-2, 1), rational(2, 1)),
    2,
)
self.assertEqual(len(sturm_sequence(self.sqrt_two_polynomial)), 3)

一つの大きなテストではなく、表現の不変条件、代入、根の個数という責務ごとに失敗箇所を切り分けられる形になっています。

高速化のための選択

構成を見せる経路と、高速に符号だけ調べる経路¶

evaluate は、これまで構成した Rational の演算を実際に使うため、数体系の積み上げを観察できます。一方、二分法やスツルム列では、多項式の 符号だけを何度も調べます。

そこで sign_at は、同じ整数比をPythonの任意精度整数を使う Fraction へ一時的に写し、ホーナー法で厳密に符号を判定します。浮動小数点の近似ではありませんが、ペアノ表現だけを通す経路でもありません。

def sign_at(self, value: object) -> int:
    point = _as_fraction(cast2r(value))
    result = Fraction(0)
    for coefficient in reversed(self.coefficients):
        result = result * point + _as_fraction(coefficient)
    return (result > 0) - (result < 0)

ここでは、構成を見せる処理と、繰り返し計算を現実的な時間で行う処理を意図的に分けています。

メソッド	優先するもの
`evaluate`	構成した数体系の演算が見えること
`sign_at`	同じ有理数の意味を保ったまま、反復判定を実用的に行うこと

中間値の定理は根が少なくとも一つあること、スツルムの定理は相異なる根の個数を保証します。プログラムの出力だけから定理そのものが証明されるわけではありません。

次へ¶

多項式と区間 (1,2) を結び付け、そこに根が一つだけあると確認できました。次章では、この条件をオブジェクトの不変条件として保存し、区間を半分ずつ狭めます。

代数的実根へ進む →