ホロライブ換字暗号のこと好きすぎ?!

概要

マドロミはかなたそトワさんが歌う楽曲である.この楽曲中には暗号らしきものが含まれており,作詞者の傾向から換字式暗号と推定される.本記事ではマドロミに含まれる暗号文の完全解読を行う.

序論

最近は僕の中でかなたそがアツい.最近は何故かわからないが忙しい気がしているので,頻繁にかなたその配信やら動画やらニュースをチェックできているわけではないのだが,最近面白いものを見つけたのでそれについて考えたいと思う.それがこちら.換字暗号である.

youtu.be

Twitterで「マドロミ 暗号」と検索しても誰のツイートも出てこないので,どうやら誰も暗号に興味が無いようである.いやわかる.かなたそとトワさんが歌ってたらみんな暗号より歌の方が気になる.

せっかく暗号があるのに誰も解かないのはもったいない上,OW-COA安全ですらない解読がお手軽な古典暗号を使ってくれたので解読していこうと思う.

前提知識

マドロミの作詞は 砂守岳央 さんが担当された.彼は実はゼロの足跡も作詞されている.ゼロの足跡は全世界の暗号オタクのわためいとを熱狂させた換字暗号を含んだ角巻わためさんの曲である.ホロライブ関連楽曲だけですでに2曲に換字暗号を含ませた砂守岳央さんはきっと暗号が好きに違いない.ゼロの足跡の暗号については以下を参照してほしい.

ouchiminh.hatenablog.com

また,マドロミの暗号文及び歌詞の全文は法律上このブログには掲載できないので,上の動画を見るなどして得て欲しい.また,マドロミに含まれている暗号文は歌として発音する必要性から子音と母音がペアになっており,さらに子音は子音に,母音は母音へと移すような暗号になっていることは注目に値する.

解読

換字暗号の解読方法の定番と言えば頻度分析だが,マドロミの暗号の解読においては頻度分析は正直役に立たなかった.そんなものよりずっと重要だったのは子音と母音の並び方と閃きであった.

最初に解読できたのは暗号文 fa ~ve ~te ~fu だった.この暗号文は4音あるが,母音も子音も3つしか含まれていなかった.さらに偶然にも fa ~ve~ te~ fu はこの楽曲のタイトルである「マドロミ」に対応する暗号文だった.一気に6文字も解読に成功し,このうち母音は3/5を解読できたので,残った母音 u,~eは消去法でも選べた.復号をi\mapsto eとしてしまうと,暗号文の ko~ fi~ to~ la~ ci ~ta ~wa ~ku に対応する平文が ?u~ ma ~ru ~  ?a ~  ?e ~ra ~ ?a ~  ?i となり,最初の3文字に「詰まる」くらいしか対応する日本語を見つけられなかったため,これを棄却し,i\mapsto u,~o\mapsto eとした.ここに母音の対応を得た.

暗号文 a i u e o
平文 a u i o e

母音が対応すればもう適当に母音だけを読んでいても何かしら思いつく. e~la~ci~fu~wa~ya~u~u~weは「マドロミ」の直後にあり,部分的に解読すると?o~?a~?u~mi~?a~?a~i~i~?oとなる.最初の四文字には「おやすみ」しか入らない気がしてくる.この時点でかなりの部分が解読できているので残りの部分の解読は益々簡単になる.

暗号文la~ci~ta~wa~kuは現時点でya~su~ra~?a~?iと解読されており,どう考えても「安らかに」に対応する.

続けて暗号文se~ya~kuは「永久に」に対応する.

あとは偏見と独断でh\mapsto h,\quad n\mapsto z,\quad b\mapsto gとすれば全て解読が完了する.

解読結果は 「おやすみ わがこよ とわに やすらかに」 「まどろみ おやすみ かわいいこ ねむれ ねむれ ねむれ ひとみをとじて まどろみ おやすみ かわいいこ ねむれ ねむれ ねむれ やすらかに」 となる.

解読に役立つツール

普通なら頻度分析が役に立つが,この暗号に限っては

  1. 暗号文のサンプルが非常に少ない
  2. (おそらく意図的に) 平文と暗号文の連接頻度に強い関連がない

などの理由から頻度分析があまり役に立たなかった.例えば,平文中に最も出てくる2文字の組み合わせ「おち」は暗号文中には一度も出てこない.故に,この暗号の解読に一番役に立ったのはvimだった.vimをただのテキストエディタと侮ることなかれ.正規表現で文字を検索できるので,できることがかなり広がる.この章ではvimがどのように役に立ったかを軽く紹介する.

本暗号で最初に解読できたのは  fa~ve~te~fu の部分だった.これをどのようにして見つけたかを説明するにはvim正規表現が欠かせない.まずは過去の日本語暗号文の解読経験から母音が判明すれば解読が非常にやりやすくなることは分かっていた.そのため,まずは母音の対応を見つけることに専念した.vimで平文歌詞をローマ字に書き下し,vimノーマルモード\cdot  a~\cdot e~\cdot e ~\cdot u に対応する可能性のある以下の文字のパターンを順に入力して対応する母音を探す.

  • /.a.i.i.u
  • /.a.i.i.e
  • /.a.i.i.o

\quad\vdots

  • /.a.o.o.i

\quad\vdots

/.a.o.o.iを入力した時点でローマ字に書き下した平文歌詞から  fa~ve~te~fu と対応に矛盾がないmadoromiを入手できる.ここで母音3文字の対応を得られるのであとは勘と経験に任せよう.

