Rokiのチラ裏

学生による学習のログ

markdown から数式のみを取り出して PNG として出力し該当箇所をパスで置換

なんだか簡単な使い捨てスクリプト程度にしようと思っていたものの, 微妙にしっかり作ってしまったのでその記録. この記事の執筆時現在, このログのように, gitbook でのビルドが失敗することがある. 発生する条件としては, 文中に大量の数式があり, それを mathjax プラグインによってレンダリングしようとすると起こるようで, 元々の問題源は gitbook が利用する mozilla 発の nunjucks によるものとされていたようだが, その修正のマージ後も未だ治ってはいない*1. 元からこの問題を認知していれば良かったものの, コンテンツ内に数式が増えてきてから発生したので対応も出来ず, 今さら gitbook からなにか他のものに移行するわけにもいかず, とりあえず応急処置として, 特定のディレクトリツリーを再帰的に検査し, markdown テキストから独自のブロックで囲まれた数式のみを取り込み(このファイルのように```mr```mrendで囲んでいる), 該当部分を PNG へ出力し, さらに画像のパスに自動で置き換えるようにした. 数式が書き加えられると, 随時それを PNG にしてそのパスで置換するようにしているので, 数式を画像にしたことによって, 作業間で特別な手間を感じずには済んでいるものの, やはり画像の読み込みは重いので, そのうち自前でホストするなど、移行を考えざるをえない状況になっている.
特に gitbook の PDF 出力機能では SVG がサポートされていないために, 出力を PNG 形式で行っている. また, gitbook IO によるプラグイン, gitbook-mathjax は目次から mathjax によって変換される数式を含んだコンテンツを開くと正常にプラグインが読み込まれず, 結果として再読み込みを行う必要があるなどの問題がある. はっきり言って, 今から数式を多分に含むような文書を gitbook によって管理することは, このような理由から全くもってお勧めできない.
今回の実装では Typescript を利用したわけだが(tumblr クライアントの開発を含めて今回で 2 回目の利用), 今更ながら全体的に関数型言語っぽく書けるリスト処理に少し感動した.

Proof of infinite geometric series

16.8 float のしくみ · ThePolitewaylearntoCPP17 の補足記事。命題「 a + ax + ax^{2} + \cdots -1 \lt x \lt 1 のとき収束し、その値は  \dfrac{a}{1-x} である。 x \leq -1, 1 \leq x のとき発散する」は、高校数学の範囲なので、特別に取り上げる必要はないと思ったが、コンテンツ内容をなるべく自分で書いた文章のみ完結できるよう、本エントリで扱うこととした。といっても、証明はとても簡単に行うことができる。


 x \neq 1 とする。無限等比級数 n 項までの和を  S_{n} としたとき  S_{n} = a + ax + ax^{2} + \cdots + ax^{n-2} + ax^{n-1} \tag{1} が得られる。両辺を  x 倍すると  xS_{n} = ax + ax^{2} + ax^{3} + \cdots + ax^{n-1} + ax^{n} \tag{2} が得られる。このとき式  (1) で 式  (2) を引くと  S_{n} - xS_{n} = a -
 ax^{n} が得られる。 S_{n} について解くと  S_{n} = \dfrac{a-ax^{n}}{1-x} \tag{3} となる。また  \displaystyle\lim_{n\to\infty}x^{n}=0\ \ (-1 \lt x \lt 1) \tag{4} がいえる(絶対値が  1 未満である値を掛け続けていくといずれ  0 になる)。

求める無限級数の値は  \displaystyle\lim_{n\to\infty}S_n であるので 式  (3), (4) を利用して  \displaystyle\lim_{n\to\infty}\dfrac{a(x^{n}-1)}{x-1} = \dfrac{-a}{x-1} = \dfrac{a}{1-x} と導出される。  \therefore 題意は示された。

Divergence of the sum of the reciprocals of the primes

