ZDDを真面目に解説する話【Daily Akari ソルバー実装編 1/3】

はじめに

前回の記事↓ で、Daily Akari というパズルの ZDD による解法をサラッと紹介しました。

kntychance.hatenablog.jp

しかし、どうにも 「ZDD とはなんぞや？」「よくわからんが難しそう」といった風に感じている人がなんとなく多そうな気がしており、これは由々しき事態だなと。

ZDD はすんごくシンプルなアイデアで色々な問題を解けてしまう激面白データ構造なので、何としてでも布教したいと思いました。

なので、アルゴリズム学徒だけでなくプログラミングにそこまで明るくない非応用数学系の人々にも届くように一から丁寧に解説していきます。ゴールがあった方がいいと思うので、Daily Akari ソルバーの作成を目標にしてしまいましょう。

目指せ、一家に一台 ZDD！

なお、全体の実装例は ↓ に置いてあります。

（本記事シリーズでは Java の方を実装していきます。）

github.com

この記事でやること

• ZDD の概要の解説
• ZDD のベースとなるクラスの解説・実装（Java JDK 24）
• ZDD の布教！！！！！！

この記事でやらないこと

• ZDD の演算の話（第2回でガッツリやる予定。ZDD の激面白ポイントはここ。）
• Daily Akari の話（第3回で入る予定。気長にやりましょう。）

1. ZDD の概要

1.1. ZDD とは

Zero-supressed Binary Dicision Diagram の略です。日本語では「ゼロサプレス（抑制？）型二分決定グラフ」になります。あまり日本語で呼ばれているところは見ないですが、提案されたのは日本の教授、湊真一先生です[1]。2025年現在は京大で教鞭をとられています。

与えられた集合 𝐸 :={𝑒0,…,𝑒𝑛−1} に対して、ZDD は部分集合族 P ⊆2𝐸 を表現します。以下では 𝐸 を 台集合、𝑒0,…,𝑒𝑛−1 を アイテム、各部分集合 𝑃 ∈P を 組合せ集合 と呼びます。

例えば、𝐸 :={𝑎,𝑏,𝑐} に対し P :={{𝑎,𝑏},{𝑎,𝑐},{𝑐}} を表す ZDD は【図 1.1-1】になります。

【図 1.1-1】　𝐸 :={𝑎,𝑏,𝑐} に対し P :={{𝑎,𝑏},{𝑎,𝑐},{𝑐}} を表す ZDD

【図 1.1-1】において、赤い辺を 0-枝、青い辺を 1-枝と呼びます。また、一番上の 𝑎 のラベルがついた節点を 根、一番下の ⊥ のラベルがついた節点を 0-終端節点、⊤ のラベルがついた節点を 1-終端節点 と呼びます。

終端節点を除く各節点からは、それぞれ0-枝と1-枝が1つずつ出ています。これが、「二分決定グラフ」と呼ばれる所以です。

ZDD において、根から1-終端節点までのパスが各組合せ集合に対応します。

すなわち、各パスについて節点から出る際に

• 1-枝を通るなら、その節点のアイテムを組合せ集合に追加する。
• 0-枝を通るなら、追加しない。

とすることで組合せ集合と対応付けます。（【図 1.1-2】）

【図 1.1-2】　根から1-終端節点までのパスが組合せ集合に対応

特に、この対応付けになぞらえると、0-終端節点のみからなる ZDD は部分集合族 {} を、 1-終端節点のみからなる ZDD は部分集合族 {∅} をそれぞれ表していることに注意してください。

以上が、ZDD の"読み方"になります。

この説明だけ聞くと「どこからそのグラフは降ってきたんですか？」という疑問が出ると思うので、次はその解説をします。

1.2. ZDD の基礎的な構成方法

前節に引き続き、台集合 𝐸 :={𝑎,𝑏,𝑐} に対し部分集合族 P :={{𝑎,𝑏},{𝑎,𝑐},{𝑐}} を表す ZDD を例にしましょう。

まず、【図 1.2-1】に示す完全二分決定木を考えます。

【図 1.2-1】　完全二分決定木

この木において、根から各終端節点までのパスは 𝐸 の全ての部分集合と対応付けられます。特に、終端節点の ⊥ と ⊤ を切り替えることで任意の部分集合族（アイテム数が 𝑛 なら、22𝑛 通り）を表現できることに注意してください。

さて、このままだと節点数が 2𝑛+1 −1 となり、𝑛 が大きいと困ってしまいます。実用上は表現したい部分集合族がべき集合に対して十分疎であることが多いので、もっと効率的に表現したいところ。

そこで、次の2つの圧縮規則（【図 1.2-2】）を導入します。

1. 削除規則：
   1-枝の先が0-終端節点であるような節点を省略する。
2. 共有規則：
   等価な節点（ラベル、0-枝の行先、1-枝の行先が同一である節点）を纏める。特に、終端節点も等価な節点とみなす。

【図 1.2-2】　削除規則と共有規則

この2つの圧縮規則を、【図 1.2-1】で示した完全二分決定木に対し（任意の手順で）繰り返し適用することで、最終的に【図 1.1-1】に示す ZDD が得られます*1。また、この圧縮規則が適用できない状態のことを既約であると言います*2。

というわけで、ZDD の構成法が得られました。

しかし、この構成法には大きな問題があります。

というのも、最初に完全二分決定木を作ってしまうと、その時点で計算量が爆発してしまいます。

目標としている Daily Akari ソルバーにおいては、台集合をマス全体として、解であるような灯りの置き方（=マスの部分集合）全体から成る ZDD を構成します。グリッドサイズは最大 24 ×14 =336 なので*3、完全二分決定木の節点数は 2336+1 −1⁢( >10100) になってしまいます。これは宇宙の原子の数より多いらしいので[2]、この構成法は物理的に不可能です。

というわけで、本節で解説した完全二分決定木を経由する構成法は、現実的ではありません。ただし、2種類の圧縮規則とアイテムの添字順（=完全二分決定木上でアイテムをどの順に置くか）によって ZDD の形状が定まっているということは、知識として知ってくおくとよいでしょう。

次項では、「じゃあ結局 ZDD ってどう作るの？」という疑問に答えていきます。

1.3. ZDD の効率的な構成方法

1.3.1. 部分集合族の分解

ZDD の各接点は、それ自体を根とする ZDD と捉えられ、対応する部分集合族が存在します。そこで、根の0-枝の先と1-枝の先がそれぞれどのような部分集合族を表しているかを考えます。

以下、台集合を 𝐸 :={𝑒0,…,𝑒𝑛−1} とします。

また、便利な記号として、部分集合族 P ∈2𝐸 とアイテム 𝑒𝑖 ∈𝐸 に対する二項演算子 ◃ を次のように定義します。

P◃𝑒𝑖:={𝑃∪{𝑒𝑖}∣𝑃∈P}(1)

すると、任意の非空な組合せ集合を少なくとも1つもつ部分集合族 P ∈2𝐸 について、

𝑒𝑖:=min⁡(⋃𝑃∈P𝑃)P0:={𝑃∣𝑒𝑖∉𝑃∈P}P1:={𝑃∖{𝑒𝑖}∣𝑒𝑖∈𝑃∈P}(2)

（アイテムの大小は添字で比較するものとします。）とおいたとき、次式が成り立ちます。

P=P0∪(P1◃𝑒𝑖)(3)

この 式(3) は、【図 1.3.1-1】に示す ZDD の抽象図をそのまま表現しています。

【図 1.3.1-1】　式(3) と ZDD

ふわふわとお絵描きしていたものが、かっちりと扱えそうになってきましたね。

1.3.2. 圧縮規則の再帰的な適用

さて、勘のいい読者は気づかれているかもしれませんが、式(3) は ZDD の再帰的な構成方法を示唆しています。

すなわち、P0 と P1 に対応した ZDD が得られていると仮定すると、そこに 𝑒𝑖 をラベルに持つ節点を追加し、枝を ｴｲﾔｯ!! とそれぞれの根に突き刺すことで、P の ZDD が完成しそうです。

この枝を ｴｲﾔｯ!! する操作に圧縮規則を適用したものが、次の手続き '_' allowed only in math mode になります。

'_' allowed only in math mode

• P1 ={} なら、P0 の根を返す。

• そうでないなら、以下の節点を作成して返す。ただし、等価な節点を既に作成済みならそれを返す。
  • ラベル： 𝑒𝑖
  • 0-枝の先： P0 の根
  • 1-枝の先： P1 の根

この手続きを用いることで、与えられた部分集合族 P に対し ZDD を構築しその根を返す手続き build⁡(P) は、下記のように再帰的に記述できます。

build⁡(P)

• P ={} なら、0-終端節点を返す。
• P ={∅} なら、1-終端節点を返す。
• 上記以外なら、式(2) により 𝑒𝑖,P0,P1 を得て、'_' allowed only in math mode を返す。

というわけで、ZDD の良い構成方法が得られました。

しかし、この build という手続きは、部分集合族 P を表現する手段が存在する前提で記述されており、結局それをどう実装するのか？という点で言えばあまり意味を成しません*4。大事なのは、このように再帰的に記述できることを理解しておくことです。

