提出用ではない。
やっぱりアルファベットが並んでいるとプログラミング言語らしい。

乱数をきちんと発生させているつもりなのに、十字の模様がうっすらと見えてしまう。
乱数生成まで自分で書かなければならないかも。

//モンテカルロ
title "モンテカルロ?"

//使用変数(円)
max_in = 1000000     ; 総ドット数(入力用
max = 0;           
in_dot = 0;            ; 内部ドット数
rest = 10000         ; await をいれる間隔
visible = 0         ; 1 で、描画過程を表示する
vx = 390 : vy = 00    ; 表示座標
r = 200                ; 半径
pi = 0                ; 求めるpi

//グラフ
graph_x = 100 : graph_y = 450
graph_w = 600 : graph_h = 150

randomize            ; 乱数初期化
rand_max = 32768    ; 乱数の最大値

//メインウィンドウ初期化
screen 0, 800, 600
//コントロールの配置
objsize 200, 25
pos 0, 25    : combox dimensions, 100, "円\n球"

pos 120, 0    : chkbox "過程描画", visible_chk

pos 0, 55    : mes "ドット数"
pos 80, 55    : input max_in, 120, 25
pos 0, 80    : mes "内部ドット数"
line 0, 100, 100, 100
pos 0, 102    : mes "全ドット数"
pos 100, 90    : mes "x volume="

//数式を読む
pos 25, 200 : picload "eq.jpg", 1

objsize 200, 45
pos 0, 350        : button "Start",*start
obj_start = stat

//軸の描画
line vx-3, vy, vx-3, vy+2*r
line vx, vy+2*r + 3, vx+2*r, vy+2*r + 3
//目盛の描画
pos vx-20, vy : mes "1"
line vx-9, vy, vx-3, vy

pos vx-20, vy+r : mes "0"
line vx-9, vy+r, vx-3, vy+r

pos vx-20, vy+r*2 : mes "-1"
line vx-9, vy+r*2, vx-3, vy+r*2
line vx, vy+r*2+3, vx, vy+r*2+9

pos vx+r-9, vy+r*2+6 : mes "0"
line vx+r, vy+r*2+3, vx+r, vy+r*2+9
pos vx+r*2-9, vy+r*2+6 : mes "1"
line vx+r*2, vy+r*2+3, vx+r*2, vy+r*2+9

//グラフ軸の描画
gosub *init_graph
stop
*start
    dimension = dimensions//これ以降dimensions はない
    max = max_in;
    visible = visible_chk
    //円描画用バッファ
    buffer 2, r*2, r*2;
    cls
    if(visible == 1): gsel 0;
    pos vx, vy
    gcopy 2, 0, 0, r*2, r*2
    //円を作る繰り返し
    gosub *init_graph;
    in_dot = 0
    for i, 0, max, 1
        //必要な値の計算
        rx = rnd(rand_max) - rand_max/2;
        ry = rnd(rand_max) - rand_max/2;
        if(dimension == 1) : rz = rnd(rand_max) - rand_max/2:else:rz = 0;
        len = sqrt(rx*rx + ry*ry + rz*rz)
        px = rx * r / (rand_max / 2);
        py = ry * r / (rand_max / 2);
        //直接描画の場合、描画位置を座標に加算
        if(visible == 1) : px += vx : py += vy ;
        //色の指定
        if (len/(32768/2) < 1){
            in_dot ++
            color 100, 100, 255
        }else{
            color 0, 0, 0
        }
        //描画
        pset px + r, py+ r
        //一定間隔で待機を行い、進行状況・円周率を表示
        if(dimension == 0): vol = 4:else : vol = 8;
        pi = (((1.0*in_dot)/(i+1))*vol)
        if( i \ rest == 0){
            prev = ginfo(3)
            gsel 0
            color 100, 100, 255
            boxf 0, 0, 100, 10
            color 0, 0, 0;
            boxf 0, 0, (i * 100) / max, 10

            color 255, 255, 255
            boxf 180, 90, 300, 115
            color 0, 0, 0
            pos 180, 90
            mes ""+pi//円周率
            gsel prev
            await 10
        }
        //グラフに打点すべきタイミング
        if(max / graph_w < 1){
            if(dimension == 0):def = 3.1416:else:def = 4.1888
            prev = ginfo(3)
            gsel 0
            color 255,0,0
            pset graph_x + i, graph_y + graph_h/2 - pi*500 + def*500
            gsel prev
        }else{
            if(i\(max / graph_w) == 0){
                if(dimension == 0):def = 3.1416:else: def = 4.1888
                prev = ginfo(3)
                gsel 0
                color 255,0,0
                pset graph_x + (1.0*i/max)*(graph_w), graph_y + graph_h/2 - pi*500 + def*500
                gsel prev
            }
        }
    next
    gsel 0
    color 0, 0, 0
    boxf 0, 0, 100, 10
    //非表示モードではまとめて円を描画
    if(visible == 0) :     pos vx, vy:  gcopy 2, 0, 0, r*2, r*2
    stop
*init_graph
    if(dimensions == 0):def = 3.1416:else:def = 4.1888
    prev = ginfo(3)
    gsel 0
    color 255, 255, 255
    boxf 0, graph_y, graph_x + graph_w, graph_y + graph_h
    color 0, 0, 0
    pos 30, graph_y + graph_h/2 : mes ""+def
    line 99, graph_y, 99, graph_y + graph_h
    line 99, graph_y + graph_h/2, 700, graph_y + graph_h/2

    pos 30, graph_y
    mes "" + (def + (1.0 * graph_h / 500))
    line 99, graph_y, graph_x + graph_w, graph_y
    gsel prev
    return ;

ベジエ曲線が四角形の中にあるのか、無いのかを知るには、どうすればよいのだろうか。

1. ベジエ曲線は、制御点の凸包の外側には飛び出さないということなので、四角形の中に凸包の領域が無ければ曲線は含まれず、有れば曲線が含まれるかもしれない。
2. 最初か最後の制御点が四角形に含まれていれば、四角形内に曲線が含まれる。
3. 曲線と四角形の辺に交点があれば、四角形内に曲線が含まれる。

１の判別：
　凸包を求める必要がある。
　グラハムスキャンという方法があるらしい。

2の判別：
　if 文を使う

3の判別：
　辺の直線の式にベジエ曲線の式を代入し、t についての方程式を解く。
　が、次数が高くなるかもしれないことになっているので、工夫がいる。

　ニュートン法という方法で 解が求められるというので、実装してみたのだが、どんなケースでも解が求まるわけではないらしい。 難しそうなのでスルーしていた bezier clipping に手を出してみる。

こちらのページとにらめっこして理解を試みる。
http://nis-lab.is.s.u-tokyo.ac.jp/nis/bezclip.html

Bezier Clipping は、最近のアルゴリズムらしい。
曲線の制御点による凸包の  y= 0 を通る辺 の x軸との交点２つをそれぞれmin, max として、曲線の当区間を切り出し、それを再帰的に行うというもの。

まず、 曲線の式を用意する

次に、直線の式を用意する

そして、解くべき式は次のようになる。

この式の各項 の f i が制御点 i の値になる。 各項での t は、 1 を 項数で均等に割ったものらしい。

この式に対して、次のように処理する。参考にしたのはこちらのページ
http://www-ui.is.s.u-tokyo.ac.jp/~kenshi/assignment/cg07/bezierclip.html

1. 凸包を得る
2. 各辺のうち、f(t)=0 を通る辺の、f(t) = 0 での t を求め、min, max とする
3. そのような辺がなければこの部分には解がない
4. min, max の範囲が、もともとの半分以上なら、解が複数あるとみて、曲線を半分に分割する。
   その両方の曲線を再帰し、両方の結果を戻す
5. そうでない場合、収束しているか調べ、していれば値を返す
6. していなければmin, max の範囲で曲線を切り出し、再帰し、結果を戻す。

自分の場合、切り取った曲線の t の変域は 0 <= t <= 1 だったので、 再帰した処理が戻ってくるとき、(max – min)倍に縮めて最小値を加算する処理を付け加えた。

//このオブジェクトをベジエクリッピングする
public float[] bezierclipping() {
    //Pointf にする t は項数で均等に分けたもの
    PointF[] pts = new PointF[points.Length];
    PointF[] convex;
    for (int i = 0; i < points.Length; i++) {
        pts[i].X = (1.0f / (points.Length-1)) * i;
        pts[i].Y = points[i];
    }

    //凸包を得る
    convex = mymath.convex_hull(pts);

    //凸包のうち、ゼロを通る辺を最小・最大値とする
   float min = float.MinValue, max = float.MaxValue;
   for (int i = 0; i < convex.Length; i++)
   {
       int j = (i + 1) % convex.Length;
       if (convex[i].Y * convex[j].Y < 0)
       {
           if (min == float.MinValue)
           {
               min = (convex[i].X * (float)Math.Abs(convex[j].Y)
                   + convex[j].X * (float)Math.Abs(convex[i].Y))
                   / ((float)Math.Abs(convex[i].Y) + (float)Math.Abs(convex[j].Y));

           }
           else if (max == float.MaxValue)
           {
               max = (convex[i].X * (float)Math.Abs(convex[j].Y)
                   + convex[j].X * (float)Math.Abs(convex[i].Y))
                   / ((float)Math.Abs(convex[i].Y) + (float)Math.Abs(convex[j].Y));
               break;
           }
       }
       else if (convex[i].Y == 0) {
           if (min == float.MinValue) min = convex[i].X;
           else if (max == float.MaxValue) { max = convex[i].X; break; }
       }
   }
    if (min == float.MinValue || max == float.MaxValue)
        return new float[] { };

    if(min > max){
        float tmp;
        tmp = min; min = max; max = tmp;
    }
    if (max - min < (float)Math.Abs(pts[pts.Length - 1].X - pts[0].X) * 0.5f)
    {
        //0.5以下に小さくなるのならば、解はひとつ
        if (max - min < 0.0001) //収束判別
        {
            return new float[] { (max - min) / 2 };
        }

        //もう一度縮めたい
        float[] tmp = cut(min, max).bezierclipping();
        for (int cnt = 0; cnt < tmp.Length; cnt++) {
            tmp[cnt] *= (max - min)/1;
            tmp[cnt] += min;
        }
        return tmp;
    }
    else 
    { 
        //0.5以上ならば解が複数ある
        bezierfunc[] cut_half = cut(0.5f);
        float[] a, b, c;
        a = cut_half[0].bezierclipping();
        b = cut_half[1].bezierclipping();

        c = new float[a.Length + b.Length];
        for (int i = 0; i < c.Length; i++) {
            c[i] = (i < a.Length) ? a[i]/2 : (b[i - a.Length]/2 + 0.5f);
        }
        return c;
    }
}

このシリーズの他記事 http://prog-city.spaces.live.com/lists/cns!E5379F11D9E8BB3F!1047/

計算が複雑になってきてきた。そこで、多項式の計算をクラスにしてみた。
四則演算と微分・積分を実装した。
ただし、 1/x の積分は出来ない。 とか、速度の問題とかいろいろとある。

カーブ半径を取得したい。
曲線の曲がり具合は、曲率というもので表されるらしく、
一般に用いる半径は曲率半径といって曲率の逆数らしい。

したがって、曲率を求めたい。が、
http://ja.wikipedia.org/wiki/%E6%9B%B2%E7%8E%87
わからない。

検索して調べていると

(x’y’’ – x’’y’ ) / √{({x’}^2 + {y’}^2)^3}

x(t), y(t) : 曲線上の点
x’(t), y’(t) : 接線の傾き

こんな式を見つけたのでそのまま計算して逆数にすると、正しいらしい結果が出てきた。

カーブ半径に合わせて円を描いてみた。

下がコードの一部。 tmp1 が上式の 分子、 tmp2 が分母。

       public float curvature( float t) {
            float tmp1 = -f_tangental_x.differentiate().GetValue(t) * f_tangental_y.GetValue(t) +
                            f_tangental_y.differentiate().GetValue(t) * f_tangental_x.GetValue(t);
            float tmp2 = (float)Math.Sqrt((double)powFI((powFI(f_tangental_x.GetValue(t), 2) 
                            + powFI(f_tangental_y.GetValue(t), 2)), 3));
            if (tmp2 == 0) return 0;
            return  (tmp1 / tmp2);
        }

differentiate は多項式クラスで微分をするメソッドで、 (x^a)’ = ax^(a-1) をすべての項に対して行う。
GetValue は多項式の変数に代入し、値を得るもの。

        //微分する
        public polynomial differentiate() {
            polynomial result = new polynomial();

            foreach (KeyValuePair<float, float> kpy in this.values)
            {
                if (kpy.Key != 0)
                    result.values[kpy.Key - 1.0f] = kpy.Key * kpy.Value;
            }

            return result;
        }

多項式クラスでは、各項を、キーを指数、値を係数とするリストに保存している。

正確さや速度に問題が出てくるかもしれないが、機能としてはこれくらいかと思う。
次はカーブをつなげたレールの設計とか。

このシリーズの他記事
http://prog-city.spaces.live.com/lists/cns!E5379F11D9E8BB3F!1047/

ビワを食べた。
家族には不人気のビワ、好きなのは自分だけらしい。
おいしかったからこれも植えてしまおうか。

トカゲを見た。 トカゲと呼んでいるがこれはカナヘビのような気がする。
http://ja.wikipedia.org/wiki/%E3%82%AB%E3%83%8A%E3%83%98%E3%83%93
 
 
飼いたいが餌は虫。しかも手ごろなデュビアというムシがなんと ゴキの仲間。
無理だ。

C# の機能を一通り勉強。 しかし ・・・ 情報が古い。
とりあえず気にせずに コーディング。

C#エッセンシャルズ 第2版 ベン アルバーリ ブラッド メリル ピーター ドレイトン