離散対数問題とは

整数 $p \ge 1, \ 0 \le a, b \lt p$ に対して, $a^{x} \equiv b \mod p$ を満たす最小の非負整数 $x$ を求める問題です.

例えば $p=40, \ a = 2, \ b = 24$ のとき,

$a^{0} \equiv 1$
$a^{1} \equiv 2$
$a^{2} \equiv 4$
$a^{3} \equiv 8$
$a^{4} \equiv 16$
$a^{5} \equiv 32$
$a^{6} \equiv 24$

より $x=6$ が解となります.

Baby-step Giant-stepアルゴリズム

底 $a$ と mod $p$ が互いに素なときに使える離散対数問題の $O(\sqrt{p})$ 解法です.

$h := \lceil \sqrt{p} \rceil$ とします. $x = ih - j \ (0 \lt i, j \le h)$ とおくと, $a^{x} \equiv b \iff a^{ih} \equiv a^{j}b$ となります.

したがって $a^{1}b, \dots, a^{h}b$ を hash table に保存しておいて, $a^{h}, a^{2h}, \dots, a^{h^{2}}$ の順に存在判定をすれば良いです.

底とmodが互いに素でない場合

$p$ の素因数分解を $p = p_1^{d_1} \cdots p_m^{d_m} \ q_1^{e_1} \cdots q_n^{e_n}$ とおきます. ただし $p_1, \dots, p_m$ は $a$ の素因数であり, $q_1, \dots, q_n$ は $a$ の素因数ではないとします.

このとき, $x \ge \lfloor \log_2 p \rfloor$ に対して $a^{x}$ は $g := p_1^{d_1} \dots p_m^{d_m}$ の倍数になります. したがって,

(1). $b$ が $g$ の倍数でないとき

$x \ge \lfloor \log_2 p \rfloor$ の範囲で $a^{x} \equiv b \mod p$ の解は存在しません.

(2). $b$ が $g$ の倍数であるとき

$x \ge \lfloor \log_2 p \rfloor$ に対して

$a^{x} \equiv b \mod p \iff a^{x} / g \equiv b / g \mod p / g$

が成立します. $a$ と $p/g = q_1^{e_1} \cdots q_n^{e_n}$ は互いに素なので, これは Baby-step Giant-step アルゴリズムで解くことができます.

全体の流れ

$0 \le x \lt \lfloor \log_2 p \rfloor$ については愚直に $a^{x} \equiv b \mod p$ かどうかを判定する. 解が見つかれば終了.
$b$ が $g = p_1^{d_1} \dots p_m^{d_m}$ の倍数でないなら, 解なし.
$a^{x} / g \equiv b / g \mod p / g$ を Baby-step Giant-step アルゴリズムで解く.

実装解説

離散対数問題(Mod-Log) | Luzhiled’s memo の実装を, 上記の説明に照らし合わせて解説します.

(※ 元のコードを一部変更しております.)

int64_t mod_log(int64_t a, int64_t b, int64_t p) {
  // g := p_1^{d_1} ... p_m^{d_m} を求める.
  int64_t g = 1;
  for(int64_t i = p; i; i /= 2) (g *= a) %= p;
  g = __gcd(g, p);

  // a^{x} が g の倍数になる手前までは愚直に計算する.
  int64_t t = 1, c = 0;
  for(; t % g; c++) {
    if(t == b) return c;
    (t *= a) %= p;
  }
  // a^{c} = t

  // b が g の倍数でないなら解なし.
  if(b % g) return -1;

  // a^{y} (t / g) = b / g mod p / g を解いて, y + c を解とする.
  t /= g;
  b /= g;
  p /= g;
  // a^{y} t = b mod p を解く.

  int64_t h = 0, gs = 1;

  // h := ceil(sqrt(p))
  // gs := a^{h}
  for(; h * h < p; h++) (gs *= a) %= p;

  // a^{ih-j} t = b mod p (0 < i, j <= h)
  // (a^{h})^{i} t = b a^{j} mod p
  // (gs)^{i} t = b a^{j} mod p

  // b a^{j} (0 < j <= h) を bs に登録.
  unordered_map< int64_t, int64_t > bs;
  for(int64_t s = 0, e = b; s < h; bs[e] = ++s) {
    (e *= a) %= p;
  }

  // (gs)^{i} t (0 < i <= h) で bs を検索. y = s - bs[e]
  for(int64_t s = 0, e = t; s < p;) {
    (e *= gs) %= p;
    s += h;
    if(bs.count(e)) return c + s - bs[e];
  }

  // 見つからなければ解なし.
  return -1;
}

参考

概要

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

シャンテン数とは

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

考える問題

清一色だけを考える。以下では萬子🀇🀈🀉🀊🀋🀌🀍🀎🀏だけを扱う。
- 複数色ある場合や字牌を含む場合でも、清一色のシャンテン数さえ計算できていれば、軽いdpで解くことができるため。
- 正確には、x面子y雀頭を完成形とみなしたときの清一色のシャンテン数を全ての $0 \le x \le 4, \ 0 \le y \le 1$ について求めておく必要があるが、これらはほぼ同様に計算できる。
4面子1雀頭の一般手だけを考える。
- 七対子や国士無双のシャンテン数は、簡単な計算で求めることができるため。
副露は考えない。
- x回副露している場合は、4-x面子1雀頭を完成形とみなしてシャンテン数を求めれば良く、ほぼ同様の計算となるため。

Notation

手牌をで表す。
- $h_i :=$ 萬子 $i$ が何枚あるか。
次の条件を満たす手牌全体の集合をとおく。
- $0 \le h_i \le 4, \quad 1 \le i \le 9$
- $\sum_{i=1}^9 h_i \le 14$
手牌に対して、4面子1雀頭を作るのに必要な牌の最小枚数をとおく。
- 厳密には、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

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

概要

$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 M$ は $x = r^{u}, r^{d+u}, \dots, r^{(L-2)d+u} \bmod M$ を根に持つ。
- $L-1$ が $M-1$ の約数のときに限る。
- $d := (M-1) / (L-1)$ とおいた。
- $r$ を $M$ の原始根とした。
- $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+j$ が $0, \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」と答えるブログ

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)$ で計算する方法を考える。

アルゴリズム

簡単のため、以降 $n$ を $2$ の冪乗とする。

$(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配列に設定します.

具体的には,

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

環境

macOS Catalina 10.15.4

手順

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のキーマッピングが有効になってしまうので, 全角の \ と | が打てない.

定義

有限集合 $X$ , (可換な)半群 $G$ について, 写像 $f:2^ X\rightarrow G$ , $g:2^ X\rightarrow G$ が $g(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 G$ が $g(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$ を計算する.