清一色のシャンテン数を01BFSで計算する

概要

  • 清一色のシャンテン数は01BFSで計算できる。
  • 機械的な計算なので、正当性が簡単に証明できる。

シャンテン数とは

  • テンパイするまでに必要な牌の最小枚数のこと。
  • 例えば 🀉🀊🀍🀎🀑🀒🀙🀚🀛🀜🀜🀀🀁 のとき、🀈🀌 を引いて🀀🀁を捨てればテンパイとなる。これより少ない枚数の入れ替えでテンパイすることはできないので、この手牌のシャンテン数は2となる。
  • この例のようにシャンテン数がわかりやすいケースがほとんどだが、清一色のような場合はシャンテン数を求めにくい。
  • 特にプログラムでシャンテン数を計算する場合、正しいコードを書くことは難しい。また、正当性を保証するのも難しい。

考える問題

  • 清一色だけを考える。以下では萬子🀇🀈🀉🀊🀋🀌🀍🀎🀏だけを扱う。

    • 複数色ある場合や字牌を含む場合でも、清一色のシャンテン数さえ計算できていれば、軽いdpで解くことができるため。
    • 正確には、x面子y雀頭を完成形とみなしたときの清一色のシャンテン数を全ての 0 \le x \le 4, \ 0 \le y \le 1 について求めておく必要があるが、これらはほぼ同様に計算できる。
  • 4面子1雀頭の一般手だけを考える。

    • 七対子国士無双のシャンテン数は、簡単な計算で求めることができるため。
  • 副露は考えない。

    • x回副露している場合は、4-x面子1雀頭を完成形とみなしてシャンテン数を求めれば良く、ほぼ同様の計算となるため。

Notation

  • 手牌を h = (h_1, \dots, h_9) で表す。
    • h_i := 萬子 i が何枚あるか。
  • 次の条件を満たす手牌全体の集合を H とおく。
    • 0 \le h_i \le 4, \quad 1 \le i \le 9
    • \sum_{i=1}^9 h_i \le 14
  • 手牌 h \in H に対して、4面子1雀頭を作るのに必要な牌の最小枚数を r(h) とおく。
    • 厳密には、4面子1雀頭を余りなく作れる手牌(以下完成形と呼ぶ)の集合を W \subseteq H とおいて、r(h) を次で定義する。

\displaystyle r(h) := \min_{w \in W} \sum_{i=1}^{9} \max(w_i - h_i, 0) \tag{1}

手牌 h のシャンテン数は r(h) - 1 となるので、r(h) を計算すれば良い。

r(h) の計算

イデア

重み付き有向グラフを以下で定める。

  • 手牌全体の集合 H を頂点集合とする。

  • 次のルールで辺を張る。

    • 手牌 h \in H から、1枚牌を追加して得られる手牌に距離0の辺を張る。

      • 厳密には、各 i=1,\dots,9 に対して h^{+i} := (h_1, \dots, h_i+1, \dots, h_9) とおき、h^{+i} \in H ならば h から h^{+i} に距離0 の辺を張る。
    • 手牌 h \in H から、1枚牌を削除して得られる手牌に距離1の辺を張る。

      • 厳密には、各 i=1,\dots,9 に対して h^{-i} := (h_1, \dots, h_{i}-1, \dots, h_9) とおき、h^{-i} \in H ならば h から h^{-i} に距離1 の辺を張る。

このとき、完成形 w \in W から h \in H の最短距離を d(w, h) とおくと、次の関係が成立する。


r(h) = \min_{w \in W} d(w, h) \tag{2}

証明

(1) より、\sum_{i=1}^{9} \max(w_i - h_i, 0) = d(w, h) を示せば良い。

これは辺の張り方から簡単に従う。(w の牌を追加・削除して h に一致させていく様子をイメージするとわかる。)

完成形 w の列挙

  • 面子は 7 + 9 通り
  • 雀頭は 9 通り

であるから、(7 + 9)^{4} \times 9 = 589824 通りを考えて、その各牌の枚数 h を計算する。