暗号文の異なる部分に着目すれば,全く違う対応を得て惑わされたかもしれない.例えば,暗号文 ye~se~nu~so に着目した場合はこれが平文「わたしと」に対応するように思えてしまう.これは異なる対応なのでいずれ間違いだと気づくが,直ちに誤りだとは分からない.時間がたってから暗号文中の ya を見つけ,それが「わ」でも「を」でもありえないことに絶望しながら対応表を書きなおす羽目になった.

C++で標数2の体

概要

 本記事では標数2の体をC++で定義する際の要点について概説する.異なる拡大による体を型で静的に区別できるような仕組みを作り,C++における自然で美しい体の実装を目指す.

動機

 夏休み中,学友がC++を学ぶために勉強会を開いた.勉強会は様々な数学的対象をC++で表現する内容だった.本記事は勉強会の一回で扱った内容の焼き直しである.

準備

 ある集合Gと演算\circ:G\times G\rightarrow Gが次の1~3の条件を満たすとき,集合Gと演算\circの組(G,\circ)を群という.また,1~4の条件を満たすとき,可換群またはアーベル群という.考えている演算が明らかなとき単にGを群と呼ぶ.

  1. 演算の結合性 : ^\forall a, b, c \in G, (a\circ b)\circ c = a\circ(b\circ c)
  2. 単位元の存在 : ^\forall a \in G, ^\exists e \in G\hspace{5pt} \text{ここで}\hspace{5pt} e\circ a = a \circ e = a
  3. 逆元の存在 : ^\forall a \in G, ^\exists x \in G \hspace{5pt} \text{ここで}\hspace{5pt} a\circ x = x\circ a = e
  4. 交換法則 : ^\forall a,b\in G, a\circ b= b\circ a

 演算\circが加法的である場合,aの逆元x-aと書き,演算\circが乗法的である場合,xa^{-1}と書く.

 ある集合Rと加法+:R\times R\rightarrow Rの組(R,+)が可換群であり,かつ以下の1~3の条件を満たすとき,R,+と乗法\ast:R\times R\rightarrow Rの組(R,+,\ast)を環という.また,1~4の条件を満たすとき,可換環という.

  1. 乗法の結合性 : ^\forall a, b, c\in R, (a\ast b)\ast c = a\ast (b\ast c)
  2. 乗法単位元の存在 : ^\forall a, \in R, ^\exists e _ \ast \in R \,\text{ここで}\, a \ast e _ \ast = e _ \ast \ast a = a
  3. 分配律 : ^\forall a, b, c \in R, (a+b)\ast c=(a\ast c)+(b\ast c)かつ a\ast(b+c)=(a\ast b)+(a\ast c)
  4. 交換法則 : ^\forall a,b\in R, a\ast b = b\ast a

 環Rの乗法単位元e _ \astについて,ne _ \ast=e _ \ast + e _ \ast + ... + e _ \ast = e _ \astとなるn>0が存在するとき,その最小の値をR標数とする.また,そのようなnが存在しないとき,R標数0とする.

イデアル

 環の部分集合Iが以下の1~3の条件を満たすときIイデアルという.1~4の条件を満たすときIを素イデアルという.

  1. 0\in I
  2. ^\forall a\in I, ^\forall r \in R, a\ast r \in I
  3. ^\forall a, b \in I, a + b\in I
  4. ^\forall r,s\in R, r\ast s \in I \text{.このとき} r\in I \text{または} s\in I

 例えば任意の整数の倍数全体の集合は,整数全体と通常の加法,乗法からなる環のイデアルとなっている.また,任意の素数の倍数全体の集合は素イデアルとなっている.可換環RイデアルIを求めることができれば,RIで割った余りからなる新たな環を定義できる.^\forall a,b\in Rについて,a-b\in Iならばa\sim bとする.この同値関係による同値類を\lbrack x\rbrack = x+Iと書くと加法 : \lbrack x \rbrack + \lbrack y \rbrack = \lbrack x + y \rbrackと乗法 : \lbrack x \rbrack \ast \lbrack y \rbrack = \lbrack x \ast y \rbrackによって剰余環R/Iが定義できる.

多項式環

 可換環や体を係数とする多項式環を定義する.Xを変数とし,可換環Rを係数とする一変数多項式環R\lbrack X \rbrackと書く.多項式 F=a_0 + a_1 X + a_2 X^{2} + ... + a_d X^{d} \in R\lbrack X \rbrackの各項のXの次数のうち最高の数を\deg F = dとし,Fの次数がdである,またはFd多項式であるという.F\deg F次の項の係数が1である多項式をモニック多項式であるという.定数でない多項式Fが,定数でない多項式の二つの積で表せないときF既約多項式という.

 ある集合\mathbb{F}可換環であり,かつ\mathbb{F}\backslash\lbrace 0 \rbraceが乗法に関して可換群をなすとき,\mathbb{F}を体という.つまり体は四則演算を自由に行える代数系である.環の素イデアルで割った剰余環は体となる.

体の代数拡大

 体の拡大とは,体\mathbb{F}\mathbb{K}\supset\mathbb{F}の組\mathbb{K}/\mathbb{F}のことをいう.\alpha \in \mathbb{K}F\in\mathbb{F}[X]の根であるとき,\alpha\mathbb{F}上代数的であるという.また,\mathbb{F}の拡大体\mathbb{K}の任意の元が代数的であるとき,\mathbb{K}/\mathbb{F}を代数拡大という.代数拡大の例に\mathbb{C}/\mathbb{R}がある.\mathbb{C}X ^ 2 + 1 \in \mathbb{R}[X]の根\sqrt{-1}\mathbb{R}に追加した体である.また,代数拡大でない拡大の例に\mathbb{R}/\mathbb{Q}がある.\pi,\mathrm{e} \in \mathbb{R}などは\mathbb{Q}上代数的でない.

 \mathbb{F}[X]の既約多項式の根\alpha\mathbb{F}に加えると\mathbb{F}を拡大できる.そのように拡大した体\mathbb{K}\supset\mathbb{F}の任意の元x\alphaa _ i \in \mathbb{F}を用いてx= a _ 0 + a _ 1\alpha + a _ 2 \alpha ^ 2 ... + a _ {\deg F - 1} \alpha ^ {\deg F - 1}と表され,\mathbb{K}\mathbb{F}のベクトル空間となる.\mathbb{F}のベクトル空間としての\mathbb{K}の次元数を拡大次数という.

