グレブナー基底を使って数独を解く

これは「FOLIO Advent Calendar 2024」18日目の記事です。

---

数独は9×9のグリッドに1から9までの数字を埋める論理パズルです。グリッドは縦横9列ずつの合計81マスに分かれ、さらに3×3の小さなブロックに区切られています。ルールは単純で、各行、各列、各3×3ブロックには数字1から9が重複せず1回ずつ入るようにします。初期状態ではいくつかのマスに数字が配置されており、その数字をヒントに空白を埋めていきます。

例えば以下のような数独の問題を考えましょう。

┏━━━┯━━━┯━━━┳━━━┯━━━┯━━━┳━━━┯━━━┯━━━┓
┃ 3 │ 9 │ 1 ┃ 4 │ 8 │   ┃ 6 │ 2 │ 7 ┃
┠───┼───┼───╂───┼───┼───╂───┼───┼───┨
┃ 2 │ 7 │ 6 ┃   │ 9 │ 1 ┃ 4 │ 8 │ 5 ┃
┠───┼───┼───╂───┼───┼───╂───┼───┼───┨
┃ 8 │ 5 │   ┃ 2 │ 7 │ 6 ┃ 3 │ 9 │   ┃
┣━━━┿━━━┿━━━╋━━━┿━━━┿━━━╋━━━┿━━━┿━━━┫
┃ 9 │ 1 │   ┃ 8 │ 5 │ 4 ┃ 2 │ 7 │ 6 ┃
┠───┼───┼───╂───┼───┼───╂───┼───┼───┨
┃ 5 │ 4 │ 8 ┃ 7 │ 6 │   ┃ 9 │ 1 │ 3 ┃
┠───┼───┼───╂───┼───┼───╂───┼───┼───┨
┃ 7 │   │ 2 ┃ 9 │ 1 │ 3 ┃   │ 5 │ 4 ┃
┣━━━┿━━━┿━━━╋━━━┿━━━┿━━━╋━━━┿━━━┿━━━┫
┃ 1 │ 3 │ 9 ┃ 5 │   │ 8 ┃   │   │ 2 ┃
┠───┼───┼───╂───┼───┼───╂───┼───┼───┨
┃   │ 2 │ 7 ┃ 1 │ 3 │ 9 ┃ 5 │ 4 │ 8 ┃
┠───┼───┼───╂───┼───┼───╂───┼───┼───┨
┃ 4 │ 8 │ 5 ┃ 6 │ 2 │   ┃ 1 │ 3 │ 9 ┃
┗━━━┷━━━┷━━━┻━━━┷━━━┷━━━┻━━━┷━━━┷━━━┛

1行目に注目すると1マスだけ空白になっており、他のマスに使われていない数字から空白マスには5が入ることが分かります。

┏━━━┯━━━┯━━━┳━━━┯━━━┯━━━┳━━━┯━━━┯━━━┓
┃ 3 │ 9 │ 1 ┃ 4 │ 8 │   ┃ 6 │ 2 │ 7 ┃
┠───┼───┼───╂───┼───┼───╂───┼───┼───┨

また1列目に注目すると1マスだけ空白になっており、他のマスに使われていない数字から空白マスには6が入ることが分かります。

┏━━━┯
┃ 3 │
┠───┼
┃ 2 │
┠───┼
┃ 8 │
┣━━━┿
┃ 9 │
┠───┼
┃ 5 │
┠───┼
┃ 7 │
┣━━━┿
┃ 1 │
┠───┼
┃   │
┠───┼
┃ 4 │
┗━━━┷

同様に左上の3×3のブロックに注目すると1マスだけ空白になっており、他のマスに使われていない数字から空白マスには4が入ることが分かります。

┏━━━┯━━━┯━━━┳
┃ 3 │ 9 │ 1 ┃
┠───┼───┼───╂
┃ 2 │ 7 │ 6 ┃
┠───┼───┼───╂
┃ 8 │ 5 │   ┃
┣━━━┿━━━┿━━━╋

数独はこのような推論を繰り返し、時に空白マスが2つ以上ある場合は数字の仮置きをして推論を進め矛盾したら初めに戻るバックトラッキングも使いながら、最終的に全ての空白マスを埋めることを目指すパズルです。

通常、数独をプログラムで解くには探索問題として考えます。

• 深さ優先探索で数独を解く #Haskell - Qiita
• Sudoku - HaskellWiki
• sudoku generation in Haskell

しかし本稿では探索問題ではなく、数独の問題を多項式の言葉に翻訳してグレブナー基底を求めることにより解くといった一風変わった方法を紹介したいと思います。

多項式の割り算

まずは簡単な自然数の割り算を思い出しましょう。7を2で割るというのは、2に掛け算した結果を7から引くと2より小さくなる数（つまり3）を商とし、7から2×3を引いた数（つまり1）を余りとするものでした。

次に一変数多項式の割り算です。本稿を通して多項式は有理数を係数に持つものを考えます。多項式 $f$ を 多項式 $g$ で割る場合、自然数の時と同様に $g$ に掛け算した結果を $f$ から引いたものと $g$ の大きさを比較したいのですが、そのためにまずは多項式の大きさを測る方法を定義する必要があります。そこで多項式 $f$ の次数 ${\rm deg}(f)$ を