素数の逆数和が発散する事(式  (1))についてメモ。
 \displaystyle\sum_{p\ prime} \dfrac{1}{p} = \dfrac{1}{2} + \dfrac{1}{3} + \dfrac{1}{5} + \dfrac{1}{7} + \dfrac{1}{11} + \dfrac{1}{13} + \dfrac{1}{17} + \cdots = \infty \tag{1} 素数の逆数和が発散することの証明には Harmonic series*1 が発散する事(式(2))を利用する \displaystyle\lim_{n\to\infty}\sum_{k=1}^n\dfrac{1}{k}= 1 + \dfrac{1}{2} + \dfrac{1}{3} + \dfrac{1}{4} + \cdots = \infty \tag{2} 式 (2) の証明には、式変形、等式変形による不等式の評価を行う方法と、積分を用いた不等式評価など*2による方法がある。まず式  (2) を証明し、その後式  (1) を証明する。

不等式による (2) の証明

他の発散系列と比較する事で Harmonic series の発散を容易に証明できる。\displaystyle 1 + \dfrac{1}{2} + \dfrac{1}{3} + \dfrac{1}{4} + \dfrac{1}{5} + \dfrac{1}{6} + \dfrac{1}{7} + \dfrac{1}{8} + \dfrac{1}{9} + \cdots \tag{3}  \displaystyle 1 + \dfrac{1}{2} + \dfrac{1}{4} + \dfrac{1}{4} + \dfrac{1}{8} + \dfrac{1}{8} + \dfrac{1}{8} + \dfrac{1}{8} + \dfrac{1}{16} + \cdots \tag{4} このとき、\displaystyle 式 (3) \gt 式 (4) である。式  (4) は次の等式  (5) の通り発散する。  (4) = 1 + (\dfrac{1}{2}) + (\dfrac{1}{4} + \dfrac{1}{4}) + (\dfrac{1}{8} + \dfrac{1}{8} + \dfrac{1}{8} + \dfrac{1}{8}) + (\dfrac{1}{16} + \cdots + \dfrac{1}{16}) + \cdots = 1 + \dfrac{1}{2} + \dfrac{1}{2} + \dfrac{1}{2} + \dfrac{1}{2} + \cdots = \infty \tag{5} 現時点で\displaystyle 式 (3) \gt 式 (4)の関係性から明白であるが、さらに、この変形を一般化する事ができ、次の不等式が得られる。 \displaystyle \sum_{n=1}^{2^{k}} \dfrac{1}{n} \geq 1 + \dfrac{k}{2} \tag{6}  (6) において  k \to \infty とすると、右辺が発散するため左辺も発散する事が言える。 \therefore 題意は示された。

積分による (2) の証明

Integral Test.svg
By Jim.belk - Own work, Public Domain, Link

Harmonic series の総和と広義積分を比較する事で、発散する事を証明できる。各図の面積の和は、 \displaystyle 1 + \dfrac{1}{2}
 + \dfrac{1}{3} + \dfrac{1}{4} + \dfrac{1}{5} + \cdots + \dfrac{1}{n} と同等である。一方で曲線  y =  \dfrac{1}{x} に着目し、 \displaystyle x \in [1, \infty) の部分の下にある面積は  \displaystyle \int_{1}^{n + 1}\dfrac{1}{x}dx であり、図から次の不等式が成り立つ。 \displaystyle \int_{1}^{n + 1}\dfrac{1}{x}dx \lt 1 + \dfrac{1}{2}
 + \dfrac{1}{3} + \dfrac{1}{4} + \dfrac{1}{5} + \cdots + \dfrac{1}{n}
ここで  \displaystyle \int_{1}^{n + 1}\dfrac{1}{x}dx = \ln{(n + 1)} であるため、次の不等式が成り立つ。 \displaystyle \ln{(n + 1)} \lt \dfrac{1}{2}
 + \dfrac{1}{3} + \dfrac{1}{4} + \dfrac{1}{5} + \cdots + \dfrac{1}{n} < 1 + \ln(n) \tag{7}
 (7) において  n \to \infty とすると、 \displaystyle \lim_{n \to \infty}\ln{(n+1)}=\infty から左辺が発散する事が言え、右辺も発散する事が言える。 \therefore題意は示された。

