Raspberry Pi 4×4台でKubernetesクラスタを構成するためにしたこと(Kubernetes v1.20 with kubeadm)

概要

Raspberry Pi 4台でkubeadmを使用したシングルコントロールプレーンクラスターの作成を参考にKubernetesクラスタを構成する際に躓いたところを中心に記録した記事です。(2021年1月3日時点)

構成

  • Raspberry Pi 4 ModelB 8GB:4台(Master 1台、Worker3台)
    • SDカード:32GB
  • kubeadm v1.20(GitVersion:v1.20.1)
  • コンテナランタイム:containerd(1.4.3)
  • ネットワークアドオン:Flannel

Raspberry Piの設定

以下、Kubernetesクラスタ構成時特有の設定ではないが、計4つのRaspberry Piに行います。

  1. OS書き込み
  2. ホスト名変更
    • ディスプレイ、キーボードを接続するなどしてログインし、raspi-configからホスト名を変更(System Optionsより)
  3. SSH設定
    • raspi-configから設定。(Interface Optionsより)
  4. kubeletのためにswapをオフにする
$ sudo dphys-swapfile swapoff
$ sudo systemctl stop dphys-swapfile
$ sudo systemctl disable dphys-swapfile
  1. cgroupを有効にするため、/boot/cmdline.txt を次のようにする(cgroup_enable以降を加筆。改行せず、1行に記載する)
console=serial0,115200 console=tty1 root=PARTUUID=eec83cbb-02 rootfstype=ext4 elevator=deadline fsck.repair=yes rootwait cgroup_enable=cpuset cgroup_memory=1 cgroup_enable=memory

kubeadmインストール

kubernetes.io より、 Rasberry Pi OSはDebian(Buster)系のため、iptables設定として

net.bridge.bridge-nf-call-ip6tables = 1
net.bridge.bridge-nf-call-iptables = 1

となるようにし、加えてnftablesバックエンドを使わないように切り替えました。

また、 /sys/class/dmi/id/product_uuid が存在しなかったため、ユニーク性の検証は省略しています。

containerdインストール

Kubernetes v1.20から、採用するコンテナランタイムの選択肢としてのDockerが非推奨となったため、今回はcontainerdを選択しました。 kubernetes.io

kubernetes.io より、Ubuntu 16.04の手順に倣ってインストールします。

ただし、Raspberry Piではadd-apt-repositoryの実行時に

aptsources.distro.NoDistroTemplateException: Error: could not find a distribution template for Raspbian/buster

のようなエラーが出ること、CPUアーキテクチャが異なること、ubuntuではなくdebianを指定するため、以下のように一部手順を変更しました。

# curl -fsSL https://download.docker.com/linux/ubuntu/gpg | apt-key add -
curl -fsSL https://download.docker.com/linux/debian/gpg | apt-key add - # GPG鍵追加

# add-apt-repository \
#    "deb [arch=amd64] https://download.docker.com/linux/ubuntu \
#   $(lsb_release -cs) \
#   stable"

 echo  "deb [arch=armhf] https://download.docker.com/linux/debian\
$(lsb_release -cs) \
 stable" > /etc/apt/sources.list.d/docker.list #リポジトリ設定追加

cgroupドライバーの設定

cgroupドライバーをcgroupfsではなくsystemdに変更します。 必須の項目ではないですが、

systemdと一緒に cgroupfs を使用するということは、2つの異なるcgroupマネージャーがあることを意味します。 コントロールグループはプロセスに割り当てられるリソースを制御するために使用されます。 単一のcgroupマネージャーは、割り当てられているリソースのビューを単純化し、 デフォルトでは使用可能なリソースと使用中のリソースについてより一貫性のあるビューになります。

ということを踏まえて、systemdへ統一します。

ここで、containerdがsystemdのcgroupドライバを使う方法について、

systemdのcgroupドライバーを使うには、/etc/containerd/config.toml内でplugins.cri.systemd_cgroup = trueを設定してください。

について具体的にどの項目を指すかというところに躓きました。

containerd config default > /etc/containerd/config.toml で得られた設定ファイルについて、systemd_cgroupという項目があるため、当初はtrueに変更すればよいように思われたのですが、実際はdeprecatedになっていました。

 [plugins."io.containerd.grpc.v1.cri"]
    disable_tcp_service = true
    stream_server_address = "127.0.0.1"
    stream_server_port = "0"
    stream_idle_timeout = "4h0m0s"
    enable_selinux = false
    selinux_category_range = 1024
    sandbox_image = "k8s.gcr.io/pause:3.2"
    stats_collect_period = 10
    systemd_cgroup = false # ここかと思いきやここではない

github.com

runtiime_type= io.containerd.runc.v2 の場合は次のように設定します。

      [plugins."io.containerd.grpc.v1.cri".containerd.runtimes]
        [plugins."io.containerd.grpc.v1.cri".containerd.runtimes.runc]
          runtime_type = "io.containerd.runc.v2"
          runtime_engine = ""
          runtime_root = ""
          privileged_without_host_devices = false
          base_runtime_spec = ""
          [plugins."io.containerd.grpc.v1.cri".containerd.runtimes.runc.options]
            SystemdCgroup = true # 追加