（ちなみに、最初からこの build という手続きで ZDD を定義してしまっても問題ありません。数学的にはその方が簡潔かもしれませんね。）

「え、じゃあ結局 ZDD ってどうやって使うの？？？」

とツッコミが入るかもしれませんが、解説する分量の都合上、次回までお待ちください。実際のところそのあたりが ZDD で一番面白い話で「おわりに」の章でも軽く触れますが、本記事ではとにかく ZDD のベースとなる部分の実装を目指します。

1.3.3. ハッシュ関数

1.3.2項 にて '_' allowed only in math mode という手続きを紹介しましたが、これを実装するうえでは等価節点をどのように保存・検索するのかという点にも着目しなくてはなりません。いわゆるデータ構造ってやつです。（ZDD も部分集合族を表現するためのデータ構造と言えます。）

計算機で扱える基本的なデータ構造として、配列があります。これは、複数のデータを1列に並べたもので、任意の非負整数 𝑖 に対し 𝑖 番目*5のデータの読み書きを高速に処理できます。なお、計算機のメモリにも限りがあるため、配列の長さとして有限な非負整数 𝑚 を決めておく必要があります。

そして、配列を流用するデータ構造として、ハッシュテーブルというものがあります。

ja.wikipedia.org

ハッシュテーブルは、保存したいデータに対し予め用意した"いい感じの関数"を使って 0  以上 𝑚 未満の非負整数 𝑖 に変換し、配列の 𝑖 番目に保存すると取り決めたものです。この"いい感じの関数"をハッシュ関数と呼びます。また、ハッシュ関数で得られた整数値のことをハッシュ値と呼びます。

（用意した関数が単射でない場合、異なるデータに対して同じ番地が充てられてしまう場合があり、これを衝突と呼びます。衝突した場合にどう保存するかは別途取り決めておく必要があるのと、そもそも衝突があまり起きないように関数を定義することが重要です。）

ハッシュテーブルに ZDD の節点を保存できれば、検索が高速にできてうれしいです。なので、サクッとハッシュ関数を作ってしまいましょう*6。

ハッシュ関数を作る上での便利な道具として、Xorshift と RollingHash という 2 つの関数を使います。具体的な中身はリンク先の解説が易しいのでそちらを参照ください。

• Xorshift⁡(𝑥);ℕ>0 ↦ℕ>0 :
  正整数値 𝑥 を何かしらの正整数値に飛ばす疑似乱数発生器
  blog.visvirial.com
  
  
• RollingHash⁡(𝑆);ℕ𝑘>0 ↦ℤ/𝑚⁢ℤ:
  任意長の正整数列 𝑆 を何かしらの法 𝑚 の整数値に飛ばすハッシュ関数
  maspypy.com
  
  

0-終端節点のハッシュ値を Xorshift⁡(𝑛 +1)、
1-終端節点のハッシュ値を Xorshift⁡(𝑛 +2) とします。

終端節点でない節点 P について、ラベルのアイテムを 𝑒𝑖、0-枝の先を P0、1-枝の先を P1 とし、P0 と P1 のハッシュ値は既に求まっていると仮定します。

この P のハッシュ値を次のように定義します：

hash⁡(P):=RollingHash(Xorshift⁡(𝑖+1),Xorshift⁡(hash⁡(P0)+1),Xorshift⁡(hash⁡(P1)+1))(4)

要するに、アイテムと枝の先にある2節点のハッシュ値から自身のハッシュ値を定めているわけですね。これでハッシュ関数の完成です。

以上で ZDD のベースとなる考え方の解説は終了です。

次章からは実装に入っていきます。

2. 実装

zddForDailyAkari パッケージに、続く3つのクラスを実装します。

なお、3つ目の ZDD_Visualizer はなくてもDailyAkari は解けるので、飛ばしても問題ありません。

2.1. RollingHash クラス

import java.util.List;
import java.util.Random;

public class RollingHash {

    public int m; // 法
    private int r; // 基数

    public RollingHash(int m) {
        this.m = m;
        Random random = new Random();
        this.r = random.nextInt(2, m);
    }

    public int hash(List<Object> list) {
        long l = r;
        long h = 0;
        for (Object e : list) {
            h += xorshift(e.hashCode() + 1) * l;
            h %= m;
            l = l * r % m;
        }
        return (int) h;
    }

    public int xorshift(int x) {
        x ^= x << 13;
        x ^= x >> 17;
        x ^= x << 15;
        return Integer.remainderUnsigned(x, m);
    }
}

式(4) に基づいたハッシュ関数を実装したクラスです。 RollingHash という名前ですが、中に Xorshift も内蔵しています。

コンストラクタは 法 m を引数に取り、コンストラクタ内で基数 r を乱択しています。これは Rolling Hashについて（survey + 研究） | maspyのHP にて解説された内容を実装したものになります。

2.2. ZDD<T>、ZDD<T>.ZDD_Node クラス

import java.util.ArrayList;
import java.util.Arrays;
import java.util.HashMap;

public class ZDD<T> {
    
    public int n;
    public ArrayList<T> items;
    public HashMap<T, Integer> idx;
    public ZDD_Node zero_terminal;
    public ZDD_Node one_terminal;
    
    public RollingHash rollingHash = new RollingHash(1401055423); // 1401055423 は素数。
    public HashMap<ZDD_Node, ZDD_Node> node_table = new HashMap<>();
    
    public ZDD(ArrayList<T> arrayList) {
        this.n = arrayList.size();
        this.items = (ArrayList<T>)arrayList.clone();
        this.idx = new HashMap<>();
        this.zero_terminal = new ZDD_Node(null, null, null);
        this.one_terminal = new ZDD_Node(null, null, null);
        zero_terminal.hash = rollingHash.xorshift(n+1);
        one_terminal.hash = rollingHash.xorshift(n+2);
        for(int i = 0; i < n; ++i) idx.put(items.get(i), i);
    }
    
    public ZDD_Node get_node(T top, ZDD_Node zero, ZDD_Node one) {
        if (one == zero_terminal) return zero;
        ZDD_Node node = new ZDD_Node(top, zero, one);
        if (node_table.containsKey(node)) {
            node = node_table.get(node);
        } else {
            node_table.put(node, node);
        }
        return node;
    }



    public class ZDD_Node {
        
        public T top;
        public ZDD_Node zero;
        public ZDD_Node one;
        public int hash;

        public ZDD_Node(T top, ZDD_Node zero, ZDD_Node one) {
            this.top = top;
            this.zero = zero;
            this.one = one;
            if (top != null) {
                this.hash = rollingHash.hash(Arrays.asList(idx.get(top), zero, one));
            }
        }
    
        @Override
        public int hashCode() {
            return this.hash;
        }

        @Override
        public boolean equals(Object other) {
            if (this == other) return true;
            if (other instanceof ZDD<?>.ZDD_Node o) {
                if (this.top == null || o.top == null) return this == o;
                return this.top.equals(o.top)
                && this.zero == o.zero
                && this.one == o.one;
            } else {
                return false;
            }
        }
    }
}

ZDD クラスと、その内部クラスとして節点を表す ZDD_Node クラスを実装しています。これは、台集合や節点のハッシュテーブルを ZDD クラスに管理してもらいたいという思想があるからです。 また、1.3.2項 ので解説した '_' allowed only in math mode も ZDD クラスのメソッドとします。

T はアイテムのクラスを表します。ZDD<T> のコンストラクタは台集合として ArrayList<T> を引数に取り、ArrayList での順に沿ってアイテムの添字を決定します。

決定した添字は HashMap<T, Integer> idx に保存しています。

ZDD_Node クラスでは hashCode() と eauals() をオーバーライドしています。

注意点として、equals() 内の

this.zero == o.zero && this.one == o.one

としている箇所を 

this.zero.equals(o.zero) && this.one.equals(o.one)

としてしまうと、equals() が終端節点まで再帰的に呼び出されてしまい計算量が爆発してしまいます。'_' allowed only in math mode で適用される共有規則を踏まえると、「==」での比較（同一のインスタンスかどうかの比較）で十分です。

同様に、終端節点まで hashCode() が再帰的に呼び出されることを回避するため、コンストラクタ内でハッシュ値を予め計算し保存しておくようにしています。

2.3. ZDD_Visualizer<T> クラス

外部ライブラリである jung-2.0.1 を利用するため、導入方法をまず紹介します。

まず、下記リンクからライブラリの圧縮ファイルをダウンロードしてください。

sourceforge.net

ダウンロードした圧縮ファイルを適当なディレクトリに解凍しておき（.jar ファイルがいっぱい入っています。）、自身の java 実行環境に合わせてクラスパスを通せばよいです。

ここでは一応、筆者の環境（vscode, Extension Pack For Java）でのパスの通し方を紹介しておきます。

1. 画面左下の ☕java Ready と表示された箇所をクリック

【図 2.3-1】　Java Ready をクリック

2. Open Project Settings を選択
   

   

   【図 2.3-2】　Open Project Settings

   
   
   
3. Libraries タブを開き、解凍した .jar ファイルを追加する。
   

   

   【図 2.3-4】　Libraries タブを開き、解凍した .jar ファイルを追加する。

   
   
   
4. Apply Settings をクリック
   

   

   【図 2.3-4】　Apply Settings をクリック

   以上で完了です。