余談

調和級数  \displaystyle\lim_{n\to\infty}\sum_{k=1}^n\dfrac{1}{k} は見ての通りゼータ関数  \displaystyle \zeta(s)=\displaystyle\sum_{n=1}^{\infty}\dfrac{1}{n^{s}} \displaystyle \zeta(1) に等しく、ゼータ関数 s \geq 2のときのみ収束する。

(1) の証明

 n 以下の素数を昇順に  \displaystyle p_{1}, p_{2}, \cdots, p_{t} とすると次の式が成り立つ。  \displaystyle \sum_{k=1}^{p_{n}} \frac{1}{k} \lt\prod_{k=1}^{n}(1+\dfrac{1}{p_{k}}+\dfrac{1}{p_k^{2}}+\dfrac{1}{p_k^{3}}+\cdots) = \prod_{k=1}^{n} \frac{1}{1-\frac{1}{p_{k}}} \tag{8} (8) においては、 \alpha_{1},...,\alpha_{n} 0 以上の自然数とし、左辺を全て展開すると、 1 以上  p_{n} 以下のいかなる自然数 p_{1}^{\alpha_{1}}p_{2}^{\alpha_{2}}\cdots p_{n}^{\alpha_{n}}素因数分解できる。また、同式の右辺を展開すると  \alpha_{1},...,\alpha_{n}
 0 以上のどのような自然数だとしても必ず  \frac{1}{p_{1}^{\alpha_{1}}p_{2}^{\alpha_{2}}\cdots p_{n}^{\alpha_{n}}} となる項がでてくるため、成り立つことが言えるのである。
次に式  (8) の右辺を等比級数の和の公式で変形すると次の式を得る。  \displaystyle\prod_{k=1}^{t}\dfrac{1}{1-\frac{1}{p_{k}}}=\prod_{k=1}^{t}(1+\dfrac{1}{p_{k-1}}) \tag{9} (8) は両辺とも正であるため、次のように不等式の対数を取ることができる。  \ln \left(\displaystyle\sum_{k=1}^{n}\dfrac{1}{k}\right) \lt \displaystyle\ln\left(\prod_{k=1}^{t}(1+\dfrac{1}{p_{k-1}})\right)\
\tag{10}

ln(x+1) と x のグラフ
 \ln{(x+1)} x のグラフ*3
また  x \gt 0 のとき、 \ln{(x+1)} \lt x である。これは、 f(x) = x - \ln{(x + 1)} f'(x) = 1-
 \dfrac{1}{x+1}としたとき、 f'(x)=1-
 \dfrac{1}{x+1}=\dfrac{x}{x+1}\gt0 であり、 f(x)x\gt0で増加関数となることがいえ、また  f(0)
 = 0 - \ln{1} = 0 であることから、 x\gt0のときf(x)\gt0が導かれる。したがって  x\gt0 のとき \ln{(x+1)}\lt xが成り立つ。また  \ln{(x+1)} x のグラフが同様に示す。

この事実から、式  (10) は次のように変形できる。  \displaystyle\sum_{k=1}^{t}\ln(1+\dfrac{1}{p_{k-1}})\
\lt \displaystyle\sum_{k=1}^{t}\dfrac{1}{p_{k-1}} \tag{11} ここで、 \displaystyle p_{k-1}\leq p_{k}-1( k - 1 番目の素数より  k 番目の素数 1 以上大きい)であるから式  (11) を変形し次の式を得る。   1+\displaystyle\sum_{k=2}^{t}\dfrac{1}{p_{k-1}}
\leq  1+\displaystyle\sum_{k=2}^{t}\dfrac{1}{p_{k-1}}\
=1+\displaystyle\sum_{k=1}^{t-1}\dfrac{1}{p_{k}}\tag{12} (10)、式  (11) の関係性から、式  (12) は次のように変形できる。  \ln{\left(\displaystyle\sum_{k=1}^{n}\dfrac{1}{k}\right)} \lt 1+\displaystyle\sum_{k=1}^{t(n)-1}\dfrac{1}{p_{k}}\tag{13} ここで  n\to\infty とすると式 (2) より左辺が発散するため、右辺も発散する。 \therefore 題意は示された。

