📐
geobucketの解説

2025/09/23に公開
 概要多くの多項式を足すことを考える。(ソート済み配列のマージと似た操作になる)

項の数の合計が N 個の場合, 最悪で O(N^2) かかってしまう。[1]

そこで, geobucket を使うと最悪でも O(N\log N) で済むため高速である。

 準備この節はざっくりした解説しかしないので詳しく知りたい人は以下の本を読むことをオススメします。

https://www.maruzen-publishing.co.jp/book/b10121907.html

 多変数多項式変数が複数(または1つ)ある多項式を多変数多項式と呼ぶ。
変数をx_1, x_2, \dots, x_nとすると, Kを係数体[2]とするn変数多項式は

\sum_{e_1, e_2, ..., e_n \in \mathbb{Z}_{\geq 0}} c_{e_1,e_2,\dots,e_n}x_1^{e_1} x_2^{e_2}\cdots x_n^{e_n}
という形で表せる。ここで、c_{e_1,e_2,\dots,e_n} \in Kを係数、x_1^{e_1} x_2^{e_2}\cdots x_n^{e_n}を単項式と呼ぶ。また、係数と単項式を合わせたものを項と呼ぶ。多変数多項式では0でない係数は有限個である必要がある。

多変数多項式の例としては x + y + z, x^2 + y^2 - 1, xyz, x^2 + 2x + 1, 0 などが挙げられる.

指数関数(e^x)や三角関数(\sin(x))などが含まれている場合は多変数多項式ではない.

 単項式順序一変数多項式では 1 < x < x^2 < x^3 < \cdotsという自然な(全)順序が考えられる。

しかし、多変数多項式の場合は同様に考えても上手く行かない。例えばxとyはどちらが大きいのかやxy^2とx^2zはどちらが大きいのかなどは自明でない。
まずは一変数多項式の(全)順序の性質を考えると(天下り的ではあるが)以下のような性質があることが分かる。
a < b \Longrightarrow \forall c, ac < bc.
単調減少列(x_1 > x_2 > x_3 > \cdots)は有限の長さで停止する[3]
1の性質は両辺に単項式を掛けた場合順序が保たれるという意味である。これはそのまま多変数多項式でも成り立つようにできる。

2の性質は何がうれしいのか分かりにくいと思われるが、アルゴリズムの停止性を証明する際に鍵となる性質になる。これもそのまま多変数多項式でも成り立つようにできる。
よって多変数多項式の単項式の(全)順序が満たすべき性質は以下となる。[4]
a < b \Longrightarrow \forall c, ac < bc.
単調減少列(x_1 > x_2 > x_3 > \cdots)は有限の長さで停止する
よく使われる単項式順序として、辞書式順序(lex)、次数付き辞書式順序(grlex)、次数付き逆辞書式順序(grevlex)や消去順序(elimination order)がある。
多項式に含まれる項の内、単項式が最も大きい項を先頭項(leading term)と呼ぶ。

また先頭項の係数を先頭係数(leading coefficient), 先頭項の単項式を先頭単項式(leading monomial)と呼ぶ。

 割り算アルゴリズム疑似コードを載せると以下のようになる。
fn division(mut f: Poly, g: &[Poly]) -> (Vec<Poly>, Poly) {
    let mut q = vec![Poly::zero(); g.len()];
    let mut r = Poly::zero();
    while !f.is_zero() {
        // fの先頭単項式を割り切るgの先頭単項式を探す
        if let Some(i) = g.iter().position(|gi| gi.lm().divides(f.lm())) {
            // 見つかった場合、商に新しい項を追加
            let new_q = f.lt() / g[i].lt();
            f -= g[i] * new_q;
            q[i] += new_q;
        } else {
            // 見つからなかったので余りに追加
            r += f.pop_lt().unwrap();
        }
    }
    (q, r)
}
この割り算アルゴリズムの f -= g[i] * new_q; がボトルネックとなる。

fは長い多項式になりやすいため、毎回足し算(マージ)を行うと時間がかかってしまう。

while ループは反復ごとにfの先頭項が小さくなるため、性質2より有限回のループで停止することが証明できる。

 geobucket本題のデータ構造について解説する。元論文は以下である。

https://www.sciencedirect.com/science/article/pii/S0747717197901760

 基本的な考え方長い多項式に短い多項式を繰り返し足しこむ場合に非効率で、同程度の長さの多項式を足すのは効率的である。よって、前者を避け、なるべく後者の状況にすることで高速化する。

また、除算を行っている途中では先頭項さえ分かっていれば良い。つまり、普通の多項式と異なるデータ構造で保持しても良い。

 詳細長さごとに多項式を保存するバケットを用意する。