導入出来たら、実装例をそのままコピペしてもらえばよいです。ZDD の本筋とは関係ないので解説は省略します。

参考にした記事 ↓（ライブラリの設計思想がなんとなくわかります）

JUNG 2: Javaライブラリ: 篠田孝祐

あとは API とにらめっこしました。

jung2 2.0 API

import javax.swing.JFrame;

import org.apache.commons.collections15.Transformer;

import java.awt.Color;
import java.awt.Dimension;
import java.awt.Font;
import java.awt.Paint;
import java.awt.Rectangle;
import java.awt.Shape;
import java.awt.geom.Arc2D;

import edu.uci.ics.jung.algorithms.layout.FRLayout;
import edu.uci.ics.jung.algorithms.layout.Layout;
import edu.uci.ics.jung.graph.DirectedSparseMultigraph;
import edu.uci.ics.jung.visualization.BasicVisualizationServer;
import edu.uci.ics.jung.visualization.decorators.ToStringLabeller;
import edu.uci.ics.jung.visualization.renderers.Renderer.VertexLabel.Position;

import java.util.ArrayDeque;
import java.util.HashSet;

public class ZDD_Visualizer<T> {

    public class Vertex{
        
        ZDD<T>.ZDD_Node node;
        String label;


        public Vertex(ZDD<T>.ZDD_Node node) {
            this.node = node;
            this.label = setLabel();
        }

        @Override
        public String toString() {
            return label;
        }

        @Override
        public boolean equals(Object other) {
            if (this == other) return true;
            if (other instanceof ZDD_Visualizer<?>.Vertex o) {
                return this.node == o.node;
            } else {
                return false;
            }
        }

        @Override
        public int hashCode() {
            return this.node.hashCode();
        }
        
        private String setLabel() {
            if (this.node == zdd.zero_terminal) return "⊥";
            if (this.node == zdd.one_terminal) return "T";
            return this.node.top.toString();
        }
    }

    public class Edge {

        Vertex vFrom;
        Vertex vTo;
        int label; // 0 or 1

        public Edge(Vertex vFrom, Vertex vTo, int label) {
            this.vFrom = vFrom;
            this.vTo = vTo;
            this.label = label;
        }
    }

    ZDD<T> zdd;

    public ZDD_Visualizer(ZDD<T> zdd) {
        this.zdd = zdd;
    }

    public void visualize(ZDD<T>.ZDD_Node root) {

        DirectedSparseMultigraph<Vertex, Edge> graph = new DirectedSparseMultigraph<>();
        Vertex vRoot = new Vertex(root);
        graph.addVertex(vRoot);

        HashSet<Vertex> foundNodes = new HashSet<>();
        ArrayDeque<Vertex> queue = new ArrayDeque<>();
        queue.add(vRoot);
        foundNodes.add(vRoot);

        while (!queue.isEmpty()) {
            Vertex vNode = queue.poll();
            if (vNode.node == zdd.zero_terminal || vNode.node == zdd.one_terminal) continue;

            Vertex vZero = new Vertex(vNode.node.zero);
            Vertex vOne = new Vertex(vNode.node.one);
            if (!foundNodes.contains(vZero)) {
                foundNodes.add(vZero);
                graph.addVertex(vZero);
                queue.add(vZero);
            }
            if (!foundNodes.contains(vOne)) {
                foundNodes.add(vOne);
                graph.addVertex(vOne);
                queue.add(vOne);
            }
            if (!graph.addEdge(new Edge(vNode, vZero, 0), vNode, vZero)) System.err.println("can't create");;
            if (!graph.addEdge(new Edge(vNode, vOne, 1), vNode, vOne)) System.err.println("can't create");;
        }

        Layout layout = new FRLayout(graph);
        BasicVisualizationServer visualizationServer = new BasicVisualizationServer(layout, new Dimension(800, 600));

        visualizationServer.getRenderContext().setVertexLabelTransformer(new ToStringLabeller<>()); // toString() で得られる文字列を節点ラベルとして表示する
        visualizationServer.getRenderer().getVertexLabelRenderer().setPosition(Position.CNTR); // 節点ラベルを節点内に表示する

        // 節点ラベルのフォントの設定
        Transformer<Vertex, Font> vertexLabelFontTransformer = new Transformer<>() {
            public Font transform(Vertex vertex) {
                return new Font(visualizationServer.getFont().getFontName(), Font.BOLD, 24);
            }
        };
        visualizationServer.getRenderContext().setVertexFontTransformer(vertexLabelFontTransformer);

        // 節点の色の設定
        Transformer<Vertex, Paint> vertexColorTransformer = new Transformer<>() {
            public Paint transform(Vertex vertex) {
                if (vertex.node == zdd.zero_terminal) return new Color(0xff6666);
                if (vertex.node == zdd.one_terminal) return new Color(0x6666ff);
                return Color.WHITE;
            }
        };
        visualizationServer.getRenderContext().setVertexFillPaintTransformer(vertexColorTransformer);
        
        // 節点の形状の設定
        Transformer<Vertex, Shape> vertexShapeTransformer = new Transformer<>() {
            public Shape transform(Vertex vertex) {
                if (vertex.node == zdd.zero_terminal || vertex.node == zdd.one_terminal) return new Rectangle(-5,-5,30,40);
                return new Arc2D.Double(-7, -7, 30, 30, 0, 360, Arc2D.OPEN);
            }
        };
        visualizationServer.getRenderContext().setVertexShapeTransformer(vertexShapeTransformer);

        // 枝・矢頭の色の設定
        Transformer<Edge, Paint> edgeColorTransformer = new Transformer<>() {
            public Paint transform(Edge edge) {
                return edge.label == 0 ? Color.RED : Color.BLUE;
            }
        };
        visualizationServer.getRenderContext().setEdgeDrawPaintTransformer(edgeColorTransformer); // 枝の色
        visualizationServer.getRenderContext().setArrowDrawPaintTransformer(edgeColorTransformer); // 矢頭の枠の色
        visualizationServer.getRenderContext().setArrowFillPaintTransformer(edgeColorTransformer); // 矢頭の塗りつぶし色

        JFrame frame = new JFrame("ZDD_visualizer");
        frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        frame.getContentPane().add(visualizationServer);
        frame.pack();
        frame.setVisible(true);
    }
}

2.4. テスト

import java.util.ArrayList;
import java.util.Arrays;

public class JungTest {
    public static void main(String[] args) {
        String[] items = {"a", "b", "c"};
        ArrayList<String> arrayList = new ArrayList<>(Arrays.asList(items));
        ZDD<String> zdd = new ZDD<>(arrayList);
        ZDD_Visualizer<String> visualizer = new ZDD_Visualizer<>(zdd);

        ZDD<String>.ZDD_Node n1 = zdd.get_node("c", zdd.zero_terminal, zdd.one_terminal);
        ZDD<String>.ZDD_Node n2 = zdd.get_node("b", n1, zdd.one_terminal);
        ZDD<String>.ZDD_Node n3 = zdd.get_node("a", n1, n2);
        visualizer.visualize(n3);
    }
}

実行結果：

【図 2.4-1】　ZDD_Visualizer による図示

おわりに+次回予告

お疲れ様でした。

今回は ZDD のベースとなる考え方として、再帰的な構築法を解説しました。

しかし、実際のところまだ ZDD の嬉しさについては全く解説していません。

ZDD のハイパー面白ポイントは部分集合族の“演算”が ZDD 上でも扱えるという点にあり、これにより様々な組合せ問題、すなわち条件を満たす組合せ集合全体からなる部分集合族を求める問題が簡単に解けてしまうのです。

そして、まさに今回目標としている Daily Akariも組合せ問題の一種で、ZDD の演算だけで解けちゃうわけです。凄い。

次回はそんな ZDD の演算を解説していきます。お楽しみに。

参考文献

1. S. Minato: "Zero-Suppressed BDDs and Their Applications", International Journal on Software Tools for Technology Transfer, Vol. 3, No. 2, pp. 156-170, Springer, May 2001.

2. 宇宙の原子の数は10^80程度 - よくある誤解と正しい知識

3. Google Chromeが採用した、擬似乱数生成アルゴリズム「xorshift」の数理 – びりあるの研究ノート

4. Rolling Hashについて（survey + 研究） | maspyのHP

5. JUNG 2: Javaライブラリ: 篠田孝祐

6. jung2 2.0 API

*1:任意の完全二分木について、圧縮規則を繰り返し適用することで、対応した ZDD が手順に依らず一意に得られます。このことは、後ほど登場する 式(2), 式(3) をもとに帰納法を用いることで示せます。

*2:既約になっていない二分決定グラフも広く ZDD と呼ぶ場合もありますが、ここでは既約なものを指して ZDD と呼ぶことにします。

*3:筆者調べ。もっと大きい日もあるかもしれません。

*4:無論、P をプログラム上で表現する別の手段があれば、この手続きにより ZDD を構築できます。

*5:配列は、0番目から数え始めるのが一般的です。

*6:ここで紹介するハッシュ関数は筆者が実装時に作成した一例であり、人々がどのように実装しているのかは知りません。また、この関数がハッシュ関数として優れているかどうかの数学的な考察も一切していません。ごめんなさい。

Daily Akari を ZDD で解く話