1 \le i \le 9 に対して 0 \le h_i \le 4 を満たすものを完成形として列挙する。

\min_{w \in W} d(w, h) の計算

辺の距離が0と1しかないので、01BFSを用いて計算することができる。

頂点数と辺の数を計算すると

  • 頂点数 |H| = 405350
  • 辺の数 |E| = 3648150 4791600

となる。O(|H| + |E|)\min_{w \in W} d(w, h) を計算することができる。

実装とテスト

https://github.com/habara-k/flush-shanten

3.2GHz CPUと8GB RAMを用いたところ、

  • python実装では42±1[sec]
  • rust実装では2.15±0.08[sec]

程度となった。

python実装

import numpy as np
from collections import deque
from copy import deepcopy
import time


# 面子のリスト
sets = [
    np.array([1, 1, 1, 0, 0, 0, 0, 0, 0]),
    np.array([0, 1, 1, 1, 0, 0, 0, 0, 0]),
    np.array([0, 0, 1, 1, 1, 0, 0, 0, 0]),
    np.array([0, 0, 0, 1, 1, 1, 0, 0, 0]),
    np.array([0, 0, 0, 0, 1, 1, 1, 0, 0]),
    np.array([0, 0, 0, 0, 0, 1, 1, 1, 0]),
    np.array([0, 0, 0, 0, 0, 0, 1, 1, 1]),
    np.array([3, 0, 0, 0, 0, 0, 0, 0, 0]),
    np.array([0, 3, 0, 0, 0, 0, 0, 0, 0]),
    np.array([0, 0, 3, 0, 0, 0, 0, 0, 0]),
    np.array([0, 0, 0, 3, 0, 0, 0, 0, 0]),
    np.array([0, 0, 0, 0, 3, 0, 0, 0, 0]),
    np.array([0, 0, 0, 0, 0, 3, 0, 0, 0]),
    np.array([0, 0, 0, 0, 0, 0, 3, 0, 0]),
    np.array([0, 0, 0, 0, 0, 0, 0, 3, 0]),
    np.array([0, 0, 0, 0, 0, 0, 0, 0, 3]),
]

# 雀頭のリスト
heads = [
    np.array([2, 0, 0, 0, 0, 0, 0, 0, 0]),
    np.array([0, 2, 0, 0, 0, 0, 0, 0, 0]),
    np.array([0, 0, 2, 0, 0, 0, 0, 0, 0]),
    np.array([0, 0, 0, 2, 0, 0, 0, 0, 0]),
    np.array([0, 0, 0, 0, 2, 0, 0, 0, 0]),
    np.array([0, 0, 0, 0, 0, 2, 0, 0, 0]),
    np.array([0, 0, 0, 0, 0, 0, 2, 0, 0]),
    np.array([0, 0, 0, 0, 0, 0, 0, 2, 0]),
    np.array([0, 0, 0, 0, 0, 0, 0, 0, 2]),
]


# hand -> int (5進数とみなす)
def encode(hand):
    ret = 0
    for e in hand:
        ret = (ret * 5) + e
    return ret


# int -> hand
def decode(hash):
    ret = np.zeros(9, dtype=int)
    for i in range(9):
        ret[8-i] = hash % 5
        hash //= 5
    return ret


# 4面子1雀頭の形を列挙
def complete_hands():
    ret = {}
    for s1 in sets:
        for s2 in sets:
            for s3 in sets:
                for s4 in sets:
                    for head in heads:
                        hand = s1 + s2 + s3 + s4 + head
                        if (hand <= 4).all():
                            ret[encode(hand)] = hand
    return ret


# 01BFSを用いてシャンテン数を計算
def bfs(W):
    dist = {}
    deq = deque()
    for hash, hand in W.items():
        dist[hash] = 0
        deq.append((hash, hand))

    while len(deq) > 0:
        hash, hand = deq.popleft()
        for k in range(9):
            if hand[k] < 4 and sum(hand) < 14:
                hand_add = deepcopy(hand)
                hand_add[k] += 1
                hash_add = encode(hand_add)
                if hash_add not in dist:
                    dist[hash_add] = dist[hash]
                    deq.appendleft((hash_add, hand_add))
            if hand[k] > 0:
                hand_sub = deepcopy(hand)
                hand_sub[k] -= 1
                hash_sub = encode(hand_sub)
                if hash_sub not in dist:
                    dist[hash_sub] = dist[hash] + 1
                    deq.append((hash_sub, hand_sub))

    return dist