総括

最後に素数の逆数和をプロット。

f:id:Rok1:20180219061945p:plain

グラフの通り、発散のスピードが非常に遅いことが窺える。素数列の生成などで篩落とす際は、十分に値が増えたとしてももかなり時間計算量を抑えられるであろうことが視覚的に分かる。プロットで用いた値は、自作ライブラリを利用して次のようにして生成した*4

#include <srook/math/primes/progression/sieve_atkin.hpp>
#include <srook/math/constants/algorithm/pow.hpp>
#include <srook/math/constants/algorithm/approximate/reciprocal.hpp>
#include <srook/string/to_string.hpp>
#include <fstream>

int main()
{
    constexpr auto max = static_cast<srook::uint32_t>(srook::math::pow(10, 6));
    std::ofstream ofs("data.csv");
    std::size_t i = 0;
    float sum = 0.0;

    srook::math::primes::progression::compute_sieve_atkin_if(max, 
            [&ofs, &i, &sum](srook::uint32_t x) {
                sum += srook::math::approximate::reciprocal(static_cast<float>(x));
                ofs << srook::to_string(i++) + "," + srook::to_string(sum) + "\n";
            },
            [&ofs, &i, &sum](srook::uint32_t) {
                ofs << i++ << "," + srook::to_string(sum) + "\n";
            });
}

グラフは次のように描いた。

#!/usr/local/bin/gnuplot
reset
set object 1 rect behind from screen 0,0 to screen 1,1 fc rgb "#333631" fillstyle solid 1.0
set term png size 800, 600 enhanced
set output 'image.png'
set grid lc rgb "white" lt 2
set border lc rgb "white"
set key tc rgb "white"
set title tc rgb "white"
set xlabel "x" tc rgb "white" font ",30"
set ylabel "y" tc rgb "white" font ",30" offset 1, 0
set size ratio 2.0/(1.0 + sqrt(5.0))
set palette rgbformulae 22, 13, -31
set style line 1 lw 2 lt rgb "#008AD4"

set xrange [0:]
set logscale x
set format x "10^{%L}"
set format y "%3.1f"

set datafile separator ","
set title "The sum of the reciprocals of the primes"
plot "data.csv" using 1:2 w l ls 1

参考

*1:調和級数。音楽の高調波という概念から派生してこのように名付けられている。

*2:英 Wikipedia Harmonic series には "There are several well-known proofs of the divergence of the harmonic series. A few of them are given below." とある通り、これら以外の証明方法がマクローリン型不等式を用いた方法などを含めて、いくつか存在する

*3:画像は gnuplot で生成した

*4:素数で単に 1 を割れば良いのだが、グラフの様子は大して変わらないので、速度の向上を図るために近似値でプロットしている。言い始めれば、浮動小数点数を使っている時点で近似であるわけだが、106 までの素数の逆数の総和となると図の通り小さい値しか出てこないので、その影響はそこまで酷くもないだろうなどと思いつつ。

zlib ラッパー

少し以前に実装したものの話題となるが、zlib のラッパーを自作ライブラリの方に実装したものの、ブログに特に書いていなかったのでそのメモ。次のように利用できる。