はじめに

Daily Akari というパズルがあります。

dailyakari.com

ルールを見た瞬間、ZDD で解けるなぁと思ったので、解きました。

ZDD が何かについては、過去記事を参照してください。(全3部)

過去記事↓

ZDDと8-Queen problemの話【卒研お勉強枠①】 - 競プロ勉強メモ

ZDDによるN-Queen problemの解法を効率化する話（卒研お勉強枠②） - 競プロ勉強メモ

ZDDでマインスイーパーを解く話（卒研お勉強枠➂（最終回）） - 競プロ勉強メモ

Daily Akari のルール

• グリッド上に「空きマス」、「ブロック」、「数字付きブロック（以下、書かれた数字が 𝑥 のとき、𝑥-ブロック と呼ぶ）」が配置された盤面が与えられる。
• 任意の空きマスには「灯り」を配置できる。灯りは上下左右方向のブロックに遮られていない空きマスを全て照らす。
• 以下の条件を満たすように灯りを配置出来ればクリア。
  1. 全ての空きマスを照らす。
  2. 全ての𝑥-ブロックについて、隣接する上下左右の空きマスに灯りをちょうど 𝑥 個配置する。
  3. 灯りは自身以外に照らされてはならない。

Daily Akari プレイ画面

クリア例

ZDDで扱える形への帰着

𝑥-ブロックは卒研で扱ったマインスイーパーと全く同じなので処理できるとして、他の2ルールについても次のように「空きマスの部分集合に対する個数制限」の形で表現できます。（下記は、以下 条件A, B, C と呼びます。）

1. 全ての空きマスについて、上下左右方向のブロックに遮られていない空きマスに灯りを少なくとも 1 つ配置する。
2. 全てのブロック（グリッド左端の壁を含む）について、右方向にあるブロックに遮られていない空きマスに灯りを高々 1 つ配置する。
3. 全てのブロック（グリッド上端の壁を含む）について、下方向にあるブロックに遮られていない空きマスに灯りを高々 1 つ配置する。

まあ実質マインスイーパーですね。（個人の感想）

空きマス全体を台集合として、べき集合を表す ZDD （0-枝と1-枝を順に繋げていくだけ）を構築しておきます。

あとは、ZDDでマインスイーパーを解く話（卒研お勉強枠➂（最終回）） - 競プロ勉強メモ で紹介した hint 演算や、それをヒント数字「以上」・「以下」に対応させた演算 '_' allowed only in math mode、'_' allowed only in math mode を適用しまくれば解けます。やったね。

ソルバーの実装例

githubに置いときます。

ちなみに

ルール説明時に貼っていた問題は唯一解らしいです。

※ 頂点の数字を9で割ったときの商と余りがグリッドの座標に対応。青辺を通る時そのマスに灯りを配置する。

ルール説明時の問題に対してソルバーが構築したZDD

 

高速化について [追記：2025/04/12]

本解法は、盤面が大きくなるほど内部で構築される ZDD の節点数が指数的に増え、計算時間もメモリも爆増します。

しかし、ZDD の特徴として、演算の順序だったり組合せだったりを適切に決めることで同じ問題でもいい感じに節点数を抑えられることが知られています*1。そういったチューニングについては、思いつき次第追記していこうかなと思います。（データ集計とかは面倒なので、体感でよくなってそうなのを書きます。）

1. 条件A について、「少なくとも 1 つ」ではなく「1 個以上 2 個以下」に制限し、'_' allowed only in math mode 演算で処理する。[追記: 2025/04/12]

*1:卒研解説記事の第2部 でも触れています。

鳩ノ巣原理の良問を懐かしむ話【ABC227-D】

はじめに

鳩が 𝑁 匹と鳩の巣が 𝑀 個あり、各鳩はいずれか一つの巣に入っているとします。このとき、 𝐍>𝐌 ならば、必ず 𝟐 匹以上の鳩が入った巣が存在します。あたりまえ体操ですね。この主張を、鳩ノ巣原理と呼びます。

この鳩ノ巣原理、競プロでたまーに出てきては猛威を振るっていくんですが、本記事では昔出てきてちょっと感心した良問を紹介します。

問題概要

扱う問題はこれ ↓ です

atcoder.jp

問題概要

• 𝑁 個の部署があり、各部署 𝑖 (1 ≤𝑖 ≤𝑁) には 𝐴𝑖 人の社員が所属している。
• 𝐾 人組のプロジェクトを作る。ここで、
  • 一つのプロジェクトに同じ部署に所属する社員が 2 人以上参加してはならない。
  • 一人が 2 つ以上のプロジェクトに参加してはならない。
• 最大でいくつのプロジェクトを作ることができるか。

• 実行時間制限: 2sec
• 1 ≤𝐾 ≤𝑁 ≤2 ×105
• 1 ≤𝐴𝑖 ≤1012 (1 ≤𝑖 ≤𝑁)

（2 ×1017人の社員を持つこの社、ヤバすぎる）

やりがちな嘘解法

単純な解法として、貪欲法が浮かびます。

例えば 𝐾 =3として、（必要ならソートをして）部署が人数の昇順に並んでいるとして、

• 部署 1 ~ 3 から 1 人ずつ取ってプロジェクトを作れるだけ作る。
• 部署 2 ~ 4 から 1 人ずつ取ってプロジェクトを作れるだけ作る。
• 部署 3 ~ 5 から 1 人ずつ取ってプロジェクトを作れるだけ作る。
• （以下繰り返し）

しかし、この解法は次のようなケースで最適解となりません。

• 𝐾 =3
• 𝐴 =[1,1,2,2]

（貪欲法だと 1 個しかプロジェクトが作れないですね。）

正しい解法

決め打ち二分探索

いつものです。

𝑃⁡(𝑥) :=𝑥 個のプロジェクトを作れる という判定問題を考え、𝑃⁡(𝑥) =𝑇⁢𝑟⁢𝑢⁢𝑒 となる最大の 𝑥 を二分探索します。

判定問題の解法

実は、任意の 𝑥 に対し次が成り立ちます：

𝑃⁡(𝑥)⇔𝑁∑𝑖=1min⁡(𝐴𝑖,𝑥)≥𝐾⁢𝑥

∵ 各部署について 𝑥 人を超える分は無視してよいことに注意して（鳩ノ巣原理！）、各部署 𝑖 (1 ≤𝑖 ≤𝑁) の社員数を min⁡(𝐴𝑖,𝑥) とおきなおしたとき、

• ( ⇒): 対偶を考えると、社員総数が足りていないので 𝑃⁡(𝑥) =𝐹⁡𝑎⁢𝑙⁢𝑠⁢𝑒 は明らか。
• ( ⇐): プロジェクトをそれぞれ 𝑝0,𝑝1,…,𝑝𝑥−1 とする。部署順に社員全員を一列に並べて（同一部署の社員は自由な順で並べてok）前から順に 0,1,…,(𝐾⁢𝑥 −1),… と番号をつけたとき、各社員 𝑗 (0 ≤𝑗 <𝐾⁢𝑥) を 𝑝𝑗%𝑥 に参加させればよい*1。

提出コード

https://atcoder.jp/contests/abc227/submissions/64164861

おわりに

いかがでしたか？

鳩ノ巣原理には、しれっと出てきて「それはそう…」と落ち込ませる破壊力があります。つらい。

余談

鳩ノ巣原理を一般化したような立ち位置にある、Ramsey理論というものがあります。気になる人は組合せ論をやりましょう（布教）。

www.maruzen-publishing.co.jp

*1:この構成法は、問題概要の定義 2 つを両方とも満たしています。確認してみましょう。

【ABC388-G 解説？】並列二分探索に裏切られた話

はじめに

前回の記事 で並列二分探索というテクニックについて触れましたが、ABC388 にてついに並列二分探索の出番が来たので早速試してきました！

この記事でやること

• G - Simultaneous Kagamimochi 2 の並列二分探索解法の解説

問題概要

問題概要

• 𝑁 個の餅がサイズの昇順に並んでおり、各餅 𝑖 のサイズは 𝐴𝑖 である。（つまり、各 𝑖 について 𝐴𝑖 ≤𝐴𝑖+1）
• 𝑄 個の質問に答える。各質問 𝑖 (1 ≤𝑖 ≤𝑄) の内容は次の通り：
  • 餅 𝐿𝑖 から 𝑅𝑖 までを用いて同時に作れる鏡餅の個数の最大値を求めよ。
  • ただし、鏡餅は2つの餅を上下に重ねたものとし、上の餅のサイズは下の餅のサイズの半分以下でなくてはならないものとする。

• 実行時間制限: 2sec
• 1 ≤𝑁 ≤2 ×105
• 1 ≤𝑄 ≤2 ×105

なお、各クエリは E問題 に相当します。

各クエリについての自然な発想として、以下のような貪欲法が浮かびます：

1. 上の餅を昇順に見る。
2. 1.の餅に対して、「まだ使っていない餅であって、下の餅として採用できる最小の餅」とマッチングさせる。

しかし、残念ながらこれは次のようなケースで壊れます。悲しい。

𝐴=[1,2,3,4]

というわけで、各クエリだけでも解法をちゃんと考えなくてはなりません。

解法解説

クエリごとの解法について