def main():
    W = complete_hands()
    shanten = bfs(W)
    assert(len(shanten) == 405350)

    # 計算結果を保存
    with open("shanten.txt", mode='w') as f:
        for hash, shanten in shanten.items():
            hand = decode(hash)
            f.write("{} {} {} {} {} {} {} {} {} {}\n".format(
                hand[0], hand[1], hand[2], hand[3], hand[4],
                hand[5], hand[6], hand[7], hand[8], shanten))


if __name__ == '__main__':
    start = time.time()
    main()
    elapsed_time = time.time() - start
    print("elapsed_time: {} [sec]".format(elapsed_time))

rust実装

use std::collections::{VecDeque, BTreeSet, BTreeMap};
use std::fs;
use std::io::Write;
use std::time::Instant;

// 4面子1雀頭の形を列挙
fn complete_hands() -> BTreeSet<Vec<i32>> {
    // 面子のリスト
    let sets = vec![
        vec![1, 1, 1, 0, 0, 0, 0, 0, 0],
        vec![0, 1, 1, 1, 0, 0, 0, 0, 0],
        vec![0, 0, 1, 1, 1, 0, 0, 0, 0],
        vec![0, 0, 0, 1, 1, 1, 0, 0, 0],
        vec![0, 0, 0, 0, 1, 1, 1, 0, 0],
        vec![0, 0, 0, 0, 0, 1, 1, 1, 0],
        vec![0, 0, 0, 0, 0, 0, 1, 1, 1],
        vec![3, 0, 0, 0, 0, 0, 0, 0, 0],
        vec![0, 3, 0, 0, 0, 0, 0, 0, 0],
        vec![0, 0, 3, 0, 0, 0, 0, 0, 0],
        vec![0, 0, 0, 3, 0, 0, 0, 0, 0],
        vec![0, 0, 0, 0, 3, 0, 0, 0, 0],
        vec![0, 0, 0, 0, 0, 3, 0, 0, 0],
        vec![0, 0, 0, 0, 0, 0, 3, 0, 0],
        vec![0, 0, 0, 0, 0, 0, 0, 3, 0],
        vec![0, 0, 0, 0, 0, 0, 0, 0, 3],
    ];

    // 雀頭のリスト
    let heads = vec![
        vec![2, 0, 0, 0, 0, 0, 0, 0, 0],
        vec![0, 2, 0, 0, 0, 0, 0, 0, 0],
        vec![0, 0, 2, 0, 0, 0, 0, 0, 0],
        vec![0, 0, 0, 2, 0, 0, 0, 0, 0],
        vec![0, 0, 0, 0, 2, 0, 0, 0, 0],
        vec![0, 0, 0, 0, 0, 2, 0, 0, 0],
        vec![0, 0, 0, 0, 0, 0, 2, 0, 0],
        vec![0, 0, 0, 0, 0, 0, 0, 2, 0],
        vec![0, 0, 0, 0, 0, 0, 0, 0, 2],
    ];

    let mut ret: BTreeSet<Vec<i32>> = BTreeSet::new();
    
    for s1 in sets.iter() {
        for s2 in sets.iter() {
            for s3 in sets.iter() {
                for s4 in sets.iter() {
                    for head in heads.iter() {
                        let hand: Vec<i32> = (0..9).map(|i| {
                            s1[i] + s2[i] + s3[i] + s4[i] + head[i]
                        }).collect();
                        if hand.iter().all(|&h| h <= 4) {
                            ret.insert(hand);
                        }
                    }
                }
            }
        }
    }

    ret
}