標数2の体のビット表現

 整数環\mathbb{Z}の部分集合2\mathbb{Z}=\{...,-4,-2,0,2,4,...\}\mathbb{Z}の素イデアルである.剰余環\mathbb{Z}/2\mathbb{Z}標数2の体\mathbb{F} _ {2}になる.\mathbb{F} _ {2}の拡大体\mathbb{F} _ {2 ^ d}の元xは,\mathbb{F} _ {2} [X]上既約なd多項式の根\alphab _ i = 0,1 \in \mathbb{F} _ 2によって以下のように表される.

 \displaystyle
x=b _ 0 + b _ 1 \alpha + b _ 2 \alpha ^ 2 + ... + b _ {d-1} \alpha ^{d - 1}

 ここで上式のb _ idビットの整数の各ビットに対応させると\mathbb{F} _ {2 ^ d}の元を2進数として表現できる.また,\alphaを変数に置き換えれば\mathbb{F} _ 2 [X]多項式も2進数として表現できる.

実装

型の定義

 型は値に意味を与える.あるビット表現が単なる整数なのか,あるいは標数2の体の元を表現しているのか,あるいは標数2の体の元を表現しているとしてどの体の元なのかを適切に表現し計算するために,異なる体には異なる型を与える必要がある.C++ではテンプレート引数に整数引数を取れるため,体を拡大する原始多項式のビット表現をテンプレート引数で与えることができる.本記事では簡単のために拡大次数がd=8, 16, 32, 64\mathbb{F} _ {2 ^ d}のみ扱うこととするが,任意次数の拡大であっても理論は同じである.

 原始多項式のビット表現Polynomialをテンプレート引数で受け取れば,異なる原始多項式による異なる拡大体を異なる型で表現することができる.また,原始多項式のビット表現では最高次数の項がビット数によって常に一定なので省略されている.例えばT = unsigned charなら,Polynomialで直接表現されているのはX ^ 7以下の項のみで,X ^ 8の項は常に係数が1である.

template<class T, std::enable_if_t<std::is_unsigned_v<T>, T> Polynomial>
class gf {
    T value_;
};

 ここに元のビット表現を格納するメンバ変数を定義し,基本的な演算を行う演算子を追加していく.

加法

 \mathbb{F} _ {2 ^ d}の加法は\mathbb{F} _ 2と同様に排他的論理和演算で表すことができる.以下の表のように\mathbb{F} _ 2の加算をビット毎 (=元の式の項ごと) に行えばよい.加法の逆演算については-1=1\in\mathbb{F} _ 2, -0 = 0 \in\mathbb{F} _ 2であるため,加法と全く同じ演算で逆演算を行える.

+ 0 1
0 0 1
1 1 0

乗法

 \mathbb{F} _ {2 ^ d}の乗法は\mathbb{F} _ 2[X]Xを原始元に置き換えた乗法と変わらないが,元の多項式基底表現の次数がdを超える元は原始多項式で割って余りをとらなければコンピューターで表現できる桁数を超えてしまう.乗法の最も素朴な実装は,筆算をするように次数の低い項から掛けていき,ビット表現が型でサポートされる最大値を上回るとき,すなわち計算途中で\alpha ^ dを加算するとき,その代わりに,原始多項式F(\alpha) = 0の最高次数の項を移行した残りの辺のビット表現を加算するものである.例えばx=\alpha, y=\alpha + 1\in \mathbb{F} _ {2 ^ 2}の乗法x\ast yを考える.ただし,原始多項式F = X ^ 2 + X + 1とする.x\ast y = \alpha \ast (\alpha + 1) = \alpha ^ 2 + \alphaだが,このままビット表現にすると101になり,2ビットで表せる数を超える.このままでは決まったサイズで表すことができないので,\alpha ^ 2を等しい値で置き換える.F(\alpha) = \alpha ^ 2 + \alpha + 1 = 0なので,\alpha ^ 2 = \alpha + 1よりx \ast y = \alpha ^ 2 + \alpha = \alpha + 1 + \alpha = 1が求まる.

template<class T, std::enable_if_t<std::is_unsigned_v<T>, T> Polynomial>
class gf {
public:
    gf() = default;
    gf(const gf&) = default;
    explicit gf(T value)
        : value_{value}
    {}
    friend gf operator+(gf a, gf b) noexcept {
        return gf{a.value_ ^ b.value_};
    }
    friend gf operator-(gf a, gf b) noexcept { return a + b; }
    friend gf operator*(gf a, gf b) noexcept {
        gf res{};
        while(b.value_) {
            res.value_ ^= (b.value_ & 1 ? a.value_ : 0);
            a.value_ = !(a.value_ & (1 << (8*sizeof(T) - 1))) ? a.value_ << 1 : ((a.value_ << 1) ^ Polynomial);
            b.value_ >>= 1;
        }
        return res;
    }

private:
    T value_;
};

カラツバ法やSchönhage–Strassen法といった乗法の高速化手法もある.

