底とmodが互いに素でない離散対数問題
離散対数問題とは
整数 に対して, を満たす最小の非負整数 を求める問題です.
例えば のとき,
より が解となります.
Baby-step Giant-stepアルゴリズム
底 と mod が互いに素なときに使える離散対数問題の 解法です.
とします. とおくと, となります.
したがって を hash table に保存しておいて, の順に存在判定をすれば良いです.
底とmodが互いに素でない場合
の素因数分解を とおきます. ただし は の素因数であり, は の素因数ではないとします.
このとき, に対して は の倍数になります. したがって,
(1). が の倍数でないとき
の範囲で の解は存在しません.
(2). が の倍数であるとき
に対して
が成立します. と は互いに素なので, これは Baby-step Giant-step アルゴリズムで解くことができます.
全体の流れ
- については愚直に かどうかを判定する. 解が見つかれば終了.
- が の倍数でないなら, 解なし.
- を 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となる。
- この例のようにシャンテン数がわかりやすいケースがほとんどだが、清一色のような場合はシャンテン数を求めにくい。
- 特にプログラムでシャンテン数を計算する場合、正しいコードを書くことは難しい。また、正当性を保証するのも難しい。
考える問題
清一色だけを考える。以下では萬子🀇🀈🀉🀊🀋🀌🀍🀎🀏だけを扱う。
4面子1雀頭の一般手だけを考える。
副露は考えない。
- x回副露している場合は、4-x面子1雀頭を完成形とみなしてシャンテン数を求めれば良く、ほぼ同様の計算となるため。
Notation
- 手牌を で表す。
- 萬子 が何枚あるか。
- 次の条件を満たす手牌全体の集合を とおく。
- 手牌 に対して、4面子1雀頭を作るのに必要な牌の最小枚数を とおく。
- 厳密には、4面子1雀頭を余りなく作れる手牌(以下完成形と呼ぶ)の集合を とおいて、 を次で定義する。
手牌 のシャンテン数は となるので、 を計算すれば良い。
の計算
アイデア
重み付き有向グラフを以下で定める。
手牌全体の集合 を頂点集合とする。
次のルールで辺を張る。
手牌 から、1枚牌を追加して得られる手牌に距離0の辺を張る。
- 厳密には、各 に対して とおき、 ならば から に距離0 の辺を張る。
手牌 から、1枚牌を削除して得られる手牌に距離1の辺を張る。
- 厳密には、各 に対して とおき、 ならば から に距離1 の辺を張る。
このとき、完成形 から の最短距離を とおくと、次の関係が成立する。
式 より、 を示せば良い。 これは辺の張り方から簡単に従う。( の牌を追加・削除して に一致させていく様子をイメージするとわかる。)証明
完成形 の列挙
- 面子は 7 + 9 通り
- 雀頭は 9 通り
であるから、 通りを考えて、その各牌の枚数 を計算する。
に対して を満たすものを完成形として列挙する。
の計算
辺の距離が0と1しかないので、01BFSを用いて計算することができる。
頂点数と辺の数を計算すると
- 頂点数 = 405350
- 辺の数 =
36481504791600
となる。 で を計算することができる。
実装とテスト
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
参考
数列比較のRolling Hashは危険だよという話
概要
- mod Rolling Hash が任意の基数で落ちる(数列比較の)問題が作れた。
- 基数を2つとってハッシュのペアで比較するか、もっと大きいmodを使おう。
Rolling Hashとは
- 文字列や整数列に対するハッシュ関数の一つ。以下では整数列のみを扱う。
- 法 と基数 を用いたときの、数列 に対するRolling Hash は次で定義される。
衝突確率
- 異なる数列が同じハッシュ値をとることを、ハッシュの衝突と呼ぶ。
法 が十分大きな素数のとき、長さ の数列 のハッシュが衝突するような基数 は高々 個しかない。
互いに異なる数列 のどこかで衝突が起こるような基数 は高々 個しかない。
法 とした場合、どんな基数 を取っても衝突してしまうケースがあるのではないか? そのような問題が構築できた。
問題の構築
アイデア
は を根に持つ。
- が の約数のときに限る。
- とおいた。
- を の原始根とした。
- は任意の整数。
証明
に対して
のとき、 となる(衝突する)。
- とおいた。
- とおいた。
証明
とおけば上の議論により直ちに結論を得る。
となるように をとれば、 が をカバーすることで でハッシュの衝突が生じる。 は原始根なので、これは任意の基数 で衝突が生じることを意味する。
仕上げ
- のとき、 程度にすれば とできる。
- をそのまま与えると、 入力サイズが 程度になってしまう。
- はほとんど で埋められた数列なので、情報量は小さい。
- クエリの問題にしてみる。各クエリで数列が編集されて、途中で出てきた数列の種類数を答えさせる。
完成した問題
- Input: 全て整数
L, Q k_1 v_1 ... k_Q v_Q
制約
Output
- 長さ の数列 ] を定義する。
- 数列 の 番目の要素を で置き換えた数列を とおく。
- 数列 は全部で何種類あるかを出力せよ。
mod Rolling Hash を落とすケース
の原始根 をとる。
が に、 が に対応している。
参考
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
使えるようになりました。
高速フーリエ変換を図で理解する
問題
次多項式 について、 での値を計算したい。
(※ ただし とおいた。)
愚直に計算するとかかるが、これをで計算する方法を考える。
アルゴリズム
簡単のため、以降 をの冪乗とする。
次多項式 を以下で定義する。
すると、 が成り立つ。
今 での の値を求めたいので、
について の値が求まっていれば良い。
実は
より、 であることがわかる。
つまり について が求まっていれば良い。
(※ 以下にの場合の図を示す。 は添字が2づつ進んで2周するため、●における の値が求まっていれば良いことがわかる。)
したがって、以下の再帰的な計算で 次多項式 に対する を求めることができる。
各呼び出しで扱う多項式は、の場合以下の図のようになる。 ただし、 から定義した2つの多項式をそれぞれ と表記した。 再帰の終端において、 になっていることがわかる。 次節の実装ではこれを利用している。
計算量はマージソートと同様の解析で となる (分割統治法)。
実装
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; } };
参考
MacBookのJISキーボードをUSに設定する
概要
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のキーマッピングが有効になってしまうので, 全角の \
と |
が打てない.
AtCoder頻出?貼るだけ高速ゼータ・メビウス変換&約数拡張
定義
有限集合 , (可換な)半群 について, 写像 , が を満たすときを考える.
(※ の形をしたものもあるが, と定義すれば前者と同じ形になるので割愛する.)
高速ゼータ変換
から を求める.
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; }
例題
の中で一番大きな数, 二番目に大きな数
と定義する.
と定めると, が成立する.
包含関係が逆になってることに注意しつつ から を求めて, あとは を1ずつ増やしながらを更新していく.
Submission #11912671 - AtCoder Regular Contest 100
高速メビウス変換
に逆元があるとき, から を求める.
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; }
約数に拡張した高速メビウス変換
自然数の集合 , (可換な)半群 について, 写像 , が を満たすとき, から を求める.
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; }
例題
が になる数列の総数
が の倍数になる数列の総数
とおくと が成立する. から を求めて, を計算する.
Submission #11899603 - AtCoder Beginner Contest 162
が になる組の総数
が の倍数になる組の総数
とおくと が成立する.
における の出現数
を計算しておく.
の中で の倍数であるものを とおくと
であり, を用いることで を で計算できる.
あとは から を求めて, を計算する.
Submission #11899671 - AtCoder Grand Contest 038