// 01BFSを用いてシャンテン数を計算
fn bfs(ws: BTreeSet<Vec<i32>>) -> BTreeMap<Vec<i32>, i32> {
    let mut dist: BTreeMap<Vec<i32>, i32> = BTreeMap::new();
    let mut deq: VecDeque<Vec<i32>> = VecDeque::new();

    for hand in ws {
        dist.insert(hand.clone(), 0);
        deq.push_back(hand);
    }

    while !deq.is_empty() {
        let hand = deq.pop_front().unwrap();
        let sum: i32 = hand.iter().sum();
        for k in 0..9 {
            if hand[k] < 4 && sum < 14 {
                let mut hand_add = hand.clone();
                hand_add[k] += 1;
                if !dist.contains_key(&hand_add) {
                    dist.insert(hand_add.clone(), dist[&hand]);
                    deq.push_front(hand_add);
                }
            }
            if hand[k] > 0 {
                let mut hand_sub = hand.clone();
                hand_sub[k] -= 1;
                if !dist.contains_key(&hand_sub) {
                    dist.insert(hand_sub.clone(), dist[&hand] + 1);
                    deq.push_back(hand_sub);
                }
            }
        }
    }

    dist
}

fn main() {
    let start = Instant::now();
    let ws = complete_hands();
    let shanten = bfs(ws);
    assert_eq!(shanten.len(), 405350);

    // 計算結果を保存
    let mut f = fs::File::create("shanten-rs.txt").unwrap();
    for (hand, sh) in shanten {
        f.write_all(format!("{} {} {} {} {} {} {} {} {} {}\n",
            hand[0], hand[1], hand[2], hand[3], hand[4],
            hand[5], hand[6], hand[7], hand[8], sh
        ).as_bytes()).unwrap();
    }
    let elapsed_time = start.elapsed().as_nanos() as f64 
        / 1_000_000_000 as f64;
    println!("elapsed_time: {} [sec]", elapsed_time);
}

また、こちらのリポジトリの実装と結果が一致することを確かめた。

https://github.com/tomohxx/shanten-number-calculator

参考

[麻雀]シャンテン数計算アルゴリズム - Qiita

シャンテン数計算アルゴリズム

数列比較のRolling Hashは危険だよという話

概要

  • 2^{61}-1 mod Rolling Hash が任意の基数で落ちる(数列比較の)問題が作れた。
  • 基数を2つとってハッシュのペアで比較するか、もっと大きいmodを使おう。

Rolling Hashとは

  • 文字列や整数列に対するハッシュ関数の一つ。以下では整数列のみを扱う。
  • M と基数 B を用いたときの、数列 s=s_0, \dots, s_{L-1} に対するRolling Hash H(s; M, B)は次で定義される。

\displaystyle H(s; M, B) := \sum_{i=0}^{L-1} B^{L-1-i}s_i \bmod M

衝突確率

  • 異なる数列が同じハッシュ値をとることを、ハッシュの衝突と呼ぶ。
  • M が十分大きな素数のとき、長さ L の数列 a, b のハッシュが衝突するような基数 B は高々 L-1 個しかない。

    証明

    
H(a; M, B) = H(b; M, B)
\iff \sum_{i=0}^{L-1} B^{L-1-i}(a_i - b_i) = 0 \pmod M

    ここで法 M素数であることから剰余類 \mathbb{Z}/M\mathbb{Z} は整域である。

    整域上の L-1多項式の根は高々 L-1 個であることから、結論を得る。

  • 互いに異なる数列 s^{0}, \dots, s^{n-1} のどこかで衝突が起こるような基数 B は高々 {}_nC_2 (L-1) 個しかない。

  • M=2^{61}-1 \approx 2.3 \times 10^{18}, n\approx 2\times10^{6}, L\approx 2\times10^{6} とした場合、どんな基数 B \in \{0, M-1\} を取っても衝突してしまうケースがあるのではないか? \implies そのような問題が構築できた。

問題の構築