$$f = a_nx^n + a_{n-1}x^{n-1} + \cdots + a_0$$

とした時に

$${\rm deg}(f) = n$$

と定義します。また $f=0$ の場合は ${\rm deg}(f)=-\infty$ とします。この次数を使って一変数多項式の割り算を以下のように定義します。多項式 $f$ と0でない多項式 $g$ に対して

$$f = qg+r$$

かつ ${\rm deg}(r) < {\rm deg}(g)$ を満たすような多項式 $q, r$ が存在し、$q$ を商、$r$を余りと呼ぶ。

例えば $f = x^3+2x^2-11x+7$ を $g = x+5$ で割ることを考えましょう。$q_1 = x^2$ という単項式（一つの項しか持たない多項式）を考えると $f - q_1g = -3x^2-11x+7$ となり $f$ より次数の低い多項式が作れます。同様に $q_2 = -3x$ とすると $f - q_1g - q_2g = 4x+7$ となり、$q_3 = 4$ とすると $f - q_1g - q_2g - q_3g = -13$ となります。よって

$$f = (q_1+q_2+q_3)g - 13 = (x^2-3x+4)g-13$$

となり $f$ を $g$ で割った商は $x^2-3x+4$ で余りは $-13$ であることが分かりました。

自然数の割り算に比べて一変数多項式の割り算は少し複雑に見えますが、実は両者はユークリッド整域という言葉を使って統一的に理解できます。しかし次の多変数多項式になるとそう簡単には行きません。

多変数多項式 $f$ を複数の多変数多項式 $g_1, g_2, \dots, g_n$ で割ることを考えましょう。いきなり割る数（多項式）が複数に増えたのでびっくりするかもしれません。実は今までの自然数や一変数多項式では複数の割る数を考えても、余りを求める計算はそれらの最大公約数による割り算に帰着できるのであまり意味がなかったのです。しかし多変数多項式の場合は単項イデアル環とならないためこのような考え方ができず、複数の割る数による割り算をちゃんと定義する必要があります。

多変数多項式の場合にも割り算を考えるために多項式の大きさを測る方法を定義する必要があります。一変数多項式の時と同様に次数の最も大きな項の次数としたいところですがそもそも多変数だと単項式の大小から定義しなければなりません。例えば以下の3変数 $x,y,z$ の単項式はどれが最も大きいでしょうか？

$$xy^3z^2,\ \ x^3y^2z,\ \ z^4$$

このような順序は単項式順序と呼ばれ複数の定義の仕方があることが知られています。例えばその一つである辞書式順序は、まず変数の順序を定め（例えば $x>y>z$）大きい変数から順番に変数の次数を比較して最初に異なる次数の大小で単項式の大小を決めるという方法です。この方法だと先程の例の単項式は降順に

$$x^3y^2z,\ \ xy^3z^2,\ \ z^4$$

と並べることができます。以下では単項式順序としてこの辞書式順序を考えることとします。

ここで多項式に関する記号をいくつか定義しておきます。

1. ${\rm LT}(f)$ は多項式 $f$ の最も大きな単項式の項を表す
2. ${\rm LC}(f)$ は ${\rm LT}(f)$ の係数部分を表す
3. ${\rm LM}(f)$ は ${\rm LT}(f)$ の変数部分を表す

定義から ${\rm LT}(f) = {\rm LC}(f){\rm LM}(f)$ となります。

例えば

$$f = 5x^3y^2z + 2xy^3z^2 + 3z^4$$

の場合 ${\rm LT}(f) = 5x^3y^2z, {\rm LC}(f) = 5, {\rm LM}(f) = x^3y^2z$ となります。

さていよいよ多変数多項式の割り算を考えましょう。多変数多項式 $f$ を0でない多変数多項式 $g_1, g_2, \dots, g_n$ で割るとは以下の条件を満たす多項式 $q_1, q_2, \dots, q_n$ と余り $r$ を見つけることとします。

1. $f = q_1g_1 + q_2g_2 + \dots + q_ng_n + r$
2. $r$ に現れるどの単項式も ${\rm LM}(g_1), {\rm LM}(g_2), \dots {\rm LM}(g_n)$ のいずれでも割り切れない

このような条件を満たす多項式は以下の手順により見つけることができます。

1. $f = f, q_1 = 0, q_2 = 0, \dots q_n = 0, r = 0$ と置く
2. ${\rm LM}(g_i)$ が ${\rm LM}(f)$ を割り切る最小の $i$ が存在すれば、$f - \frac{{\rm LT}(f)}{{\rm LT}(g_i)}g_i$ を新たに $f$ とし、$q_i + \frac{{\rm LT}(f)}{{\rm LT}(g_i)}$ を新たに $q_i$ とする
3. 2の条件を満たすものが存在しなければ、$f - {\rm LT}(f)$ を新たに $f$ とし、$r + {\rm LT}(f)$ を新たに $r$ とする
4. 2,3 を $f = 0$ となるまで繰り返す

例えば $f = x^3y^2+xy+x$ を $g_1 = y^2+1, g_2 = xy+1$ で割ることを考えます。まず ${\rm LM}(f)$ は ${\rm LM}(g_1)$ で割り切れるので、$f$ を新たに