作る鏡餅の個数を決め打ち二分探索すればよいです。

すなわち、𝑘 個の鏡餅を作れるか？という判定問題を考えると、小さい餅上位 𝑘 個と大きい餅上位 𝑘 個について、小さい順にマッチングさせられるか を見ればよいです。

つまり、各判定問題は O⁡(𝑘) ⊆O⁡(𝑁) で解けます。

※ 正当性の説明については E問題の公式解説 を参照してください

判定問題を並列に解く

前節での決め打ち二分探索を各クエリごとに行うと、全体 O⁡(𝑄⁢𝑁⁢log⁡𝑁) となり間に合いません。

前回の記事 でも解説した通り、並列二分探索は 𝑄 個の判定問題を並列して解くことで高速化します。

各クエリ 𝑖 ごとに 𝑘𝑖 個の鏡餅を作れるか？という判定問題を考え、並列に解く方法を考えます。

上側に置く餅を昇順に走査します。

このとき、前節での判定問題の解法によると、各クエリについて、マッチング相手として見るべき餅は1つだけです。

また、餅がサイズの昇順に並んでいることに注意すると、クエリを「今見ている上側の餅のindexからマッチング相手のindexまでの距離」*1で昇順にソートしてある時、マッチング相手としてよいかの判定結果に単調性が現れます。

赤と青の矢印が各クエリのマッチング相手。判定結果に単調性がある。

一度 False と判定されたクエリはその時点で 𝑘𝑖 個の鏡餅が作れないことが確定するので、以降は考えなくてよいです。

一度 False となったクエリは見なくてよい

以上を踏まえると、

1. 上側の餅を走査しつつ、
2. イベントソートによって判定するクエリを「今見ている上側の餅のindexからマッチング相手のindexまでの距離」の昇順に持ち、
3. 「持っているクエリのうち先頭にあるもののマッチング相手が False になるならそのクエリの判定結果を False として捨てる」を while文 で処理

の要領で、𝑄 個の判定問題を並列して解くことができます。

クエリの持ち方としては、最小ヒープを使えばよいです。イベント終了の処理は 遅延評価 で行えばよいです。

計算量は O⁡((𝑁 +𝑄)⁢log⁡𝑁⁢log⁡𝑄) です。

実装例

TLEしました笑。

atcoder.jp

おわりに

いかがでしたか？

TL が 2sec な時点で嫌な予感がしていたのですが、案の定でした。悲しい…

PyPyで通る気は今のところしていないのですが、高速な言語だと通ったりするんですかね。

最後までお読みいただきありがとうございました。

PyPyでもちゃんと通る解法（2024/1/15 追記）

各餅に対し「初めてサイズが2倍以上になる餅までの距離（存在しない場合は∞）」を並べた配列を二分探索や尺取法などで作成しておき、これを Range Max が O⁡(1) で取れるデータ構造（Sparse Table 等）に載せておくことで、各クエリの決め打ち二分探索が O⁡(log⁡𝑁) で解けるようです*2。何故気が付かなかったのか…。

atcoder.jp

*1:これは上側の餅に寄らず 𝑘𝑖 に対して一意に定まります。

*2:一応、Range Max でセグ木を使って余分な log をつけても通ります。log は定数。

典型を組み合わせて解くのが楽しかった話【ABC233-Ex 解説】

はじめに

問題の解説をしつつ、競プロ楽しいよ～という布教をする回です。

特に、競プロを未経験の方に、競プロにおいて「典型*1を組み合わせて問題を解く」とはどのようなものなのかを実感してもらえたらと思います。

解説する問題はこいつ ↓ です。とりあえず読んでみてください。

atcoder.jp

※ 問題文中にあるマンハッタン距離とは、縦横でのみ移動できるとしたときの2点間の距離のことです。厳密には、2点 (𝑥1,𝑦1),(𝑥2,𝑦2) 間のマンハッタン距離は |𝑥1 −𝑥2| +|𝑦1 −𝑦2| と表されます*2。

一見するとかなりシンプルな問題ですが、実は出題されているコンテスト（AtCoder Beginner Contest）の中ではかなりの難問に分類されており、典型知識を多重に組み合わせて解く必要があります。

どういった点が難しいのか、どのように解いていくのかを、競プロ関連の前提知識は使わないように（無理なところは適宜端折って）解説していくので、是非最後まで読んで「解けた！」が届く快感を味わってください。

また、できる限り丁寧に解説したつもりではありますが、問題自体が難しいので「解けた！」が万人には届けられないかもしれません。それでも、途中まで読んだ中で「なるほど！」が届けばなと思います。

何が難しいのか

改めて、問題の概要を確認します。

問題概要

• 二次元座標平面上に 𝑁 個の相異なる点が与えられる。
• 𝑄 個の質問に答える。各質問 𝑖 (1 ≤𝑖 ≤𝑄) の内容は次の通り：
  • マンハッタン距離を用いるとして、座標 (𝑋𝑖,𝑌𝑖) から 𝐾𝑖 番目に近い点までの距離は？

• 実行時間制限: 7sec
• 1 ≤𝑁 ≤105
• 1 ≤𝑄 ≤105
• 与えられる各座標の x, y はいずれも 0 以上 105 以下

この問題、答えを求めること自体はそう難しくはありません。

各質問 𝑖 について、以下のような手順で質問に答えることができます：

1. dist⁡[𝑣] :=点⁢𝑣⁢か⁢ら⁡(𝑋𝑖,𝑌𝑖)⁢ま⁢で⁢の⁢マ⁢ン⁢ハ⁢ッ⁢タ⁢ン⁢距⁢離 とする配列 dist を作成する。
2. dist を昇順にソートする。
3. dist の 𝐾𝑖 番目が答え。

しかし、残念ながらこの解法で実装したコードを提出しても、TLE（Time Limit Exceeded / 実行時間制限超過）となり、不正解となります。プログラムの実行時間が長すぎるのでダメ！と怒られるわけですね。

この問題の実行時間制限は 7sec と設定されています。対して、プログラムの実行時間はどのくらいになるのでしょうか。

上記解法においてボトルネック（= 最も実行時間への影響が大きい箇所）となっているのは 2. のソート部分で、計算量は Θ⁡(𝑁⁢log⁡𝑁) である（≒ 𝑁⁢log⁡𝑁 の定数倍程度のステップ数を要する）ことが知られています*3。さらに、𝑄 回質問することを踏まえると、全体の計算量は Θ⁡(𝑄⁢𝑁⁢log⁡𝑁) になります。

一般的なPCが1秒間に処理できるステップ数はおよそ 109 回程度と言われています。つまり、実行時間はざ～っくりと見積もると

105×105×log⁡105109>50

という感じで、ガバガバ議論にはなりますが 7sec には到底収まりそうにないことが分かります。

更に言うと、上記解法の 1. だけ見ても計算量は Θ⁡(𝑁⁢𝑄) で、この時点で既に 7sec は厳しいです。つまりは、「各質問ごとに各点を見る」という当たり前の処理すらさせてもらえない時間制限になっているのです。

シンプルな見た目をしていますが、実はかなりの無茶ぶりを要求されているということがご理解いただけたでしょうか？

いよいよここからは、競プロの典型知識を用いてこの問題を解きほぐしていきます！

典型1: K番目の取得は決め打ち二分探索

早速ですが、問題の言い換えを行います。

𝑖 番目の質問

• 座標 (𝑋𝑖,𝑌𝑖) から 𝐾𝑖 番目に近い点までの距離は？

は、次と同値です。

• 「(𝑋𝑖,𝑌𝑖) からの距離が 𝑑 未満であるような点の個数は 𝐾𝑖 個未満」を満たす、最大の 𝑑 は？

ちょっと長いので、判定問題として

𝑃𝑖⁡(𝑑):=(𝑋𝑖,𝑌𝑖)⁢か⁢ら⁢の⁢マ⁢ン⁢ハ⁢ッ⁢タ⁢ン⁢距⁢離⁢が⁢𝑑⁢未⁢満⁢で⁢あ⁢る⁢よ⁢う⁢な⁢点⁢の⁢個⁢数⁢は⁢𝐾𝑖⁢個⁢未⁢満⁢で⁢あ⁢る⁢か？

とおいてあげると、最終的に求めたいものは 𝑃𝑖⁡(𝑑) =𝑇⁢𝑟⁢𝑢⁢𝑒 となる最大の 𝑑 となります。

さて、ここからが重要です。

𝑃𝑖⁡(𝑑) は、𝑑 についての単調性、すなわち

• 𝑑1 <𝑑2 なる 𝑑1,𝑑2 について、𝑃𝑖⁡(𝑑1) =𝐹⁡𝑎⁢𝑙⁢𝑠⁢𝑒 ⇒𝑃𝑖⁡(𝑑2) =𝐹⁡𝑎⁢𝑙⁢𝑠⁢𝑒

という性質を持ちます。（𝑑 を大きく取れば距離 𝑑 未満の範囲にある点の個数は単調増加するので、当然です。）

この性質を用いれば、次のスライドのように範囲を半分ずつ絞っていくことで、判定問題 𝑃𝑖⁡(𝑑) を高々 log2⁢(𝑑⁢の⁢範⁢囲) 回程度解くだけで所望の 𝑑 を求めることができます。いわゆる二分探索というやつです。