イデア

  •  x^{L-1} - r^{u(L-1)} = 0 \pmod Mx = r^{u}, r^{d+u}, \dots, r^{(L-2)d+u} \bmod M を根に持つ。

    •  L-1M-1 の約数のときに限る。
    • d := (M-1) / (L-1) とおいた。
    • rM の原始根とした。
    • u は任意の整数。
      証明

      k=0,\dots, L-2 に対して 
(r^{kd+u})^{L-1} = (r^{M-1})^{k}r^{u(L-1)} = 1^{k}r^{u(L-1)} = r^{u(L-1)} \pmod M

  • B = r^{in+j+kd} \bmod M, k=0, \dots, L-2 のとき、H(a^{i}; M, B) = H(b^{j}; M, B) となる(衝突する)。

    • a^{i} := [r^{-in(L-1)} \bmod M, 0, \dots, 0], \ 0 \le i \lt n とおいた。
    • b^{j} := [0, \dots, 0, r^{j(L-1)} \bmod M], \ 0 \le j \lt n とおいた。
      証明

      
H(a^{i}; M, B) = H(b^{j}; M, B)
\iff B^{L-1} - r^{(in+j)(L-1)} = 0 \pmod M

      in+j=u とおけば上の議論により直ちに結論を得る。

n^{2} \ge d となるように n をとれば、in+j0, \dots, d-1 をカバーすることで B=r^{0}, r^{1}, \dots, r^{M-2} \bmod M でハッシュの衝突が生じる。r は原始根なので、これは任意の基数 B で衝突が生じることを意味する。

仕上げ

  • M=2^{61}-1 \approx 2.3\times10^{18} のとき、 L \approx 2\times10^{6}, n \approx 2\times10^{6} 程度にすれば n^{2} \ge d := (M-1)/(L-1) とできる。
  • a^{i}, b^{j} をそのまま与えると、 入力サイズが 4\times10^{12} 程度になってしまう。
  • a^{i}, b^{j} はほとんど 0 で埋められた数列なので、情報量は小さい。
  • クエリの問題にしてみる。各クエリで数列が編集されて、途中で出てきた数列の種類数を答えさせる。

完成した問題

  • Input: 全て整数
L, Q
k_1 v_1
...
k_Q v_Q
  • 制約

    • 0 \le L \lt 2\times10^{6}
    • 0 \le Q \lt 3\times10^{6}
    •  0 \le k_t \lt L, \quad t = 1, \dots, Q
    •  0 \le v_t \lt 2^{61}-1, \quad t = 1, \dots, Q
  • Output

    • 長さL の数列 s^{0} := [0, \dots, 0] を定義する。
    • 数列 s^{t-1}k_t 番目の要素を v_t で置き換えた数列を s^{t} とおく。
    • 数列 s^{0}, \dots, s^{Q} は全部で何種類あるかを出力せよ。

2^{61}-1 mod Rolling Hash を落とすケース

M = 2^{61}-1 の原始根 r := 37 をとる。

  • L := 1998910
  • Q := 2n + 1 := 2 \times 1074035 + 1
  • k_t := 0, v_t := r^{-(t-1)n(L-1)} \bmod M, \quad t=1, \dots, n
  • k_t := 0, v_t := 0, \quad t = n+1
  • k_t := L-1, v_t := r^{(t-n-2)(L-1)} \bmod M, \quad t=n+2, \dots, 2n+1

s^{1}, \dots, s^{n}a^{0}, \dots, a^{n-1} に、 s^{n+2}, \dots, s^{2n+1}b^{0}, \dots, b^{n-1} に対応している。

参考

安全で爆速なRollingHashの話 - Qiita

http://hos.ac/blog/#blog0003

Rolling Hashを殺す話

ハッシュの衝突 - あなたは嘘つきですかと聞かれたら「YES」と答えるブログ

clion-live-templates-generator をリリースした

TL;DR

競技プログラミングのライブラリから、CLionのスニペット機能である「ライブラリテンプレート」の設定ファイルを生成してくれるパッケージを公開しました。 github.com

使い方

ローカルで使う

既にライブラリテンプレートを設定している場合は、上書きする前にバックアップを取っておいてください。