$$\begin{matrix} f - \frac{{\rm LT}(f)}{{\rm LT}(g_1)}g_1 &=& x^3y^2+xy+x - \frac{x^3y^2}{y^2}(y^2+1) \\ &=& x^3y^2+xy+x - x^3(y^2+1) \\ &=& x^3y^2+xy+x - x^3y^2-x^3 \\ &=& -x^3+xy+x \\ \end{matrix}$$

と置き $q_1$ を新たに

$$\begin{matrix} q_1 + \frac{{\rm LT}(f)}{{\rm LT}(g_1)} &=& 0 + \frac{x^3y^2}{y^2} \\ &=& x^3 \\ \end{matrix}$$

と置きます。

次に ${\rm LM}(f) = x^3$ は ${\rm LM}(g_1), {\rm LM}(g_2)$ いずれでも割り切れないので、$f$ を新たに

$$f - {\rm LT}(f) = xy+x$$

と置き、$r$ を新たに

$$r + {\rm LT}(f) = 0 + (-x^3) = -x^3$$

と置きます。

次に ${\rm LM}(f) = xy$ は ${\rm LM}(g_2)$ で割り切れるので、$f$ を新たに

$$\begin{matrix} f - \frac{{\rm LT}(f)}{{\rm LT}(g_2)}g_2 &=& xy+x - \frac{xy}{xy}(xy+1) \\ &=& xy+x - (xy+1) \\ &=& x-1 \\ \end{matrix}$$

と置き $q_2$ を新たに

$$\begin{matrix} q_2 + \frac{{\rm LT}(f)}{{\rm LT}(g_2)} &=& 0 + \frac{xy}{xy} \\ &=& 1 \\ \end{matrix}$$

と置きます。

次に ${\rm LM}(f) = x$ は ${\rm LM}(g_1), {\rm LM}(g_2)$ いずれでも割り切れないので、$f$ を新たに

$$f - {\rm LT}(f) = -1$$

と置き、$r$ を新たに

$$r + {\rm LT}(f) = -x^3 + x$$

と置きます。

再び ${\rm LM}(f) = 1 (= x^0y^0)$ は ${\rm LM}(g_1), {\rm LM}(g_2)$ いずれでも割り切れないので、$f$ を新たに

$$f - {\rm LT}(f) = 0$$

と置き、$r$ を新たに

$$r + {\rm LT}(f) = -x^3 + x + -1$$

と置きます。$f = 0$ となったのでここで終了です。

以上により $f$ の $g_1, g_2$ による割り算は

$$f = q_1g_1 + q_2g_2 + r = x^3(y^2+1) + (xy+1) + (-x^3 + x + -1)$$

となることが分かりました。

この手続きを（いちいち書き下すのは大変なので…）プログラムで実装してみましょう。

まずは単項式の実装です。単項式の表現には vector-sized の Vector を使うことにします。Vector が持つ型レベルの長さは変数全体の数を表し、要素の値は0以上の整数にしてその変数の何乗を考えているか（$x^n$の$n$）を表すことにします。Vector の Ord のインスタンスは標準で辞書式順序になっているので特に追加で実装する必要はありません。

import GHC.TypeLits (KnownNat)

import qualified Data.Vector.Sized as VS

-- | n変数の単項式
type Monomial n = VS.Vector n Int

-- | 定数項1を表す単項式
one :: KnownNat n => Monomial n
one = VS.replicate 0

-- | n番目の変数を表す単項式
var :: KnownNat n => Int -> Monomial n
var n = VS.generate (\k -> if fromIntegral k == n then 1 else 0)

-- | 単項式が割れるかどうか判定する
divisible :: Monomial n -> Monomial n -> Bool
divisible f g = VS.and $ VS.zipWith (<=) f g

-- | 単項式の割り算
divM :: Monomial n -> Monomial n -> Monomial n
divM = VS.zipWith (-)

-- | 単項式の最小公倍数
lcm' :: Monomial n -> Monomial n -> Monomial n
lcm' = VS.zipWith max

次にこの単項式を使って多項式を実装します。多項式は単項式（と係数の組）の集合なのでリストで表してもよいですが、後の計算で最高次の単項式を取り出す処理が多用されるためリストを毎回ソートすると効率が悪いので、今回は単項式の順序により構成された木構造（containersのMap）を使用します。

import Data.Ratio (Ratio)

import Data.Map (Map)
import qualified Data.Map as Map

-- | 有理数
type Q = Ratio Int

-- | 係数付きの単項式
type Term n = (Q, Monomial n)

-- | n変数の多項式
-- | 単項式が降順に並ぶようにする
type Polynomial n = Map (Monomial n) Q

-- | 0を表す多項式
zero :: KnownNat n => Polynomial n
zero = Map.singleton one 0

-- | 単項式と多項式の掛け算
times :: KnownNat n => Term n -> Polynomial n -> Polynomial n
times (0, _) = const zero
times (c, m) = Map.map (*c) . Map.mapKeysMonotonic (VS.zipWith (+) m)

-- | 多項式の足し算
add :: KnownNat n => Polynomial n -> Polynomial n -> Polynomial n
add p1 p2 =
  let p = Map.filter (/=0) $ Map.unionWith (+) p1 p2
   in if Map.null p then zero else p