テストは次のように、フリードキュメントである赤毛のアンのテキストを圧縮/解凍をし、それぞれ gzip との実行結果で差異がないかをチェックすることで行なっている。 1 種類のフリードキュメントに対する圧縮, 解凍ではテストが不十分であるように感じたのと, 単一のサーバーに対して負荷がかかってしまうことを避けるために, 候補からランダムでフリードキュメントを選択し, それを利用してテストするようにした.

最後にこれはお遊びであるが、当然 git の blob オブジェクトも同ラッパーで解凍できる(意味はない)。次は普通にgit cat-fileによって blob オブジェクトを解凍/閲覧している例。

$ git clone https://github.com/falgon/git_playground.git
$ cd git_playground
$ git cat-file -p HEAD | grep tree | sed -e s/tree// | xargs git cat-file -p  | cut -d ' ' -f 3 | cut -f1 | xargs git cat-file -p | cut -d ' ' -f 3 | cut -f1  | xargs git cat-file -p
hello

似た機能を同ラッパーで次のように実装できる。blob オブジェクトは、コード内のコメントにある通りある一定のフォーマットのヘッダが先頭にあるので、それをスキップする。

次のように実行できる。

$ g++-7 -std=c++1z -Wall -Wextra -pedantic -lz blob_decompress.cpp -o blob_decompress
$ git clone https://github.com/falgon/git_playground.git
$ cd git_playground
$ ../blob_decompress .git/objects/ce/013625030ba8dba906f756967f9e9ca394464a
hello

TMP におけるバインドとその利用

私の観測範囲内ではあまり有名で無いようなので発信。バインドとは C++ において Callable を満たすものに対し引数の束縛を行い、束縛済みの新たな Callable オブジェクトを生成する事を広義に言うが、これを TMP で行うにはどうすれば良いか。つまり例えば次のようにするにはどうすれば良いか。

typedef bind<std::is_same, int> type;
type::type<int>::value; // == std::is_same<int, int>::value

C++98 から次期標準規格である C++20 の間で、テンプレートテンプレートパラメータを typedef することはできない。

template <template <class...> class F, class... Ts>
struct bind { typedef ??? type; };

しかしこれは意外と単純で、次のようにすれば書ける。

template <template <class...> class F, class... Ts>
struct bind {
    template <class... Us>
    struct type : std::bool_constant<F<Ts..., Us...>::value> {};
};

TMP とは純粋関数型プログラミングであるので、バインドのような機能は役立つ場面が多い。例えば、与えられたポリシーに従う型でグループ化するといった処理には、次のようにして利用するととても実装が楽だ。

これは次のように動作する。

using srook::tmpl::vt::packer;
struct X {};
struct Y { Y(const X&); };

static_assert(
    std::is_same<
        srook::tmpl::vt::group_by_t<std::is_convertible, packer<int, int, X, Y>>,
        packer<packer<int, int>, packer<X, Y>>
    >::value
);

この実装で利用されている take_whileD と drop_whileD は、バインドされたメタ関数を受け付けるメタ関数であるが、基本的には haskell のそれと同等である。尚、この group_by も haskell の それと同等である。さらに、ソート等の処理を記述する際には、指定された順序の逆順を表せる機能があると実装が楽であるがこれも同様にしてバインド機能で達成できる。以下はその利用によるクイックソートの実装である。

この実装では、placeholder の機能を利用して引数の入力をコントロールしている。これは、C++ 標準ライブラリに属するstd::bindを意識して作られたメタ機能であり、次の通り同じようにして動作する。

#include <type_traits>
#include <srook/tmpl/vt/bind.hpp>

template <class T1, class T2, class T3>
struct F : std::conjunction<std::is_same<T1, double>, std::is_same<T2, int>, std::is_same<T3, char>> {};

typedef srook::tmpl::vt::bind<F, srook::tmpl::vt::placeholders::_1, int, srook::tmpl::vt::placeholders::_2> bind_type1;
static_assert(bind_type1::type<double, char>::value);