$ pip install clion-live-templates-generator
$ lt-generate -d <YOUR_LIBRARY_DIR>
$ cp -i C_C__.xml ~/Library/Application\ Support/JetBrains/CLion2020.1/templates/C_C__.xml

最後の行は、OSによって移動先のフォルダを変えてください。詳細はリポジトリのREADME.mdで確認してください。

Github Actions で自動生成する

Workflow file の一例。

name: CI

on: push

jobs:
  verify:
    runs-on: ubuntu-latest

    steps:
    - uses: actions/checkout@v1

    - name: Set up Python
      uses: actions/setup-python@v1

    - name: Install dependencies
      run: pip3 install -U clion-live-templates-generator
    - name: Create Live-templates file
      run: lt-generate -d lib

    - name: Upload Live-templates file
      uses: actions/upload-artifact@master
      with:
        name: C_C__.xml
        path: C_C__.xml
      if: always()

使ってみた

自分のライブラリで試してみます。

GitHub - habara-k/procon-library: 競技プログラミング用のライブラリ

  • ファイル生成
$ pip install clion-live-templates-generator
$ git clone https://github.com/habara-k/procon-library.git
$ lt-generate -d procon-library
$ cp -i C_C__.xml ~/Library/Application\ Support/JetBrains/CLion2020.1/templates/C_C__.xml
  • CLionを再起動
  • cmd-j

f:id:habara_k:20200719200921p:plain

使えるようになりました。

高速フーリエ変換を図で理解する

問題

(n-1)多項式 f(x) := a_0 + a_1x + \ldots + a _ {n-1}x ^{n-1} について、 x = \zeta_n^ i \ (i = 0, 1, \ldots, n-1) での値を計算したい。

(※ ただし\zeta_n := \exp\left(\frac{2\pi\sqrt{-1}}{n}\right) とおいた。)

愚直に計算するとO(n^ 2)かかるが、これをO(n\log n)で計算する方法を考える。

アルゴリズム

簡単のため、以降 n2の冪乗とする。

(n/2-1)多項式 f_0(x), \ f_1(x) を以下で定義する。

  • f_0(x) := a_0 + a_2x + \ldots + a _ {n-2}x^ {n/2 - 1}
  • f_1(x) := a_1 + a_3x + \ldots + a _ {n-1}x^ {n/2 - 1}

すると、f(x) = f_0(x^ 2) + x \ f_1(x^ 2) が成り立つ。

x = \zeta_n^ i \ (0 \leq i \lt n) での f(x)の値を求めたいので、

0 \leq i \lt n について f_0(\zeta_n^ {2i}), \ f_1(\zeta_n^ {2i})の値が求まっていれば良い。

実は

  • \zeta_n^ {2i} = \zeta_ {n/2}^ i \ \ (0 \leq i \lt n)
  • \zeta _ {n/2}^ i = \zeta _ {n/2}^ {i + n/2} \ \ (0 \leq i \lt n/2)

より、\zeta_n^ {2i} = \zeta _ {n/2}^ {i \ \mathrm{mod} \ n/2} \ \ (0 \leq i \lt n) であることがわかる。

つまり 0 \leq i \lt n/2 について f_0(\zeta _ {n/2}^ i), \ f_1(\zeta _ {n/2}^ i) が求まっていれば良い。

(※ 以下にn=8の場合の図を示す。x = \zeta_n^ {2i} は添字が2づつ進んで2周するため、●における f_0(x), \ f_1(x)の値が求まっていれば良いことがわかる。) f:id:habara_k:20200520234237p:plain

