Jul 26, 2008

Maximaでドラゴンカーブ

MathematicaとかMapleとかの数式処理ソフトは、何であんなに高いのだろう。2008/7/26現在、個人向けライセンス価格はMathematicaは40万円超え、Mapleも30万円超えだ。教育機関にも政府機関にも属していない、日曜数学ファンの私にはとても手が出せない。
会社でも研究職とかでなく、ソフトウェア開発の付帯業務をやってるので、数学は一切使っておらず、買ってもらえる見込みは無い。バグ収束曲線が云々で必要とか言っても、そんな統計処理ならExcelで十分と言われてしまうであろう。
従って、MathematicaやMapleに触れる機会が無い。

そんな折、Maximaというフリーの強力な数式処理ソフトがあることを知ったので、今回ヤボ用があって、インストールして使ってみた。
とりあえず、テーラー展開させてみる。

 KURO-BOX> maxima

Maxima 5.10.0 http://maxima.sourceforge.net
Using Lisp GNU Common Lisp (GCL) GCL 2.6.7 (aka GCL)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
This is a development version of Maxima. The function bug_report()
provides bug reporting information.
(%i1) taylor(sin(x),x,0,10);
                          3    5      7       9
                         x    x      x       x
(%o1)/T/             x - -- + --- - ---- + ------ + . . .
                         6    120   5040   362880
(%i2) taylor(sqrt(r^2-y^2),y,0,10);
                    2      4      6         8       10
                   y      y      y       5 y     7 y
(%o2)/T/       r - --- - ---- - ----- - ------ - ------ + . . .
                   2 r      3       5        7        9
                         8 r    16 r    128 r    256 r
(%i3)
まあ言いたいことは十分わかるが、表示がなんか残念だ。
そこで、満足のいく表示を得るためにTeXmacsを使ってみる。以下の画面はMaxima 5.11.0とTeXmacs 1.0.6.11を使って実際に表示されたものである。%i??という青い部分が入力した部分だ。(表示環境はFreeBSD 6.3+Xorg 1.4.0)

積分もしてみる。

e^(-x^2)を-∞から+∞まで積分すると√πになるのは、手計算でやると面積分にして極座標変換してと技を尽くして偶然解ける感じの有名な例題だが、Maximaは区分求積による近似値計算ではなく、あっさり解析的に計算したようだ。

せっかくTeXmacsを使ってるので、見た目が派手な連分数を表示してみる。

cf(continued fraction)の出力は連分数のリスト表記になっていて、それをcfdisrepすると連分数表記になる。それをratで実数にすると(無限の連分数でないので元々実数だが)約分され、floatで小数にすると元の√2に限りなく近いことがわかる。
蛇足だが、cfdisrepしないと実数や小数に変換できない。

√3についてもやってみる。

方程式を解いてみる。solveで方程式を解析的に解ける。

解けない場合は、xの解にxが混ざったりする。(例は省略)

allrootsを使うと、多項式の方程式を数値的に計算できる。

指数の方程式を渡すと文句を言われる。

find_rootを使うと、(おそらくニュートン法による)近似計算で、解が1つだけ含まれる範囲の解を計算できる。

x^3=3^xのe~100の間の解は当然3なのだが、誤差が出た。


さて、せっかくフリーソフトなので、参考資料も無料のものを、と思って、図書館でMaximaの本を探したら、無かった。その代わり、あるMathematicaの本で、わずか十数行のコードでドラゴンカーブ(ドラゴン曲線)を描くMathematicaのコードを見つけた。マンデルブロー集合とかペアノ曲線とかのフラクタル図形の簡単なやつだ。Maximaの練習とMathematicaとの比較を兼ねて、その本を参考にMaximaで似たようなコードを書いてみた。

dragon(segments_):= block([p, q, r, s],
 p: map(
  lambda([x], [
   [x[1][1], x[1][1] + (x[1][2] - x[1][1]) * (1-%i)/2, x[1][2]],
   [x[2][1], x[2][1] + (x[2][2] - x[2][1]) * (1+%i)/2, x[2][2]]
  ]), segments_),
 q: map(
  lambda([x], [
   {[x[1][1], x[1][2]], [x[1][2], x[1][3]]},
   {[x[2][1], x[2][2]], [x[2][2], x[2][3]]}
  ]), p),
 r: map(lambda([x], listify(x)), flatten(q)),
 s: map(lambda([x], [realpart(x), imagpart(x)]), flatten(r)),
 plot2d([discrete, float(s)],
  [gnuplot_preamble, "set xrange [-1.2:1.2];set yrange [-1.2:1.2]"]),
 r
)$