typedef srook::tmpl::vt::bind<F, srook::tmpl::vt::placeholders::_2, int, srook::tmpl::vt::placeholders::_1> bind_type2;
static_assert(bind_type2::type<char, double>::value);

その他にも、リポジトリにて Variadic templates に対する高度なリスト処理が可能なメタ関数、高階関数などをいくつか公開している。

maximal length sequence

M系列に関する学習メモ*1

import Test.HUnit
import System.IO

mulcon :: Int -> Int
mulcon 0 = 1
mulcon n = (a * mulcon(n - 1) + b) `mod` m
    where
        a = 3
        b = 0
        m = 7

mc :: [Int] -> [Int]
mc = map mulcon

main :: IO (Counts, Int)
main = runTestText (putTextToHandle stderr False) tests
    where
        x = mc [0 .. 11]
        tests = TestList [ "mc test" ~: x ~?= [1, 3, 2, 6, 4, 5, 1, 3, 2, 6, 4, 5] ]

これは線形合同法*2(漸化式  \displaystyle X_{n} = ax_{n - 1} + b \pmod{m})の単純な実装であるが、結果からも分かるように、周期は  m に強く依存するため*3、大きな値  m を用いることで長い周期の乱数を得る事ができる。しかしそれに応じて計算処理の膨大化の他、既知の様々な問題*4がある。M系列は  \bigoplus を Bitwise XOR とした時、次の高次線形漸化式で発生する 1 ビット数列を言い、線形合同法の欠点を克服する。

 \displaystyle \large \begin{equation} x_{n} = x_{n - p} \bigoplus x_{n - q}  \ \ \ \ (p \gt q) \end{equation}

各項の値は 1 ビット、すなわち 0 か 1 で事前に  \displaystyle x_{n-1},\ x_{n-2},\ \cdots,\ x_{n-p} に対して初期値を与えておく。次のコードは、 x_{0} = 1,\ x_{1} = 0,\ x_{2} = 0,\ p=3,q=1 とした場合(話を簡単にするため整数をビットをみなしている)のM系列の実装である。加えて、得られたM系列の先頭から  p ビットを 1 ビットずつずらしてリスト内にリスト化していく関数と、リスト内のリストがその自身を内包するリスト内で唯一のリストであるかを検査する関数を実装した。

import Test.HUnit
import System.IO
import Data.Bits

p :: Int
p = 3

mseq :: Int -> Int
mseq 0 = 1
mseq 1 = 0
mseq 2 = mseq 1
mseq n
  | p > q = (mseq (n - p) `xor` mseq (n - q)) :: Int
  where
    q = 1

invmseq :: [Int] -> [Int]
invmseq [] = []
invmseq (x:xs)
  | x >= 0 = mseq x : invmseq xs

isUniqueImpl :: [Int] -> [[Int]] -> Bool
isUniqueImpl x (y:z)
  | x == y = False
  | otherwise = isUniqueImpl x z
isUniqueImpl x y = [x] /= y

isUnique' :: [[Int]] -> Bool
isUnique' [] = True
isUnique' [_] = True
isUnique' (x:xs)
  | not (isUniqueImpl x xs) = False
  | otherwise = isUnique' xs

shiftpPattern :: [Int] -> [[Int]]
shiftpPattern [] = []
shiftpPattern xs
  | length xs >= p = take p xs : shiftpPattern (drop 1 xs)
  | otherwise = []


main :: IO (Counts, Int)
main = runTestText (putTextToHandle stderr False) tests
  where
    x = invmseq [0 .. 8]
    shiftp = shiftpPattern x
    isunique = isUnique' $shiftp
    tests = TestList
        [ 
            "mseq" ~: x ~?= [1, 0, 0, 1, 1, 1, 0, 1, 0],
            "shift pattern" ~: shiftp ~?= [[1, 0, 0], [0, 0, 1], [0, 1, 1], [1, 1, 1], [1, 1, 0], [1, 0, 1], [0, 1, 0]],
            "unique check" ~: isunique ~?= True
        ]