したがって、以下の再帰的な計算で (n-1)多項式 f(x) = \sum_{i=0}^ {n-1} a_i x^ i に対する f(\zeta_n^ i) \ (0 \leq i \lt n) を求めることができる。

  • n = 1のとき、f(x) = a_0 より f(\zeta_1^ 0) = a_0を返す。
  • そうでないとき、

    1. f_0(x) := \sum _ {i=0}^ {n/2-1} a _ {2i} x^ i, \ f_1(x) := \sum _ {i=0}^ {n/2-1} a _ {2i+1} x^ i を定義する。

    2. (n/2-1)多項式 f_0(x), \ f_1(x) に対する f_0(\zeta _ {n/2}^ i), \ f_1(\zeta _ {n/2}^ i) \ \ (0 \leq i \lt n/2) をそれぞれ再帰的に計算する。

    3. 0 \leq i \lt n について f(\zeta_n^ i) = f_0(\zeta _ {n/2}^ {i \ \mathrm{mod} \ (n/2)}) + \zeta _ n^ i \ f_1(\zeta _ {n/2}^ {i \ \mathrm{mod} \ (n/2)}) を計算して返す。

各呼び出しで扱う多項式は、n=8の場合以下の図のようになる。 ただし、f _ {*} から定義した2つの多項式をそれぞれ f _ {0*}, \ f _ {1*} と表記した。 f:id:habara_k:20200520182929p:plain 再帰の終端において、f _ {*}(x) = a _ {*} になっていることがわかる。 次節の実装ではこれを利用している。

計算量はマージソートと同様の解析で O(n\log n)となる (分割統治法)。

実装

using Complex = std::complex<double>;

struct FFT {

    std::vector<Complex> a_;
    int n;

    FFT(const std::vector<Complex>& a) : a_(a), n(1) {
        // n を2の冪乗にする
        while (n < a.size()) n <<= 1;
        a_.resize(n);
    };

    std::vector<Complex> solve() {
        return fft(0, 0);
    }

    std::vector<Complex> fft(int d, int bit) {
        int sz = n >> d;
        if (sz == 1) return {a_[bit]};

        auto f0 = fft(d+1, bit);
        auto f1 = fft(d+1, bit | 1<<d);

        std::vector<Complex> f(sz);

        for (int i = 0; i < sz; ++i) {
            Complex x = std::polar(1.0, 2*M_PI / sz * i);
            f[i] = f0[i % (sz / 2)] + x * f1[i % (sz / 2)];
        }
        return f;
    }
};

参考

http://compro.tsutajiro.com/archive/fft.pdf

C - 高速フーリエ変換

MacBookのJISキーボードをUSに設定する

概要

MacBookのJIS配列のキーボードをUS配列に設定します.

具体的には,

  • 英数入力時, 記号だけUS配列に対応させる.
  • `~ は対応するキーがないので, deleteの一つ左にあるキーに割り当てる.

環境

手順

1. 入力ソースの設定

システム環境設定 -> キーボード -> 入力ソース に移動して,

  • 日本語の平仮名入力以外をOFFにする.

  • オランダを追加する

入力ソースの設定
入力ソースの設定

これで英数入力時だけ, 記号がUS配列になることを確認してください. このままだと `~ が入力できないので次に進みます.

2. Karabiner-Elements を入れる.

Karabiner-Elements をダウンロード・インストールしてください.

3. `~ を割り当てる.

Karabiner-ElementsのSimple modification を開いて,

  • Add item を押す.

  • From key に international3 を, To key に grave_accent_and_tilde (`) を設定する.

キーの割り当て
キーの割り当て