なお、0 ≤𝑑 ≤2 ×105 と評価できるので、𝑃𝑖⁡(𝑑) を解く回数は高々 log2⁢(2 ×105) <20 回程度です。

このように、何かしらの数値を求める問題に対して単調性のある判定問題に置き換えて二分探索を適用する手法を、決め打ち二分探索と呼びます。

とはいえ、まだ問題はほとんど解決していません。

というのも、決め打った 𝑑 に対して判定問題 𝑃𝑖⁡(𝑑) が高速に解けないようでは、判定問題に置き換えた意味がないからです。

ここからは、この判定問題を解く方法を掘り下げていきます。

典型2: マンハッタン距離は45°回転

(𝑋𝑖,𝑌𝑖)からマンハッタン距離が 𝑑 未満である領域を図示すると、次のような菱形になります。𝑃𝑖⁡(𝑑) は、この領域に含まれる点が 𝐾𝑖 個未満かどうかを判定する問題に相当します。

(𝑋𝑖,𝑌𝑖) からのマンハッタン距離が 𝑑 未満の領域（境界線は含めない）

ただ、この菱形というのは色々と扱いづらいです。

そこで、以下のような座標変換を行ってみます：

𝑢:=𝑥+𝑦𝑣:=𝑥−𝑦

各位置ベクトルに対して45°の回転行列をかけたあと √2 倍すると理解してもよいでしょう。この変換により先ほどの領域は、u-v座標系で表すと次の通り正方形になります。

u-v座標系における、x-y座標系で𝑃𝑖⁡(𝑑)の対象だった領域

これにより、正方形区画に点が何個含まれるかを求める問題に帰着することができました。これが解けるかどうかは別として、少なくとも菱形で考えるよりはなんとなく見通しがよさそうですよね。

このように、マンハッタン距離に関する問題はとりあえず45°回転してから考えるというのが典型となっています。

以降はこのu-v座標系で議論を進めていきます。

なお、座標変換により

• 0 ≤𝑢 ≤2 ×105
• −105 ≤𝑣 ≤105

となっています。

典型3: 平面走査

「正方形区画内の点の個数」を、もう少しだけ扱いやすい形で表現します。

座標平面上の各格子点を二次元グリッドのマス目に置き換え、点が存在するなら 1, しないなら 0 を書き込みます。これにより、二次元グリッド上で正方形区画内の総和を求める問題として定式化できます。

グリッドに置き換える

なお、𝑢,𝑣 の範囲を思い出すと、グリッドサイズが縦横いずれも 2 ×105 となっていることに注意しましょう。（つまり、プログラム上で明にこのグリッドを二次元配列で持つことはできません*4。悲しい。）

さて、このような大きなグリッドに関する問題に対しては、平面走査と呼ばれる典型的なアプローチがあります。まずは、次のスライドをご覧ください。

何がしたいのかは一旦置いておくとして、行ごとに走査しながら配列に値を加算していく様子が伺えます。このように、何かしらのデータを持ち、行（あるいは列）ごとにデータを更新しながら平面を走査することを、総称して平面走査と呼びます。

急にわけのわからない操作をし始めたぞコイツ…？と思われているかもしれませんが、実はこの操作に超重要な秘密が隠されています！

↓この緑部分の総和（4）から…

↓この黄色部分の総和（1）を引くと…

求めたかった正方形区画の総和（4-1 = 3）が得られます！

この理由は単純で、矩形和の引き算に対応しているからです。

2つの矩形和の差が、求めたかった正方形の総和に対応する。

ということで、正方形区画の総和は次の要領で求められます：

1. 平面走査を行う。その途中で、次の二つの値を別途保存しておく。
   • 黄色部分の総和
   • 緑部分の総和
2. 上記で保存した2数の差が答え。

計算量を考えてみましょう。グリッドの横幅を 𝑊, 𝑑⁢の⁢範⁢囲 を 𝐷 とおくと、

• 平面走査については、長さ 𝑊 の配列を用意したうえで更新が 𝑁 回発生するので、Θ⁢(𝑊 +𝑁)
  • 黄色部分、緑部分の総和を求めるのは Θ⁡(𝑑) で、これは平面走査の 𝑊 があるためボトルネックにならない。
• → 判定問題 𝑃𝑖⁡(𝑑) 1回あたり Θ⁢(𝑊 +𝑁)
• → 全体 Θ⁡(𝑄⁢(𝑊 +𝑁)⁢log⁡𝐷)

あれ…？？？？？？全然よくなってない…（´・ω・｀）

典型4: 並列二分探索

平面走査などというトンデモびっくりアプローチまで取り出して頑張ったのに、計算量は何も改善されていませんでした。果たして、今までの努力は何だったのでしょうか…。

ここで、全く新しい視点を導入します。

今までは質問1つごとにどう解くか？ということだけを考えていましたが、思い切って 𝑄 個の質問を一斉に解くことを考えてみましょう。

より正確には、𝑄 個の質問に対して二分探索を同時に実施し、𝑄 個の判定問題 𝑃𝑖⁡(𝑑𝑖) (1 ≤𝑖 ≤𝑄) を1回の平面走査中にまとめて解くことにします。

今までの議論と同様に考えると、グリッド上に 𝑄 個の正方形区画が与えられ、それぞれの総和を求める問題を解くことになります。実はこれ、正方形が何個だろうとやることは全く変りません。すなわち、平面走査の途中で各正方形区画ごとに対応した黄色部分の総和と緑部分の総和を保存し、最後にそれぞれ差をとればよいです。

ありがとう、平面走査…！！！！

計算量を考えます。

• 平面走査については、長さ 𝑊 の配列を用意したうえで更新が 𝑁 回発生するので、Θ⁢(𝑊 +𝑁)
  • 黄色部分、緑部分の総和を求めるのは各 𝑖 ごとに Θ⁡(𝑑𝑖) で、まとめると Θ(𝑄⁢𝐷)
• → 判定問題 𝑃𝑖⁡(𝑑𝑖) (1 ≤𝑖 ≤𝑄) をまとめて解く処理1回あたり Θ(𝑊 +𝑁 +𝑄⁢𝐷)
• → 全体 Θ((𝑊 +𝑁 +𝑄⁢𝐷)log⁡𝐷)

惜しい！！

ここにきて、配列の区間和を計算する箇所がボトルネックになってしまいました。

典型5: Segment Tree

いよいよ最後の典型です。

平面走査中に行いたい処理を改めて書き下すと、次の2つです：

• 1点更新
• 区間和取得

平面走査のデータを通常の配列で持ってしまうと、どうしても区間和取得に区間長分の計算量がかかってしまいます。

そこで登場するのが、競プロにおいて典型中の典型、キング・オブ・典型データ構造である Segment Tree です*5。

Segment Tree でデータを持つと、上記の2操作がいずれも Θ⁡(log⁡𝑊)（𝑊: グリッドの横幅）で処理できます。

詳細は端折りますが、ざっくり概要だけ解説します。

Segment Tree は図のような完全二分木の構造をしており、葉の部分には平面走査のデータ、その他の節点には左右の子の和が保存されています。ここで、完全二分木の高さは Θ⁡(log⁡(𝑊)) 、節点数は Θ⁡(𝑊) になることに注意してください。

Segment Tree

このようにデータを持つと、1点更新は図のように Θ⁡(log⁡𝑊) 個の節点の更新で済みます。

1点更新

また、任意の区間和についても参照する節点数が Θ⁡(log⁡𝑊) で抑えられることが示せます*6。

区間和の取得

というわけで、Segment Tree を用いれば1点更新と区間和の取得の両方を高速に処理できます。

Segment Tree を用いた場合の計算量は、次の通りです：

• 平面走査については、長さ 𝑊 の配列を用意したうえで更新が 𝑁 回発生するので、Θ⁢(𝑊 +𝑁)
  • 黄色部分、緑部分の総和を求めるのは各 𝑖 ごとに Θ(log⁡𝑊) で、まとめると Θ(𝑄⁢log⁡𝑊)
• → 判定問題 𝑃𝑖⁡(𝑑𝑖) (1 ≤𝑖 ≤𝑄) をまとめて解く処理1回あたり Θ(𝑊 +𝑁 +𝑄⁢log⁡𝑊)
• → 全体 Θ((𝑊 +𝑁 +𝑄⁢log⁡𝑊)log⁡𝐷)

長い長い戦いの末、ついに2乗の項を完全に消すことができました。

ステップ数の計算はもう面倒なのでしませんが、実行時間はちゃんとこの解法で制限に収まります。

以上、解説を終わります。お疲れさまでした。

おわりに

典型に典型を重ねて考察を進める雰囲気を、少しは味わっていただけたでしょうか？

個人的には、自身の知りうる典型知識を如何にうまく組み合わせて解法に到達するかを考えるのが、競プロ（というか、AtCoder Beginner Contest）の醍醐味のように思います。

この記事を読んで、内容を全部理解できたかどうかは正直どうでもいいです。（そもそも初学者が解く問題の難易度じゃないので。）

それよりも、なんとなく「面白そう！」と思った方は、間違いなく競プロ向いてると思うので、是非この機会に始めてみてはいかがでしょうか。