加えて、kubeletからもsystemdを使用するよう、 /etc/default/kubelet を作成し、

KUBELET_EXTRA_ARGS="--cgroup-driver=systemd"

を記載しました。この作業は計4台全ノードで必要なものです。

クラスタ作成

マスタノードでクラスタを作成し、残りのノードで作成したクラスタに参加します。 まず、マスタノードで kubeadm init としてクラスタを作成します。

Flannelを使うために、kubeadm init 時に --pod-network-cidr=10.244.0.0/16 を渡す必要がありますが、設定項目を保存するため、あらかじめ kubeadm config init-defaults > ~/kubeadm_config.yaml としてから kubeadm init --config kubeadm_config.yaml とすることとしました。デフォルト設定との差分は以下のとおりです。

<   advertiseAddress: 1.2.3.4
---
>   advertiseAddress: <IP address> # masterノードのIPアドレス。worker側でkubeadm joinする際にアクセスするためのものも兼ねる

<   criSocket: /var/run/dockershim.sock
---
>   criSocket: "unix:///run/containerd/containerd.sock"

また、pod-network-cidrは以下の箇所で指定しました。

networking:
  dnsDomain: cluster.local
  serviceSubnet: 10.96.0.0/12
  podSubnet: 10.244.0.0/16 # 追加

最後に、 kubeadm init --config kubeadm_config.yamlクラスタ作成を試みると、以下のエラーで失敗しました。

[ERROR Mem]: the system RAM (1 MB) is less than the minimum 1700 MB

これはkubeadmのメモリ容量計算式のバグで、Raspberry Pi 4 8GBのような32bitシステムかつ4GB以上のメモリを持つシステムでのみ起きる(メモリサイズを誤判定する)挙動のようです。

github.com

そこで、(上記GitHub Issueでも示されているように)

kubeadm init --config kubeadm_config.yaml --ignore-preflight-errors=Mem

とすることで作成できました。
※ Issueに対応するPull Requestがマージ済みのため、後のkubeadmでは --ignore-preflight-errors=Mem は不要になる

最後に、クラスタ作成時に表示された kubeadm join コマンド(トークンを含む)を他の3ノードでも実行して完了です。

結果

kubectl get pod -n kube-system としてapiserverなどの各種コンポーネントが確認でき、実際に kubectl run としてPodが動作したため、ひとまず構築できたようです。
Raspberry Piなのでやはり kubectl get pod のようなリクエストでも重く、時間がかかるようですが、Kubernetesの実験のために使っていきたいところです。

マクドナルドのクーポン番号のヒント問題を手計算で解いてみた

概要

2020年7月31日に突如公開された以下の問題が興味深かったので解いてみました。

今回はプログラムを作成したりWolfram|Alphaのような高性能な計算機(?)を使わずに挑戦してみました。

使う道具:

方針

  •  \pi に床関数、天井関数を施した値を評価する
  • 無限級数の計算を頑張る

 \pi に床関数、天井関数を施した値を評価する

まず、問題中の意味深長な床関数、天井関数を評価します。特筆すべきなのは、  \lfloor \sqrt{{\lfloor \pi \rfloor }^{\pi}} \rfloor \lfloor \sqrt{{\pi}^{\pi}} \rfloor です。

 \lfloor \sqrt{{\lfloor \pi \rfloor }^{\pi}} \rfloor の値

 \lfloor \pi \rfloor = 3より、  \lfloor \sqrt{3^{\pi}} \rfloor = \lfloor  3^{\frac{\pi}{2}} \rfloor が得られます。 次に、値を見積もるために \log_{10} {3^{\frac{\pi}{2}}} の値を考えます。


\displaystyle{ \log_{10} {3^{\frac{\pi}{2}}}  =  \frac{\pi}{2} \log_{10} 3}


\displaystyle{\frac{\pi}{2} \log_{10} 3 \simeq 3.141592 \times 0.4771 \div 2 = 0.7494}

となります。

一方、  \log_{10} 5 = 1 - \log_{10} 2 = 0.6990  \log_{10} 6 = \log_{10} 2 + \log_{10} 3 = 0.7781 であることから、
  \displaystyle{\log_{10} 5 \lt \frac{\pi}{2} \log_{10} 3 \lt \log_{10} 6} を満たすため、
 \displaystyle{ 5 \lt 3^{\frac{\pi}{2}} \lt 6  }、つまり  \lfloor \sqrt{3^{\pi}} \rfloor = \lfloor  3^{\frac{\pi}{2}} \rfloor = 5 が得られます。

 \lfloor \sqrt{{\pi}^{\pi}} \rfloor の値

前節と同様に、 \lfloor \sqrt{{\pi}^{\pi}} \rfloor の値を見積もるために  \log_{10} {\pi^{\frac{\pi}{2}}} の値を考えます。 今、  3.14 \lt \pi \lt 3.15 であることから  {3.132}^{\frac{\pi}{2}} \lt {\pi}^{\frac{\pi}{2}} \lt {3.15}^{\frac{\pi}{2}} を満たすため、  \log_{10} {3.14^{\frac{\pi}{2}}} の値と  \log_{10} {3.15^{\frac{\pi}{2}}} の値を考え、整数部分の大小を見積もります。

 \displaystyle{ \log_{10} {3.14^{\frac{\pi}{2}}}} について