-- | 多項式の引き算
minus :: KnownNat n => Polynomial n -> Polynomial n -> Polynomial n
minus p1 p2 = p1 `add` Map.map negate p2

また多項式を表示する際に分かりやすく整形して表示する関数も用意しておきましょう。

-- | 単項式を表示する関数を受け取って多項式を表示する関数
showPoly :: KnownNat n => (Monomial n -> String) -> Polynomial n -> String
showPoly showMono p =
  let ((t0, c0) : ts) = Map.toDescList p
   in showCoef c0 t0 ++ showMono t0 ++ concat (map showTerm ts)
  where
    showCoef c t
      | c == 1 = if t == one then "1" else ""
      | c == -1 = if t == one then "" else "-"
      | denominator c == 1 = show (numerator c)
      | otherwise = show (numerator c) ++ "/" ++ show (denominator c)
    showTerm (t, c) = (if c > 0 then " + " ++ showCoef c t else " - " ++ showCoef (-c) t) ++ showMono t

-- | 単項式を表示する関数（変数は26文字まで x, y, z, w, a, b, c, ... と表す）
showMonoXYZ :: KnownNat n => Monomial n -> String
showMonoXYZ = concatMap (\(x, n) -> if n == 0 then "" else if n == 1 then [x] else x : '^' : show n) . zip ("xyzw" ++ ['a'..'v']) . VS.toList

-- | 単項式を表示する関数（変数は x_0, x_1, x_2, x_3, ... と表す）
showMonoXN :: KnownNat n => Monomial n -> String
showMonoXN = concatMap (\(i, n) -> if n == 0 then "" else "x_" ++ show i ++ (if n == 1 then "" else '^' : show n)) . zip [0::Integer ..] . VS.toList

-- | 多項式を表示する関数（変数は26文字まで x, y, z, w, a, b, c, ... と表す）
showPolyXYZ :: KnownNat n => Polynomial n -> String
showPolyXYZ = showPoly showMonoXYZ

-- | 多項式を表示する関数（変数は x_0, x_1, x_2, x_3, ... と表す）
showPolyXN :: KnownNat n => Polynomial n -> String
showPolyXN = showPoly showMonoXN

ここまで準備すれば多項式の割り算は以下のように実装できます。