逆に、「難しそう…こわ…近寄らんとこ…」と感じた方であっても、「はじめに」で述べた通りこの記事の中に少しでも「なるほど！」と高揚感を得られた箇所があれば、それが競プロを楽しむモチベーションになっているということだけでもを知っていただければ幸いです。

それでは。

参考文献

Editorial - AtCoder Beginner Contest 233

マンハッタン距離と45度回転 - ま

CODE THANKS FESTIVAL 2017 H - Union Sets （並列二分探索解法） - ARMERIA

*1:「よく使うテクニック」といった意味です。

*2:ちなみに、一般的に「距離」と呼ばれる三平方の定理のアレは「ユークリッド距離」と呼びます。

*3:教科書の最初の方に載ってるような内容なので、このくらいは認めさせてください…。

*4:余談ですが、仮にこのグリッドが明に持てるサイズ （106 程度）に収まっていれば、二次元累積和と呼ばれるテクニックを用いることで正方形区画内の総和を定数時間で求められるため、全体 Θ⁢(𝑁 +𝑄⁢log⁡(𝑑⁢の⁢範⁢囲)) で解くことができます。

*5:諸説あります

*6:各深さにおいて、参照する節点が高々 2 個で済むからです。3つ以上になると、親を代わりに参照できるはずですよね。

【Trie木】ABC377G - Edit to Match の図説

この記事でやること

ABC377G-Edit to Match の図説です。やってることは公式解説と同じです。

解法

便宜上 𝑆0 を空文字列として、𝑆0,…,𝑆3 まで処理済みで 𝑆4 を push したTrie木を考えます。（Trie木を知らない方は、参考文献に記載した記事などを参照してください。）

𝑆4 を push したTrie木

ここで、𝑆4 に対する解は min0≤𝑖<4⁢{dist⁡(𝑆4,𝑆𝑖)} になります。上図の場合、𝑆4 から 𝑆0 または 𝑆2 が最短で、答えは 3 になります。

これは、各ノード 𝑣 について 𝑣.𝑑⁢𝑝 :=𝑣⁢を⁢根⁢と⁢す⁢る⁢部⁢分⁢木⁢内⁢で⁢の⁢min0≤𝑖<4⁢{dist⁡(𝑣,𝑆𝑖)} を記録しておくことで高速に解くことができます。

𝑣.𝑑⁢𝑝、部分木内に 𝑆 が存在しない場合は ∞

Trie木の根から 𝑆4 のノードまでのパスを 𝑃 としたとき、求める解は min𝑣∈𝑃⁢{𝑣.𝑑⁢𝑝 +dist⁡(𝑣,𝑆4)} となります。

また、𝑣.𝑑⁢𝑝 は push の帰りがけにchminすることで簡単に更新できますね。

帰りがけにchmin

計算量は、各文字列 𝑆𝑖 に対し O⁡(|𝑆𝑖|)*1 になります。

提出コード

Python (PyPy 3.10-v7.3.12)

Submission #59196277 - TOYOTA SYSTEMS Programming Contest 2024（AtCoder Beginner Contest 377)

おわりに

公式解説、こういう挿絵とかあったらもっとわかりやすいのにね…。

参考文献

トライ木(Trie木) の解説と実装【接頭辞(prefix) を利用したデータ構造】 | アルゴリズムロジック

*1:Trie木の枝をhashmapで持った場合の期待計算量。配列で持つなら文字種数σ倍が付きます。

Aho-Corasick法をしっかりと理解する話

はじめに

𝑁 個のキーワード 𝑆𝑘⁡(0 ≤𝑘 <𝑁) が定められている環境下で、入力テキスト 𝑇 からそれらの出現位置を全列挙するクエリを考えます。

文字列の一致判定はローリングハッシュなどを用いて定数時間でできるにしても、含まれている場所が分からない以上は、テキスト1つあたり O⁡(|𝑇|⁢𝑁) 程度の計算量になってしまいそうです。

Aho-Corasick法は、キーワードに対してO⁡(∑|𝑆𝑘|) で工夫したTrie木を作る*1ことで、クエリあたりO⁢(|𝑇| +∑|𝑆𝑘| +出⁢現⁢回⁢数) で処理できます。強い。*2

本記事は、Aho-Corasick法の考え方やざっくりとしたアルゴリズム、計算量評価の方法をまとめたものです。これを読めば、Aho-Corasick法をしっかりと理解でき、自力で実装できるようになる…かも？

この記事でやること

Aho-Corasick法をできる限り丁寧に解説します。

また、最後にPython(PyPy)での実装例を載せてあるので、ご自由にお使いください。

Aho-Corasick法の解説

Trie木を効率よく辿る

まずキーワードでTrie木を作っておきます。

【図1】𝑆 =[𝑎,𝑎⁢𝑏⁢𝑏⁢𝑐,𝑏⁢𝑎,𝑏⁢𝑏⁢𝑐⁢𝑎,𝑐⁢𝑏⁢𝑎] のTrie木

このTrie木において、𝑇 の 𝑖 文字目（以後、𝑇⁡[𝑖] と表記、0-indexedとする。）を左端とするキーワードを見つけたい場合、

• Trie木の根に"駒"を置き、𝑗 :=𝑖 とする。
• 駒のいるノードに文字 𝑇⁡[𝑗] の枝が存在する限り、その先のノードに駒を進めて 𝑗 をインクリメントすることを繰り返す。

とすればよいです。

具体例として、𝑇 =𝑎⁢𝑏⁢𝑏⁢𝑐⁢𝑏⁢𝑎⁢𝑐 として考えてみましょう。

【図2-1】 𝑇 =𝑎⁢𝑏⁢𝑏⁢𝑐⁢𝑏⁢𝑎⁢𝑐 に対し、左端を 𝑇⁡[0] とするキーワードを見つけたい。

【図2-2】𝑎⁢𝑏⁢𝑏⁢𝑐⁢𝑏⁢𝑎⁢𝑐、キーワード 𝑎,𝑎⁢𝑏⁢𝑏⁢𝑐 を発見！

しかし、これを各左端について行うと 𝑂⁡(|𝑇|⁢max⁡(|𝑆𝑘|)) になるので、まだまだ嬉しくありません。

Aho-Corasick法では、この左端の切り替えを効率よく行います。

試しに、【図2-2】の続きを考えてみましょう。

今、𝑇⁡[4]⁢( =𝑏) の枝が駒のいるノードに存在しないため、駒を進めるループを脱出した状態です。なので、次は駒を根に戻し 𝑗 :=1 として再度駒を進めるループに入ります。

このループは 𝑎⁢𝑏⁢𝑏⁢𝑐⁢𝑏⁢𝑎⁢𝑐 と駒を進めますが、この 𝑏⁢𝑏⁢𝑐 と進める部分は、前回のループで 𝑎⁢𝑏⁢𝑏⁢𝑐 と進めたことを思い出すと無駄な動きをしているような気がします。どうせなら根に戻らず直接 𝑏⁢𝑏⁢𝑐 のノードに移したいですね。

【図3】駒を 𝑎⁢𝑏⁢𝑏⁢𝑐 から直接 𝑏⁢𝑏⁢𝑐 に移したい

これを一般的に述べます。

便宜上、根から駒がいるノードまでのパスにより得られるテキストを読み込んでいるテキストと呼ぶことにします。

ループが止まった後の駒の移動先は、「読み込んでいるテキストの接尾辞（自身は除く）であって、キーワードの接頭辞として含まれる最長のもの」を表すノードとすればよいです。（なので、𝑎⁢𝑏⁢𝑏⁢𝑐⁢𝑏⁢𝑎⁢𝑐 で止まったあと駒は 𝑐 に移動します。）

Trie木の構造上、根以外のどのノードに駒があったとしても対応した上記の移動先のノードが一意に定まります。そこで、各ノード 𝑛⁢𝑜⁢𝑑⁢𝑒 について対応する移動先のノードを 𝑛⁢𝑜⁢𝑑⁢𝑒.fail と表すことにしましょう。根以外の全てのノードに対して fail が既に求まっているものと仮定すると、ざっくり次のような手順で全てのキーワードを列挙できます：

• 駒をTrie木の根に置き、𝑖 =0,…,|𝑇| −1 に対して、順に以下を行う:
  
  • 駒がいるノードが根でなく、かつ文字 𝑇⁡[𝑖] の枝が存在しない限り、駒を fail に移すことを繰り返す。
  • 駒がいるノードに 𝑇⁡[𝑖] の枝が存在するなら、その先のノードに移す。
  • 読み込んでいるテキストの接尾辞として含まれるキーワードを列挙する。（※）

この手順で全てのキーワードが過不足なく列挙できていることは、（※）の操作によって 𝑇⁡[𝑖] を末尾とするキーワードが全て列挙されていることから分かります。

列挙の方法としては、 fail を根に到達するまで辿り、キーワードそのものになっているものを抽出すればよいです。さらに、fail を適切に経路圧縮して得られる「読み込んでいるテキストの接尾辞（自身は除く）として含まれる、最長のキーワード」を表すノードへの辺を用意しておけば、計算量を出現回数で抑えられます。*3