常用対数表によれば、 \log_{10} {{3.14}^{\frac{\pi}{2}}} = \frac{\pi}{2} \log_{10} {3.14} > \frac{3.141592}{2} \times 0.4969 = 0.7805で、  \log_{10} 6より大きいため、 6 \lt 3.14^{\frac{\pi}{2}}を満たします。

 \displaystyle{ \log_{10} {3.15^{\frac{\pi}{2}}}} について

  \log_{10} {3.15^{\frac{\pi}{2}}} = \frac{\pi}{2}\log_{10} 3.15 =\frac{\pi}{2}\log_{10} {\frac{315}{100}} = \frac{\pi}{2} \log_{10} 315 - \pi
 \frac{\pi}{2} \log_{10} 315 - \pi = \frac{\pi}{2} \log_{10} \left(  5\times  3^2 \times 7 \right) - \pi = \frac{\pi}{2} \left( \log_{10} 5 + 2 \log_{10} 3 + \log_{10} 7 \right) - \pi

次に、やや大きい値を使って上から抑えます。1
 \frac{\pi}{2} \left( \log_{10} 5 + 2 \log_{10} 3 + \log_{10} 7 \right) - \pi  \lt \frac{3.1416}{2} \left( 0.6990 + 0.9544 + 0.85 \right) = 0.838

ここで、 \log_{10} 7 = 0.8451 (常用対数表より)であることから、   \log_{10} {3.15^{\frac{\pi}{2}}} \lt \log_{10} 7 を満たすため、   3.15^{\frac{\pi}{2}} \lt 7が言えます。

したがって、
 6 \lt {3.14}^{\frac{\pi}{2}} \lt {\pi}^{\frac{\pi}{2}} \lt {3.15}^{\frac{\pi}{2}} \lt 7 より、  6 \lt {\pi}^{\frac{\pi}{2}} \lt 7 、つまり  \lfloor \sqrt{{\pi}^{\pi}} \rfloor = 6 となります。

無限級数の計算

前節までの結果を使って問題の式を整理すると、

 \displaystyle{ 12 \pi^2 \left( \sum_{\tau = 0}^{\infty}\left( \frac{5}{\left(3\tau + 1 \right)^2}  - \frac{4}{\left(6\tau + 1 \right)^2}\right) \right)^{-1}}

となります。 特に無限和の部分に着目すると、

 \displaystyle{ \sum_{\tau = 0}^{\infty}\left( \frac{5}{\left(3\tau + 1 \right)^2}  - \frac{4}{\left(6\tau + 1 \right)^2}\right) =  \sum_{\tau = 0}^{\infty}\left( \frac{1}{\left(3\tau + 1 \right)^2} + 4\left(\frac{1}{\left(3\tau + 1 \right)^2} - \frac{1}{\left(6\tau + 1 \right)^2}\right)\right) }

と変形できます。 ここで、

 \displaystyle{ \left(\frac{1}{\left(3\tau + 1 \right)^2} - \frac{1}{\left(6\tau + 1 \right)^2}\right) }

について、

f:id:je6bmq:20200801141349p:plain
tauが偶数の場合に打ち消しあう

 \tau = 2k  k自然数 )の場合に  \tau = k の項と打ち消しあう項が存在し、残る項の分母が 4^2,10^2,16^2, \ldotsと続くことから、 2

この項は

 \displaystyle{ \frac{1}{\left(6\tau + 4 \right)^2} =\frac{1}{4\left(3\tau + 2 \right)^2} }

となります。

したがって、

 \displaystyle{ \sum_{\tau = 0}^{\infty}\left( \frac{1}{\left(3\tau + 1 \right)^2} + 4\left(\frac{1}{\left(3\tau + 1 \right)^2} - \frac{1}{\left(6\tau + 1 \right)^2}\right)\right) = \sum_{\tau = 0}^{\infty}\left( \frac{1}{\left(3\tau + 1 \right)^2} + \frac{1}{\left(3\tau + 2 \right)^2}\right) }

が得られます。 これは自然数の逆数の2乗の和( \zeta \left(2\right) = \frac{\pi^{2}}{6})と3の倍数の逆数の2乗の和との差なので。。。といったところで作問者の@mathlavaさんから同様の内容の投稿がなされていました。

まとめ

常用対数表とゼータ関数は偉大。


  1. 上から抑えるために  \log_{10} 7 \simeq 0.84 \lt 0.85 \pi \lt 3.1416 を使っています

  2. 無限和の順序交換を行っていいことを仮定しています

10進数→2進数変換のアルゴリズムで出力を逆に読む理由

概要

手計算で10進数を2進数に手計算で変換するアルゴリズムとして,筆算を連鎖したあと,剰余を逆順に並べて出力とするものがありますが,なぜ逆順でうまくいくのか考えたことがなかったので,メモがてらまとめてみます.

