はじめに

この記事はTUT Advent Calendar 2018 三日目の記事（になる予定だった記事）です．遅れてしまい申し訳ありません．

二日目はT_T3_M5さんの「けんきうのおはなし」でした．生き物はどうやってものを「見て」いるか？という観点でしたが，「人は見たいものを見たいように見る」みたいな言葉を思い出しました．

自己紹介

今年は二つのTUT合同でのAdCということで，開催おめでとうございます．
（そういえば，確かフィンランドの方にもTUTってありましたよね．）

はじめまして（？），「むしぱん」です．東京じゃない方の人間です．プログラミング，数学が好きです．
よく使うプログラミング言語はRust，Juliaです．最近は全然書いてないですが，Kotlin，Elixirも好きです．

私も「けんきうのはなし」をするか？と思いましたが用意に身元を特定できそうなので避けようと思います．

私はプログラミングと数学が好きなので，この記事では数学とプログラミング，両方に関わる話題として，（自然数の）最大公約数を求めるアルゴリズム「ユークリッドの互除法」を多項式に拡張し，（一般化して）Rustで実装した話をします．

ユークリッドの互除法（自然数）

そもそも，一般に言われるユークリッドの互除法は自然数の最大公約数（Greatest Common Divisor）を求めるアルゴリズムです．

自然数 ${a,b}$ に対して，以下の式が成り立つとします．つまり， ${a}$ を ${b}$ で割ったときの商が ${q}$ で，あまりが ${r}$ ということです．

${a= bq + r}$

ここで，最大公約数を求める関数を ${\mathrm{gcd}(x,y)}$ とおくと，以下の式が成り立ちます．

${\mathrm{gcd}(a,b)=\mathrm{gcd}(b,r)}$

ここで重要なのは，常に ${r\leq b}$ が成立する点です．上の式を繰り返し適用することで最大公約数を求めます．このとき， ${r\leq b}$ であることから，徐々に ${a,b}$ の値が小さくなり，最終的にあまり ${r=0}$ となります． あまりが0となったときの ${b}$ の値が最大公約数 となります．

例えば， ${a=28,b=6}$ の例を考えます．最終的に4割る2を計算することとなり，あまりが0であることから，最大公約数は2となります．

f:id:je6bmq:20181204000324p:plain — ユークリッドの互除法による最大公約数の導出（ ${a=28,b=6}$ ）

これをプログラムとして記述すると，次の様に記述できます．
Rustにおける実装例ですが，他の言語でも同様の実装になるかと思います．
（ ${a \leq b}$ だとしても，一度 ${q=0}$ となって大小関係が逆転するのみなので，大小関係のチェックは不要です）

fn gcd(a: u32, b: u32) -> u32 {
     if b == 0 { a } else { gcd(b, a % b) }
}

多項式における最大公約数って？

この記事の目的は，この「ユークリッドの互除法」を多項式に応用することでした．今回は各係数が実数の， ${x}$ についての多項式を考えます．
例えば， ${ x^2+1}$ や， ${x^3 +\sqrt{2}x + 1}$ が挙げられます．

この「（ ${x}$ についての多項式）における最大公約数」を考える前に，一般的な（自然数における）最大公約数とはどういう数か考えます．

28と6の最大公約数，280と588の最大公約数( ${\mathrm{gcd}(288,588)=28}$ )を例に，両整数の素因数分解と最大公約数に着目すると，実は， 最大公約数は共通の素因数の積 になります．

f:id:je6bmq:20181204004508p:plain — 最大公約数は，共通の素因数の積になる

多項式における最大公約数も同様に， 各多項式を素因数分解（に対応する操作をした場合の）共通の素因数（に対応するもの）の積 である，とします．

では，この多項式における素因数分解（に対応する操作）と素因数（に対応するもの）は何なのでしょう．
そもそも，自然数における素因数分解は，最大限細かく素数の積で表す操作（分解）でした．そして，素数は，1と自身以外を約数に持たない数，すなわち，1と自身の積以外に分解できない数でした．

多項式における「分解」として，一番最初に学校で習う操作があります．「因数分解」です．因数分解は，ある多項式を最大限細かく別の多項式の積で表す操作でした．
実はこの因数分解したときの各多項式（例： ${x^2-1=(x+1)(x-1)}$ における ${(x+1)}$ と ${(x-1)}$ ）には「既約多項式」という名前がついています．

そして，自然数における素数と素因数分解が，多項式における既約多項式と因数分解に対応しています．（厳密には対象にしている集合が整数か，自然数か，実数か，によってある多項式が既約多項式であるかどうかは変わってきますが…）

したがって，多項式における最大公約数は，両多項式で共通の既約多項式の積となります．

f:id:je6bmq:20181204012259p:plain — 自然数と多項式の「分解」の対応

ユークリッドの互除法（多項式）

では，自然数におけるユークリッドの互除法を多項式に拡張します．つまり，これによって多項式の最大公約数（共通の既約多項式の積）が求められればよいということです．

多項式 ${f(x),g(x)}$ に対して，以下の式が成り立つとします．つまり， ${f(x)}$ を ${g(x)}$ で割ったときの商が ${q(x)}$ で，あまりが ${r(x)}$ ということです．

${f(x)= g(x)q(x) + r(x)}$

自然数同様，最大公約数を求める関数を ${\mathrm{gcd}(f(x),g(x))}$ とおくと，以下の式が成り立ちます．