shiftpが示しているように、先頭から 1 ビットずつずらした  p ビットの数列は、それぞれが他の数列と一切異なる数列、つまり  p ビットで表現できるビットパターンを網羅している事が分かる。しかしそれは全てではなく、全てのビットパターンが 0 となる数列は発生されない*5。よって、M系列の周期  T T=2^{p}-1 で求まる事が言える。
この時選ばれる  p q の値は次の式が既約多項式となるようにして選ばれる。

 \large x^{p} + x^{q} + 1 \ \ \ \ (p \gt q)

この時  x の係数はガロア*6  GF(2) であり 0 か 1 を取る。
次にこれまでで用いた  p=3, q=1 が、上記の条件を満たす事を背理的に証明する。次の方程式   x^{3} + x + 1 = 0 \tag{1}有理数解を持つと過程する。有理数解であるならば、互いに素な p, q を用いて  \displaystyle x = \dfrac{p}{q} とする事ができるため、\displaystyle \dfrac{p^{3}}{q^{3}}+\dfrac{p}{q}+1=0 p^{3}+pq^{2}+q^{3}=0 p(p^{2}+q^{2})=-q^{3} であり、これは  q p の倍数である事を意味する。 p, q は互いに素であるため  p = 1 であり、\displaystyle \dfrac{1}{q^{3}}+\dfrac{1}{q+1}=0 1+q^{2}+q^{3}=0 q^{2}(q+1)=-1 となるが、左辺は  q(q+1) を因数に持つため偶数、右辺は奇数となり矛盾する。よって、方程式  (1)有理数解を持たず、既約である事が言える。
[※ 追記 もっと単純化できた...
 f(x) = x^{3} + x + 1 とする。 f(x) は 三次式であるため  f(x) が位数 2 のガロア GF(2) 上可約であるとすると  f(x)=(x-a)(x^{2}+bx+c) となる  a, b, c \in GF(2) が存在するはずである。すなわち  f(a) = 0 である。しかし  GF(2) = {0, 1} f(0) \neq 0, f(1) \neq 0 であるためそれに矛盾する。よって  f(x) は既約である。
-- 追記ここまで]

蛇足

有名なメルセンヌ・ツイスタ(MT)はこのM系列を発展させたものであり、Twisted General Feedback Shift Register(TGFSR) がはじめに発案された。これは、次の漸化式  x_{n+p} = x_{n+q} + x_{n}A\ \ \ \ (p \gt q \gt 0) によって発生する。 A は Twister と呼ばれる  GF(2) の元を要素とする  \omega \times \omega\ \ \ \ (x_{n} \in GF(2)^{\omega}) の行列である。TGFSR はM系列と同様  2^{p\omega} - 1 の最大周期を持つ、すなわちM系列よりも周期が長く、かつ高速な演算となるような  A を選択できることから、M系列よりも優れる。MT はこの TGFSR を改良したものであり、次の漸化式  x_{n+p}=x_{n+q}+x_{n+1}B+x_{n}C\ \ \ \ (p \gt q \gt 0) \tag{2} によって発生する。式  (2) による最大周期が  2^{19937} - 1 である時、特に MT19937 と呼ばれる*7

*1:コードと数式との対称性からプログラミング言語 Haskell を利用しているが、筆者はプログラミング言語 Haskell に関して初心者であるので、より良い方法などあればコメントして頂けると幸いである

*2:c を 0 としているので厳密には乗算合同法である

*3:この場合周期は 6

*4:下位ビットのランダム性が低い、多次元で規則的な分布となるなど。線形合同法が一次の漸化式であるため多次元疎結晶構造となることに起因する。

*5:原始多項式である事に起因する。

*6:有限体。有限個の元から成る有限集合(すなわち体の公理を満たす)を言う。位数は素数の冪乗となる。

*7:TMP によるコンパイル時 MT の実装などを以前に実装している。テストはこちら