https://ymos-hobby-programing.hatenablog.com/entry/2021/07/17/015738

それでは今日もいきます。

というか、日付変わってるから、実は1日飛ばしちゃったことになるけど気にしない。

本日のお題

漸化式を計算するPiet

まず作ったのが、漸化式を計算するプログラム。

漸化式、何それ、おいしいの？

漸化式は数列 $\{x_n\}$ を作るための規則となる数式です。

もっとも簡単な例で言うと、次のようなものです。

$x_{n+1}=ax_n+b$

この式を使って、n番目の数値 $x_n$ からn+1番目の数値 $x_{n+1}$ を求めます。

これを繰り返すと、数列 $\{x_n\}$ が得られます。やったね！

自動生成

例によって「GridPietGenerator」で生成します。まず、処理フローファイル「reccur.txt」は次のようになります。

inn inn inn      #[a,b,x
:reccur          #[a,b,x
push3 push2 roll #[b,x,a
dup              #[b,x,a,a
push4 push3 roll #[x,a,a,b
dup              #[x,a,a,b,b
push5 push3 roll #[a,b,b,x,a
mul add          #[a,b,a*x+b(=>x)
:print
outn pop pop     #[
end

あらシンプル。

これは、 $a$ 、 $b$ 、 $x_n$ の順に値を入力すると、

上式から次の値 $x_{n+1}$ を算出し、出力してくれます。

Pietに変換するとこんな感じ。

漸化式計算プログラム１

npietで実行すると、出力は次のようになります。

./npiet reccur.ppm
? 2
? 1
? 3
7

手計算するとわかりますが、確かにあってますね。

また、トレース画像はこんな感じになります。

漸化式計算プログラム１のトレース画像

やはり単純ですね。

漸化式を次々と計算するPiet

さっきは $x_n$ から $x_{n+1}$ を計算するという、

1回分の計算しかしませんでした。

ですが、せっかくなので、数列 $\{x_n\}$ の変化をもう少し追ってみたいところ。

そこで、先ほどのプログラムに改良を加え、

漸化式を続けて計算していき、数列 $\{x_n\}$ を出力してみたいと思います。

処理フローファイル「reccur2.txt」は以下の通りです。

inn inn inn inn  #[n,a,b,x
:reccur_loop     #[n,a,b,x
dup outn         #[n,a,b,x ==> x
push4 dup mul push2 mul outc #[n|a,b,x ==>' '
:reccur          #[...|a,b,x
push3 push2 roll #[...|b,x,a
dup              #[...|b,x,a,a
push4 push3 roll #[...|x,a,a,b
dup              #[...|x,a,a,b,b
push5 push3 roll #[...|a,b,b,x,a
mul add          #[...|a,b,a*x+b(=>x)
:reccur_loop_jud #[n|a,b,x
push4 push3 roll #[a,b,x,n
push1 sub        #[a,b,x,n-1(=>n)
dup if::end      #[a,b,x,n
push4 push1 roll #[n,a,b,x
goto:reccur_loop #[n,a,b,x
:end
pop pop pop pop  #[
end

5-11行目は、1つ目のプログラムの一部を拝借しました。

このように、作ったプログラムの一部を、他のプログラムに使い回すことも可能です。

この時に注意することは、使いまわした部分の最初と最後で、

スタックの状態がきちんと管理されていることです。

今回の場合は、初めも終わりもスタックトップから順に

$x_n$ 、 $b$ 、 $a$ が順に格納されており、

これを前提として他の部分を実装しました。

Pietソースコード

GridPietGeneratorで生成された画像は次のようになります。

漸化式を次々と計算するプログラム

これをnpietで実行したときの結果は次のようになります。

./npiet reccur2.ppm
? 10
? 2
? 1
? 3
3 7 15 31 63 127 255 511 1023 2047

はじめにPietから入力を要求されるので、入力します。

入力すべき値は、上から順に、添字 $n$ の最大値 $N$ 、 $a$ 、 $b$ 、及び、数列 $\{x_n\}$ の初期値 $x_0$ です。

つまり、このプログラムは、冒頭の漸化式と初期値 $x_0$ で生成される数列 $\{x_n\}$ を、

$N$ 番目の値まで計算して表示します。

つまり、表示されるのは、 $x_0$ 、 $x_1$ 、・・・、 $x_{N-1}$ の $N$ 個の値です。

上の出力結果を見ると、なんとなくあっている気はします。

トレース画像は次のようになります。