{-# LANGUAGE TupleSections #-}

import Data.Functor ((<&>))

-- | List の Zipper から List を得る
integrate :: ([a], [a]) -> [a]
integrate ([], ys) = ys
integrate (x:xs, ys) = integrate (xs, x:ys)

-- | 多項式の割り算
divP :: KnownNat n => Polynomial n -> [Polynomial n] -> ([Polynomial n], Polynomial n)
divP f0 fs0 = go1 f0 (fs0 <&> (,zero)) zero
  where
  -- f が 0 になるまで繰り返す処理
  go1 f fgs r =
    let (f', fgs', r') = go2 f ([], fgs) r
     in if f' == zero then (map snd fgs', r') else go1 f' fgs' r'
  -- 割り算における場合分け②のケース
  go2 f (fgs', fg@(fj, gj):fgs) r =
    let (t, c) = Map.findMax f
        (tj, cj) = Map.findMax fj
        (t', c') = (t `divM` tj, c / cj)
     in if divisible tj t
          then (f `minus` ((c', t') `times` fj), integrate (fgs', (fj, gj `add` Map.singleton t' c'):fgs), r)
          else go2 f (fg:fgs', fgs) r
  -- 割り算における場合分け③のケース
  go2 f (gs', []) r =
    let (t, c) = Map.findMax f
        p = Map.singleton t c
     in (f `minus` p, integrate (gs', []), r `add` p)

実装した多項式の割り算 divP を使って先程の例の計算結果が再現できるか試してみましょう。

{-# LANGUAGE DataKinds #-}

example1 = do
  let vec ls = fromJust $ VS.fromList ls
      f  = Map.fromList [(vec [3, 2], 1), (vec [1, 1], 1), (vec [1, 0] :: Monomial 2, 1)] -- f = x^3y^2 + xy + x
      g1 = Map.fromList [(vec [0, 2], 1), (vec [0, 0], 1)] -- g_1 = y^2 + 1
      g2 = Map.fromList [(vec [1, 1], 1), (vec [0, 0], 1)] -- g_2 = xy + 1
      (qs, r) = f `divP` [g1, g2]
  putStrLn $ "q1 = " ++ showPolyXYZ (qs !! 0)
  putStrLn $ "q2 = " ++ showPolyXYZ (qs !! 1)
  putStrLn $ "r  = " ++ showPolyXYZ r

> example1
q1 = x^3
q2 = 1
r  = -x^3 + x - 1

うまく動いてますね👏

ところで多項式の割り算の計算手順は割る数の順番に依存していました。その影響を見るために先程の例で割る数の順番を入れ替えた場合の計算結果を見てみましょう。

example2 = do
  let vec ls = fromJust $ VS.fromList ls
      f  = Map.fromList [(vec [3, 2], 1), (vec [1, 1], 1), (vec [1, 0] :: Monomial 2, 1)] -- f = x^3y^2 + xy + x
      g1 = Map.fromList [(vec [0, 2], 1), (vec [0, 0], 1)] -- g_1 = y^2 + 1
      g2 = Map.fromList [(vec [1, 1], 1), (vec [0, 0], 1)] -- g_2 = xy + 1
      (qs, r) = f `divP` [g2, g1]  -- 割る数の順番を入れ替えている
  putStrLn $ "q1 = " ++ showPolyXYZ (qs !! 0)
  putStrLn $ "q2 = " ++ showPolyXYZ (qs !! 1)
  putStrLn $ "r  = " ++ showPolyXYZ r

> example2
q1 = x^2y - x + 1
q2 = 0
r  = 2x - 1

先程とは異なる計算結果になってしまいました。余りの値の計算結果が計算順序に依存して異なってしまうのはあまり望ましい性質ではありません。

グレブナー基底

余りの値が一意に定まるような多項式の割り算を考えることはできないでしょうか？実は割る数の組を等価なものに変形することで割り算の余りが一意になるような組が見つけられることが知られています。このような組のことを グレブナー基底 と呼びます。ここで多項式 $F = \{f_1, f_2, \dots f_n\}$ と $G = \{g_1, g_2, \dots, g_m\}$ が等価であるとは、連立方程式 $\{f_1=0, f_2=0, \dots f_n=0\}$ と $\{g_1=0, g_2=0, \dots g_m=0\}$ の解の集合が等しい（これは $F$ が生成するイデアルと $G$ が生成するイデアルが等しいことに同じ）という意味です。

与えられた多項式の集合 $F=\{f_1, f_2, \dots, f_n\}$ からグレブナー基底を生成する方法としてブッフベルガーアルゴリズムという手法があります。これは以下のような手順で構成されます。

1. $G = F$ と置く
2. $G$ から2つの多項式 $f_i, f_j$ を取り出す
3. $f_i$ と $f_j$ の最高次の単項式 ${\rm LM}(f_i)$ と ${\rm LM}(f_j)$の最小公倍数を $a_{ij}$ として $S_{ij} = \frac{a_{ij}}{{\rm LT}(f_i)}f_i - \frac{a_{ij}}{{\rm LT}(f_j)}f_j$ と計算する（これは S多項式 と呼ばれる）
4. $S_{ij}$ を $G$ で割って余りを計算し、余りが0でなければ $G$ を新たに $G\cup\{S_{ij}\}$ と置く
5. 2~4 を $G$ の任意の2つの元に対するS多項式が0になるまで繰り返す

この手続きによって得られる多項式の集合 $G$ はグレブナー基底になっています。またこの手続きが有限回で停止することはディクソンの補題によって保証されています。

このアルゴリズムをプログラムで実装すると以下のようになります。多項式の集合の要素の数は計算の中で変化していくので IntMap を使って可変にできるように番号と多項式の対応付けを管理しています。

import qualified Data.IntMap as IntMap

-- | S多項式の計算
sPolynomial :: KnownNat n => Polynomial n -> Polynomial n -> Polynomial n
sPolynomial fi fj =
  let (gi, ai) = Map.findMax fi
      (gj, aj) = Map.findMax fj
      gij = lcm' gi gj
   in ((recip ai, gij `divM` gi) `times` fi) `minus` ((recip aj, gij `divM` gj) `times` fj)

-- | Buchberger アルゴリズム
buchberger :: KnownNat n => [Polynomial n] -> [Polynomial n]
buchberger gs =
  let gs' = IntMap.fromList $ zip [0..] gs
      n = length gs
      pairs = [(i, j) | i <- [0..n-1], j <- [i+1..n-1]]
   in go gs' n pairs
  where
  go gs' _ [] = IntMap.elems gs'
  go gs' n ((i,j):pairs) =
    let fi = gs' IntMap.! i
        fj = gs' IntMap.! j
        s = sPolynomial fi fj
        (_, r) = divP s (IntMap.elems gs')
     in if r == zero
          then go gs' n pairs
          else go (IntMap.insert n r gs') (n+1) (pairs ++ ([0..n-1] <&> (,n)))

実は余りの一意性を与えるという条件だけではグレブナー基底は一意には定まりません。等価なグレブナー基底の中で唯一定まる多項式の組を考えるために、極小グレブナー基底と被約グレブナー基底という概念を考えます。

まず極小グレブナー基底は全ての多項式がモニック（最高次係数が1）であり、任意の元がそれ以外の多項式で割り切れないようなものです。これはグレブナー基底から以下のように計算できます。

-- | グレブナー基底から極小グレブナー基底を求める
minimize :: KnownNat n => [Polynomial n] -> [Polynomial n]
minimize gs0 = go ([], map monicize gs0)
  where
  monicize f =
    let (_, c) = Map.findMax f
     in (recip c, one) `times` f
  go gs@(_, []) = integrate gs
  go (gs', g:gs) =
    let (_, r) = g `divP` integrate (gs', gs)
     in go (if r == zero then gs' else g:gs', gs)

次に被約グレブナー基底は極小グレブナー基底であって、任意の元に現れるどの単項式もそれ以外の多項式で割り切れないようなものです。これは極小グレブナー基底から以下のように計算できます。

-- | 極小グレブナー基底から被約グレブナー基底を求める
reduce :: KnownNat n => [Polynomial n] -> [Polynomial n]
reduce gs0 = go ([], gs0)
  where
  go gs@(_, []) = integrate gs
  go (gs', g:gs) =
    let (_, r) = g `divP` integrate (gs', gs)
     in go (if r == zero then gs' else r:gs', gs)

この被約グレブナー基底はイデアルと項順序により一意に定まることが知られています。

先程の例の被約グレブナー基底を求めてみましょう。

example3 = do
  let vec ls = fromJust $ VS.fromList ls
      g1 = Map.fromList [(vec [0, 2], 1), (vec [0, 0] :: Monomial 2, 1)] -- g_1 = y^2 + 1
      g2 = Map.fromList [(vec [1, 1], 1), (vec [0, 0], 1)] -- g_2 = xy + 1
      gs = reduce . minimize $ buchberger [g1, g2]
  mapM_ (putStrLn . showPolyXYZ) gs

> example3
y^2 + 1
x - y

式自体は変わっていますが得られた被約グレブナー基底が定める解（$=0$とした時の答え）と元の多項式の組が定める解は変わっていません。得られた被約グレブナー基底による割り算の余りが順序によらないことを確認してみましょう。

example4 = do
  let vec ls = fromJust $ VS.fromList ls
      f  = Map.fromList [(vec [3, 2], 1), (vec [1, 1], 1), (vec [1, 0] :: Monomial 2, 1)] -- f = x^3y^2 + xy + x
      g1 = Map.fromList [(vec [0, 2], 1), (vec [0, 0], 1)] -- g_1 = y^2 + 1
      g2 = Map.fromList [(vec [1, 0], 1), (vec [0, 1], -1)] -- g_2 = x - y
      (_, r1) = f `divP` [g1, g2]
      (_, r2) = f `divP` [g2, g1]  -- 割る数の順番を入れ替えている
  putStrLn $ "r1 = " ++ showPolyXYZ r1
  putStrLn $ "r2 = " ++ showPolyXYZ r2

> example 4
r1 = 2y - 1
r2 = 2y - 1

どちらも同じ結果になっていますね👏

グレブナー基底の重要な性質の一つに消去定理があります。消去定理から消去順序という適切な項順序を選べば（辞書式順序もその一つ）、もし多項式の集合が唯一の解を持つならば各変数からその解を引いた一次式 $x_i - a_i$ がグレブナー基底に現れるということが導かれます。つまりグレブナー基底を求めることにより連立方程式の解を計算することができるのです。

数独を解く

いよいよグレブナー基底を使って数独を解くことを考えましょう。そのためにまずは数独のルールを多項式の言葉に翻訳していきます。

まず数独の各マスに以下のように番号を割り振り、$i$ 番目のマスに入る数を変数 $x_i$ で表します。

┏━━━━┯━━━━┯━━━━┳━━━━┯━━━━┯━━━━┳━━━━┯━━━━┯━━━━┓
┃  0 │  1 │  2 ┃  3 │  4 │  5 ┃  6 │  7 │  8 ┃
┠────┼────┼────╂────┼────┼────╂────┼────┼────┨
┃  9 │ 10 │ 11 ┃ 12 │ 13 │ 14 ┃ 15 │ 16 │ 17 ┃
┠────┼────┼────╂────┼────┼────╂────┼────┼────┨
┃ 18 │ 19 │ 20 ┃ 21 │ 22 │ 23 ┃ 24 │ 25 │ 26 ┃
┣━━━━┿━━━━┿━━━━╋━━━━┿━━━━┿━━━━╋━━━━┿━━━━┿━━━━┫
┃ 27 │ 28 │ 29 ┃ 30 │ 31 │ 32 ┃ 33 │ 34 │ 35 ┃
┠────┼────┼────╂────┼────┼────╂────┼────┼────┨
┃ 36 │ 37 │ 38 ┃ 39 │ 40 │ 41 ┃ 42 │ 43 │ 44 ┃
┠────┼────┼────╂────┼────┼────╂────┼────┼────┨
┃ 45 │ 46 │ 47 ┃ 48 │ 49 │ 50 ┃ 51 │ 52 │ 53 ┃
┣━━━━┿━━━━┿━━━━╋━━━━┿━━━━┿━━━━╋━━━━┿━━━━┿━━━━┫
┃ 54 │ 55 │ 56 ┃ 57 │ 58 │ 59 ┃ 60 │ 61 │ 62 ┃
┠────┼────┼────╂────┼────┼────╂────┼────┼────┨
┃ 63 │ 64 │ 65 ┃ 66 │ 67 │ 68 ┃ 69 │ 70 │ 71 ┃
┠────┼────┼────╂────┼────┼────╂────┼────┼────┨
┃ 72 │ 73 │ 74 ┃ 75 │ 76 │ 77 ┃ 78 │ 79 │ 80 ┃
┗━━━━┷━━━━┷━━━━┻━━━━┷━━━━┷━━━━┻━━━━┷━━━━┷━━━━┛

まず各マスに1~9の数字しか入らないというルールは以下のような多項式の解に限るという条件で表すことができます。

$$F_i = (x_i - 1)(x_i - 2)\cdots(x_i-9)$$

このような条件式を81個全てのマスについて考えます。

次に異なるマスに異なる数字が入るというルールは以下のような多項式の解に限るという条件で表すことができます。

$$G_{ij} = \frac{F_i - F_j}{x_i - x_j}$$

このような条件式を縦列と横列と3×3の小さなブロック全てにおいて $i<j$ を満たす変数のペアで考えます。

以上で数独のルールの表現は終わりました。最後に与えられた問題を一次式に翻訳していきます。例えば与えられた問題の左上のマスが以下のようになっている場合

┏━━━┯━━━┯━━━┳
┃ 3 │ 9 │ 1 ┃
┠───┼───┼───╂
┃ 2 │ 7 │ 6 ┃
┠───┼───┼───╂
┃ 8 │ 5 │   ┃
┣━━━┿━━━┿━━━╋

この部分は

$$\begin{matrix} x_0 - 3, & x_1 - 9, & x_2 - 1 \\ x_9 - 2, & x_{10} - 7, & x_{11} - 6 \\ x_{18} - 8, & x_{19} - 5 & \\ \end{matrix}$$

という連立多項式の解に限るという条件で表すことができます。

以上で数独のルールと具体的な問題を連立方程式で表すことができました。数独の問題は答えが一意に定まるので、これらの多項式の組からグレブナー基底を求めれば連立方程式の解が求まるはずです。

実際にプログラムで数独の問題を解いてみましょう。まずは各マスに1~9の数字しか入らないというルールの実装です。

-- | i番目の変数には1~9のいずれかが入るというルールを表す多項式
sudokuF :: KnownNat n => Int -> Polynomial n
sudokuF i = foldr (\j f -> ((1, var i) `times` f) `minus` ((j, one) `times` f)) (Map.singleton one 1) [1..9]

-- | 各マスには1~9のいずれかが入るというルールを表す多項式
sudokuFAll :: KnownNat n => [Polynomial n]
sudokuFAll = [sudokuF i | i <- [0..9*9-1]]

次に異なるマスに異なる数字が入るというルールの実装です。

-- | i番目のマスとj番目のマスには異なる数字が入るというルールを表す多項式
sudokuG :: KnownNat n => Int -> Int -> Polynomial n
sudokuG i j = head . fst $ (sudokuF i `minus` sudokuF j) `divP` [Map.singleton (var i) 1 `minus` Map.singleton (var j) 1]

-- | 各列中のマスには全て異なる数字が入るというルールを表す多項式
sudokuColumns :: KnownNat n => [Polynomial n]
sudokuColumns = [sudokuG (i+j*9) (i+(j+k)*9) | i <- [0..9-1], j <- [0..9-2], k <- [1..9-j-1]]

-- | 各行中のマスには全て異なる数字が入るというルールを表す多項式
sudokuRows :: KnownNat n => [Polynomial n]
sudokuRows = [sudokuG (i*9+j) (i*9+j+k) | i <- [0..9-1], j <- [0..9-2], k <- [1..9-j-1]]

-- | 各セル中のマスには全て異なる数字が入るというルールを表す多項式
sudokuCells :: KnownNat n => [Polynomial n]
sudokuCells = [sudokuG (i+k*27+l*3) (j+k*27+l*3) | (i, j) <- pairs, k <- [0..2], l <- [0..2]]
  where
  cells = [0, 1, 2, 9, 10, 11, 18, 19, 20]
  pairs = [(i, j) | i <- cells, j <- cells, i < j]

最後に数独の問題として冒頭の例を考えましょう。

┏━━━┯━━━┯━━━┳━━━┯━━━┯━━━┳━━━┯━━━┯━━━┓
┃ 3 │ 9 │ 1 ┃ 4 │ 8 │   ┃ 6 │ 2 │ 7 ┃
┠───┼───┼───╂───┼───┼───╂───┼───┼───┨
┃ 2 │ 7 │ 6 ┃   │ 9 │ 1 ┃ 4 │ 8 │ 5 ┃
┠───┼───┼───╂───┼───┼───╂───┼───┼───┨
┃ 8 │ 5 │   ┃ 2 │ 7 │ 6 ┃ 3 │ 9 │   ┃
┣━━━┿━━━┿━━━╋━━━┿━━━┿━━━╋━━━┿━━━┿━━━┫
┃ 9 │ 1 │   ┃ 8 │ 5 │ 4 ┃ 2 │ 7 │ 6 ┃
┠───┼───┼───╂───┼───┼───╂───┼───┼───┨
┃ 5 │ 4 │ 8 ┃ 7 │ 6 │   ┃ 9 │ 1 │ 3 ┃
┠───┼───┼───╂───┼───┼───╂───┼───┼───┨
┃ 7 │   │ 2 ┃ 9 │ 1 │ 3 ┃   │ 5 │ 4 ┃
┣━━━┿━━━┿━━━╋━━━┿━━━┿━━━╋━━━┿━━━┿━━━┫
┃ 1 │ 3 │ 9 ┃ 5 │   │ 8 ┃   │   │ 2 ┃
┠───┼───┼───╂───┼───┼───╂───┼───┼───┨
┃   │ 2 │ 7 ┃ 1 │ 3 │ 9 ┃ 5 │ 4 │ 8 ┃
┠───┼───┼───╂───┼───┼───╂───┼───┼───┨
┃ 4 │ 8 │ 5 ┃ 6 │ 2 │   ┃ 1 │ 3 │ 9 ┃
┗━━━┷━━━┷━━━┻━━━┷━━━┷━━━┻━━━┷━━━┷━━━┛

これを実装すると以下のようになります。空白マスは 0 で表します。

import Data.Maybe (mapMaybe)

problem :: KnownNat n => [Polynomial n]
problem = mapMaybe interpret $ zip [0..] $ concat [
  [3, 9, 1, 4, 8, 0, 6, 2, 7],
  [2, 7, 6, 0, 9, 1, 4, 8, 5],
  [8, 5, 0, 2, 7, 6, 3, 9, 0],
  [9, 1, 0, 8, 5, 4, 2, 7, 6],
  [5, 4, 8, 7, 6, 0, 9, 1, 3],
  [7, 0, 2, 9, 1, 3, 0, 5, 4],
  [1, 3, 9, 5, 0, 8, 0, 0, 2],
  [0, 2, 7, 1, 3, 9, 5, 4, 8],
  [4, 8, 5, 6, 2, 0, 1, 3, 9]
  ]
  where
  interpret (_, 0) = Nothing
  interpret (i, n) = Just $ Map.singleton (var i) 1 `minus` Map.singleton one n

以上で定義した全ての多項式のグレブナー基底を求めてみましょう。

example5 = do
  let sudokuRules = sudokuFAll <> sudokuColumns <> sudokuRows <> sudokuCells :: [Polynomial 81]
      answers = reduce $ minimize $ buchberger (problem <> sudokuRules)
   in mapM_ (putStrLn . showPolyXN) answers

x_0 - 3
x_1 - 9
x_2 - 1
x_3 - 4
x_4 - 8
x_6 - 6
x_7 - 2
x_8 - 7
x_9 - 2
x_10 - 7
x_11 - 6
x_13 - 9
x_14 - 1
x_15 - 4
x_16 - 8
x_17 - 5
x_18 - 8
x_19 - 5
x_21 - 2
x_22 - 7
x_23 - 6
x_24 - 3
x_25 - 9
x_27 - 9
x_28 - 1
x_30 - 8
x_31 - 5
x_32 - 4
x_33 - 2
x_34 - 7
x_35 - 6
x_36 - 5
x_37 - 4
x_38 - 8
x_39 - 7
x_40 - 6
x_42 - 9
x_43 - 1
x_44 - 3
x_45 - 7
x_47 - 2
x_48 - 9
x_49 - 1
x_50 - 3
x_52 - 5
x_53 - 4
x_54 - 1
x_55 - 3
x_56 - 9
x_57 - 5
x_59 - 8
x_62 - 2
x_64 - 2
x_65 - 7
x_66 - 1
x_67 - 3
x_68 - 9
x_69 - 5
x_70 - 4
x_71 - 8
x_72 - 4
x_73 - 8
x_74 - 5
x_75 - 6
x_76 - 2
x_78 - 1
x_79 - 3
x_80 - 9
x_63 - 6
x_46 - 6
x_12 - 3
x_58 - 4
x_61 - 6
x_26 - 1
x_5 - 5
x_20 - 4
x_29 - 3
x_41 - 2
x_51 - 8
x_60 - 7
x_77 - 7

2時間ちょっとかけて結果が出力されました…！空白マスに対応する変数だけを抜き出すと以下のようになっており、ちゃんと問題を解けていることがわかります👏

x_63 - 6
x_46 - 6
x_12 - 3
x_58 - 4
x_61 - 6
x_26 - 1
x_5 - 5
x_20 - 4
x_29 - 3
x_41 - 2
x_51 - 8
x_60 - 7
x_77 - 7

まとめ

本稿では数独を多項式の言葉に翻訳してグレブナー基底を解く方法とその実装を紹介しました。多項式の言葉に変換できるパズルは数独だけでなく、ラテン方格 をベースにした様々なパズルに応用が可能です（例えばカックロとかサムナンプレとか）。グレブナー基底を使ってパズルを解く話の文献としては以下を辿るのがオススメです。

• 割り算アルゴリズムとグレブナー基底 〜連立方程式でナンプレを解こう！〜（ロング版）
• グレブナー基底計算プログラムはどの程度のパズルを解くことができるか？
• Sudokus and Gröbner Bases: Not Only a Divertimento

また本稿ではグレブナー基底にまつわる事実を何も証明せずに用いましたが理論的な側面が気になる方はこちらの本が分かりやすかったのでオススメです。

• グレブナー基底とその応用

本稿のプログラムを実装するにあたっては上記の本と以下の文献を参考にしています。

• Computing Gröbner Bases in Python with Buchberger’s Algorithm

グレブナー基底という代数幾何学的な手法が数独という身近なパズルを解くのに使えるのはとても面白い考え方ですが、こちらの文献にもあるようにこの方法は問題の難易度が上がると答えが出るまでに時間がかかってしまい、高速な方法というわけではありません。また本稿での実装は理論から素直に実装したものであり効率が良いものでもありません。グレブナー基底を高速に求める方法はF5アルゴリズムを始め様々な工夫が研究されています。また数独を解くにあたってはブーリアングレブナー基底というブール多項式環上のグレブナー基底を用いる方法もあるようです。

挿絵イラスト: Loose Drawing