概要としては,

  • ビットマスクとビットシフト
  • 数式

の両面から納得することができました.

アルゴリズム

10進数表記の数  aを2進数表記するとして,この記事で対象にするアルゴリズムは以下のとおりです.

  1.  aを2で割り,剰余と商a'を記録する
  2. a a' を代入し, aが2を下回るまで,これを繰り返す
  3. 剰余を逆順に左から並べる

 a=29 の場合を例に図解すると以下のようになります.

f:id:je6bmq:20190507204803p:plain
2進数変換の例(29の場合)

Python(3.6.2)のコード例は次のようになります.(実行例

bits_list = []
while a > 0:
    bits_list.append(a % 2)
    a = a // 2

bits_list.reverse()

この記事では,割り算の方向に対して出力の読み上げが逆方向に行われる理由について考えていきます.

準備

2進数で n桁の数は,式のように表されます.

 \sum_{i=0}^{n-1}  d_i  \times 2^i

 d_i は, i番目の桁の値を表し, d_i=0または d_i=1です. 例えば, a=29の場合は d_0=1 d_1=0 d_2=1 d_3=1 d_4=1となります.

もちろん,2進数,10進数は同じ数の違う表現というだけなので,

 2 \times 10^1 + 9 \times 10^0 = 1 \times  2^4 + 1 \times  2^3 + 1 \times  2^2 + 0 \times  2^1 + 1 \times  2^0

が成り立ちます.
(左辺:10進数表現,右辺:2進数表現)

ビットマスクとビットシフト

このアルゴリズムの要点は,2で割った場合の商と剰余を使い続けることにあります.ここで,2で割ることと2で割った場合の剰余をビット演算で表現することを考えます. このとき,2で割ることは 右に1bitシフトすること に,2で割った場合の剰余は 0b1との論理積 に対応します.

2で割った場合の剰余について

右シフトで商が得られることを用いると,ある数bについて2で割った剰余は a - (a>>1)<<1で求められます.
(正の数において,シフト時に足りないbitは0が補われるため)
ここで,最下桁以外の各桁の値はaと(a>>1)<<1の間で必ず等しくなるので,結果として0b1との論理積と同義です.

剰余を求めて記録し,2で割る操作は,まるで だるま落としのように,最下層の桁を取り出し,桁を落とす操作に対応します. したがって,最下桁が最初に取り出され,最上位桁が最後まで残るため,逆順に出力されます.

f:id:je6bmq:20190507215330p:plain
最下桁から順番に取り出される

ビット演算を用いて前述のコードを書き直すと,次のようになります.(実行例

bits_list = []
while a > 0:
    bits_list.append(a & 0x1)
    a = a >> 1

bits_list.reverse()

数式

2で割った場合の商と剰余を求めることは,数式から2をくくりだして,2 \times 商 + 剰余 の形式にすることと同義です. そこで,

 \sum_{i=0}^{n-1} d_i \times  2^i

を変形すると,

 =d_{n-1} \times 2^{n-1} + d_{n-2} \times 2^{n-2} + \cdots  + d_1 \times 2^1 + d_0 \times 2^0  =2 \times \left(d_{n-1} \times 2^{n-2} + d_{n-2} \times 2^{n-3} + \cdots + d_1 \times 2^0 \right) + d_0 \times 2^0

となります.ここで,剰余  d_0 \times 2^0を無視して(アルゴリズム上は記録することになりますが), この商の部分に着目すると,再び2でくくることができます. したがって,  =2 \times \left( 2 \times \left (d_{n-1} \times 2^{n-3} + d_{n-2} \times 2^{n-4} + \cdots + d_2 \times 2^0 \right)+ d_1 \times 2^0 \right) + d_0 \times 2^0

となります.

このように,剰余の項を無視して2でくくりだす,を繰り返すことができ,これが 最下桁の値から 剰余として分離することに相当するため,出力を逆順に読むことになります.

まとめ

10進数→2進数変換アルゴリズムの出力を逆順に読み上げる理由について,

意味付けを行いました.

補足

個人的には再帰で書くほうが好みです.(実行例

# start
def convert_to_bin(n, acc_list=[]):
    return reversed(acc_list) if n == 0 else convert_to_bin(n // 2, acc_list + [n % 2]) 
# end

a = 29 
result = "".join([str(b) for b in convert_to_bin(a)])

print(result)

ユークリッドの互除法を多項式に拡張して,(Rustで)実装する

はじめに

この記事はTUT Advent Calendar 2018 三日目の記事(になる予定だった記事)です. 遅れてしまい申し訳ありません.

二日目はT_T3_M5さんの「けんきうのおはなし」でした. 生き物はどうやってものを「見て」いるか?という観点でしたが, 「人は見たいものを見たいように見る」みたいな言葉を思い出しました.

自己紹介

今年は二つのTUT合同でのAdCということで,開催おめでとうございます.
(そういえば,確かフィンランドの方にもTUTってありましたよね.)

はじめまして(?),「むしぱん」です.東京じゃない方の人間です. プログラミング,数学が好きです.
よく使うプログラミング言語はRust,Juliaです.最近は全然書いてないですが,Kotlin,Elixirも好きです.

私も「けんきうのはなし」をするか?と思いましたが用意に身元を特定できそうなので避けようと思います.

私はプログラミングと数学が好きなので,この記事では数学とプログラミング,両方に関わる話題として, (自然数の)最大公約数を求めるアルゴリズムユークリッドの互除法」を多項式に拡張し,(一般化して)Rustで実装した話をします.

ユークリッドの互除法自然数

そもそも,一般に言われるユークリッドの互除法自然数の最大公約数(Greatest Common Divisor)を求めるアルゴリズムです.

自然数{a,b}に対して,以下の式が成り立つとします.つまり, {a}{b}で割ったときの商が{q}で,あまりが{r} ということです.

{a= bq + r}

ここで,最大公約数を求める関数を{\mathrm{gcd}(x,y)}とおくと,以下の式が成り立ちます.

{\mathrm{gcd}(a,b)=\mathrm{gcd}(b,r)}

ここで重要なのは,常に{r\leq b}が成立する点です. 上の式を繰り返し適用することで最大公約数を求めます.このとき,{r\leq b}であることから,徐々に{a,b}の値が小さくなり,最終的にあまり{r=0}となります. あまりが0となったときの{b}の値が最大公約数 となります.

例えば,{a=28,b=6}の例を考えます.最終的に4割る2を計算することとなり,あまりが0であることから,最大公約数は2となります.

f:id:je6bmq:20181204000324p:plain
ユークリッドの互除法による最大公約数の導出({a=28,b=6}

これをプログラムとして記述すると,次の様に記述できます.
Rustにおける実装例ですが,他の言語でも同様の実装になるかと思います.
{a \leq b}だとしても,一度{q=0}となって大小関係が逆転するのみなので,大小関係のチェックは不要です)

fn gcd(a: u32, b: u32) -> u32 {
     if b == 0 { a } else { gcd(b, a % b) }
}

多項式における最大公約数って?

この記事の目的は,この「ユークリッドの互除法」を多項式に応用することでした.今回は各係数が実数の,{x} についての多項式を考えます.
例えば,{ x^2+1}や,{x^3 +\sqrt{2}x + 1}が挙げられます.

この「({x}についての多項式)における最大公約数」を考える前に,一般的な(自然数における)最大公約数とはどういう数か考えます.

28と6の最大公約数,280と588の最大公約数({\mathrm{gcd}(288,588)=28})を例に, 両整数の素因数分解と最大公約数に着目すると, 実は, 最大公約数は共通の素因数の積 になります.

f:id:je6bmq:20181204004508p:plain
最大公約数は,共通の素因数の積になる

多項式における最大公約数も同様に, 多項式素因数分解(に対応する操作をした場合の)共通の素因数(に対応するもの)の積 である,とします.

では,この多項式における素因数分解(に対応する操作)と素因数(に対応するもの)は何なのでしょう.
そもそも,自然数における素因数分解は,最大限細かく素数の積で表す操作(分解)でした.そして,素数は,1と自身以外を約数に持たない数,すなわち,1と自身の積以外に分解できない数でした.

多項式における「分解」として,一番最初に学校で習う操作があります.「因数分解」です.因数分解は,ある多項式を最大限細かく別の多項式の積で表す操作でした.
実はこの因数分解したときの各多項式(例:{x^2-1=(x+1)(x-1)}における{(x+1)}{(x-1)})には「既約多項式」という名前がついています.

そして,自然数における素数素因数分解が,多項式における既約多項式因数分解に対応しています. (厳密には対象にしている集合が整数か,自然数か,実数か,によってある多項式が既約多項式であるかどうかは変わってきますが…)

したがって,多項式における最大公約数は,多項式で共通の既約多項式の積となります.

f:id:je6bmq:20181204012259p:plain
自然数多項式の「分解」の対応

ユークリッドの互除法多項式

では,自然数におけるユークリッドの互除法多項式に拡張します.つまり,これによって多項式の最大公約数(共通の既約多項式の積)が求められればよいということです.

多項式{f(x),g(x)}に対して,以下の式が成り立つとします.つまり, {f(x)}{g(x)}で割ったときの商が{q(x)}で,あまりが{r(x)} ということです.

{f(x)= g(x)q(x) + r(x)}

自然数同様,最大公約数を求める関数を{\mathrm{gcd}(f(x),g(x))}とおくと,以下の式が成り立ちます.

{\mathrm{gcd}(f(x),g(x))=\mathrm{gcd}(g(x),r(x))}

例えば,次の二つの多項式について,互除法を適用します.各々の因数分解から,最大公約数は{(x-1)(x^2+1)=x^3-x^2+x-1}であることがわかります.

 {f(x) =2 x^5 - 7 x^4 + 4 x^3 - 4 x^2 + 2 x + 3= (x-1)(x-3)(2x+1)(x^2+1)}

 {g(x)  = x^4 -1 = (x-1)(x+1)(x^2+1)}

では,互除法を使って最大公約数が求めます.

f:id:je6bmq:20181204145738p:plain
多項式の最大公約数を互助法によって求める例

互除法によると,最大公約数は{4x^3-4x^2+4x-4}だそうです.
あれ?{(x-1)(x^2+1)=x^3-x^2+x-1}ではない?と思うかもしれませんが実は{4x^3-4x^2+4x-4=4(x-1)(x^2+1)}であるため,両多項式とも{4x^3-4x^2+4x-4}によって割り切れるのです.

このように,多項式では最大公約数の定数倍の多項式ではすべて割り切れてしまうため,実際は既約多項式の最高次の係数を1に正規化して扱うことがあります.
例えば,先ほどの互除法の例の最後の商,{\frac{1}{4}x+\frac{1}{4}}から{\frac{1}{4}}をくくり出して,{\frac{1}{4}(x+1)}と扱えば,定数倍の差を吸収できます.

最大公約数関数を多項式に拡張,実装する(Rust)

ユークリッドの互除法を行う上で自然数多項式共通の操作は四則演算と剰余演算,余りがゼロであるという判定です. つまり,四則演算と剰余演算ができて,かつゼロが定義されていれば互除法が適用できることが示唆されます.

そこで今回は 新たに多項式を表すPolynomial型を定義し,Polynomial型を実装し,ゼロの概念と等号判定を導入する ことで自然数多項式両方に適用できる最大公約数関数を実装します. 詳細の実装の説明は省略して,概要のみ説明します.詳細はこちらのGistを参照ください.(思いの外,数式パートが長引いたので…)

Rustでは,トレイト(他の言語でいうインターフェースのようなもの)を用いて関数の型制約を表現できるので,次の様に実装できます.

fn gcd<T>(a: T, b: T) -> T
where
    T: Add<Output = T>
        + Sub<Output = T>
        + Mul<Output = T>
        + Rem<Output = T>
        + Zero
        + PartialEq
        + Clone,
{
    if b == <T as Zero>::get_zero() {
        a
    } else {
        gcd(b.clone(), a % b)
    }
}

ここで,型TはAddトレイト,Subトレイト,Mulトレイト,Remトレイト,Zeroトレイト,PartialEqトレイトを実装している必要があることを意味します.そして,各トレイトは次の意味を持ちます.

  • Addは加算可能,かつ計算結果の型がT
  • Subは減算可能,かつ計算結果の型がT
  • Mulは乗算可能,かつ計算結果の型がT
  • Remは剰余演算可能,かつ計算結果の型がT
  • PartialEqは == によって同値関係が確かめられること

上記のトレイトは標準ライブラリに搭載されているので,自然数(符号なし整数型)には当然実装されています.したがって,新たにPolynomial型にも実装することとしました.

新たに定義したのはZeroトレイトです.Zeroトレイトは,名前の通りその型における ゼロ を定義したいため,次の様な実装にしました.
符号なし整数(ここでは32bitのみ対象)では0をそのまま返し,Polynomial型においては0次の項の値が0の多項式を返します.

impl Zero for u32 {
    fn get_zero() -> u32 {
        0
    }
}
impl Zero for Polynomial {
    fn get_zero() -> Polynomial {
        Polynomial::new(&vec![0f32], false)
    }
}

実行結果

一般化したgcd関数はもちろん32bit符号なし整数型にもそのまま適用可能ですが,ここでは先ほどの二つの多項式

 {f(x) =2 x^5 - 7 x^4 + 4 x^3 - 4 x^2 + 2 x + 3= (x-1)(x-3)(2x+1)(x^2+1)}

 {g(x)  = x^4 -1 = (x-1)(x+1)(x^2+1)}

の最大公約数を計算した結果を示します.ここで,^ (ハット)の後ろの数字が指数を表します.

実際に,先ほどと同様に最大公約数を求めることができています.

f:id:je6bmq:20181204153247p:plain
一般化したgcd関数(Rust)の実行結果

まとめ

この記事では,最大公約数を求めるアルゴリズムユークリッドの互除法自然数から多項式に拡張して,かつプログラミング言語Rustで実装できることを確認しました. 少なくともRustではトレイトで型の制約を表現できるので, 数学における概念と満たすべき性質 がうまく対応付けられて嬉しいと思っています.

多項式が既約多項式かどうか,といったような概念は有名な理論である「 ガロア理論 」とも関わりが強く, 私がこの実装を試みたきっかけにもなっているので,興味のある方は調べてみてはいかがでしょう.(というか私も教えてほしい.)

本来の四日目の担当はyanteneさんで,「何書くかまだ何も考えてない!遅れたらごめんなさい!」とのことです. この記事自体が遅れてしまったので何も言えませんが期待しています.

Juliaの構造体とジェネリクスについていろいろ勘違いしていた話

Juliaについて

最近、数値計算用のプログラミング言語、Juliaにはまっています。
Juliaは、動的型付けでありながら、ユーザ(プログラマ側)が型を明示したり、関数単位のJIT(Just-In-Time)コンパイルに伴う最適化によって高速に実行可能であるため、注目されています。

Juliaにおいても、多言語のジェネリクスのように静的型付け言語にあるようなサブタイプ、スーパータイプを指定することによる型の制約を設けることが可能です。 (2018年5月29日時点、Juliaのバージョンは0.6) 今回はJuliaの構造体とジェネリクス(?)に少々はまったので記しています。

コンストラクタにおける勘違い

まず、以下のようなコードを考えます。
ここで、TはAbstractFloatのサブタイプであることを意味しており、 Hoge{T}(_a::T) = new(_a) はコンストラクタの定義を表しています。 また、AbstractFloatは64bit浮動小数型であるFloat64、32bit浮動小数型であるFloat32のスーパータイプとなっています。

struct Hoge{T<: AbstractFloat}
      a::T
      Hoge(_a::T) = new(_a)
end

このコードを実行すると、型パラメータと制約を明記するよう、警告されます。 Juliaでは型パラメータTをそのまま利用できるわけではないようです。これは 構造体の定義外にもコンストラクタを定義できることから、構文を統一されるためのもの だと推測しています。
(Juliaのドキュメントでは、コンストラクタの定義位置が構造体定義の内外に応じて、Innner、Outerと区別されている)

WARNING: deprecated syntax "inner constructor Hoge(...) around C:\ProjectYukito\julia\famiconlike_sound_spectrum.jl:21".
Use "Hoge{T}(...) where T" instead.

そこで、次のように変更します。これにより、警告が無くなりました。

struct Hoge{T<:AbstractFloat}
    a::T
    Hoge{U}(_a::U) where U<:AbstractFloat = new(_a)
end

ジェネリクスにおける勘違い

次に、Hoge型の変数を宣言することを考えます。

hoge = Hoge(0.5)

一見問題ないように見えますが、 このコードは、以下のエラーのため実行できません。

ERROR: LoadError: MethodError: Cannot `convert` an object of type Float64 to an object of type Hoge
This may have arisen from a call to the constructor Hoge(...),
since type constructors fall back to convert methods.

=コンストラクタがconvert関数として扱われるが、convert関数がFloat64→Hogeの変換に対応していない(?)

コンストラクタが引数の型に対応していない場合、convert関数でエラーを出力することは型変換に関するJuliaのドキュメントにも記されています。

Conversion and Promotion · The Julia Language

問題は、 Float64に対してHogeのコンストラクタが対応していない点 です。
Float64はAbstractFloatのサブタイプであるはずなので、エラーが出力されることが解せませんでした。そこで、Float64であることを明示してみます。

hoge = Hoge{Float64}(0.5)

これにより、初めて正常に実行されました。 ちなみに、Float64の代わりにAbstractFloatのサブタイプであるFloat32を指定すると、同様に cannnot convert ... のエラーが出力されます。一方で、AbstractFloatのサブタイプでない型(UInt32など)を指定すると、型制約に伴うエラーが出力されます。

これらのことから、 Juliaでは、引数の型に応じて推論の結果型パラメータの値が定まるのではなく、型を明示する必要がある という結論に至っています。

まとめ

普段Rustばかり書いているせいか、Juliaのジェネリクスについて混乱してしまいました。

  • 内部コンストラクタ(Innner-Constructor)において、構造体の定義に用いる型パラメータは直接適用されない
  • 変数の宣言時、型パラメータは明示する必要がある(引数の型推論から確定するわけではない?)

WindowsでSource Code Pro をビルドしようと思うも一筋縄ではいかなかった話

最近、エディタをSpacemacsに乗り換えました。 Spacemacsのデフォルト設定では、フォントがSource Code Proなので、Source Code Proが入っていない場合は起動後にWarningが出ます。
おそらく代替フォントで描画されているかと思いますが、特に日本語の描画が大変厳しいようで、お世辞にも綺麗とは言えない状況でした。

そこで、Source Code Proを導入することとしました。

Source Code Pro(GitHub)にもあるように、 リリースページからダウンロードすればよい のですが、 github.com 見事にそのリンクの情報を見逃した私は、ソースコードからビルドすることにしました。

基本的にはREADMEにしたがって行えばよいのですが、私の環境(Windows10)ではbuild.cmdの実行時にエラーが発生しました。

syntax error at "}" missing ";" 

当該リポジトリのIssueにも同様の事例が報告されていませんでしたので、困惑していたのですが、別リポジトリのIssueを参考に解決することができました。 つまり、familyTables.feaの2行目のセミコロンを削除する ことで対処できました。

github.com

いやわざわざビルドしなくていいのに、と言われればそれまでなのですが。。。

Rustで数値シミュレーションをする際に使用しているライブラリ(クレート)の紹介

この記事はTUT Advent Calendar 2017の18日目の記事です。

概要

普段私はRustを使ってシミュレータを作成しているのですが、そこでお世話になっているライブラリ(クレート)を簡単に紹介します。
今回紹介するライブラリ(用途)は以下の2つです。

  • rayon (データ並列化)
  • statrs (様々な確率分布に従う乱数生成)

rayon

すでにいくつかの記事で紹介されていますが、rayonによって主に以下の処理を簡潔に行うことができます。両者とも、Work stealingと呼ばれる方法で実現されているそうです。

  • スレッドによる並列処理
  • 再帰関数の分割処理

スレッドによる並列処理

rayonでは、Parallel Iteratorの導入によって通常のイテレータに対する関数(mapなど)の処理を並列化することができます。D言語でいうstd.parallelismに似ているかもしれません。 ただし、Rustのシステムによって並列時の安全性は保証されています。つまり、標準ライブラリのスレッド同様、Sendトレイト、Syncトレイトの制約によってコンパイル時に検査が行われます。

iter関数やinto_iter関数が実装されている型は基本的にParallel Iteratorに変換することができ、以下のような対応関係を持ちます。

コード例は以下のとおりです。

extern crate rayon;
use rayon::prelude::*;

fn main() {
    (0..100)
    .into_par_iter()
    .map(|i| i * 2)
    .for_each(|x| println!("{}",x));
// 200未満の偶数が出力される
}

大変便利ですが、Parallel Iteratorに対してcollect::<Vec<_>>()のようなことはできず、 あらかじめ格納先のベクタをミュータブルとして宣言したのち、collect_into関数を用いて格納する必要があります。

コード例は以下のとおりです。

extern crate rayon;
use rayon::prelude::*;

fn main() {
    let mut even_vec = vec![0;100];
    (0..100)
    .into_par_iter()
    .map(|i| i * 2)
    .collect_into(&mut even_vec);
// 200未満の偶数がeven_vecに格納される
}

並列処理時のスレッド数は自動的に決定されますが、スレッド数を能動的に設定したい場合は、num_threads関数によって設定可能です。
また、実行時にスレッド数を設定したい場合は、環境変数RAYON_NUM_THREADSの値を設定してから実行することで対応できます。

再帰関数の分割処理

join関数によって、再帰関数を効率的に並列処理します。
rayonのGithubリポジトリの例が簡潔であるため、引用します。 コード例は、並列化したクイックソートのようです。

fn quick_sort<T:PartialOrd+Send>(v: &mut [T]) {
    if v.len() <= 1 {
        return;
    }

    let mid = partition(v);
    let (lo, hi) = v.split_at_mut(mid);
    rayon::join(|| quick_sort(lo), || quick_sort(hi));
}

リポジトリのREADMEによれば、

Note though that calling join is very different from just spawning two threads in terms of performance. 
This is because join does not guarantee that the two closures will run in parallel. 
If all of your CPUs are already busy with other work, Rayon will instead opt to run them sequentially. 
The call to join is designed to have very low overhead in that case, so that you can safely call it even with very small workloads (as in the example above).

とのことですので、join関数の引数のクロージャを、可能なら並列処理し、CPUの負荷が高い場合は逐次的に処理するよう、うまくスケジューリングされる(?)ようです。

statrs

statrsはRustにおける統計的な数値計算のためのクレートで、特に私はdistributionモジュールを頻繁に利用しています。distributionモジュールを利用することで、各種確率分布に従う乱数を容易に生成することができます。
ただし、rust statrsで検索すると検索エンジンが変に気を利かせて別の検索ワードを推奨してくるので注意です。

離散型分布、連続型分布問わず様々な分布に対応しており、以下の分布ごとに構造体が定義されています(非常に多い)。また、分布の期待値、分散を返す関数が各分布に実装されており、離散型確率分布に対しては確率質量関数の値、連続型確率分布に対しては確率密度関数の値を返す関数が実装されています。

さらに、sample関数(rand::Rngトレイトを実装する乱数生成器を引数)の実装によって、各分布に従う乱数を生成可能です。 rand::Rngトレイトを実装している乱数生成器であればよいため、StdRngやThreadRng、XorShiftなど様々なものに対応しています。

  • 離散型確率分布
    • ベルヌーイ分布
    • 二項分布
    • カテゴリカル分布(Categorical Distribution)
    • 離散一様分布
    • 幾何分布
    • 超幾何分布
    • 多項分布
    • ポワソン分布
  • 連続型確率分布
    • ベータ分布
    • コーシー分布
    • カイ分布
    • カイ2乗分布
    • ディリクレ分布
    • アーラン分布
    • 指数分布
    • F(Fisher-Snedecor)分布
    • ガンマ分布
    • 逆ガンマ分布
    • 対数正規分布
    • 正規分布
    • パレート分布
    • t分布
    • 三角分布
    • 連続一様分布
    • ワイブル分布

本記事ではdistributionモジュールにのみ着目しましたが、他のモジュールにてベータ関数やエラー関数などの複雑な関数の値を計算する関数も用意されているようです。

まとめ

rayonやstatrsは、数値シミュレーションをするうえで、大変有用です。並列処理により計算時間を削減することができ、randクレートでは対応していない様々な確率分布に従う乱数生成を行うことで様々なシミュレーションが可能です。
また、Rustはコンパイラが大変優秀なので、コンパイル時の検査のおかげでシミュレーション実行時のエラーがなく、重宝しています。その一方で、コンパイルを通すのは大変ですが。

私の身の回りでもRustが流行ればいいな、と思うこの頃です。