漸化式を次々と計算するプログラムのトレース画像

ループが入ったので、多少処理が複雑になりました。

疑似乱数生成

さて、なぜこんな小難しい話をしたかというと、

Pietで乱数を発生させたかったからです。

乱数発生のアルゴリズムは色々ありますが、実装が簡単そうな「線形合同法」を実装します。

線形合同法では、次の漸化式を使って乱数を発生させるそうです。

$x_{n+1}=(ax_n+b)~\mathrm{mod}~m$

ただし、 $m>a$ 、 $m>b$ 、 $a>0$ 、 $b\geq{0}$ だそうです。

このように、乱数発生の際に漸化式をつかうため、

あらかじめ似たような形の漸化式を計算するプログラムを作っておいたのです。

Pietソースコード生成

それでは処理フローファイル「rand1.txt」をば。

inn inn inn inn inn  #[m,n,a,b,x
:rand_loop           #[m,n,a,b,x
push5 push4 roll     #[n,a,b,x,m
dup                  #[n,a,b,x,m,m
push6 push1 roll     #[m,n,a,b,x,m
mod                  #[m,n,a,b,x%m(=>x)
:print               #[m,n,a,b,x
dup outn             #[m,n,a,b,x ==> x
push4 dup mul push2 mul outc #[m,n|a,b,x ==>' '
:reccur              #[...|a,b,x
push3 push2 roll     #[...|b,x,a
dup                  #[...|b,x,a,a
push4 push3 roll     #[...|x,a,a,b
dup                  #[...|x,a,a,b,b
push5 push3 roll     #[...|a,b,b,x,a
mul add              #[...|a,b,a*x+b(=>x)
:rand_loop_jud       #[n|a,b,x
push4 push3 roll     #[a,b,x,n
push1 sub            #[a,b,x,n-1(=>n)
dup if::end          #[a,b,x,n
push4 push1 roll     #[n,a,b,x
goto:rand_loop       #[n,a,b,x
:end
pop pop pop pop  #[
end

これは、先程の「漸化式の値を次々計算する」プログラムにちょっと付け足しをしたものです。

線形合同法では、漸化式計算ののちに剰余計算（mod）が入りますので、

剰余計算するときの「割る値 $m$ 」を入力に追加しました。

また、剰余計算の処理を3-6行目に追加しています。

こうしてできたPietソースコードがこれです。

乱数を発生させるPietソースコード

npietで動かしてみます。まず処理結果はこのようになります。

./npiet rand1.ppm
? 24
? 100
? 13
? 5
? 1
1 18 23 16 21 14 19 12 17 10 15 8 13 6 11 4 9 2 7 0 5 22 3 20 1 18 23 16 21 14 19 12 17 10 15 8 13 6 11 4 9 2 7 0 5 22 3 20 1 18 23 16 21 14 19 12 17 10 15 8 13 6 11 4 9 2 7 0 5 22 3 20 1 18 23 16 21 14 19 12 17 10 15 8 13 6 11 4 9 2 7 0 5 22 3 20 1 18 23 16

入力した値は、上から順に、「割る値 $m$ 」、添字 $n$ の最大値 $N$ 、 $a$ 、 $b$ 、及び、数列 $\{x_n\}$ の初期値 $x_0$ です。

具体的な値はWikipedia先生の仰せの通りにし、初期値1として100個出力しました。

線形合同法で得られる乱数列 $\{x_n\}$ は、特定の条件が整うと、最大周期が $m$ に一致するそうで、

今回はそのようなパラメータを指定しました。

確かに周期24（ $=m$ ）の数列になっていますが、その間の数値はランダムに変化しているように見えます。

（このように「周期を持つ」などの性質から、線形合同法で得られる数列は、厳密には「乱数列」ではないようです。）

ちなみにトレース画像は次のようになります。

乱数を生成するPietプログラムのトレース画像

今回は全体的に少し控えめな複雑さでした。

終わりに

さらっと書こうと思っていましたが、けっこう重かった。。。

Pietプログラミングしている時間より、ブログ書いている時間の方が長いミステリー。

眠気もすごいけど、もやもやもすごいです。

というのも、なぜ乱数生成プログラムを作ろうと思ったのか、思い出せず・・・。

何か別の用途に使おうとしたはずですが、なんだったかが思い出せません。

そんなモヤモヤに目をつむりつつ、

ブログを途切れさせないためのエネルギーを蓄えるべく、

寝ます。

では。