${\mathrm{gcd}(f(x),g(x))=\mathrm{gcd}(g(x),r(x))}$

例えば，次の二つの多項式について，互除法を適用します．各々の因数分解から，最大公約数は ${(x-1)(x^2+1)=x^3-x^2+x-1}$ であることがわかります．

${f(x) =2 x^5 - 7 x^4 + 4 x^3 - 4 x^2 + 2 x + 3= (x-1)(x-3)(2x+1)(x^2+1)}$

${g(x) = x^4 -1 = (x-1)(x+1)(x^2+1)}$

では，互除法を使って最大公約数が求めます．

f:id:je6bmq:20181204145738p:plain — 多項式の最大公約数を互助法によって求める例

互除法によると，最大公約数は ${4x^3-4x^2+4x-4}$ だそうです．
あれ？ ${(x-1)(x^2+1)=x^3-x^2+x-1}$ ではない？と思うかもしれませんが実は ${4x^3-4x^2+4x-4=4(x-1)(x^2+1)}$ であるため，両多項式とも ${4x^3-4x^2+4x-4}$ によって割り切れるのです．

このように，多項式では最大公約数の定数倍の多項式ではすべて割り切れてしまうため，実際は既約多項式の最高次の係数を1に正規化して扱うことがあります．
例えば，先ほどの互除法の例の最後の商， ${\frac{1}{4}x+\frac{1}{4}}$ から ${\frac{1}{4}}$ をくくり出して， ${\frac{1}{4}(x+1)}$ と扱えば，定数倍の差を吸収できます．

最大公約数関数を多項式に拡張，実装する（Rust）

ユークリッドの互除法を行う上で自然数，多項式共通の操作は四則演算と剰余演算，余りがゼロであるという判定です．つまり，四則演算と剰余演算ができて，かつゼロが定義されていれば互除法が適用できることが示唆されます．

そこで今回は 新たに多項式を表すPolynomial型を定義し，Polynomial型を実装し，ゼロの概念と等号判定を導入する ことで自然数，多項式両方に適用できる最大公約数関数を実装します．詳細の実装の説明は省略して，概要のみ説明します．詳細はこちらのGistを参照ください．（思いの外，数式パートが長引いたので…）

Rustでは，トレイト（他の言語でいうインターフェースのようなもの）を用いて関数の型制約を表現できるので，次の様に実装できます．

fn gcd<T>(a: T, b: T) -> T
where
    T: Add<Output = T>
        + Sub<Output = T>
        + Mul<Output = T>
        + Rem<Output = T>
        + Zero
        + PartialEq
        + Clone,
{
    if b == <T as Zero>::get_zero() {
        a
    } else {
        gcd(b.clone(), a % b)
    }
}

ここで，型TはAddトレイト，Subトレイト，Mulトレイト，Remトレイト，Zeroトレイト，PartialEqトレイトを実装している必要があることを意味します．そして，各トレイトは次の意味を持ちます．

Addは加算可能，かつ計算結果の型がT
Subは減算可能，かつ計算結果の型がT
Mulは乗算可能，かつ計算結果の型がT
Remは剰余演算可能，かつ計算結果の型がT
PartialEqは == によって同値関係が確かめられること

上記のトレイトは標準ライブラリに搭載されているので，自然数（符号なし整数型）には当然実装されています．したがって，新たにPolynomial型にも実装することとしました．

新たに定義したのはZeroトレイトです．Zeroトレイトは，名前の通りその型におけるゼロを定義したいため，次の様な実装にしました．
符号なし整数（ここでは32bitのみ対象）では0をそのまま返し，Polynomial型においては0次の項の値が0の多項式を返します．

impl Zero for u32 {
    fn get_zero() -> u32 {
        0
    }
}
impl Zero for Polynomial {
    fn get_zero() -> Polynomial {
        Polynomial::new(&vec![0f32], false)
    }
}

実行結果

一般化したgcd関数はもちろん32bit符号なし整数型にもそのまま適用可能ですが，ここでは先ほどの二つの多項式

${f(x) =2 x^5 - 7 x^4 + 4 x^3 - 4 x^2 + 2 x + 3= (x-1)(x-3)(2x+1)(x^2+1)}$

${g(x) = x^4 -1 = (x-1)(x+1)(x^2+1)}$

の最大公約数を計算した結果を示します．ここで，^ （ハット）の後ろの数字が指数を表します．

実際に，先ほどと同様に最大公約数を求めることができています．

f:id:je6bmq:20181204153247p:plain — 一般化したgcd関数（Rust）の実行結果

まとめ

この記事では，最大公約数を求めるアルゴリズム，ユークリッドの互除法を自然数から多項式に拡張して，かつプログラミング言語Rustで実装できることを確認しました．少なくともRustではトレイトで型の制約を表現できるので， 数学における概念と満たすべき性質 がうまく対応付けられて嬉しいと思っています．

多項式が既約多項式かどうか，といったような概念は有名な理論である「 ガロア理論 」とも関わりが強く，私がこの実装を試みたきっかけにもなっているので，興味のある方は調べてみてはいかがでしょう．（というか私も教えてほしい．）

本来の四日目の担当はyanteneさんで，「何書くかまだ何も考えてない！遅れたらごめんなさい！」とのことです．この記事自体が遅れてしまったので何も言えませんが期待しています．

je6bmq’s diary

ユークリッドの互除法を多項式に拡張して，（Rustで）実装する