p: [[[0, 1/2-%i/2], [1/2-%i/2, 1]]]$
for i:1 thru 10 do p:dragon(p)$
gnuplotが使える環境でこれを実行すると、以下のようなドラゴンカーブが描かれる。TeXmacsは必要ない。順に、L字型に配置した2本の線から始めて、全体を90°回転させて先端に繋げる、というのを1回、2回、…10回繰り返したものである。
アルゴリズムとしては、全ての線をL字に折り曲げる、というのを繰り返せばドラゴンカーブと同じ形になるというのを使っている。必要最小限にコードを簡略化して難解になったし、同様のアルゴリズムの説明は多くあると思うので、ここではコードの説明は省く。

See more ...

Posted at 22:09 in 数学 | WriteBacks (12)
WriteBacks

Barton のflatten.lisp を改良してMathematicaのように使えるようにしてみました。 こんな感じ (%i4) load("newflatten.lisp"); (%o4) newflatten.lisp (%i5) nflatten([[[a,b,c],d],[e,f]],1); (%o5) [[a, b, c], d, e, f] (%i6) nflatten([[[a,b,c],d],[e,f]]); /*default level 3*/ (%o6) [a, b, c, d, e, f] (%i7) nflatten([[[a,b,c],d],[e,f]],2); (%o7) [a, b, c, d, e, f] 関数定義はいかのとおりです。 (defun $nflatten (e &optional (n 5)) (setq e (ratdisrep e)) (cond ((or ($atom e) ($subvarp e)(or (member ($inpart e 0) (list '&^ '&=)))) e) (t (let ((op (multiple-value-list (get-op-and-arg e)))) (setq e (cadr op)) (setq op (car op)) (setq e (mapcar #'(lambda (x) (flatten-op x op n)) e)) (setq e (reduce #'append e)) (cond ((and (consp (car op)) (eq (caar op) 'mqapply)) (append op e)) (t `(,op ,@e))))))) (defun flatten-op (e op nlev) (let ((e-op) (e-arg)) (setq e-op (multiple-value-list (get-op-and-arg e))) (setq e-arg (cadr e-op)) (setq e-op (car e-op)) (cond ((and (>= nlev 1)(equal e-op op)) (mapcan #'(lambda (x) (flatten-op x op (1- nlev))) e-arg)) (t (list e)))))

Posted by GO_MAXIMA at 01/09/2009 09:47:22 PM

ありがとうございます。 最初のコードは、安全でない{}(set)とmap1つが無くなって、少しシンプルになりました。 dragon(segments_):= block([p, q, r], p: map( lambda([x], [ [x[1][1], x[1][1] + (x[1][2] - x[1][1]) * (1-%i)/2, x[1][2]], [x[2][1], x[2][1] + (x[2][2] - x[2][1]) * (1+%i)/2, x[2][2]] ]), segments_), q: nflatten(map( lambda([x], [ [[x[1][1], x[1][2]], [x[1][2], x[1][3]]], [[x[2][1], x[2][2]], [x[2][2], x[2][3]]] ]), p), 1), r: map(lambda([x], [realpart(x), imagpart(x)]), nflatten(q,2)), plot2d([discrete, float(r)], [gnuplot_preamble, "set xrange [-1.2:1.2];set yrange [-1.2:1.2]"]), q )$ p: [[[0, 1/2-%i/2], [1/2-%i/2, 1]]]$ for i:1 thru 5 do p:dragon(p)$ なお、ご提示頂いたnflatten()を動かすには、flatten.lispにあるget-op-and-arg()を明示的に定義する必要があるようです。 当方の環境では、load(newflatten);の前にload(flatten);すると動きました。

Posted by ynomura at 01/10/2009 02:06:28 PM

連分数で漂着致しました。 http://f.hatena.ne.jp/mathnb/20100528021239 (内部リンク:20100528021239.gif ) この 廣大の函数fが    ◆  いい加減法 (と命名します);    x^2=7 3倍し;3x^2=3*7   8*xを(いい加減)加え 3x^2+8*x=3*7+8*x x*(3*x+8)=8*x+21 から 生まれた。なんて 信じる 学習者は 世界に 存在しない。 授業で いい加減法で 導出される方 は 存在しそう(嗚呼)......◆ ★★ 廣大の函数f の導出過程を ご教示ください★★ (f の 導出にこそ 意味が在ると 考えます ので)  ---------------------------------------------------------------------            また  http://f.hatena.ne.jp/mathnb/20100528021239     に倣い 例えば Sqrt[3], Sqrt[109], Sqrt[263], Sqrt[431], Sqrt[601], Sqrt[773], Sqrt[971], Sqrt[1153]      等のそれぞれについて 廣大の函数f に相当する函数の導出を、 遊び心で、お願い致します; f(Sqrt[3])=Sqrt[3](不動点) f[x]= f(Sqrt[109])=Sqrt[109](不動点) f[x]=

Posted by gb at 06/05/2010 11:14:56 PM

連分数とは関係ないですね。 問題を解いてみればわかりますが、Nを問題で言う所の7として、f(Sqrt[N])=Sqrt[N]でf(x) > xでf'(x) > 0であれば何でもいいのです。(収束が速いか遅いかだけの違い、件の広大の入試問題からはそれ以上の条件は読み取れませんでした) f(x)の形を(bx+c)/(ax+b)に限定し、a,b,cが全て整数、b^2 - a^2 * N = 1という制約を付けると、それは整数問題なので、総当たりで計算するしか無いと思います。 ちなみに、 N=3: a=1, b=2, c=3 (→ f(1)=1.667, f(f(1))=1.732) N=263: a=8579, b=139128, c=2256277 N=109, 431, 601等については10^6までに適当な整数の組み合わせはありませんでした。 N=107ならa=93, b=962, c=9951 がありました。 b^2 - Na^2 = 1という制約を外すと、条件はc=aNとb^2>N*a^2だけになりますので、例えば広大の問題の8を9以上に変えても良いことになります。 N=109: a=3, b=32, c=109*3 N=431: a=3, b=63, c=431*3 N=601: a=3, b=74, c=601*3 (以下略) いかがでしょうか。

Posted by ynomura at 06/06/2010 10:48:45 PM

同じ文面で たずねた ところ 或る方から 連分数     [2,1,1,1,2+x]=2+1/(1+1/(1+(1/(1+1/(2+x)))))=(8*x + 21)/(3*x + 8) が広島大学で出題されたf(x)の正体。 と 返事をいただきました。 連分数数と大いに関係が在ります。 Sqrt[61]等で f を 導出願います。

Posted by gb at 06/07/2010 11:09:00 AM

面白い情報ありがとうございます。 そういう連分数で書ける式に限定する条件、またはそれに等価な条件は、件の広大の問題には含まれていないと思いますが、まあ形からして元ネタはそれなんでしょうね。 Maximaでやってみました。 f(N) := block([y], y:cf(sqrt(N)), y[length(y)] : x+floor(sqrt(N)), ratsimp(cfdisrep(y))); f(7); →(8*x+21)/(3*x+8) f(61); →(29718*x+232105)/(3805*x+29718) f(109); →(8890182*x+92816225)/(851525*x+8890182)

Posted by ynomura at 06/07/2010 09:08:04 PM

再読し 気付き  追記します;http://www.wolframalpha.com/input/?i=e%5E%28-x%5E2%29

Posted by gb at 06/08/2010 01:37:19 AM

>Maximaでやってみました。多謝!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!mathematica に ヤラセるには如何に記述すればよいのでしょうか?ご教示下さい。ContinuedFraction[Sqrt[N] 云々....

Posted by gb at 06/08/2010 01:51:54 AM

再読し 追記; http://www.wolframalpha.com/input/?i=z%5E17%3D1 0の原像 f^(-1)(0) のみでなく fによる像 f[D] etc mathematicaで叶いmath.

Posted by gb at 06/08/2010 01:59:34 AM

Mathematicaだと、こんな感じでヤラセられると思います。F[n_] := Module[{y = Flatten[ContinuedFraction[Sqrt[n]]]}, y[[Length[y]]] = x + Floor[Sqrt[n]]; Simplify[FromContinuedFraction[y]] ]F[7]F[61]いつの間にか、MathematicaのWebインターフェイスができてたんですね。見よう見まねで、一応そのWolframAlphaでもやってみました。http://www.wolframalpha.com/input/?i=Simplify%5BFromContinuedFraction%5BReplacePart%5BFlatten%5BContinuedFraction%5BSqrt%5B%23%5D%5D%5D%2C+-1+-%3E+x+%2B+Floor%5BSqrt%5B%23%5D%5D%5D%5D%26+%2F%40+%7B7%2C41%2C61%7D%5D残念ながら1行1式でしか書けないようですので、逐次処理のプログラミングに慣れた者には結構厳しいですね。

Posted by ynomura at 06/08/2010 07:43:42 PM

In[62]:= Table[Prime[k], {k, 4, 19}] Out[62]= {7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67} In[64]:= {F[Prime[4]], F[Prime[13]], F[Prime[18]]} Out[64]= {(21 + 8*x)/(8 + 3*x), (205 + 32*x)/(32 + 5*x), (232105 + 29718*x)/(29718 + 3805*x)} WolframAlphaは更に或点で展開表示していますが「なんでだろー」

Posted by gb at 06/09/2010 01:31:55 AM

b^2 - a^2 * N = 1を満たす整数a,bを求める問題は、「ペル方程式」と呼ばれる、色々な所で出てくる方程式だそうです。上記の広大の問題は、やはり連分数と直接の関係は無く、「ペル方程式」を応用して作られた問題だと考えられます。 また、言うまでもなく、この問題はニュートン法の問題として作られたものであり、導出過程と出題意図とは関係ありません。 f(x)の形を(bx+c)/(ax+b)に限定し(平方根を分数式で近似することを考える、分子を(cx+d)としても結局b=cになる)、a,b,cが全て整数、f(√N)=√Nとすると、このfを繰り返し適用する時に収束が速くなるには、x=√Nの近傍でf'(x)=(b^2-aN)/(ax+b)^2の絶対値がなるべく小さい方が良い(大きいほどf(x)が√Nから遠ざかる)ので、b^2 - a^2 * N = ±1という制約がつきます。従って、広大の問題のa,bは「ペル方程式」の最小解のことだと思います。(連分数でも求まるが、連分数を使わなくても求まる) f(x)を繰り返し適用すると、分母にペル方程式の解が次々に出てきますが、これもペル方程式の特性を利用したものと考えられます。 (参考)271828さんの滑り台Logへのリンク ・ペル方程式と共役の考え無理数と連分数(追記あり) 残念なことに、gbさんのコメントは、同一文面を多方面に送信していたマルチポストだったようです。何人かの方が個別に同じ問題に取り組んでおられました。無駄な労力を使わされたようで、あまり気分の良いものではありません。

Posted by ynomura at 08/10/2010 05:27:06 PM

Jul 20, 2008

整数演算のみで楕円を描く

前のエントリーでブレゼンハムのアルゴリズムを使って乗除を使わずに整数演算だけで円を描いてみたが、同じ原理を応用して、乗除を使わずに整数演算で楕円を描けないだろうか、と考えてみた。
Javaのデモページへのリンク
ソースコード

横幅が2a、縦幅が2bの楕円を考え、円の時と同様に、(a,0)や(0,b)から始めて、楕円の外側か内側かを判定してポインターを1ずつ縦横斜めに動かしながら、点を描いていく。第2象限~第4象限(90°~360°)はx,yの符号を変えながら第1象限(0°~90°)と同じ曲線を描けば良いので、第1象限だけを考える。

楕円の式は(x/a)2+(y/b)2=1であるが、整数演算のために割り算を避けて、b2x2+a2y2=a2b2とする。
まず、(x,y)=(a,0)からyを+1していって、xは(x,y)と(x-1,y)の中点が実際の楕円の外側になったら-1する(前のエントリー参照)。それを、xの動きがyの動きより大きくならない所、すなわち接線の傾きが135°になる所まで繰り返す。
次に、(x,y)=(0,b)からxを+1していって、同様に(x,y)と(x,y-1)の中点が楕円の外側になったら-1する、というのを、接線の傾きが135°になる所まで繰り返す。

式で書くと、(x,y)と(x-1,y)の中点は(x-1/2,y)なので、楕円の式に当てはめて、b2(x-1/2)2+a2y2 > a2b2なら中点が外側ということになる。E(x,y)を左辺-右辺とすると、E(a,0)=b2(-a+1/4)であり、E(x,y)はyが+1するとa2(2y+1)増え、xが-1するとb2(2x-2)減る。このE(x,y)の1回の増減分をそれぞれdEdy, dEdxと置くと、dEdyはyが+1すると2a2増え、dEdxはxが-1すると2b2減る。a2やb2は先に計算しておけばいいので、これで乗算を使わずに点を動かしていけることがわかる。
(x,y)と(x,y-1)の中点については、今の話のxとyを置き換えれば良い。

ということで書いたのが、サンプルコードのdraw1()である。コード上では上記のa,bはWx,Wyとしている。whileループの中では乗算は使っていない。
デモでは、赤線で描いているのがそれである。

さて、デモでは、比較のために青線でjava.awt.Graphics#drawArc()で楕円を描いて、その上にdraw1()で同じ大きさの楕円を赤線で描いているが、Java VMにも依るだろうが、青線と赤線が結構ずれている。原因として1つ間違いないのは、X軸、Y軸をy=0, x=0に取って上下対象、左右対称にしているので、楕円の横幅や縦幅が偶数の時も奇数しか書いていないことだ。-1<=y<=1で書くと縦幅は3、-2<=y<=2で書くと縦幅は5である。実は前のエントリーのブレゼンハムの円描画のアルゴリズムでも偶数幅の円が描けないという同じ問題があるのだが、コードの簡潔さを重視して敢えて触れなかった(巷のサンプルコードでも大概無視してるし)。
今回はI/F仕様をjava.awt.Graphics#drawArc()に合わせて楕円のバウンディングボックス(外接矩形)のサイズを入力するようにしているので、それに見合った動作をさせたいものである。

横幅が偶数の場合、-a < x < a-1の範囲で描くとすると、横方向の対称軸はx=-0.5となる。x=0とx=-1とでY座標が同じになるのだが、対称軸がx=0の楕円のx=0.5, x=-0.5の時のyをX方向に-0.5ずらして描くと、yの最大点、最小点が無くなってしまう。実際、a,bをそのままにして(x,y)=(a-1/2,0)や(0,b-1/2)からプロットすると、楕円の上下左右の端が十分に膨らまず、少し四角に近くなってしまった。
そもそも、バウンディングボックスを中心に考えると、x=-a, x=a-1の点はそれぞれ必ずxの最小点と最大点なのであり、中心軸がx=-0.5なので、X軸方向の半径がa-1/2の楕円を描かないといけないのであり、半径がaの円を平行移動する方針ではきれいな楕円にならないようだ。

そこで、曲線描画中は、楕円の横幅や縦幅が偶数の時は径を1/2減らして計算する方向で考えてみた。2aや2bが偶数の場合、楕円の式は(b-1/2)2(x-1/2)2+(a-1/2)2y2 = (a-1/2)2(b-1/2)2となり、(x,y)=(a-1,0)から始めてxを-1、yを+1していく時、(x-1,y)と(x,y)の中点が楕円の外側かどうかを判定する式は(b-1/2)2(x-1)2+(a-1/2)2y2-(a-1/2)2(b-1/2)2 > 0 となる。以下、左辺を4倍して考える。左辺の初期値はy=0なので(2b-1)2(-x+1/4)であり、yが+1すると(2a-1)2(2y+1)増え、xが-1すると(2b-1)2(2x-1)減る。以下略。結局、辛うじて乗算を使わずに点を動かしていける。

ということで書いたのが、サンプルコードのdraw2()である。コード上では上記の2a,2bがWx,Wyに対応する。
デモでは、比較のために青線でjava.awt.Graphics#drawArc()で楕円を描いて、その上にdraw2()で同じ大きさの楕円を緑色の線で描いている。JRE 6.0だと、緑線と青線はほとんど差が無い。

See more ...

Posted at 23:16 in PC一般 | WriteBacks (2)
WriteBacks

初めまして、自作OSを作っている者です。 現在、GUI周りをちょこちょこと実装していまして、その中で楕円の描画にこちらのアルゴリズム(draw2)を使わせてもらいたいのですが、許可を頂けないでしょうか? ちなみに、OS本体は近いうちにGPLv2で公開しようと考えております。 どうか宜しくお願い致します。[link]

Posted by Liva at 12/28/2010 07:37:45 PM

全然大丈夫です。ご自由にお使いください。

Posted by ynomura at 12/29/2010 06:17:36 PM

Jul 03, 2008

乗除を使わずに円を描く

三角関数や平方根はもちろん、×や÷をも使わずに円を描くアルゴリズムがあることを知った。通称ブレゼンハムのアルゴリズム(Bresenham's circle algorithm)と呼ばれるらしい。

さぞかし恐ろしく複雑な反復演算をするんだろう、乗除を使った方が速いというオチではないか、と思いながらウェブ上で調べてみると、どのサンプルコードを見てもとてもシンプルだった。ここ「伝説のお茶の間」)とかここ(Wikipedia(英語))とかに図解入りの説明とサンプルコードがある。これは速そうだ。

それらのページの図を見て基本的な原理はすぐに理解できたが、説明とソースコードが理解できなかったので、自分で考えてみた。
Javaのデモページへのリンク
ソースコード

基本的な原理は、例えば格子上のXY平面(X座標とY座標は整数のみ)で(x,y)=(r,0)の点から45°分の円弧を描く場合、同じY座標でX軸方向に複数の点を描く必要は無いので、yを1ずつ動かしながら、適当な時にxを1つ動かして、点を描いていけば良い。45°分の円弧が描ければ、そのX軸対称、Y軸対称や90°回転を使って360°分の円弧が作れる。

yを+1していく時、いつxを-1すれば良いだろうか。
円の式はx2+y2=r2であるが、x,y,rを整数に限ると、もちろんこの等式は成立せず、なるべく誤差の小さいxを選んでいくことになる。ここでは左辺と右辺の差e(x,y)=x2+y2-r2を誤差として考える。
仮にyを固定してxを-1すると、誤差は(x-1)2+y2-r2になるので、誤差の差分は2x-1である。つまり、誤差eが2x-1より大きい時は、xを-1してもまだ(x,y)は円の外側なので、xを-1する必要がある。
yが1増えるとeが2y+1増えるので、eに2y+1を足していって2x-1を超えればxを-1すれば良い、と考えられる。そう考えてとりあえず書いたのがサンプルコードのmethod1()である。

なるほど、こういう手があったか。
要領が解ったので、もう少し真剣に考えてみる。

上の方法では点が外側に行き過ぎないようにしてるだけなので、円の内側に点を取った方が円弧に近い場合も、円弧より遠い外側の点を取ってしまう。
そこで、次に(x,y)と(x-1,y)の中点が円弧の内側にあるか外側にあるかで点を選ぶことが考えられる。そのためには、円弧とその中点との誤差e(x,y)=(x-1/2)2-y2-r2が正か負かを判定すれば良い。
このeはyが+1すれば2y+1増え、xが-1すれば2x-2減る。(r,0)から始めると初期値が-r+1/4なので、eは整数にならない。従って、e=0になることはなく、e>0かe<0である。yを+1していって、e>0の時は中点が円の外側にあるから、内側の点(x-1,y)を取れば(xを-1すれば)良く、逆にe<0の時は外側の点(x,y)を取れば(xをそのままにすれば)良い。
ということで書いたのが、サンプルコードのmethod2()である。
なお、このコードは上記のWikipediaにあるコード(通称Bresenham's circle algorithm)と意味的にはほぼ同じであり、全く同じ円を描くことを確認した。

MSXを使っていた頃、マシン語で書いてもSIN()やCOS()や乗除が遅いのに対して、BASICのCIRCLE文の円てどうやってこんなに速く描いてるんだろう、と不思議に思っていたものだが、内部でこういうのが動いてたんだなと思うと、感慨深い。

See more ...

Posted at 22:08 in PC一般 | WriteBacks (0)
WriteBacks