これでdeleteの一つ左のキーに`~ が割りあてられていることを確認すれば完成です.

妥協点

日本語入力のときにもKarabinerのキーマッピングが有効になってしまうので, 全角の \| が打てない.

AtCoder頻出?貼るだけ高速ゼータ・メビウス変換&約数拡張

定義

有限集合 X, (可換な)半群 G について, 写像 f:2^ X\rightarrow G , g:2^ X\rightarrow Gg(S) = \sum_{S \subseteq T} f(T) を満たすときを考える.

(※ g(S) = \sum_{T \subseteq S} f(T) の形をしたものもあるが, \bar{g}(S) := g(2^ X - S)\ , \bar{f}(T) := f(2^ X - T) と定義すれば前者と同じ形になるので割愛する.)

高速ゼータ変換

f から g を求める. O(|X|2^ {|X|})

template<typename T>
vector<T> fast_zeta_transform(const vector<T>& f) {
    int n = f.size();
    vector<T> g = f;
    for (int i = 0; (1 << i) < n; ++i) {
        for (int j = 0; j < n; ++j) {
            if (!(j >> i & 1)) {
                g[j] += g[j | (1<<i)];
            }
        }
    }
    return g;
}

例題

atcoder.jp

G := (\mathbb{N},\mathbb{N}),

 (a, b), (c, d) \in G, (a, b) + (c, d) := ((a,b,c,d) の中で一番大きな数, 二番目に大きな数)

と定義する.

f(S) := (A_S, 0)

g(S) := \max \{(A_i, A_j) | i \neq j, i, j \subseteq S \}

と定めると, g(S) = \sum_{T \subseteq S} f(T) が成立する.

包含関係が逆になってることに注意しつつ f から g を求めて, あとは S を1ずつ増やしながら\maxを更新していく.

Submission #11912671 - AtCoder Regular Contest 100

高速メビウス変換

G に逆元があるとき, g から f を求める. O(|X|2^ {|X|})

template<typename T>
vector<T> fast_mobius_transform(const vector<T>& g) {
    int n = g.size();
    vector<T> f = g;
    for (int i = 0; (1 << i) < n; ++i) {
        for (int j = 0; j < n; ++j) {
            if (!(j >> i & 1)) {
                f[j] -= f[j | (1<<i)];
            }
        }
    }
    return f;
}

約数に拡張した高速メビウス変換

自然数の集合 M = \{1, 2, \ldots, N\} , (可換な)半群 G について, 写像 f:M \rightarrow G , g:M \rightarrow Gg(d) = \sum_{d | n, n \leq N} f(n) を満たすとき, g から f を求める. O(N\log\log N)

template<typename T>
vector<T> fast_mobius_transform_prime(const vector<T>& g) {
    int n = g.size();
    vector<T> f = g;
    vector<bool> sieve(n, true);
    for (int p = 2; p < n; ++p) {
        if (!sieve[p]) continue;
        for (int i = 1; i * p < n; ++i) {
            sieve[i * p] = false;
            f[i] -= f[i * p];
        }
    }
    return f;
}

例題

atcoder.jp

M := \{1, \ldots , K\}

f(d) := \mathrm{gcd}(A_1, \ldots, A_N)d になる数列の総数

g(d) := \mathrm{gcd}(A_1, \ldots, A_N)d の倍数になる数列の総数

とおくと g(d) = \sum_{d | n, n \leq K} f(n) が成立する. g(d) = \left(\left\lfloor\frac{K}{d}\right\rfloor\right)^ N から f(d) を求めて, \sum_d f(d) \cdot d を計算する.

Submission #11899603 - AtCoder Beginner Contest 162

atcoder.jp

M := \{1, \ldots , \max(A)\}

f(d) := \mathrm{gcd}(A_i, A_j)d になる組の総数

g(d) := \mathrm{gcd}(A_i, A_j)d の倍数になる組の総数

とおくと g(d) = \sum_{d | n, n \leq \max(A)} f(n) が成立する.

C(a) := A_1, \ldots, A_N における a の出現数

を計算しておく.

A_1, \ldots, A_N の中で d の倍数であるものを B_1, \ldots, B_k とおくと


g(d) = \sum_{i < j} B_i B_j = \frac{1}{2} \left\{ \left( \sum_{l=1}^{k}B_l \right)^ 2 - \sum_{l=1}^{k}B_l ^ 2 \right\}

であり, C を用いることで g(d)O(\max(A)\log(\max(A))) で計算できる.

あとは g(d) から f(d) を求めて, \sum_d \frac{f(d)}{d} を計算する.

Submission #11899671 - AtCoder Grand Contest 038

参考

ARC100-E 「Or Plus Max」 (700) 高速ゼータ変換解法 - Mister雑記

高速ゼータ・メビウス変換 - Mister雑記

高速ゼータ変換の約数版 - noshi91のメモ

ゼータ変換・メビウス変換を理解する - Qiita