バケット0には長さ1の多項式を、バケット1には長さ2,3の多項式を、バケット2には長さ4~7の多項式をと長さの上限が約2倍になるようにする。(長さの上限が幾何級数(geometric series)であることからGeobucketと名前が付いている)

こうすると実用上バケットが64もあれば十分である。(項の数が2^{64}個の多項式を保持するメモリも、処理する時間も現実的でないため)
なお元論文では2倍にしていくだけではなく、d倍(dは2以上の整数)にしていく手法として述べられている。

元論文ではd=4が勧められているが、d=2が採用されていることが多い気がする。

 実装疑似コードは以下のようになる。
struct Geobucket {
    bucket: [Poly; 64],
}
impl Default for Geobucket {
    fn default() -> Self {
        Self {
            bucket: [Poly::zero(); 64],
        }
    }
}
#[auto_impl_ops::auto_ops]
impl AddAssign<Poly> for Geobucket
{
    fn add_assign(&mut self, mut other: Poly) {
        if other.is_zero() {
            return;
        }
        debug_assert!(!other.is_empty());
        let s = other.len().ilog2() as usize;
        for (i, b) in self.bucket.iter_mut().enumerate().skip(s) {
            if b.is_zero() {
                // bucket[i] が空なら保存して終わり
                *b = other;
                return;
            } else {
                // bucket[i] が空でない場合は足す
                let mut p = Poly::ZERO;
                core::mem::swap(&mut p, b);
                // 長さの比率は(多くの場合)2未満
                other += p;
                // 現在のバケットの長さ上限に引っかかっていないか確認
                if other.len().ilog2() as usize == i {
                    // 収まっているので bucket[i] に保存して終わり
                    *b = other;
                    return;
                }
            }
            // 1つ上のバケットを試す
        }
        panic!("polynomial is too long")
    }
}
impl Geobucket {
    // 1つの多項式へ変換
    pub fn into_polynomial(self) -> Poly
    {
        // 短い方から全部足す
        self.bucket.into_iter().sum()
    }
    // 先頭項を探す
    fn lt_aux(&mut self) -> Option<usize>
    where
        M: Clone + Ord,
        K: Clone + Zero + AddAssign<K>,
        for<'x> &'x K: Add<Output = K>,
    {
        let mut max = None;
        let n = self.bucket.len();
        for i in 0..n {
            if self.bucket[i].is_zero() {
                continue;
            } else if max.is_none() {
                max = Some(i);
                continue;
            }
            let max_i = max.unwrap();
            let p = self.bucket[i].lm().unwrap();
            let q = self.bucket[max_i].lm().unwrap();
            use core::cmp::Ordering::*;
            match p.cmp(q) {
                Greater => {
                    // 古い先頭係数が0だった場合その項を捨てる
                    if self.bucket[max_i].lc().unwrap().is_zero() {
                        self.bucket[max_i].pop_lt();
                    }
                    // より大きい項が見つかったので置き換え
                    max = Some(i);
                }
                Equal => {
                    // 同じ項が見つかったので係数を足し合わせる
                    let (c, _) = self.bucket[i].pop_lt().unwrap();
                    *self.bucket[max_i].lc_mut().unwrap() += c;
                }
                Less => (),
            }
        }
        max
    }
    pub fn lt(&mut self) -> Option<usize>
    {
        while let Some(i) = self.lt_aux() {
            let lc = self.bucket[i].lc().unwrap();
            // 係数を足し合わせたので0になっているかもしれない
            if lc.is_zero() {
                // 0の場合は捨てて、もう一度探す
                // 繰り返し回数は最大で保持している項数の半分になる(打消しが起きるには少なくとも2項必要で、打消しが起きた項は取り除かれるため)
                self.bucket[i].pop_lt();
            } else {
                return Some(i);
            }
        }
        None
    }
}

 応用グレブナー基底(Gröbner basis)を求めるブッフベルガー(Buchberger)アルゴリズムの中で鍵となるのがS多項式の剰余を求める処理になります。よって、geobucketを使って割り算アルゴリズムを高速化することでグレブナー基底を高速に求められるようになります。

脚注
単項式をN個左結合で足していくことを考えよ ↩︎
平たく言えば足し算、引き算、掛け算、(0以外による)割り算が自由にできる集合 ↩︎
整列順序 ↩︎
実際このように定義すると上手くいくことが知られている ↩︎
概要

準備

多変数多項式

単項式順序

割り算アルゴリズム

geobucket

基本的な考え方

詳細

実装

応用

Discussion
