まずは乱数一つの分布です。サンプル数は100と少なめにしたので、荒い結果になっています。
次に5,10,20です。

実際に試してみると、納得できます。
乱数でなくても任意の関数f(x)でも結果は同じになるようです。

わかりやすくするため、a=-1 固定でbの値をだんだん大きくしてみました。
b = -0.4

右の曲線がだんだん膨らみ、左に小さな島ができます。

b = 0

その島がだんだん大きくなり、

b = 0.3849…

もとの曲線のくっつきます。判別式D=4*a^3+27*b^2がD=0

b = 0.5

そして最後に曲線に島が飲み込まれていきます。
なんか生物的で面白いです。

これだけでなく、楕円加算という曲線と直線の3つの交点(接線の場合は2点)を加算と見立てる演算があります。
これをExcelで演算してみました。(加算を34回順次行い、演算ごとに結果を描画することによって曲線上の点であることを確認)

B列:x, C列:y, D列:直線の傾きL、E列:演算誤差(ゼロになること)の確認(y^2-x^3+x)
L=(y2-y2)/(x2-x2) ※はじめの一点は接線の傾きなので微分でもとめる(L=(3*x^2+a)/(2*y) )

P(x1,y1), Q(x2,y2), R(x3, y3)として
P + Q = R

x3=L^2-x1-x2
y3=L*(x1-x3)-y1
※最初だけ、P + P = Q (x1=x2,y=1=y2とする)

P,Q,R,Sの点について、下記サイトで描画してみました。
https://www.desmos.com/calculator/zvix1zzdjm
加算結果はY軸対象の点に置き換えるのがルールのようです。

楕円曲線について、ちょっとは理解が深まりました。

1/f Noise

決定論的なのにランダムと区別がつかない数列をつくることができるので、毎回同じ振る舞いをする乱数?としてアルゴリズム作曲などに使ったこともあります。(擬似乱数もシードが同じなら同様の振る舞いをしますが、プラットホームが違っても同じ値を使いたいときとか、乱数としての性能は専用ロジックにはかなわないが、演算プロセスをわかりやすいものにしたいときなど有用)
決定論とランダム、最近話題のユヴァル・ノア・ハラリ氏の「ホモ・デウス」の中でもたびたび出てくるキーワードですが、現在でも必須のアイテムではないでしょうか。

0<=a<=1: 0に収束
1<a<=2: 1-1/aに収束
2<a<=3: 1-1/aに振動して収束
3<a<=1+sqrt(6) : 周期2
1+sqrt(6)<a<=3.5699456.. : 周期4,8,16..
3.5699456..<a<=4: 収束も周期性もない複雑な挙動

一部をExcelでテスト
A1: 定数a
B1: 初期値X{0}
B2: =$A$1*(1-B1)*B1
B3:=$A$1*(1-B2)*B2
B4..以下同様 データ数1000

a=2.5

a=3.5

aを面白くなってくる3.0から0.1ステップで4.0まで3Dで表示してみました。
Frequency関数でデータ配列上と同様1000、区間配列は0.01ステップで100で計算(0から1)

3.6 あたりから複雑になり、4では全ての区間で値が存在して両端が多いですがほぼランダムと同じくらいの頻度となります。
ウィキペディアにある図と比べるといい感じになっていることがわかると思います。

https://ja.wikipedia.org/wiki/%E3%83%AD%E3%82%B8%E3%82%B9%E3%83%86%E3%82%A3%E3%83%83%E3%82%AF%E5%86%99%E5%83%8F

地味で実用的な内容でした。

4ケースだけだったので、手作業でソルバを一つずつ設定、計算して使っています。(Dのケース)

セルの計算(Aのケース)

左の列は、支点、第一関節、第二関節の水平に対する角度。
ソルバの条件を変えると違う動きをします。
グラフをアームに見立てるのも面白いと思いました。

PrintWriter writer = createWriter("noise.csv");
background(0);
size(512, 512);

float n1 = 0;
float n2 = 0;
float n3 = 0;