ここで、次のスライドショーを見てもらうとわかる通り、読み込んでいるテキストは駒の移動に合わせて常に右方向に動きます。つまり、尺取法の要領で駒の移動回数が高々 2⁢|𝑇| 回に抑えられるわけですね。

ということで、あとは各ノードに対する fail を前計算するパートを倒せばよいです。

fail の前計算

fail の定義を思い出すと、fail の移動先のノードはTrie木において自身より浅い位置にあります。そこで、𝑑 ( ≥1) について深さ 𝑑 −1 以下の全てのノードの fail が求まっていると仮定して、深さ 𝑑 のノード 𝑞 について 𝑞.fail を求め方を考えます。また、便宜上 𝑟⁢𝑜⁢𝑜⁢𝑡.fail :=𝑟⁢𝑜⁢𝑜⁢𝑡 としておきます。

𝑞 の親を 𝑝 とします。また、根から 𝑝 までのパスによるテキストを 𝑤、𝑝 から文字 𝑐 の枝を通って 𝑞 に到達するものとします。

𝑞.fail が指すノードは、「𝑤⁢𝑐 の接尾辞（自身は除く）であって、キーワードの接頭辞として含まれる最長のもの」を表すノードでした。これは必ず 𝑤 の接尾辞に 𝑐 をつけた形をしている点に注意すると、次の手順で 𝑞.fail を求めることができます：

• 𝑛⁢𝑜⁢𝑑⁢𝑒 :=𝑝 とする。
• 𝑛⁢𝑜⁢𝑑⁢𝑒 が根でない かつ 𝑛⁢𝑜⁢𝑑⁢𝑒.fail が 𝑐 の枝を持たない限り、𝑛⁢𝑜⁢𝑑⁢𝑒 :=𝑛⁢𝑜⁢𝑑⁢𝑒.fail とすることを繰り返す。
• 𝑛⁢𝑜⁢𝑑⁢𝑒 が 𝑐 の枝を持つなら 𝑞.fail :=𝑛⁢𝑜⁢𝑑⁢𝑒⁢の⁢枝⁢𝑐⁢の⁢先⁢の⁢ノ⁢ー⁢ド とする。そうでないなら、𝑞.fail :=𝑛⁢𝑜⁢𝑑⁢𝑒 とする*4。

そして、根から始めるBFSにて 𝑝 を見るときに、各子ノード 𝑞 に対して、上記の処理をしたのちにキューに入れていくことで、全てのノードの fail を求めることができます。

計算量を考えてみましょう。BFS部分は O⁡(∑|𝑆𝑘|) でいいとして、fail を求めるループが問題です。実は、各キーワード 𝑆𝑘 について、𝑆𝑘 の各文字のノードによるループの合計回数が高々 2⁢|𝑆𝑘| 回で抑えられます（※※）。よって、前計算全体で O⁡(∑|𝑆𝑘|) となります。

（※※）の証明

キーワード 𝑆 について、「𝑆 の接尾辞（自身は除く）であって、キーワードの接頭辞として含まれる最長のもの」の長さを 𝑙0 とし、𝑆 の末尾 max⁡(𝑙0,1) 文字による文字列を 𝑆(0) とします。次に、𝑆 から末尾の 𝑆(0) を削った部分文字列対して同様の操作を行い 𝑙1 と 𝑆(1) を得ます。以下これを繰り返して 𝑙0,…,𝑙𝑚−1、𝑆(0),…,𝑆(𝑚−1) が得られたところで 𝑆 =𝑆(𝑚−1)…𝑆(0) となったとします。また、𝑙𝑚 :=0、𝑆(𝑚) :=（空⁢文⁢字⁢列） とします。

このとき、各 𝑖 =𝑚 −1,…,0 について、𝑆(𝑖) の各文字のノードの　fail を求めるループは高々 |𝑆(𝑖+1)| +|𝑆(𝑖)| 回で抑えられます。なぜなら、

• 𝑆(𝑖)⁡[0] のノードについては、親ノードの fail 移動先の深さが 𝑙𝑖+1 であることに注意すると、高々 |𝑆(𝑖+1)| 回のループで求まる。
• 𝑆(𝑖)⁡[𝑗] (1 ≤𝑗 <|𝑆(𝑖)|) のノードについては、親ノードの fail 移動先のノードに必ず 𝑆(𝑖)⁡[𝑗] の枝が存在するため、ループは1回しか発生しない。

となるからです。  ◻

競プロへの応用

冒頭で示したような列挙クエリは、計算量的に競プロではおそらく問われません。

そこで、列挙の代わりに各キーワードの出現回数を数え上げる問題を考えてみましょう。

単純なのは、各キーワードがそれぞれ何回出現したかを表す配列を用意して、前述の（※）を利用してインクリメントしていくという方法です。しかし、これでは計算量が変わらず出現回数に依存してしまいます。

解決法としては、各ノードについて（※）が始まるノードになった回数をカウントしておき、最後に木dpの要領で、回数を fail に沿って配りつつキーワードの末尾ノードに到達するごとにその値を配列に加算していけばよいです。

問題例 ↓

atcoder.jp

Pythonでの実装例

• PyPy3.10-v7.3.12 にて動作確認済。

from collections import deque
sys.setrecursionlimit(200000) # enum() を使わないならコメントアウトしてOK。
class AhoCorasick_Node:
    def __init__(self, a=None) -> None:
        self.a = a
        self.son = {}
        self.fail = self
        self.output = None
        self.nxt_output_node = self
        self.dp = 0
    
    def find_nxt(self):
        if self.nxt_output_node is self:
            return self
        if self.nxt_output_node.output is None:
            self.nxt_output_node = self.nxt_output_node.find_nxt()
        return self.nxt_output_node

class AhoCorasick:
    def __init__(self, words) -> None:
        self.words = words
        self.head = AhoCorasick_Node()
        self._build_Trie()
        self._set_fail()

    def _build_Trie(self):
        for w_idx, word in enumerate(self.words):
            node = self.head
            for a in word:
                if a not in node.son:
                    node.son[a] = AhoCorasick_Node(a)
                node = node.son[a]
            node.output = w_idx

    def _set_fail(self):
        dq = deque([self.head])
        while dq:
            now = dq.pop()
            for a, nxt in now.son.items():
                node = now.fail
                while a not in node.son and node.fail is not node:
                    node = node.fail
                if a in node.son and node.son[a] is not nxt:
                    nxt.fail = node.son[a]
                else:
                    nxt.fail = self.head
                nxt.nxt_output_node = nxt.fail
                dq.appendleft(nxt)
    
    def enum(self, text):
        node = self.head
        ret = [[] for _ in range(len(self.words))]
        for i, a in enumerate(text):
            while a not in node.son and node.fail is not node:
                node = node.fail
            if a in node.son:
                node = node.son[a]

            if node.output is not None:
                w_idx = node.output
                ret[w_idx].append(i+1-len(self.words[w_idx]))
            output_node = node
            while output_node.find_nxt() is not self.head:
                output_node = output_node.nxt_output_node
                w_idx = output_node.output
                ret[w_idx].append(i+1-len(self.words[w_idx]))

        return ret
    
    def cnt(self, text):
        tps = []
        dq = deque([self.head])
        while dq:
            now = dq.pop()
            now.dp = 0
            tps.append(now)
            for nxt in now.son.values():
                dq.appendleft(nxt)

        node = self.head
        ret = [0]*len(self.words)
        for a in text:
            while a not in node.son and node.fail is not node:
                node = node.fail
            if a in node.son:
                node = node.son[a]
            node.dp += 1

        for node in reversed(tps[1:]):
            if node.output is not None:
                w_idx = node.output
                ret[w_idx] += node.dp
            node.fail.dp += node.dp
        return ret

使用例↓

Submission #56078430 - Toyota Programming Contest 2024#7（AtCoder Beginner Contest 362）

ちなみに

これは驚愕の事実なのですが、Aho-Corasickは 「えいほこらしっく」であって、「あほこらしっく」ではないそうです。

エイホ–コラシック法 - Wikipedia

参考文献

1. Aho Corasick 法 - naoyaのはてなダイアリー

2. Alfred V. Aho and Margaret J. Corasick. 1975. Efficient string matching: an aid to bibliographic search. Commun. ACM 18, 6 (June 1975), 333–340.

*1:正確に述べると、これは内部のTrie木でハッシュテーブルを用いた場合の期待計算量になっています。配列でTrie木を作る場合は文字種数σ倍かかるかわりに最悪計算量になります。どちらで実装するかはお好みでどうぞ。

*2:なお 𝑇 =𝑎⁢𝑎⁢𝑎⁢𝑎,𝑆 =[𝑎,𝑎⁢𝑎,𝑎⁢𝑎⁢𝑎,𝑎⁢𝑎⁢𝑎⁢𝑎] のようなケースを考えれば出現回数を壊せるので、競プロ目線では出現回数に依存しないクエリを考えることが多そうです。

*3:経路圧縮のために余分に潜る操作は、高々Trie木のノード数回しか発生しないため、Trie木構築時の計算量に押しつけて"ならす"ことができます。

*4:この場合での 𝑛⁢𝑜⁢𝑑⁢𝑒 は必ず根になります。そこで、根について本来存在しない文字の枝を全て自身へのループとしておけば、ここの場合分けをなくすこともできます。