除法

 除法は割る数の逆元を掛けることで実現する.有限群Gでは任意の元a\in Gについてa ^ {|G|-1}=a ^ {-1}となる.体は零元を除いて乗法に関して群をなすので,\mathbb{F} _ {2 ^ d}においても任意の元の2 ^ d - 2乗はその元の逆元となる.a ^ nを愚直に計算してもよいが,大きな体では逆元計算に時間がかかるようになる.べき乗は次数の2進展開をもちいて高速化できる.例えばa ^ {25}ならば,25を2進展開し1 1001とし, a ^ {25} = a ^ 1 \ast a ^ 8 \ast a ^ {16}となる.8乗はa ^ {2 ^ {2 ^ {2}}}と計算すれば掛け算は3回で済む.C++で書けば以下のようになる.

template<class T>
constexpr T pow(T a, unsigned int e) {
    T r{1};
    while (e) {
        if (e & 1) {
            r *= a;
            --e;
        }
        else {
            a *= a;
            e >>= 1;
        }
    }
    return r;
}

 最終的なソースコードは以下のようになる.

#include <type_traits>
#include <limits>

template<class T>
constexpr T pow(T a, unsigned int e) {
{
    T r{1};
    while (e) {
        if (e & 1) {
            r *= a;
            --e;
        }
        else {
            a *= a;
            e >>= 1;
        }
    }
    return r;
}

template<class T, std::enable_if_t<std::is_unsigned_v<T>, T> Polynomial>
class gf {
public:
    gf() = default;
    gf(const gf&) = default;
    explicit gf(T value)
        : value_{value}
    {}
    constexpr gf inv() const noexcept {
        return pow(*this, std::numeric_limits<T>::max() - 1);
    }
    constexpr explicit operator T() const noexcept { return value_; }
    friend gf operator+(gf a, gf b) noexcept {
        return gf{a.value_ ^ b.value_};
    }
    friend gf operator-(gf a, gf b) noexcept { return a + b; }
    friend gf operator*(gf a, gf b) noexcept {
        gf res{};
        while(b.value_) {
            res.value_ ^= (b.value_ & 1 ? a.value_ : 0);
            a.value_ = !(a.value_ & (1 << (8*sizeof(T) - 1))) ? a.value_ << 1 : ((a.value_ << 1) ^ Polynomial);
            b.value_ >>= 1;
        }
        return res;
    }
    friend gf operator/(gf a, gf b) noexcept {
        return a * b.inv();
    }
private:
    T value_;
};

羊さんの換字式暗号

(一度非公開にした記事の再投稿です)

ここで扱う内容はすべて解読済みで,鍵も本人から公開されている*1が,もし鍵が公開されない場合でも現実的に解けるのかどうかを検証する.

概要

角巻わためさんのオリジナル曲『ゼロの足跡』には単一換字式暗号で暗号化された文が含まれている.単一換字式暗号では平文と暗号文の文字は一対一で対応する.また,ある言語で書かれた文章に含まれる文字の頻度は,どの十分長い文章でも一定になるため,単一換字式暗号の解読には頻度分析が有効である.しかし,ゼロの足跡に含まれる暗号文は24音しかなく,頻度分析の有効性に疑問が残る.本記事では非常に短い暗号文を用いた単一換字式暗号系の暗号文単独攻撃について考察する.

モチベーション

私の専攻は暗号理論で,普段は楕円曲線暗号などをやっているが,たまには原始に還って古典暗号を考えるのもなかなか趣があると思った.

前提知識

単一換字式暗号

単一換字式暗号とは文字集合\Sigma=\lbrace \text{a,b,c,...,x,y,z}\rbrace についての写像\Phi:\Sigma \rightarrow \Sigmaで定義される.

\Phi全単射写像でなければならない.全単射写像\Phiには逆写像\Phi^{-1}が定義される. 単一換字式暗号の暗号化処理は^\forall\sigma\in\Sigmaに対して\Phi(\sigma)である. 復号処理は\Phi^{-1}を用いて\Phi^{-1}(\sigma)と表せる.

単一換字式暗号の鍵空間のサイズは|\Sigma|!である.

頻度分析

ある言語で書かれた文章に含まれる文字の出現頻度の分布は,文章が十分長ければほかの文章でも同じような分布が現れる.文字から文字へ移す全単射写像で文字を置換するシーザー暗号や単一換字式暗号では,暗号文中の文字の出現頻度の分布は,平文の文字の出現頻度の分布と似た形になるため,頻度分析が解読の助けとなることがある.

ゼロの足跡の暗号

ゼロの足跡の暗号文は(tu tu ke fa o i ni so si ku tu tu lu)(tu tu fa fa o e o li te fo ci)であることが明かされている.

暗号の素敵な歌唱部分はここの冒頭で聴くことができるが,フルで聴くとなお良い.

解読

前述のとおりゼロの足跡の暗号文は(tu tu ke fa o i ni so si ku tu tu lu)(tu tu fa fa o e o li te fo ci)である. 24音からなる短い暗号文だが,子音と母音の組もしくは母音単独でトークンになっていることから,平文は日本語であり,暗号化写像\Phi\text{子音}\rightarrow\text{子音}, \text{母音}\rightarrow\text{母音}であると推測される.また,日本語に含まれない子音\lbrace \text{c, j, l, q, v, x} \rbraceは平文にも含まれていないと推測できるため,\Phi^{-1}の終域で\lbrace \text{c, j, l, q, v, x} \rbraceを考える必要がなくなる.この時点で,鍵空間のなかで考えるべき鍵の数は(21-6)!\cdot 5!=156920924160000 \lt 2^{48}となる.

暗号文は歌詞の一部なので,文字の出現頻度は平文の歌詞と似る可能性が高い.以下に平文歌詞と暗号文のそれぞれの文字の出現頻度を示す.

平文と暗号文のそれぞれの文字の頻度:平文歌詞は読みを平仮名に変換したうえで頻度を分析した

暗号文において最も頻度の高いtuは暗号文中ではtu tuと必ず二つ連続している.平文歌詞で特に出現頻度が高い子音と母音の組は\lbrace\text{き, こ, た, み, な, と, に, し, ...}\rbraceであり,且つ平文歌詞中で連続して出現するものは\lbrace\text{こ, す, で}\rbraceのみである.「ここ」という並びは平文中に4回含まれている.このことから,tu tuに対応する平文が「ここ」で,\Phi^{-1}\text{t}\mapsto\text{k}, \text{u}\mapsto\text{o}である可能性が高い.これが正しいとして,この時点で解読できていない文字は子音が14個,母音が4個であり, 2^{41}\gt14!\cdot 4!=2092278988800 通りの鍵が可能性としてあり得る.

暗号文中で連続しているほかのトークンにfaがある.平文歌詞中で連続して出現する「こ」以外の文字は「す」と「で」である.faが「す」であると考えるとfa o ifa o eといった暗号文が「す<母音><母音>」に対応することになるが,そのような並びは平文歌詞中には登場しないことからfaが「す」に対応するかどうかは分からない.一方,平文歌詞中に「で<母音><母音>」という並びは「であう」として頻出するため,faは「で」,oは「あ」である可能性が高い.ofaの頻度からみてもそれぞれ平文の「あ」,「で」に対応することは不思議ではない.この時点で解読できていない文字は子音が13個,母音が2個であり,2^{34}\gt 13!\cdot 2!=12454041600 通りの鍵があり得る.

現在までの解読が正しいとすれば途中までの復号は次のようになる.「(こ こ ke で あ i ni sa si ko こ こ lo)(こ こ で で あ eli keci)」ここまでで解読できていない暗号文中の子音は\lbrace \text{k,n,s,l,c} \rbraceの5つ,母音は\lbrace\text{e,i}\rbraceの2つまで減っていた.復号結果の候補は2^{19}\gt P(13,5)\cdot 2!=308880 通りとなる.

暗号文中のkeは「ここ」と「出会」に挟まれており非常に高い確率で格助詞と対応する.平文歌詞に出現する文字のうち,格助詞として挙げられるものは\lbrace \text{と, に, が, で, の, を} \rbraceであるが,ここからすでに判明している母音と子音が含まれたものを除くと\lbrace \text{に} \rbraceしか残らない.よって\text{k}\mapsto\text{n}, \text{e}\mapsto\text{i}, \text{i}\mapsto\text{u}とわかる.

現在までの復号結果は次のようになる.「(こ こ に で あ う nu sa su の こ こ lo)(こ こ で で あ い あ lu き だ cu)」現時点での復号結果の候補は2^{15}\gt P(12,4)=11880 通りある.これ以上の復号は頻度分析や日本語の文法からは期待できず,あとは経験と勘で解読する.鍵と完全な復号結果は以下から閲覧できる.

まとめ

24音という非常に短い暗号文でも同一の世界観の下で書かれた平文の標本があれば,頻度分析が一定程度解読の助けになることを確認した.実際に頻度の高いtuが同じく頻度の高い「こ」と対応している.また,同音の連続や母音の連続といった,暗号文中の特徴的な文字の並びも解読においては重要であることも確認した.

これらの重要な要素を十分に活用すれば非常に短い暗号文からでも平文に関する情報を多く得られることが判明した.攻撃者の経験と勘次第では完全な解読も可能だと思われる.


マドロミの暗号も解読したので是非見てくれよな. I also decoded a cipher in "Madoromi."

ouchiminh.hatenablog.com

【記号摂動】幾何計算で退化に対処したい

幾何学計算は、誤差が全くないような整数同士の演算であっても特定の状況で正しく計算ができなくなることで私のような一般学生を困らせる。 二次元の凸包を作るときには三角形の符号付面積が0になる状況、二次元のドロネー図を得るときには三角形の外接円周上に三角形の頂点以外の点がのっている状況などが、正しく計算できない状況ということになる。 このような状況を退化といい、退化の状況は幾何アルゴリズムによって千差万別だ。

以下計算には誤差がないと仮定して記述する。*1

個別に対処するのはしんどい

凸包を作るときに三角形の符号付面積が0になる状況が退化と呼ばれるというのは先ほど述べた通りだ。この状況では何らかの特別な処理を行わなければ計算が破綻しプログラムが無限にループするか、異常終了するか、あるいはとんでもない結果が出力される。ここではGrahamのアルゴリズムを例に退化への対処方法を説明しよう。Grahamのアルゴリズムは凸包を作るアルゴリズムで、実装も大変簡単である。2次元の点集合p_0, p_1, ..., p_{n-1}から凸包を求めるGrahamのアルゴリズムは以下のような手順で構成される。

  1. 点のx座標の小さい点を求め、それをq_0とする。
  2. q_0からの偏角で、他の点をソートする。ソートした点の並びをq_1, q_2, ..., q_{n-1}とする。
  3. スタックを用意し、q_0, q_1, q_2を入れる。
  4. i=3, ..., n-1に対して次のことを繰り返す。
    • スタックの最後から2番目の点uと最後の点vに対して、u, v, q_iが時計回りなら (=符号付面積が負なら) vをスタックからpopし、もう一度同じ操作を行う。
    • そうでなければq_iをスタックに追加する。

このアルゴリズムの1,2,4のステップで退化が生じ得る。そして各ステップごとに退化への対処方法が異なる。

ステップ1で生じ得る退化とは、x座標が最小である点が複数存在する場合だ。この場合は、y座標も考慮してx座標とy座標が共に最小となる点をq_0とする。

ステップ2ではq_0とある点を結ぶ直線上に他の点が乗っている (偏角が同じ) 場合だ。この場合はq_0からの距離を考慮して近い方を削除すればよい。

ステップ4ではu, v, q_iの符号付面積が0になる (3点が直線上に並んでいる) 場合だ。この場合はvをスタックからpopすればよいが、上の二つの例外処理が完全でない場合はスタックの容量が2未満になる可能性がある。

そしてすべてのステップで共通して問題になるのが全く同じ場所に複数の点があるという状況だ。これに対処するためには同じ場所にある点を一つのみ残して残りを削除するという処理が必要になるだろう。

このように一つのアルゴリズムでも様々な退化が生じ、そのたびに違った例外処理を書く必要があるというのは、腱鞘炎を助長する悪しき作業だ。絶対にやりたくないという強い気持ちを持とう。

アルゴリズムに手を加えず退化に対処する記号摂動

概要

私が知っているすべてのアルゴリズムに於いて、退化というのは何らかの計算によって求まった値が別の値と完全に等しいことが原因で起きている。ゆえにすべての退化は入力の値をほんの少しだけずらせば解決する。ただ、実際に値を変更したうえで計算を行うと、それによってまた別の図形が退化しているかもしれないので、無限小\varepsilonを導入して、各点に固有な無限小を足した値を入力としよう。すなわち今までp_i=(x_i, y_i)であったものが p_i = (x_i+\varepsilon ^i, y_i+\varepsilon ^{M_i}) (ただしM_iは十分大きな正の整数) となり、各座標は多項式で表されるというわけだ。

このようにして作った多項式同士の加法と乗法の演算結果は多項式になる。よって多項式全体の集合は加法と乗法に閉じており、加法に対して逆元を持つ環であることが明らかにわかる。では、この多項式環に順序関係を導入しよう。\varepsilonも所詮は定数なので特に難しいことを考えなくても自然に順序関係を導くことができる。

ここに多項式


F=f_0+f_1 \cdot \varepsilon + f_2 \cdot \varepsilon  ^2 + ... + f_n \cdot \varepsilon ^n

がある。この時、Fの符号を、f_0, f_1, ..., f_nの順に見ていった時に最初に見つかる0でないf_iの符号としよう。これはつまり、f_{i+1} \cdot \varepsilon ^{i+1}以降の項がどうであれ、f_i \cdot \varepsilon ^iの項の符号を反転させるほど影響はないという意味である。もともと\varepsilonは無限小なので、これは自然な定義だ。このことを利用すれば二つの多項式F, Gの順序関係は、F-Gの符号がマイナスの時、F \lt Gであることが導ける。

ここで具体的な計算をしてみよう。例えば3つの点p_0=(0+\varepsilon,0+\varepsilon ^{M_1}),p_1=(1+\varepsilon ^2,1+\varepsilon ^{M_2}),p_2=(2+\varepsilon ^3,2+\varepsilon ^{M_3})がどのような位置関係にあるかを計算して求めたいとする。コンピュータでは3点で構成される三角形の符号付面積を計算してそれを確かめる。



  \rm{signed\_area}(p_0, p_1, p_2) = \left|
    \begin{array}{ccc}
      x_{p_0} & y_{p_0} & 1 \\
      x_{p_1} & y_{p_1} & 1 \\
      x_{p_2} & y_{p_2} & 1 \\
    \end{array}
  \right|

無限小摂動させた入力に対して計算を行ったものを\varepsilonの指数ごとに展開すると0 -1 \cdot \varepsilon ^1 + 2 \cdot \varepsilon ^2+...となる。この時定数項は0なので普通に計算すれば退化していたが、\varepsilon ^1の項の係数が-1なので、この多項式全体の符号はマイナスということになる。結果p_0,p_1,p_2を結ぶ線は右に折れ曲がっているという判定になり、退化は解消されている。

実装

ここまでの話は計算幾何学の本を買えば必ず書いてあるような一般的なことだ。本を読めば点をちょっと動かせば退化が解消できるのは分かったが、では結局それをどのように実装するのかという問題がまとわりつく。

今回私は、入力データを多項式として扱う専用のmetric型を作り、その型に対して演算を定義した。

github.com

例えば符号付面積を求めたいときならば以下のようにして扱う。

    using namespace ouchi::geometry;
    ouchi::math::fl_matrix<metric<int>, 3,3> sa = {
        // metricのコンストラクタには座標, εの指数の順に指定する。
        metric<int>(0, 1ull), metric<int>(0, 20ull), 1,
        metric<int>(1, 2ull), metric<int>(1, 400ull), 1,
        metric<int>(2, 3ull), metric<int>(2, 8000ull), 1
    };
    auto signedarea = det(sa);
    assert(signedarea < 0);

アルゴリズムを記述したコードから個別の退化への対処が消え、一般の場合のみについて記述したコードで正しく動くようになる。また、一般の処理のみを施した既存のコードへの変更は、演算対象の型を変更するだけでよいので、大変少ない変更で済む。

記号摂動法は自動で退化に対処する非常に優れた手法だが、計算時間の増大は問題だ。浮動小数点加速のように退化が生じた時だけ整数帰着と記号摂動により退化を解消する手法もあるが、実装が面倒くさそうなので今回は扱わない。

結果

以下に120 \times 120の二次元空間内に一様に配置した10000点の集合に対してGrahamのアルゴリズムによって凸包を得た時記号摂動の有無による出力の違いが分かるよう、それぞれの場合の出力を画像で示す。例外処理がない場合は完全に壊れてしまっているが、記号摂動法によって退化を解消した場合の出力ではきちんと凸法になっていることが分かる。

f:id:ouchiminh:20200604200755p:plainf:id:ouchiminh:20200604200930p:plain
左:例外処理なし 右:記号摂動


参考文献 - 杉原厚吉.『計算幾何学』 朝倉書店 (ISBN978-4-254-11681-6)

*1:有限桁のコンピュータでも絶対値が十分小さな整数なら加算、減算、乗算には誤差がない。幾何計算を誤差なく行うためには整数帰着法を用いて、入力を十分な桁数の整数に変換してから計算する。

C++で行列のなんやかんやを実装した話

これは最近実装したやつの解説です。主に自分が忘れたときのための備忘録として書いていますから他の人が読むにはあんまり親切ではないかもしれません。

静的なのと動的なの

C++のテンプレートはとても強力なので大抵のことはテンプレートで解決できる。行列を使って何かを計算したいとき、行列のサイズ(型)*1が静的でよい場合と動的であってほしい場合がある。例えばシャミアの秘密分散法などなどで連立一次方程式を解く場合には行列のサイズが動的であると便利だ。しかし、ゲームなどでオブジェクトの描画位置を変えたいというときには3次の正方行列でよいし、ヒープにメモリを確保する必要もない。だからといって静的なサイズの行列と動的なサイズの行列で全く別のクラスにしてしまうのは面倒だし、どちらの行列も同じように扱えるほうがよいに決まっている。なので作った。

github.com

作り方

detection idiomでSizeを制約

まずは定義だけ見てみましょい。

template<class, class, class = void>
class basic_matrix;

template<class T, class Size>
class basic_matrix<T, Size, std::enable_if_t<is_size_spec_v<Size>>> {
    ...
};

Sizeというテンプレート引数にはこちらで用意したサイズ指定子以外を指定されたくない。そのため、detection idiomで、サイズ指定子以外が指定されたbasic_matrixが実体化されないようになっている。

is_size_spec_v<Size>というのがSizeがサイズ指定子かどうかを判定してbool値を返すメタ関数だ。C++ではテンプレートの実体化の優先順位が決まっていて、完全特殊化>部分特殊化>プライマリーテンプレートの順に実体化が試みられる。is_size_spec_v<Size>trueに評価されるとき、部分特殊化がwell-formedになるのでそれが実体化される。またis_size_spec_v<Size>falseになるときはプライマリーテンプレートが実体化される (が、定義がないのでコンパイルエラーになる)。

これでbasic_matrixの実体化については終わり。次は便利な関数を生やしながらそいつらを静的なサイズのやつと動的なサイズの奴の両方に適用できるようにしていく。

メンバ関数とゆかいな仲間たち

サイズ指定子

まずは先ほど紹介しそびれたサイズ指定子について話そう。サイズ指定子は行列のサイズが静的か動的かを選ぶときに指定する。サイズが静的な場合は、実際のサイズの情報もサイズ指定子が持つ。また、basic_matrixは各成分を保存するコンテナにサイズ指定子で宣言されるコンテナ型を使う。静的な場合はC++17の時点でconstexprにふるまうstd::array、動的な場合はstd::vectorを使う。

// 固定長
template<size_t Row, size_t Column>
struct fixed_length final : detail::size_base {
    static constexpr size_t row_size = Row;
    static constexpr size_t column_size = Column;
    template<class T>
    using container_type = std::array<T, row_size * column_size>;
};
// 可変長
struct variable_length final : detail::size_base {
    template<class T>
    using container_type = std::vector<T>;
};

演算可能性

行列はそのサイズによって定義される演算が異なる。たとえば違うサイズの行列には和は定義されないし、左辺の列の数と右辺の行の数が異なれば積は定義されない。ユーザーがある演算を試みた時、それがコンパイル時に演算可能か不可能かが分かるのならばコンパイル時にお知らせしたいと思うのが人間だ。ただし、basic_matrixは動的にサイズを変更できる場合もあるので演算可能性はtrue, falseの二値では表すことも出来ない。よって以下のような三値の条件値を定義した。

enum class condvalue {
    no, yes, maybe
};

例えば、basic_matrix<T, fixed_length<2, 2>>basic_matrix<T, fixed_length<2, 2>>の和を求めたいときはcondvalue::yes

basic_matrix<T, fixed_length<2, 2>>basic_matrix<T, fixed_length<2, 3>>の和を求めたいときはcondvalue::no

basic_matrix<T, fixed_length<2, 2>>basic_matrix<T, variable_length>の和を求めたいときはcondvalue::maybeを返すようなメタ関数を書けば、コンパイル時にエラーを予期したり、コンパイル時エラーにしたりできる。

例えば和

以上を踏まえて、戯れに和を計算する演算子を定義しよう。

template<class T, class Size>
class basic_matrix<T, Size, std::enable_if_t<is_size_spec_v<Size>>> {
    ...
    ...
    template<class U, class S>
    [[nodiscard]]
    friend constexpr auto operator+(const basic_matrix& a, const basic_matrix<U, S>& b)
        noexcept(is_fixed_length_v<Size>&& is_fixed_length_v<S>)
        ->std::enable_if_t<(detail::add_possibility_v<Size, S> > detail::condvalue::no), basic_matrix<std::common_type_t<T, U>, detail::add_possibility_t<S, Size>>>
    {
        basic_matrix<std::common_type_t<T, U>, detail::add_possibility_t<S, Size>> res;
        // 分岐はコンパイル時だが、副文の実行は実行時
        if constexpr (is_variable_length_v<detail::add_possibility_t<Size, S>>) {
            if (a.size() != b.size()) throw std::domain_error("two matrixes that have different size cannot be added each other.");
            res.resize(a.size().first, a.size().second);
        }
        basic_matrix::add(a, b, res);
        return std::move(res);
    }

};

まずはadd_possibility_vというメタ関数だが、これは先に述べた和が定義できるかどうかをno, yes, maybeの三値で返してくれるメタ関数だ。 ここではそのメタ関数で右辺を制約している。戻り値の型を後置する関数宣言構文なのでぱっと見わかりにくいかもしれないが、std::enable_if_tの第一引数には和が定義できる条件になっている。また、動的なサイズの行列の場合、constexpr if文で分岐し、動的に和が定義できるかチェックしている。

このように各演算ごとに定義できるかどうかを三値で返すメタ関数を書きまくれば他の演算についても同じ手法で演算子メンバ関数を書くことができる。おそらくこれが一番楽な方法だと思う。

まとめ

やっぱりC++はなんでもできる。

*1:以降型と書くとC++の型と紛らわしいため行列の型については行列のサイズと書く

定数回のループを高速化する with fold expression

久しぶりにネタになりそうなことが見つかったのだ!更新してないからか、読む人はめちゃめちゃ少ないのだけれど、ありがたいことに0人ではないので、せっかく見つかったネタを自分だけの知見とせずにアウトプットしたいと思う。

またC++の話だけれど、ライブラリを作っててちょっと高速化できた書き方を紹介するよ。

本題

forwhiledo-while goto いろいろあります 反復文。何を使ってもいいが静的に反復する回数が決まるならそれらは必要ないのだ!例えばforを書く。

constexpr int cnt = 10;
for(auto i = 0; i < cnt; ++i) { some_work(); }

上のように書いた時、cntが定数なのにiを一回ごとに1ずつ加算してcntと比較してi < cntならば最初に戻るという一連の動作をする。おかしい。どう考えても無駄だ。forの中身をそのまま10回書いた方がどう考えても速い。

some_work();
some_work();
some_work();
some_work();
some_work();

some_work();
some_work();
some_work();
some_work();
some_work();

まあだがこの書き方は誰が見てもバカなので、モダンで洗練されていてなんでも出来てかっこいいC++には似合わない。確かこういう時に使えそうな機能が…… ほら、あったじゃないか。C++17で追加されたモダンで洗練されていてかっこいいfold expressionが!!!!

template<size_t ...S>
void do_some_work(std::index_sequence<S...>) {
    ((S, some_work()), ...);
}
do_some_work(std::make_index_sequence<10>{});

速くてかっこいいC++にはこういう書き方が似合う*1

速くなるか半信半疑の人のためにwandboxへのリンクを貼っておく。

wandbox.org

*1:諸説あります

久しぶりに大学の課題で本気を出した。

情報工学科で出る定番の課題といえば、ソートアルゴリズムの実装だ (多分) (他の大学を知らない) 。 普通の人間は愚直に実行時に配列などをソートしてしまうかもしれない。私も高校時代にこの課題が出ていたら普通に実行時にソートして喜んでいただろう。 だが、ここは大学。生半可なソートを実装して喜んでいるようではレベルが低い (ハードウェアに近いという意味ではない) 。

ヒープソートを実装して実行時間を計る課題が出たので、私はコンパイル時にstd::integer_sequenceヒープソートした。これくらいやらなければ実行時間で他の学生に圧倒的な差をつけて勝つことはできない。ソートする整数列が長くなれば、恐ろしいほどメモリを食うが、一応メモリが無限にあればコンパイルできるようになっている。(余談だが、コンパイル時ソートを実装するのは2回目だ。以前にコンパイル時挿入ソートを実装したことがある)

そもそもヒープソートって?

ヒープソートはO(n log n)で要素をソートできるアルゴリズムだ。常に最小値または最大値がヒープのトップに来るように完全2分木を作り、末尾と先頭を入れ替えながらヒープを縮小していくと最後には綺麗に並んでいるというアレだ。アルゴリズムクイックリファレンスを読めば書いてあることなので知らない人は読んでくれ。

ヒープソートを実装する

まずヒープを作るには以下の操作が必要だ。

  • 整数列のi番目の親と子の添え字を求める
  • 条件を満たすとき、あるノードをその親とスワップする
  • 条件を満たすとき、あるノードをそのどちらかの子とスワップする

しかし、std::integer_sequenceの整数列のN番目とM番目の整数をスワップするのは困難を極めるか不可能なので、結果的にスワップしたように見える動作をするために、std::integer_sequenceのN番目の整数をIにする関数writeを書く。 また、そのIを求めるためにstd::integer_sequenceのN番目の整数を求める関数atも書く。

そのようにしてできたwriteatswapを作れば、条件を満たすときにswapしてくれるような関数swap_ifを書くと、いよいよヒープの本質を実装できるようになる。

実装したのがこちらだ。昨日書いたものだが、もうすでにどうして動いているのか分からない。もっといい書き方があれば教えてくれ!

github.com

サンプルコード

int main() {
    std::integer_sequence<int, 1, 3, 2, 8, 10, 0> seq;
    static_assert(std::is_same_v<
        decltype(heapsort<std::less<>>(seq)),
        std::integer_sequence<int, 0, 1, 2, 3, 8, 10>;
    >);
}