for(int x = 0 ; x < 512 ; x += 1){
  float p1 = noise(n1);
  float p2 = noise(n2);
  float p3 = noise(n3);

  n1 += 2;
  n2 += 0.1;
  n3 += 0.003;
  
  writer.println(p1 + "," + p2 + "," + p3);
 
  noStroke();
  
  fill(255,0,125);
  ellipse(x, p1 * height, 5, 5);
  fill(0,255,0);
  ellipse(x, p2 * height, 5, 5);
  fill(0,255, 255);
  ellipse(x, p3 * height, 5, 5);
}

writer.flush();
writer.close();

noiseのパラメータn1,n2,n3の増加量で、滑らかさを変えることができます。でたらめ、自然ぽい、その中間、といった意図で、パラメータを調節してみました。
そして3つの乱数列をcsvファイルの出力して、Excelで読み込みます。

座標系が逆なので、Excelに読み込むと上下逆になります。

前回と同様、パワースペクトルをとりました。
1/fにそっているように見えます。

またこのバーリンノイズは、2次元、3次元にも適用できるところが面白いです。

float d = 100;
size(512, 512);
noStroke();
background(0);

loadPixels();
for(int y = 0; y < height; y++){
  for(int x = 0; x < width; x++){
    float g = noise(x / d, y / d) * 255;
    pixels[y * width + x] = color(g);
  }
}
updatePixels();

実行するたびに、変化します。

これを自前だとどうやってつくるのだろう、と興味が広がります。

ふとしたことからこの本を見つけ、当時理解したつもりでいた、ホワイトノイズの1/2階積分を実際にテストしてみようと思いました。
前回、1/2階積分の計算方法は、この本を参考にしています。(同じではないですが・・)

上記、まず512個のランダム数(A列)を順にたしたもの(積分)をブラウンノイズ(B列)、それを前回の数値フィルタで1/2階積分したものが(E列)を示します。(係数をソルバーで決めるときに最小二乗法で直線をフィットしますが、元のデータができるだけ直線になるように係数1を決め、それから係数2をあてはめたりして調整します。) 値をとる振幅が小さすぎたので、(A)(B)を5倍にしたものを(H)(I)としました。
ランダムとブラウン運動との中間に、1/fゆらぎがあるように見えます。

次に、パワースペクルトを求めるため(H)(B)(I)をそれぞれフーリエ解析したものが、次の図の(E)(H)(K)、それぞれ複素数の絶対値IMABS()をとったのが、(F)(I)(L)になります。表示対象は、対称性があるので半分の256サンプルにします。

ブロットされたものを見ると、ランダムは水平に分布(1/f^0)し、半積分されたものは、右下に傾いて分布、積分されたものは、さらに傾いて分布しているように見えます。1/f,1/f^2の傾きと比較していみると、積分と半積分、どちらが近いかといえば、積分は1/f^2、半積分は1/fになります。

もっとデータ数をふやせば、精度もあがるかもしれません。

まず半積分をどのように実現するのか、これに悩みました。ScilabやOcatveに便利な関数はないかと調べましたが、見当たりませんでした。
そこで下の2つの式を使ってExcelでやる方法を試してみました。

1). 前のデータを足していくので積分を表現できる
y(t) = ay(t-1) + x(t)

2). 前のデータとの差分を得るので微分を表現できる
y(t) = x(t) + bx(t-1)

上の2つをつないで、1/2回の積分を表現し、これを2回すれば1回の積分と同じになることを見込んでいます。

※Excelシートの数値
A列　: 元データ
B列、D列 : 1)
C列、E列 : 2)
F列　: A列の積分(E列と近いことを確認)
B25 : 係数a
C25 : 係数b
E25 : =LINEST(E2:E21) ここが1になるように、ゾルバーで係数a,bを求める
L列　: スーパークライン(曲がり具合確認用。微分で傾きを調整している)

EとFが近いということから、1/2回積分2回で、積分が実現しているように見える。(直線が半積分で一旦曲線になってさらに半積分でまた直線に戻る)

今度は、この係数を使い、元データを定数から2次曲線に変えてみました。

うまくいっているように、見える。
この方法が正しいかどうかわかりませんが、詳しくはまた次回に。