Dec 29, 2008

統計学復習メモ6: χ^2検定で独立性の検定

以前のメモではカイ2乗検定の例としてばらつきの検定(適合度の検定)をやってみたが、もう1つのカイ2乗検定の使われ方として、独立性の検定がある。これもやってみる。

プログラマー35歳定年説というのがある。あるプロジェクトにおいて、
35歳未満で、受けたバグ指摘件数が20件未満のプログラマーが45人、
35歳未満で、受けたバグ指摘件数が20件以上のプログラマーが15人、
35歳以上で、受けたバグ指摘件数が20件未満のプログラマーが25人、
35歳以上で、受けたバグ指摘件数が20件以上のプログラマーが15人、
だったとする。表にすると、

O(a,b)B<20B≧20
age<354515
age≧352515
(Bはバグ票数)である。暗雲が立ちこめてきたが、有意水準を95%として、35歳以上であることとバグ票数20件以上であることに相関はあるだろうか。

まず、age≧35とB≧20が独立の場合の理論値を計算する。age<35が60人、age≧35が40人、B<20が70人、B≧20が30人なので、これだけから考えると、B<20の70人は60:40でage<35とage≧35に分かれるので、それぞれ42人、28人である。同様に、B≧20の30人は18人、12人に分かれる。

E(a,b)B<20B≧20
age<354218
age≧352812
帰無仮説を、相関がない(独立である)と取ると、観測値が理論値から大きくばらついていれば棄却されるので、ピアソンのカイ2乗検定が使える。χ^2分布に従う、ピアソンの検定統計量
T=¥sum_{i=1}^n¥frac{(O_i-E_i)^2}{E_i}
を計算すると、T=(45-42)^2/42 + (25-28)^2/28 + (15-18)^2/18 + (15-12)^2/12≒1.79となる。このTは、4つの理論値を固定すると、4つの観測値のどれか1つが決まると定まるので、自由度は1である。自由度が1のχ^2分布の左側が95%になるχ^2の値は3.84なので、今回のTはばらついていない方(χ^2が小さい=理論値の通りに近い方)の95%に入っており、仮説は棄却されず、相関があるとは言えないことになる。

これは、B≧20かどうかで分けたのが明暗を分けた可能性があり、もしB≧25で分けると

O(a,b)B<25B≧25
age<35555
age≧353010
だったりした場合はT=5.23でOuchとなるので、あまりいい例ではないかも知れないが、そこまでは敢えて考えないことにする。


同じ考え方で、2つの属性で2x2に分類したデータから2つの属性の独立性を検定する方法を定式化してみる。1か2かを取る2つの属性をa,bとし、母数がnの、それぞれの属性の組み合わせに属する標本数の観測値を

O(a,b)b=1b=2
a=1x11x12
a=2x21x22
(x11+x12+x21+x22=n)とする。a,bを独立とした場合の期待値は、b=1について(a=1):(a=2)=(x11+x12):(x21+x22)、a=1について(b=1):(b=2)=(x11+x21):(x12+x22)…となるので、
nE(a,b)b=1b=2
a=1(x11+x12)(x11+x21)(x11+x12)(x12+x22)
a=2(x11+x21)(x21+x22)(x12+x22)(x21+x22)
の定数倍である。これらを全て足すと(x11+x12+x21+x22)^2=n^2なので、これらが期待値のn倍であり、1/nすれば期待値になることがわかる。
E(a,b)b=1b=2
a=1(x11+x12)(x11+x21)/n(x11+x12)(x12+x22)/n
a=2(x11+x21)(x21+x22)/n(x12+x22)(x21+x22)/n
これを、
O(a,b)b=1b=2sum
a=1x11x12a1
a=2x21x22a2
sumb1b2n
という風に行和、列和を取って整理すると、
E(a,b)b=1b=2
a=1a1b1/na1b2/n
a=2a2b1/na2b2/n
となる。検定統計量Tは、各セルの(O-E)^2/Eを合計したものなので、
T=¥sum_{i=1}^{2}¥sum_{j=1}^{2}¥frac{(x_{ij}-¥frac{a_ib_j}{n})^2}{¥frac{a_ib_j}{n}}
と書ける。これを自由度1のχ^2分布のχ^2の値として見れば良いわけである。

aがp個の値、bがq個の値を取る時も同じことが言えるので、
T=¥sum_{i=1}^{p}¥sum_{j=1}^{q}¥frac{(x_{ij}-¥frac{a_ib_j}{n})^2}{¥frac{a_ib_j}{n}}
(aiはa=iの標本の和、bjはb=jの標本の和)
となる。自由度は、aについては(p-1)、bについては(q-1)なので、全体としては(p-1)(q-1)である。

以上を踏まえて、もう1つ例題をやってみる。原理はわかったので、手計算チックな計算は卒業して、次はMaximaに組み込みのカイ2乗検定関数を使おう、と思って関数を探すと、手元のMaximaには入ってなかった。ExcelにはCHI2TEST()という名前でデフォルトで入ってるのだが、Maximaにはデフォルトでは入っていないらしい。上記の計算だけなら大して複雑じゃないので、自分で関数を作ってみる。

chi2test(x_,row_,col_):=block([i,j,n,a,b,chi2],
    a:map(lambda([i],sum(x_[i][j],j,1,col_)),create_list(i,i,1,row_)),
    b:map(lambda([j],sum(x_[i][j],i,1,row_)),create_list(j,j,1,col_)),
    n:sum(a[i],i,1,row_),
    chi2:sum(sum((x_[i][j]-a[i]*b[j]/n)^2/(a[i]*b[j]/n),j,1,col_),i,1,row_),
    print("chi^2=",chi2),
    print("95% confidence at",quantile_chi2(0.95,(row_-1)*(col_-1))),
    print("90% confidence at",quantile_chi2(0.90,(row_-1)*(col_-1))),
    chi2
)$
試しに
chi2test([[45,25],[15,15]],2,2),numer;
(numerは結果が分数になるのを避けるため)とやってみると、
chi^2= 1.785714285714286
95% confidence at 3.841458820694166
90% confidence at 2.705543454095465
と出力される。有意水準90%でも、2つの属性が独立である仮説は棄却されないことになる。

次の例は、バグの原因を設計ミス要因かコーディングミス要因かに分けて、さらにそれをプログラマーの年齢層で分けて数えてみたという架空のデータである。

O(a,b)設計ミスコーディングミス
30歳以上35歳未満176195
35歳以上40歳未満85139
40歳以上45歳未満90106
年齢層とバグの原因との間に相関はあるだろうか。

chi2test([[176,195],[85,139],[90,106]],3,2),numer;
これを実行すると、次の出力が得られる。
chi^2= 5.35085981290286
95% confidence at 5.991464547107982
90% confidence at 4.605170185988094
従って、有意水準を95%とすると、年齢層とバグの原因が独立であるという仮説は棄却されない。
有意水準を95%とすると、である。

See more ...

Posted at 22:53 in 数学 | WriteBacks (0)
WriteBacks

Dec 16, 2008

統計学復習メモ6: χ2検定で独立性の検定

以前のメモではカイ2乗検定の例としてばらつきの検定(適合度の検定)をやってみたが、もう1つのカイ2乗検定の使われ方として、独立性の検定がある。これもやってみる。

プログラマー35歳定年説というのがある。あるプロジェクトにおいて、
35歳未満で、受けたバグ指摘件数が20件未満のプログラマーが45人、
35歳未満で、受けたバグ指摘件数が20件以上のプログラマーが15人、
35歳以上で、受けたバグ指摘件数が20件未満のプログラマーが25人、
35歳以上で、受けたバグ指摘件数が20件以上のプログラマーが15人、
だったとする。表にすると、

O(a,b)B<20B≧20
age<354515
age≧352515
(Bはバグ票数)である。暗雲が立ちこめてきたが、信頼度を95%として、35歳以上であることとバグ票数20件以上であることに相関はあるだろうか。

まず、age≧35とB≧20が独立の場合の理論値を計算する。age<35が60人、age≧35が40人、B<20が70人、B≧20が30人なので、これだけから考えると、B<20の70人は60:40でage<35とage≧35に分かれるので、それぞれ42人、28人である。同様に、B≧20の30人は18人、12人に分かれる。

E(a,b)B<20B≧20
age<354218
age≧352812
帰無仮説を、相関がない(独立である)と取ると、観測値が理論値から大きくばらついていれば棄却されるので、ピアソンのカイ2乗検定が使える。χ2分布に従う、ピアソンの検定統計量
T=¥sum_{i=1}^n¥frac{(O_i-E_i)^2}{E_i}
を計算すると、T=(45-42)2/42 + (25-28)2/28 + (15-18)2/18 + (15-12)2/12≒1.79となる。このTは、4つの理論値を固定すると、4つの観測値のどれか1つが決まると定まるので、自由度は1である。自由度が1のχ2分布の左側が95%になるχ2の値は3.84なので、今回のTはばらついていない方(χ2が小さい=理論値の通りに近い方)の95%に入っており、仮説は棄却されず、相関があるとは言えないことになる。

これは、B≧20かどうかで分けたのが明暗を分けた可能性があり、もしB≧25で分けると

O(a,b)B<25B≧25
age<35555
age≧353010
だったりした場合はT=5.23でOuchとなるので、あまりいい例ではないかも知れないが、そこまでは敢えて考えないことにする。


同じ考え方で、2つの属性で2x2に分類したデータから2つの属性の独立性を検定する方法を定式化してみる。1か2かを取る2つの属性をa,bとし、母数がnの、それぞれの属性の組み合わせに属する標本数の観測値を

O(a,b)b=1b=2
a=1x11x12
a=2x21x22
(x11+x12+x21+x22=n)とする。a,bを独立とした場合の期待値は、b=1について(a=1):(a=2)=(x11+x12):(x21+x22)、a=1について(b=1):(b=2)=(x11+x21):(x12+x22)…となるので、
nE(a,b)b=1b=2
a=1(x11+x12)(x11+x21)(x11+x12)(x12+x22)
a=2(x11+x21)(x21+x22)(x12+x22)(x21+x22)
の定数倍である。これらを全て足すと(x11+x12+x21+x22)2=n2なので、これらが期待値のn倍であり、1/nすれば期待値になることがわかる。
E(a,b)b=1b=2
a=1(x11+x12)(x11+x21)/n(x11+x12)(x12+x22)/n
a=2(x11+x21)(x21+x22)/n(x12+x22)(x21+x22)/n
これを、
O(a,b)b=1b=2sum
a=1x11x12a1
a=2x21x22a2
sumb1b2n
という風に行和、列和を取って整理すると、
E(a,b)b=1b=2
a=1a1b1/na1b2/n
a=2a2b1/na2b2/n
となる。検定統計量Tは、各セルの(O-E)2/Eを合計したものなので、
T=¥sum_{i=1}^{2}¥sum_{j=1}^{2}¥frac{(x_{ij}-¥frac{a_ib_j}{n})^2}{¥frac{a_ib_j}{n}}
と書ける。これを自由度1のχ2分布のχ2の値として見れば良いわけである。

aがp個の値、bがq個の値を取る時も同じことが言えるので、
T=¥sum_{i=1}^{p}¥sum_{j=1}^{q}¥frac{(x_{ij}-¥frac{a_ib_j}{n})^2}{¥frac{a_ib_j}{n}}
(aiはa=iの標本の和、bjはb=jの標本の和)
となる。自由度は、aについては(p-1)、bについては(q-1)なので、全体としては(p-1)(q-1)である。

以上を踏まえて、もう1つ例題をやってみる。原理はわかったので、手計算チックな計算は卒業して、次はMaximaに組み込みのカイ2乗検定関数を使おう、と思って関数を探すと、手元のMaximaには入ってなかった。ExcelにはCHI2TEST()という名前でデフォルトで入ってるのだが、Maximaにはデフォルトでは入っていないらしい。上記の計算だけなら大して複雑じゃないので、自分で関数を作ってみる。

chi2test(x_,row_,col_):=block([i,j,n,a,b,chi2],
    a:map(lambda([i],sum(x_[i][j],j,1,col_)),create_list(i,i,1,row_)),
    b:map(lambda([j],sum(x_[i][j],i,1,row_)),create_list(j,j,1,col_)),
    n:sum(a[i],i,1,row_),
    chi2:sum(sum((x_[i][j]-a[i]*b[j]/n)^2/(a[i]*b[j]/n),j,1,col_),i,1,row_),
    print("chi^2=",chi2),
    print("95% confidence at",quantile_chi2(0.95,(row_-1)*(col_-1))),
    print("90% confidence at",quantile_chi2(0.90,(row_-1)*(col_-1))),
    chi2
)$
試しに
chi2test([[45,25],[15,15]],2,2),numer;
(numerは結果が分数になるのを避けるため)とやってみると、
chi^2= 1.785714285714286
95% confidence at 3.841458820694166
90% confidence at 2.705543454095465
と出力される。信頼度90%でも、2つの属性が独立である仮説は棄却されないことになる。

次の例は、バグの原因を設計ミス要因かコーディングミス要因かに分けて、さらにそれをプログラマーの年齢層で分けて数えてみたという架空のデータである。

O(a,b)設計ミスコーディングミス
30歳以上35歳未満176195
35歳以上40歳未満85139
40歳以上45歳未満90106
年齢層とバグの原因との間に相関はあるだろうか。

chi2test([[176,195],[85,139],[90,106]],3,2),numer;
これを実行すると、次の出力が得られる。
chi^2= 5.35085981290286
95% confidence at 5.991464547107982
90% confidence at 4.605170185988094
従って、信頼度を95%とすると、年齢層とバグの原因が独立であるという仮説は棄却されない。
信頼度を95%とすると、である。

See more ...

Posted at 21:50 in 数学 | WriteBacks (0)
WriteBacks

Dec 07, 2008

「〜ますでしょうか」

仕事でたくさんのメールを書いていると、日本語的にというか敬語的にというか、どう書くのが適切かわからなくなることがしばしばある。元々、筆者は国語が苦手なのである。中でも、「よろしいでしょうか」「合ってますでしょうか」などの「〜でしょうか」は長きに渡り決まり文句のように使い続けているのだが、書く度に違和感がある。

なぜ「よろしいでしょうか」に違和感があるかというと、語尾を肯定文にかえると「よろしいです」であり、デアル調に変えると「よろしいである」だからである。ファミ・コン語の「よろしかったでしょうか」は論外だが、「よろしいですか」より「よろしいでしょうか」の方が柔らかい気がするから使ってるのだが、「よろしいですか」でも問題は同じである。では「よろしいですか」の正しい敬語は何なのだろうか?

「合ってますでしょうか」については、「合ってますか」より柔らかい気がするから使ってるのだが、「合ってますですか」と考えると敬語が2重だし、これで既におかしいが、肯定文にして口語体に変えると「合ってるです」「合ってるだ」だ。

「合ってますでしょうか」などと書く時に時々思い出すことがある。
ある将棋の解説で、驚きをユーモラスに表す為に、解説者がわざと「これは…恐ろしいことに…寄ったですね〜。」、その後「これは…詰んだですね〜。」と言ったのを聞き、笑いながらもその表現にとても共感したことがある。その局面は、私にとっては「寄りましたね」→「詰みましたね」と表現されるより、「寄ったですね」→「詰んだですね」と表現されるべきものだったのである(関係ないが、「〜べき。」で終わる表現も違和感がある。「〜べし。」とすべし。じゃないのか?)。解説者の横の聞き手が「詰んだですか。」と返したかどうかははっきりと覚えていないが、その時に限り「詰んだですか」はありだったと思っている。
なぜその状況において「詰んだですね」しか有り得ないと思ったのだろう、私の言語感覚がおかしいからだろうか、とずっと考えている。「詰んだですね」の口語体は「詰んだだね」だから、おかしいことはおかしい。今の所、「これは『詰んだ』ですね」という伝聞調だったのだろうと思っているが、まだ自分で納得できていない。
とりあえず、そんな私の言語感覚でも「詰んだですね」は稀にしか正しくないのだから、「よろしいでしょうか」「合ってますでしょうか」も基本的に正しくないと思う。

閑話休題、「よろしい」の敬語は「よろしゅうございます」と読んだことがあるが、するとその疑問形は「よろしゅうございますか」だと思うが、私の仕事関係で実際に使われているのは見た記憶が無い。「合ってますでしょうか」は「合っておりますか」が正しいと記憶しているが、ちょっときつい気がする。「合っておりますでしょうか」だと元の問題が解決してないし、「合ってございますか」はやはり見た記憶が無い。「合ってます?」は結構見かけるが、また別な違和感があるので、私には書けない。

私には何年も前から気になっていて、おそらくその辺りには敏感になっているにも関わらず、未だ解決していないということは、少なくとも私が読んできたメールには私が納得する答えが無かったということだろう。敬語に無頓着な人以外は「よろしい」の疑問形や「合ってます」の疑問形を使わない表現方法を身に付けているのだろうか。

明日からまた「合ってますでしょうか」と書くんだろうな。気を緩めていると「合っておりますでございますでしょうか」等に進化してしまうかも知れない。

See more ...

Posted at 23:03 in 雑記 | WriteBacks (2)
WriteBacks

初めておじゃまいたします。本日URLを教えていただいたので早速遊びに来ました。 このエントリー、笑い転げながら読みましたよ~ (^o^) 私も改まった文章を書く際にはいつも違和感と戦っているので 似たようなことを考えている人がいるとわかってうれしいです♪ 「合ってますか?」を柔らかく書きたい場合の決まり文句として 最近は「認識合いますでしょうか?」を愛用していますが、 いい表現に出会うたびに変えるので自分の中でも流行り廃れがありますね~。

Posted by ゆきねこ at 06/24/2010 12:03:19 AM

早速の御チェケありがとうございます。 認識合いますでしたようで良かったでした。

Posted by ynomura at 06/24/2010 10:51:43 PM

Dec 06, 2008

統計学復習メモ5: χ2分布とχ2検定

以前のエントリーにて実際にやってみたが、カイ2乗分布とカイ2乗検定に軽く触れてみると、χ2の定義が違うので混乱した。カイ2乗分布で定義されるχ2
¥chi^2=¥sum_{i=1}^n ¥left(¥frac{x_i-¥mu_i}{¥sigma_i}¥right)^2
(xiは標本値、μiは平均、σiは標準偏差)であるのに対して、カイ2乗検定で使われるχ2
¥chi^2=¥sum_{i=1}^n¥frac{(O_i-E_i)^2}{E_i}
(Oは観測値(observed)、Eは期待値(expected)または理論値)である。何でそれらをχじゃなくわざわざχ2と置くんだ、ということも少し引っ掛かるが、それは分散のσ2の例もあるので、値の意味的にきっと何か2乗的、2乗系、2乗fulなんだろうと思って通り過ぎるとして、2つのχ2の式は関連がありそうに見えて何と何が対応するのかよくわからない。Σ内の分母は上のが分散なのに対して、下のは期待値である。分散で割るのと期待値で割るのは違うであろう。分母だけ見るからそう思うんだろうかと思ってΣ内全体を見ても、上のはお馴染みの(xが正規分布N(μ,σ2)に従う時の)標準正規分布N(0,1)に従う値の2乗であるし、下のは平均からの差の2乗を平均で割ってて、平均の大きさに比例しそうな、2乗的には見えない、奇妙な値である。まるで1つ期待値で割り忘れたかのような形をしている。

結論から書くと、上の2つの式は、何か対応がありそうな形はしているが、真正面からパーツ毎に対応させてようと考えるのはやめた方が良さそうだ。
落ち着いてWikipediaを眺めていると、「ピアソンのカイ2乗検定」という言葉が目に入った。カイ2乗検定というのも色々ある中で、以前のエントリーに書いた検定の例は、最もポピュラーな、ピアソンのカイ2乗検定というやつらしい。数学的に説明するのは大変難しいが、検定統計量
T=¥sum_{i=1}^n¥frac{(O_i-E_i)^2}{E_i}
がカイ2乗分布に従うことがピアソンによって示されているので、カイ2乗分布を使って、この検定統計量がある値以上になる確率を知ることができるということだ。
Tは、各事象の実際の発生回数と期待値とのずれ(ばらつき)が大きいほど大きくなるので、観測されたTが十分低い確率でしか起こらないかどうか、すなわちそのばらつきが十分低い確率でしか起こらないかどうかを、このピアソンの方法で知ることができる。
逆に、実際の発生回数を真実、期待値を仮説とすると、Tが十分大きければ、もしその仮説が正しいとするとかなり稀なことが実際には起こったということになり、むしろその仮説が間違っていると考える方が合理的であるので、ある仮説を確率的に誤りと見なす(仮説を棄却する)のに使えるということである。

N(0,1)に従う標本の2乗和と、カイ2乗検定の検定統計量の両方がχ2と定義されるから、学ぶ人が混乱するんだ、と思うのは筆者だけだろうか。きっと何か合理的な理由があるんだろうが、それならその理由を先に知っておかないと混乱の元になると思うが、筆者には未だに見つけられていない。

See more ...

Posted at 23:41 in 数学 | WriteBacks (0)
WriteBacks

Dec 04, 2008

Windows XPでのかな入力のトラブル記

今年の春、会社で異動になり、転勤先で新しいパソコンを買ってもらった。仕事柄、メールをよく打つようになった。秋には仕事にも慣れてきて、横の席の人がメールを打っている音が耳に入ってくるようになった。すごい速さでキーボードをヒットしている。画面上に日本語が何文字か増える間にも、やたらキーボードを叩いている気がした。
すると、自分のタイピングも気になり出した。何せローマ字入力は指が完全に覚えているし、文字数に対して打鍵数が多い。「させて頂きます」は"SASETEITADAKIMASU"、「申し訳ございません」は"MOUSIWAKEGOZAIMASENN"、「よろしくお願い致します」は"YOROSHIKUONEGAIITASHIMASU"である。我ながらよくこんなにたくさんキーを打ってるものだ、と思った。1つのメールを書く間の打鍵量は、数えてみたらすごいことになるのではないか。
そんなことを意識すると、キーボードのJIS配列のかなキーに目の焦点が合うようになってきた。9割以上の人がローマ字入力で日本語を打つこの時代に、まだあったのか。

よく見ると、我が家のパソコンのキーボードにも、アルファベットキーや数字キーを中心に、キーの1つ1つにひらがなが刻印されている。数字キーには「ぬふあうえおやゆよわ」、アルファベットキーの1段目はUキーから左に読んでいくと「なんかすいてた」である。

20年以上前に買ったMSX2パソコンにも、JIS配列が書かれていた。「12345」が「アイウエオ」、「QWERT」が「カキクケコ」であるアイウエオ配列も書かれていたので、JIS配列で打つ訳が無かった。
大学で触ったパソコンには、ローマ字かな入力があり、実質それが標準だったので、私もローマ字入力で日本語を打ち始めた。同時に電子メールも使い始めたので、指が完全にローマ字入力になった。
以後、キーボードにJIS配列の刻印があっても、まず意識することは無かった。

1日の日本語入力量が増えて、ますますキーボードのひらがなの刻印にピントが合うことが増えた。子音、母音、子音、母音、子音、子音、母音、子音、子音、子音、母音。日本語入力中もローマ字入力のために頭の中はアルファベットが流れ続けている。指が押し続けているアルファベットにはひらがなが書いてある。

私には、かな入力を試さないことは不可能だった。

家で練習を始めてから1ヶ月ちょっと経ったが、時間に余裕がある時に会社でかな入力を試し始めると、やはり打つのに時間がかかる。指が配列を十分に覚えていないことも、英数字や記号を打つのが面倒なのも要因だが、致命的な問題が、時々勝手にかな入力モードから英数入力モードに切り替わってしまうことだ。
会社で使っているパソコンはWindows XPで、日本語入力に関しては何も設定を変更していない。ALT+ひらがなキーを押してローマ字入力とかな入力を切り替えているだけである。それでいて、漢字変換しながらひらがなを入力し続けているだけで、ふと画面を見ると「スケジュール4a30p」と表示されていたりする。間違って英数入力モードに切り替えてしまったのかと思って、反射的に「半角・全角/漢字」を押してもう一度「うちあわせ」と打ち込むと、画面には「4a30p」と出る。
突然かな入力モードから全角の英数入力モードに切り替わるのである。しかも、言語バーの表示は「A」ではなく「あ」のままなので、いつ英数入力モードに切り替わったのかは打ってみないとわからない。よく起こる時は1分に1回くらい起こる。何度も何度もそんなことが起こると、時間のロスになるだけでなく、平常心を保てなくなる。使い物にならない。

1つのパターンとして、Internet Explorerの入力フォームで、かな入力して、スペースを何度か叩いて漢字変換の候補枠を出して、カタカナを選んで確定して、数字キーに割り当たっているひらがなを打ち込もうとすると、100%の確率で発生することが判った。「スケジュール」を候補窓で確定して入力し、「う」キーを押すと全角英数に切り替わってしまうのだ。これを防ぐ方法は、候補窓でカタカナを確定した後は「半角・全角/漢字」を2回押してリセットすることしか見つけていない。しかも、この現象が起こるのはその手順だけではない。
不思議なことに、ほとんど同じ状況を家のPCで再現しても、家のPCでは全く起こらない。OSはどちらもWindows XP、日本語入力FEPはどちらもデフォルトのMS-IMEである。

Webで調べ回った限りでは、正確な原因には辿り着けなかったが、同様の症状に遭った人は少なからず居るようで、MS-IMEのバグというのが定説のようであり、確実にこれを回避する方法はただ1つ、コントロールパネルの「地域と言語のオプション」で「詳細なテキスト サービスをオフにする」を有効にすることなのだそうである。

確かに、その設定をすることにより、上記の問題は全く起こらなくなった。しかし、この設定をすると、言語バーが貧弱になり、半透明にならなくなったり、タスクバーに入れても出しても見えなくなったりして、とても悲しくなる。なんともはや、という感じである。JIS配列は至る所のキーボードに刻印されていても、既に見捨てられているのだろうか。

参考:http://www.atmarkit.co.jp/fwin2k/win2ktips/630ctfmon/ctfmon.html
本件のバグに関する直接的な記述は無いが、色々バグがあるようだ。

See more ...

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

リカバリーかけると 直るらしいですよ。

Posted by もちっこ at 04/23/2009 12:02:15 PM

そういう問題ではないようです。

Posted by ynomura at 04/23/2009 09:10:06 PM

Nov 30, 2008

統計学復習メモ4: (思考実験)学力偏差値と合格率

統計解析の入門書を初めから読み返していると、偏差値の説明がちょろっとだけ書かれていた。ああ、あの高校受験や大学受験の模試で出てきたやつだな、と反射的に思った。
筆者が受験生だったのは遠い昔で、入試というものはほとんど意識することは無くなったのだが、それでも私にとって「偏差値」と言えば入試である。というか、これまでの人生で入試関連以外で50が中心の「偏差値」が使われているのを見た記憶が無い。IQテストなんかは点数の分布が正規分布に従っていることを前提にしてるらしいので、これも広い意味では「偏差値」だろうが、一般的な定義では平均値の偏差値が50で、IQの平均点は100なので、狭い意味では偏差値とは言わないと思う。
従って、偏差値という文字を見ると、受験生時代に見た模試の結果の偏差値しか思い浮かばない。

そういえば、偏差値が60とか70とかってどういうことなんだろう?
我ながら貴様本当に理系だったのか?という感じだが、少なくとも私は授業で習わなかったし、試験には出てこなかった。
定義的には何のことはない、偏差値50が平均で、60が+σ、70が+2σだ。正規分布表を見ると、60で大体上位16%、70で上位2.3%らしい。やっぱり、だから何?という感じである。

私は高校時代、偏差値の意味を知らなかったので、模試の偏差値の値を意識することがなく、自分の学力偏差値がどれくらいだったのか覚えていない。しかし、共通一次から名前を変えて間もないセンター試験の点数は覚えていたので、ウェブ上で過去の平均点と最近の標準偏差を調べて(過去の標準偏差は見つからなかった)概算してみると、58くらいだった。自分が受けた大学のは覚えてない(意識すると凹むから)が、東大の偏差値が70だったのをかすかに覚えている。
仮に当時の私の学力の偏差値が58だったとして、仮にもしセンター試験で偏差値70の点を取ったら東大に合格したとすると、私の東大合格率はどれくらいだっただろうか?

そもそもどうやったらそういうのが計算できるのかがさっばりわからなかったので、統計解析の本を眺めながら考えてみた。

上の仮定と命題にもかなり無理があるが、さらに無理な仮定を重ねる。
800点満点で、受験者の点数が正規分布し、平均点が480、標準偏差が120とする。私が必ず偏差値58の点を取ると必ず落ちるので、私の点数も正規分布するとする。標準偏差が不明だが、感覚的な理由で40だったとする。支離滅裂になってきたが、続ける。つまり母集団の点数がN(480,120^2)に従い、私の点数がN(偏差値58の点,40^2)に従い、私の点数が偏差値70の点を超えると東大合格だった、という仮定である。偏差値58とか70とかはN(480,120^2)における話なので、偏差値xの点は480+120*(x-50)/10である。
Maximaに

cdf_normal(480+120*(70-50)/10, 480+120*(58-50)/10, 40),numer;
(私の点数の分布における、偏差値70の点までの累積密度)と入力すると、0.99984という答えが出たので、私の東大合格率は0.00016、すなわち0.016%、6250中1回の確率ということになる。毎日東大を受験しても17年に1回しか合格しない確率だ。

もし私が1浪して偏差値を5上げたとすると、

cdf_normal(480+120*(70-50)/10, 480+120*(63-50)/10, 40),numer;
は0.982と出てくるので、合格率は1.8%である。

…命題と仮定に無理があり過ぎて、話にならないのは、言うまでもない。
とりあえず、こういう計算にはいろんなデータが必要になることがわかった。

See more ...

Posted at 21:32 in 数学 | WriteBacks (0)
WriteBacks

Nov 29, 2008

統計学復習メモ3: 視聴率調査と標本数

TVの視聴率ってどうやって測ってるんだろう?というのはよく聞かれる問いであるし、モニターとして選ばれたいくらかの数の一般家庭に視聴率計測用の機器を取り付けて測っているのだ、というのはよく言われる答えである。そのモニターの数は、一例では関東・関西で600世帯、その他の地域は200世帯、とも言われる。

たったそれだけ?ある番組を20人が見てれば視聴率10%?学校の1クラスより少ない。それではたまたまその番組が好きな人が2人多めや少なめにモニターに入ってれば、すぐ1%くらい変わってしまうじゃないか、そんな値を信用できるのか、そんな値に意味はあるのか?
…とか書くと、貴様本当に大学を理系で卒業したのか?と問われてしまいそうであるが、それでもやはり反射的には疑問に思ってしまうことを自白したい。

このいわゆる標本調査、標本がいくらぐらいだとどれくらいの信頼性があるのだろうか。また、誤差をある程度に抑えるには、どれくらいの標本数が必要なのだろうか。

これを考えるのにも、区間推定の式が使える。母比率の推定区間は[標本比率±t*√(標本比率*(1−標本比率)÷標本数)]で、このtは標準正規分布に従う、と教えられている。数式で書くと、
I \in p\pm t\sqrt{\frac{p(1-p)}{n}} …(1)
(Iは推定区間の下限と上限、pは標本比率、nは標本数)である。信頼度を90%とするとt=1.645であり、n=200で、ある番組を見てたのが20人だった場合、p=0.1ということなので、母比率の推定区間は0.1±1.645*√(0.1*0.9/200) ≒ 0.1±0.035、すなわち6.5%〜13.5%となる。結構大雑把だ。信頼度を95%とするとt=1.96なので、5.8%〜14.2%ともっと幅が広くなる。
この式の±以降の部分が推定区間の幅なので、標本比率を固定すると、誤差をどれくらいにするためにはどれくらいの標本数が必要かを求めることができる。最大誤差をeとすると、
e=t¥sqrt{¥frac{p(1-p)}{n}}
なので、これをnについて解くと、
n=¥left(¥frac{t}{e}¥right)^2p(1-p)
である。例えば、視聴率10%の近辺で、信頼度を90%で視聴率の誤差を1%以内にするためには、n=(1.645/0.01)2*0.1*0.9=2435と求まるので、2400人以上のモニターが必要ということになる。誤差が2%で約600人、3%で約270人である。視聴率が20%近辺だと、さらにハードルが上がり、約4300人で誤差1%、約1100人で誤差2%、約500人で誤差3%となる。視聴率が50%に近づくほど、必要なモニターが増えることがわかる。

ところで、最初の母比率の推定区間の式(1)は、どうしてそうなるのだろうか。

以前のエントリーで、正規分布に従う値の母平均の推定区間が[標本平均±t*√(標本分散÷標本数)]で与えられる、と書いた。数式で書くと、
I \in \bar{X}\pm t\sqrt{\frac{s^2}{n}}
(Iは推定区間の下限と上限、Xは標本、s2はXの不偏分散、nは標本数)という感じであるが、この平方根部分の中身は要するに標本平均の分散ということであるので、ちょっと乱暴であるが、母平均に限らず、標本から求められる推定値からの推定区間が[推定値±t*√(推定値の分散)]で与えられる、と理解していいと思う。
なお、tはt分布に従うが、自由度が大きい場合はt分布が標準正規分布に近似できることになっている。

この推定値が比率の場合、その比率の分散が何になるかが問題だが、これは二項分布の定理から求められる。発生確率がpの場合、n回の試行で発生する回数Xの期待値はE(X)=np、分散はV(X)=np(1-p)だと二項分布の定理で確定しているので、n個の標本中の確率pで存在する何かの個数Xの期待値もE(X)=npであり、分散もV(X)=np(1-p)であることから、存在確率X/nの期待値はE(X/n)=p、分散はV(X/n)=p(1-p)/nだとわかる。これを使って、[推定値±t*√(推定値の分散)]=E(X/n)±t√V(X/n)=(1)となる訳だ。

See more ...

Posted at 21:44 in 数学 | WriteBacks (0)
WriteBacks

Nov 16, 2008

統計学復習メモ2: t分布と不偏分散

t分布を、正規分布に従う標本の平均が従う分布である、という感じに捉えて、なんとなく理解した気になると、1つ疑問が湧いた。正規分布に従う複数の値の和はやはり正規分布に従うのではなかったか?それなら、値の和を変数の数で割ったもの、つまり平均もt分布でなく正規分布に従うはずではないのか?
x1,x2をそれぞれ正規分布N(μ112),N(μ222)に従う確率変数とすると、x1+x2はN(μ121222)に従う、つまり和の平均は平均の和、和の分散も分散の和となる、というのを正規分布の性質の1つとして教えられている。x1+x2が正規分布に従うなら(x1+x2)/2も正規分布(この場合N(\frac{\mu_1+\mu_2}{2},\frac{\sigma_1^2+\sigma_2^2}{4}))に従うはずである。
というか、正規分布N(μ,σ2)に従うn個の標本の平均はN(μ,σ2/n)に従う、つまりn個の平均を取ると分散が1/nになる、という定理を覚えていないと、t分布の式を読むのが難しいと思う。それからしてもやはり標本の平均は正規分布に従うはずである。

そんなことで混乱したのは私だけだろうか?

混乱の原因は、標本平均、標本分散と母平均、毋分散をごっちゃにしたことだった。母平均をμ、毋分散をσ2、標本をx1...xn、標本平均をm=\bar{x}=\frac{1}{n}\sum_{i=1}^{n}x_i、標本分散をs_n^2=\frac{1}{n}\sum_{i=1}^{n}(x_i-\bar{x})^2とすると、mはあくまでN(μ,σ2/n)に従うのであって、N(m,Sn2/n)に従うのではないのである。つまり、μやσが不明な場合はmが従う正規分布も不明なのである。
そこで、mやSn2からμの範囲を推定するのにt分布が出てくる。
(※本エントリーではシステムの都合上、標本平均をmと書いたりxバー(xの上にバー)と書いたりするが、両者は同じものである)

…という感じに混乱を解決するのにも私には結構時間がかかったのだが、その間に私の理解を妨げたものの1つが、不偏分散の概念である。
標本分散が
s_n^2=\frac{1}{n}\sum_{i=1}^{n}(x_i-\bar{x})^2
であるのに対して、母平均の区間推定に必要になる不偏分散は
s_{n-1}^2=\frac{1}{n-1}\sum_{i=1}^{n}(x_i-\bar{x})^2
と教わる。nで割るかn-1で割るかの違いだけである。なぜn-1で割るのかというと、「自由度がn-1だから」と教わる。そんなのは丸覚えすれば引っかかる必要の無いことであるが、私の脳はそこで「???」となり先に進まなくなった。

「推定値の不偏分散の自由度が(n-1)である」理由は、ばらつきの指標の元となる
Σ(xi-m)2
のmにm=1/n Σxiという縛りがあるため、mとx1...xn-1が決まるとxnが決まるので、この2乗和がn-1個分しか動かないから、と説明されることが多い。極端な例としては、n=2の時、本当の平均μがわかってれば
Σ(xi-μ)2
はx1の分もx2の分もばらつきに応じて大きくなるのに対して、
\sum_{i=1}^{n}(x_i-m)^2 = \sum_{i=1}^{n}(x_i-\frac{x_1+x_2}{2})^2
は計算すると
(x1-x2)2/4
となり、標本が2個あるのに、その間の差1個分しかばらつきの値として反映されない。同様の理由により、標本をn個としてもn-1個分しかばらつきの値として加算されないのである。

その説明によって、私は一瞬目から鱗が落ちた気がしたが、まだ納得できなかった。もし標本の数を母集団の大きさと同じにすれば、つまり母集団の全てをサンプルとして拾えば、mが本当の平均になるから、自由度がn-1でも分散は(1/n)Σ(xn-m)2ではないか?
…ここから先はまだ疑問が解けていないが、おそらく推定値の分散と実際のデータの分散は同じ尺度では測れないということのような予感がする。

ところで、2008/11/16時点の日本語のWikipediaの分散の項に、不偏分散の計算、というのが載っている。曰く、
\frac{1}{n} \sum_{i=1}^{n}(x_i - \mu)^2

= \frac{1}{n}\sum_{i=1}^{n}(x_i - \bar{x})^2+\frac{1}{n}\frac{1}{n} \sum_{i=1}^{n}(x_i-\mu)^2
になるので、
\frac{1}{n} \sum_{i=1}^{n}(x_i - \mu)^2= \frac{1}{n-1}\sum_{i=1}^{n}(x_i - \bar{x})^2
となるとのことである。これだ!と思って計算を追ってみると、途中の
\frac{1}{n} \sum_{i=1}^{n}(x_i - \mu)^2
= \frac{1}{n}\sum_{i=1}^{n}(x_i - \bar{x})^2+\frac{1}{n^2} \{ \sum_{i=1}^{n}(x_i-\mu)^2 +2\sum_{i\neq j}^{n}(x_i-\mu)(x_j-\mu) \}
までは確認できたが、その次の
= \frac{1}{n}\sum_{i=1}^{n}(x_i - \bar{x})^2+\frac{1}{n}\frac{1}{n} \sum_{i=1}^{n}(x_i-\mu)^2
はわからなかった。というか、少なくともn=2として展開してみる限り、そうはならない。=が≒の間違いなのかも知れないが、よくわからない。

結局、私は次のように理解した。上のWikipediaの計算と同じ方針で、x_i=(x_i-\bar{x})+\bar{x}として両辺の分散を取ると、分散の加法性から
\sigma^2=Var(X-\bar{x})+\frac{\sigma^2}{n}
となる。Var(...)の部分は標本分散なので、これを毋分散σ2について解くと、
\sigma^2 = \frac{n}{n-1}Var(X-\bar{x}) = \frac{n}{n-1}\frac{1}{n}\sum_{i=1}^{n}(x_i - \bar{x})^2
= \frac{1}{n-1}\sum_{i=1}^{n}(x_i - \bar{x})^2
となる。これが標本から推測される毋分散であり、毋分散と同じ尺度で使える不偏分散ということである。

See more ...

Posted at 21:44 in 数学 | WriteBacks (0)
WriteBacks

Nov 09, 2008

統計学復習メモ1: t分布とχ2分布

高校時代は確率が得意中の得意だった私が、大学の統計学の講義にはすぐついていけなくなった。高3の時は、授業で順列と組み合わせの公式を覚えた途端、ほとんど全ての問題が難なく解けた。周りには理解につまずく人がいたが、私にはなぜ理解できないのかがわからなかった。大学2年の時に受けた統計学の講義も、半年くらいは割と楽についていってたが、ある時点から突然理解できなくなった。それは、平坦な道に突如現れた高い壁だった。どちらかというと興味があった科目なので、教科書をじっくりと何度も読んで理解しようとしたが、特に数式や日本語が難しい訳でもないのに、全然頭に入らなかった。なぜ理解できないのかがわからない不思議な感覚だった。
「推定・検定」という名前の関門だった。

結局、大学で統計学の講義を3回くらい取って、毎回1から学んだが、いつも確率分布までは理解できるのに、「推定・検定」で手も足も出ず挫折した。3回目の統計学で悟ったのは、t分布とカイ2乗分布が障害ということだった。

それから幾年月、それ以来統計学に接する機会が無かったが、推定・検定のことがずっと気にかかっていた。大学時代に最も興味があった科目の1つでありながら全く私を寄せ付けず、私に最も挫折感を味わわせた推定・検定を克服しない限り、死んでも死に切れない。今更そんなものを勉強して何の役に立つのかと言われようが関係ないのである。

という訳で先々月ぐらいからリベンジを試みているが、やはり相性が悪い。何となく少し理解した気になっても、2週間もするとすっかり忘れてしまい、先に進まない。後半が推定・検定でそれで終わりの入門書チックな統計解析の本を買って読んでるが、後半のみもう4周目くらいである。私以外のこれを読む人には前半と後半、確率分布と推定・検定の間に落差は無いのだろうか。

無駄話はこの辺にして、勉強したことをメモする。

t分布とは、xを平均μ、分散σ2の正規分布N(μ,σ2)に従う標本、(xの上にバー)を標本平均、nを標本の数、sを標準偏差(不偏分散の平方根)とした時、
t=¥frac{¥overline{x}-¥mu}{s/¥sqrt{n}}
というtが従う分布であり、
χ2分布とは、
¥chi^2=¥sum_{i=1}^n ¥left(¥frac{x_i-¥mu_i}{¥sigma_i}¥right)^2
というχ2が従う分布である、というようなことがどこにでも書いてある。
入門書で無ければ、さらにこれらの分布の確率密度関数が、ガンマ関数を使った容赦ない鮮烈な式で書かれる。
これが推定・検定の本題の前に出てくるから厄介だ。

どちらもにわかにはイメージが掴みにくい数式であり、なぜこういうものが必要かを説明される前に、まずこれを覚えて、と言われるのは結構辛いと思う。
それらの分布の意味を理解する前に、使い道を1つくらい知った方が、最初の理解の足掛かりになるし、トータルとして速く理解できるのではないだろうか。

という訳で、t分布の使用例として、母平均の区間推定を考えてみる。
大勢の人が受けたテストの点が正規分布に従うとする。その中の10人の点数(標本)が
70, 60, 40, 80, 40, 60, 30, 80, 50, 90
だとすると、全体の平均点は95%の確率でどの範囲に入ると推測できるだろうか。

ここで、その範囲が(標本平均±t*√(標本分散÷標本数))で求められるというのが、区間推定のありがたい教えである。

上の標本の場合、平均が60、分散が400、自由度(n-1)=9で信頼区間が95%となるtが2.262なので、推定区間は60±2.262√(400/10) ≒ 60±14.3ということになる。標本の平均は60点だが、95%の確率で全体の平均がこの範囲に入るということだ。

標本平均の分散が(標本分散÷標本数)になるというのは定理としてあるので、t分布のtの値というのは、標本の標準偏差を何倍すれば推定区間の幅になるのかを教えてくれるものと考えられる。

χ2分布の使用例として、ばらつきの片側検定を考えてみる。
プログラマーがそれぞれ2,3,3,4,8人の5つのソフトハウスにソフトウェアの開発を委託してるとし、あるプロジェクトにおいてそれぞれのソフトハウスが発生させたバグの数が15,18,14,19,34だったとする。バグの数がプログラマーの数に比例するという仮定は間違ってる方の5%に入る(95%の確率で間違ってる)だろうか。

ここで、((測定値−期待値)の2乗÷期待値)の和、という値がχ2分布に従うので、χ2分布を見ればその発生確率がわかる、というのが、χ2検定のマジカルかつミラクルな教えである。

上の標本の場合、バグが計100件でプログラマーが計20人なので、1人当たり5件、ソフトハウス別では10,15,15,20,40件となるのが期待値である。従って、測定値のχ2の値は、
(15-10)2/10 + (18-15)2/15 + (14-15)2/15 + (19-20)2/20 + (34-40)2/40 ≒ 4.1
である。自由度(5-1)=4で棄却域(有意水準)が5%となる、χ2分布におけるχ2の値は9.49なので、測定値のχ2の値はそれより小さく、95%の確率で起こり得る範囲内であることがわかる。
従って、バグの数がプログラマーの数に比例するという仮定は間違ってるとは言えない。

χ2分布のχ2の値というのは、同時に起こらない複数の事象の発生回数と各々の回数の期待値がわかる場合に、その期待結果とのずれがどれくらいの確率で起こるものかを調べるのに使える、発生回数のばらつきの尺度だと考えられる。

See more ...

Posted at 19:33 in 数学 | WriteBacks (0)
WriteBacks

Nov 02, 2008

MacでJava3Dを使う時の注意

自作のJava3DアプレットがMacで動かないことに気付いて調べた結果、現時点でMac OS X向けに公式にリリースされているJava3Dのバージョンが1.4.2と古いことが原因だと判明した。

APIリファレンスにあるから使っていた、javax.vecmath.Vector3dとかのgetX(),getY(),getZ()がバージョン1.5以降にしかないのだ。
こういうメソッドは最初に無いものをなぜ後から追加するのだろう?最初から無いならずっと無いままでいいと思うのだが。互換性を考慮するのが面倒で仕方ない。そんなにメンバアクセスはメソッド呼び出しでできないと気が済まない人がいるのだろうか。

ともかく、仕方ないのでgetX()とかを使わないように書き直した。
3Dカッチンアプレット
 ソースコード
J3DTest2アプレット
 ソースコード

See more ...

Posted at 21:39 in Java | WriteBacks (0)
WriteBacks

Oct 26, 2008

Blosxom環境構築(2/2)

今回このblogのBlosxom版を構築するために行ったことをまとめる。

●Blosxom本体
SourceForge.netのBlosxom projectから最新版を入手した。
なお、同ページは"The Unofficial Blosxom User Group"では"SF.net Project"と書かれてリンクされている。

●プラグイン
Blosxomは本体が小規模に作られているので、欲しい機能のプラグインを入れる所から始まる。
今回はblosxom starter kit(bsk)とBlosxom Plugin Registryを参考にプラグインを選び、"The Unofficial Blosxom User Group"のV2 Plugin Registoryから最新版を取得していった。
以下、インストールしたプラグイン名と、インストールに伴って改造した点を記す。
・archives
月毎のアーカイブページと、それらへの階層状のリンクを生成する。
特に改造無し。
・back_and_forth
個別エントリーのページから、更新日付が1つ前のページ、1つ後のページへのリンクを生成するプラグイン。
bskの"«前 | Main | 後»"という並びが気に入ったので、bskのを使わせてもらった。(V2 Plugin RegistoryのはMainの部分が出せない)
・categories
カテゴリー毎のアーカイブページと、それらへの階層状のリンクを生成するプラグイン。
ディレクトリ名と表示文字列の対応($friendly_name)を埋めたのと、'all entries'を'全てのエントリー'に書き換えたのみ。
※report_dir_start呼び出しがコメントアウトされてるバージョンだと、'all entries'を表示するにはコメントアウトを解除する必要あり。
・entries_index
エントリーの作成日時を保存する。これが無いとエントリーの作成日時がエントリーファイルの更新日時になるので、FTPでwebサーバーにアップロードするとファイルの日時がアップロード時刻になり、エントリーの作成日時が失われてしまう。
特に改造無し。(インデックスファイルは別途作成して一緒にアップロードした)
・entry_title
個別エントリーページのタイトル(<TITLE>タグの中身)をそのエントリーのタイトルにする。
"ブログ名 セパレータ タイトル"という並びが気に入ったので、Kyo Nagashimaさんのバージョンを使用。特に改造無し。
・find
ブログ内をテキスト検索して、マッチするエントリーだけを出力するプラグイン。
特に改造無し。
・geek
(出力例)。いらないけど消すには惜しいので残しておく。
・seemore
複数エントリーが表示されるページで、各エントリーを丸ごと表示するのでなく、部分的に表示して、残りの部分へのリンクを付けるようにする。
"See more ...."を「続きを読む」に変えかけたが、"See more ...."が格好よくて気に入ったので、改造無し。
・titles_index
複数エントリーが表示されるページで、各エントリーをタイトルのみで表示するようにする。
seemoreで表示するトップページ以外では表示エントリー数の制限を緩和するため、blosxom本体のgenerate{}の"# Stories"の所の

my $ne = $num_entries;
my $ne = $currentdir ? 100 : $num_entries;
に変更。
・writeback
後述。

●Writebackスパム対策
Blosxom公式サイトのWritebackプラグインを使うと、スパムの餌食になる。対策として、以下の5つのプラグインを見つけた。
(1) WritebackPlus
(2) Feedback
(3) spam_blocker
(4) hail2u.netのWriteback
(5) Writeback with Akismet
まずデファクトスタンダードと思われる(1)のWritebackPlusを試そうとしたが、その中で使われているMT-blacklistが既に配布されておらず(MovableType 3.2で他の方法が使われるようになったため)、スパム対策にならなかった。
ネット上で流し読みする限り(2)はマイナーで、(3)と(5)は導入が面倒くさそうで、(4)の方法で十分効果がありそうなので、今回はWritebackPlusに(4)の方法を適用することにした。
WritebackPlusの改造結果
(WritebackPlusはVersion 0.2.1aを使用)
改造点(コメントで"Y.Nomura"とある部分)
・250行目付近の$pathを修正
・257行目付近、(4)のWritebackの100~131行目をコピー
・305行目付近、blacklist無しでも動くようにコメントアウト
・375行目付近、Net::SMTP無しでも動くようにコメントアウト(Hi-Hoのサーバーに無かったため)

●MovableTypeのエントリーのインポート
MovableTypeの管理ページの「エントリーの書き出し」を使って全エントリーを1つのテキストファイルに出力し、Perlスクリプトをガリガリ書いてBlosxom用のファイルに変換した。

●Hi-Ho依存の設定
参考までに、Hi-HoのレンタルWebサーバーでBlosxomを動かすために設定した、blosxom.cgi内の設定項目の一部を挙げておく。(一部伏字)
$dataurlは、テンプレートHTML(フレーバー)内からの他ファイル(スタイルシートや画像)参照で必要になる絶対パス。フレーバー内の相対パスは特殊なCGIディレクトリからの相対パスになり、データディレクトリに辿り着けない。

# Where are this blog's entries kept?
$basedir = "$ENV{'HOME'}/html/*************";
$datadir = "$basedir/entries";

# What's my preferred base URL for this blog (leave blank for automatic)?
$url = "http://www.sam.hi-ho.ne.jp/cgi-bin/user/ynomura/blosxom.cgi";
$dataurl = "http://www.sam.hi-ho.ne.jp/~ynomura/blosxom";

# Where are my plugins kept?
$plugin_dir = "$basedir/plugins";


See more ...

WriteBacks

Oct 25, 2008

Blosxom環境構築(1/2)

先月、このweblogが動く自宅サーバーへのアクセス手段でお世話になってる、無料のDDNSサービス、ieServer.Netが不調で、しばしば外からここにアクセスできなくなっていた。
当初から大変ありがたく使わせてもらっているが、無料のサービスなので十分な設備で運営されるとは限らないし、いつサービス終了となるかもわからない。大してアクセス数がある訳でも無いが、せっかくネチネチと綴ってきた現役のweblogなので、まだ消滅させるには惜しい。そこで、予備のweblog環境を検討することにした。

このサーバーのblogは、MovableTypeを使って記事を静的なHTMLに変換して作っており、一度HTMLに変換すれば、WebサーバーにMovableTypeが無くても閲覧可能である(但しコメント/トラックバックの受信やコンテンツ内検索等はできない)。従って、MovableTypeを使って作成した静的なHTMLをプロバイダーが提供するWebサーバーにアップロードするだけで、公開はできる。
…と考えて作業を始めると、1つのMovableType環境で複数のURL向けに静的なHTMLを出力するのは容易でないことがわかった。MovableTypeの設定ファイルを2つ用意して手で切り替えるか、データベース内に全データを丸ごと2つ持たないといけないのだ。あまりにもすっきりしないし、不便である。

そこで、プロバイダーのWebサーバー側にブログシステムを構築することを考えた。
世の中にはMovableType以外にNucleus、tDiary、華式、blosxom等様々なブログツールや日記ツールがあるが、私が使用しているプロバイダーであるHi-HoのレンタルWebサーバーは持ち込みのプログラムを動かす制約が厳しく、ほとんど素のPerlのCGIしか動かせない。従って、多数のPerlモジュールを必要とするMovableTypeも、PHPで動く華式も、SQLデータベースを必要とするNucleusも動かない。
そんな中、blosxomは動いた。

実はこのblogを始める以前に、Kyo Nagashimaさんのblosxom starter kit(以下bsk)を使ってHi-Ho上でblosxomを動かして、少しだけblogしていた。そのbskをインストールした4年前、Hi-Hoのサーバー上で動かせたフリーウェアのブログ/日記ツールはblosxomだけだったが、現在やはり動かせるのはblosxomだけのようである。

Blosxomの公式サイトのページは更新が止まっており、そのサイトにあるblosxom本体も5年前からバージョンアップされていない。ユーザーも減っているようで、ユーザー自身が「なぜ今更blosxomを使っているかというと…」と自嘲気味に言い訳しているのを見かける。ユーザーがいなければノウハウを公開する人もいない。もうblosxomを使うのは厳しいか…と思いかけた所で、"The Unofficial Blosxom User Group"というサイトを見つけた。どうやら今年に入っても有志によってバージョンアップされ、プラグインも追加されていってるようだ。
そして、日本でもInfoseekの無料Webレンタルサーバーでblosxomしている人がいた。InfoseekでもHi-Hoでも動くのなら、レンタルウェブサーバー間のポータビリティは十分期待できる。一旦システムを構築したら、無料のレンタルサービスででもその環境を入れれば全く同じブログを使い続けられるのだ。エコロジーだ。

ひと安心してblosxomをバージョンアップしようとすると、以前にインストールしたblosxomが動かなくなっていた。Perlスクリプトの途中でエラーになってたようで、サーバーのPerlの仕様が変わったのか?とげんなりしたが、仕方なく解析すると、spamを受けまくってwriteback(コメント、トラックバック)ファイルが巨大になってたのが原因だと判明した。インストールしていたのはbsk 1.02そのままの環境で、当時spamは全く無かったので、spam対策はしていなかった。しかも、パスワードを設定していたはずなのだが、有効にしていたwikiedish(ブラウザ上でエントリーを編集するためのプラグイン)を使って、全てのエントリーがspamによって書き換えられていた。
辛うじてbsk 1.02を初期状態に戻して動かすことはできたが、blosxom本体だけをバージョンアップする方法と、spam対策されたWritebackPlusを組み込む方法がわからなかったので、bskを参考にして1からblosxom環境を構築することにした。

See more ...

WriteBacks

Maximaでマンデルブロー図形

お盆なので、特にお盆とは関係ないが、ついでに定番のマンデルブロー集合も描いておく。

f(c,z):=c-z^2$
density(x,y):= block([c,z,d,i],
c:x+y*%i,
z:c,
d:0,
for i:0 while (abs(z)<2 and d<40) do (z:f(c,z),d:d+1),
d
)$
plot3d(density,[x,-0.5,1.5],[y,-1.5,1.5],[grid,100,100],
[gnuplot_pm3d,true],[gnuplot_preamble,"set view map;unset surface;set xrange [-0.5:1.5];set yrange [-1.5:1.5];"]);

See more ...

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

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:50 in PC一般 | WriteBacks (0)
WriteBacks

Maximaで微分方程式

せっかくMaxima使い始めたので、慣れるためにさらにMaximaのマニュアルを読み進めて遊んでると、Maximaは微分方程式も解けることがわかった。難しそうな響きを擁して実際難しい、物理学の基本ツールである。一瞬気が進まなかったが、大したものだと思ったので、敬意を表して試してみた。
すると、微分方程式が見えた気分になるようなグラフが出せることがわかり、面白かったので、ここに書き記すことにした。

まず、コマンドラインで微分させてみる。diff()が微分演算だ。

例によってTeXmacsを使って出力を整形している。(%i??)に続く青い部分が入力した部分、(%o??)に続く黒い部分がMaximaの出力だ。diff()の2つ目の引数は何で微分するか、3つ目の引数(省略可)は何回微分するかを意味する。

次に、desolve()を使ってシンプルな微分方程式\frac{d^2}{dx^2}y+\frac{d}{dx}y+y=0を解かせてみる。

desolve()の2つ目の引数は、どの関数について解くかを示す。diffの前の'(アポストロフィー)は、そのLisp関数を評価(実行&展開)せず、そのまま渡すことを指示するものである(desolve()の場合はこれを省略して評価してしまっても良いようだが、後述のode2()では不可)。y(x)の(x)はyがxの関数であることを示し、これを省略するとdy/dxが0になってしまう。
atvalue()を使って、y,y'の初期値を与えてから再度解かせてみる。

入力行の行末に;の代わりに$を付けると、計算結果を出力しない。
念のため、元の方程式に当てはまるかどうか、検算してみる。

%は直前の計算結果を表す。%o7は見るからに0なのだが、なぜか計算してくれなかったので、ratsimp()で整理させた。

常微分方程式(未知関数が1つだけの微分方程式)はode2()を使って解く方法もある。ode2()を使う場合、初期値は一般解の計算結果に対してic1(),ic2(),bc2()などを使って与える。同じ方程式を解かせてみる。

ode2()の2つ目の引数は常微分方程式の未知関数、3つ目の引数はその関数が従属する変数だ。上の場合は、引数からyがxの関数であることがわかるので、yをy(x)と書かなくても良い。

desolve()とode2()の能力を比較するため、解くには少し難しそうな微分方程式xy+\frac{d}{dx}y=0を解かせてみる。

ode2()は難なく解いた。念のため検算しておく。

それに対して、同じものをdesolve()に解かせると、下のような奇妙な式になった。

おそらく、ラプラス変換してx,yについて解いてラプラス逆変換したもの、という意味なのだろうが、これでは使えないだろう。これでも使える解なのかどうかを一応確認してみる。

やはり解けていない、使えない解であることがわかる。常微分方程式を解くのはdesolve()よりode2()の方が優れているようだ。

しかし、desolve()は複数の未知関数を含む連立微分方程式も解ける。
\left\{ \begin{array}{lll} \frac{d}{dx}f(x) & = & \frac{d}{dx}g(x)+\sin{x} \\ \frac{d^2}{dx^2}g(x) & = & \frac{d}{dx}f(x)-\cos{x} \end{array} \right.

e1とe2が方程式、f(x)とg(x)が未知関数だ。desolve()には、上のように、方程式と未知関数をそれぞれリスト形式で与える。

次に、plotdf()を使って微分方程式のグラフを出してみる。plotdf()に与える微分方程式は、dy/dx=F(x,y)または[dx/dt, dy/dt]=[F(x,y), G(x,y)]に限定されている。

load("plotdf");
plotdf(x);
これを実行すると、下のような画面が開く。(FreeBSD 6.xではTcl/Tkをインストールしていても「wishが無い」というエラーが出たが、wishという名のwish?.?へのリンクを作ればOKだった)

たくさんある矢印は、その地点(x,y)におけるdy/dxの向きと大きさを表している。このグラフ上のどこかをクリックすると、その座標を初期値とするグラフが描かれる。下の図は(0,0)辺りをクリックしたもの。

さらに色々な所をクリックする。

次の2つの図は、

plotdf(sin(x));
を実行したものである。

xとyを含む式も試してみる。

plotdf(sin(x)+cos(y));

放物線を横にした、x=y^2を出そうとした。両辺をxで微分すると1=2y(dy/dx)なので、dy/dx=1/(2y)である。よりシンプルに、dy/dx=1/yで試そうとして、

plotdf(1/y);
とすると、原因がよくわからないが、Tclのエラーになった。しかし、次のように小さなオフセットを与えると、何か描画された。2つ目以降の引数に色々なオプションを指定しているが、これらが無くてもエラーにはならない。
plotdf(10^-10+1/y,[trajectory_at,0.5,1],[nsteps,200],[xradius,3],[yradius,6],[xcenter,3],[ycenter,0]);

trajectory_atで初期値を(0.5,1)として描画したものだが、左側が何かおかしい。さらに何ヶ所かクリックしてみると、

y=0の近辺でおかしくなるようだ。1/(y+⊿y)-1/yの計算でもしているのだろうか。

x=cos(t), y=sin(t)のグラフは円になるので、[dx/dt, dy/dt]=[-y,x] (=[-sin(t), cos(t)])のグラフも円になる(図は省略)。それよりもっと場を複雑にさせようと考えて、dx/dtとdy/dtを周期関数にしてみた。次のグラフは[dx/dt, dy/dt]=[sin(y),sin(x)]である。

plotdf([sin(y),sin(x)]);

色を変えながら初期値をクリックした。セル状のグラフが無限に続く。
dx/dtとdy/dtを入れ替えると、次のようなグラフになる。
plotdf([sin(x),sin(y)]);

こんなのも作ってみた。

plotdf([
10^-10 - sin(atan(y/x))*cos(%pi/4*sin(8*atan(y/x))) - cos(atan(y/x))*sin(%pi/4*sin(8*atan(y/x))),
10^-10 - sin(atan(y/x))*sin(%pi/4*sin(8*atan(y/x))) + cos(atan(y/x))*cos(%pi/4*sin(8*atan(y/x)))
]);

[dx/dt, dy/dt]=[F(x,y), G(x,y)]で表せない(dx/dtにtが入るとか)方程式も、desolve()で解けるものなら
、いつものplot2d()でグラフにできる。
\left\{\begin{array}{l}\frac{dx}{dt}=3x-7y\\ \frac{dy}{dt}=5x-4y \end{array}\right.
(この方程式は[dx/dt, dy/dt]=[F(x,y), G(x,y)]で表せるが、plotdf()ではおかしなグラフになる)

atvalue(x(t),t=0,1);
atvalue(y(t),t=0,1);
Z:desolve([diff(x(t),t)=3*x(t)-7*y(t), diff(y(t),t)=5*x(t)-4*y(t)],[x(t),y(t)]);
plot2d([parametric, x(t), y(t), [t,0,6]],[nticks,200]),Z;

アステロイド曲線\left\{\begin{array}{l} x''(t)=2y'(t)-3x(t)\\ y''(t)=-2x'(t)-3y(t) \end{array}\right.(実はx=cos(t)^3, y=sin(t)^3から逆算したもの)を描いてみる。

atvalue(x(t),t=0,1);
atvalue(y(t),t=0,0);
atvalue(diff(x(t),t),t=0,0);
atvalue(diff(y(t),t),t=0,0);
Z: desolve([
diff(x(t),t,2)=2*diff(y(t),t)-3*x(t),
diff(y(t),t,2)=-2*diff(x(t),t)-3*y(t)
],[x(t),y(t)]);
plot2d([parametric,x(t),y(t), [t,0,2*%pi]],[nticks,200]),Z;

[parametric, ...]の部分を、それらのリストにすることにより、複数のグラフを1枚に描ける。

eq: [
diff(x(t),t,2)=2*diff(y(t),t)-3*x(t),
diff(y(t),t,2)=-2*diff(x(t),t)-3*y(t)
];
atvalue(x(t), t=0, 1);
atvalue(y(t), t=0, 0);
atvalue(diff(x(t), t), t=0, 0);
atvalue(diff(y(t), t), t=0, 0);
Z1: desolve(eq, [x(t), y(t)]);
atvalue(x(t), t=0,2);
Z2: desolve(eq, [x(t), y(t)]);
atvalue(x(t), t=0,3);
Z3: desolve(eq, [x(t), y(t)]);
L1: [parametric, x(t), y(t), [t, 0, 2*%pi]], Z1;
L2: [parametric, x(t), y(t), [t, 0, 2*%pi]], Z2;
L3: [parametric, x(t), y(t), [t, 0, 2*%pi]], Z3;
plot2d([L1,L2,L3], [nticks,200]);

See more ...

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

円盤の衝突計算

弾性衝突円盤アプレットの衝突計算について、悲しいことに式を求めるのにかなり時間を費やしたので、忘れないように記録する。これが無いと、自分でもソースコードを理解できなくなってしまう。

まず直線上の衝突について考える。物体1,2の質量をm1,m2、衝突前の速度をv1,v2、衝突後の速度をv1',v2'とすると、完全弾性衝突の場合、運動量保存の法則より
m_1v_1+m_2v_2=m_1v'_1+m_2v'_2
-\frac{v'_1-v'_2}{v_1-v_2}=1
が成り立つ。速度変化に注目してこれを整理すると、
v'_1-v_1=\frac{-2m_2(v_1-v_2)}{m_1+m_2}
v'_2-v_2=\frac{-2m_1(v_2-v_1)}{m_1+m_2}
が得られる。これは、
v_{tmp}=\frac{m_1v_1+m_2v_2}{m_1+m_2}---(1)
と置くと、
v'_1-v_1=2(v_{tmp}-v_1)---(2)
v'_2-v_2=2(v_{tmp}-v_2)---(3)
と整理できる。

次に平面上の衝突について考える。
ball_hit.png
図のように円盤と円盤が衝突する場合の速度変化については、円盤の中心同士を結ぶ線上の方向の
速度が衝突によって変化し、衝突方向でない速度成分は変化しないので、衝突方向の速度V1r、V2rについて、上記の式(1)~(3)により速度の変化を求めれば良い。つまり、 Pを円盤の中心間を結ぶベクトル、V1,V2を円盤1,2の速度ベクトル、v1r,v2rを円盤1,2のP方向の速度、V1,v1rの衝突後の速度をそれぞれV1',v1r'とすると、円盤の速度変化は
\vec{V'_1}-\vec{V_1} = \frac{\vec{P}}{|\vec{P}|}(v'_{1r}-v_{1r})
である。

v1rの計算については、θをPとV1の間の角度とすると、
v_{1r} = |\vec{V_1}|\cos{\theta}
であり、PとV1との内積は
\vec{V_1}\cdot\vec{P} = |\vec{V_1}||\vec{P}|\cos{\theta}
だから、
v_{1r} = \frac{\vec{V_1}\cdot\vec{P}}{|\vec{P}|}
で求められる。同様に
v_{2r} = \frac{\vec{V_2}\cdot\vec{P}}{|\vec{P}|}
となる。

直線の場合の式(1)と同様に
v_{tmp}=\frac{m_1v_{1r}+m_2v_{2r}}{m_1+m_2} \left (=\frac{1}{|\vec{P}|}\frac{m_1<br>(\vec{V_1}\cdot\vec{P})+m_2(\vec{V_2}\cdot\vec{P})}{m_1+m_2}\right )
と置くと、
v'_{1r}-v_{1r}=2\frac{\vec{P}}{|\vec{P}|}(v_{tmp}-v_{1r})
となる。vtmpとv1rの両方に1/|P|が入っているので、計算上の便宜のため、
v_{1r} = \vec{V_1}\cdot\vec{P}
v_{2r} = \vec{V_2}\cdot\vec{P}
と置き直すと、
v'_{1r}-v_{1r}=2\frac{\vec{P}}{|\vec{P}|^2}(v_{tmp}-v_{1r})
v'_{2r}-v_{2r}=2\frac{\vec{P}}{|\vec{P}|^2}(v_{tmp}-v_{2r})
が得られる。2次元のベクトルをX軸成分とY軸成分に分けて
\vec{P}=(P_x,P_y), \vec{V_n}=(V_{nx},V_{ny})
とスカラー表記にすると、
|\vec{P}|^2 = P_x^2+P_y^2, \vec{V_n}\cdot\vec{P}=V_{nx}P_x+V_{ny}P_y
なので、円盤同士の衝突による速度変化は、
\vec{V'_1}-\vec{V_1} = \left (\frac{2P_x}{P_x^2+P_y^2}(v_{tmp}-v_{1r}),\frac{2P_y}{P_x^2+P_y^2}(v_{tmp}-v_{1r}) \right )
\vec{V'_2}-\vec{V_2} = \left (\frac{2P_x}{P_x^2+P_y^2}(v_{tmp}-v_{2r}),\frac{2P_y}{P_x^2+P_y^2}(v_{tmp}-v_{2r}) \right )
v_{tmp}=\frac{m_1v_{1r}+m_2v_{2r}}{m_1+m_2}
v_{1r} = V_{1x}P_x+V_{1y}P_y
v_{2r} = V_{2x}P_x+V_{2y}P_y
で求められることがわかる。

See more ...

Posted at 22:49 in 雑記 | WriteBacks (0)
WriteBacks

Sep 25, 2008

3Dカッチンアプレット

作ろうとしていたものが何となくそれらしく動き出したので、一旦ここに貼っておく。
3Dカッチンアプレットのページへのリンク
ソースコード

一応、衝突したら弾き合うのだが、わかりづらい。どうしたらわかり易くなるのだろう。
衝突したら赤く光るチェックボックスを付けてみたが、やはりわかりづらい。
3Dだと仕方ないのだろうか。

See more ...

Posted at 22:06 in Java | WriteBacks (0)
WriteBacks

Sep 14, 2008

Java3Dアプレット試作

Java3Dを触り始めたので、何か作ってみようと思って作り始めたのだが、想定したものが意外に難しいことがわかってきて、挫折気味である。しかし、その途中段階でできたテストアプリが、自己満足に浸れるくらいには興味深く動いたので、一旦区切りをつけて、公開することにする。

何らかのアプレットのページへのリンク
・ソースコード
 J3DTest2.java
 TimerBehavior.java
 TimerBehaviorProcess.java

Java3D(OpenGL)のイベントのスケジューラーはAWTやSwingのそれとは別であり、例えばAWTやSwingのパーツに対してrepaint()してもCanvas3Dは再描画されない。Java3Dのイベント処理モデルに合わせるには、Behavior(javax.media.j3d.Behavior)というクラスを継承するクラスでイベントを処理する必要がある。
Java3Dでアニメーションするには、その理由により、java.util.Timerやjavax.swing.Timerを使ってAWTやSwingで処理するだけではうまくいかず、Canvas3D、Bounds(BoundingBox,BoundingSphere etc.)に関連付けたBehaviorに処理を渡さないといけない。Java3Dのアニメーションをさせるだけなら、java.util.Timerやjavax.swing.Timerを使うよりBehaviorを使う方が賢明だと思われる。
とりあえず、javax.media.j3d.WakeupOnElapsedTimeのを使って周期的にBehaviorを動かせるようなので、今回のアプリでは、それを使ったTimerBehaviorというのを定義している。
一定の速度の移動や一定周期の回転や一定の拡大率の拡大など、変化量が一定なアニメーションなら、なんとかInterpolatorというBehaviorを使えばシンプルである。今回できたアプレットもInterpolatorを使うと美しく書けるのかも知れないが、拡張性を求めて敬遠した。

See more ...

Posted at 00:59 in Java | WriteBacks (0)
WriteBacks

Sep 06, 2008

Java3Dのアプレットを公開する方法

Java3Dがインストールされていない環境でもJavaアプレットでJava3Dを動かす方法があることを知ったので、試しに簡単なアプリを公開する。

Java3Dアプレットのページへのリンク
ソースコード

コードは、Java3Dのデモパッケージに入っているHelloUniverse.javaに少し味をつけただけのものである。この3次元世界にあり得ない動きなので、ちょっと気持ち悪いかもしれない。

Canvas3dのインスタンスを貼り付けたAppletのクラスを、普通にHTMLの<applet>タグで公開すると、Java3Dのランタイムがインストールされていない環境では実行できない。しかし、org.jdesktop.applet.util.JNLPAppletLauncherを使って起動するように公開すると、(メジャーな環境では)Java3Dのランタイムがインストールされていなくても実行可能になる。

より具体的な方法は、Java3Dがインストールされた環境で動くJavaアプレット(Canvas3dのインスタンスを貼り付けたAppletのクラス)をJAR形式で用意し、AppletLauncherのページにある通りに、下記のようなHTMLで起動を指示する。少々長いが、最低限の設定なら、触る所は
・width/heightの部分
・アプレットの指定部分(下記で言うと、2ヶ所あるHello4dUniverseの部分)
・subapplet.displaynameの部分
くらいで良いと思う。

<applet code="org.jdesktop.applet.util.JNLPAppletLauncher2"
  width=400 height=400
  archive="Hello4dUniverse.jar,
               http://download.java.net/media/applet-launcher/applet-launcher.jar,
               http://download.java.net/media/java3d/webstart/release/j3d/latest/j3dcore.jar,
               http://download.java.net/media/java3d/webstart/release/j3d/latest/j3dutils.jar,
               http://download.java.net/media/jogl/builds/archive/jsr-231-webstart-current/jogl.jar,
               http://download.java.net/media/gluegen/webstart/gluegen-rt.jar,
               http://download.java.net/media/java3d/webstart/release/vecmath/latest/vecmath.jar">
  <param name="codebase_lookup" value="false">
  <param name="subapplet.classname" value="Hello4dUniverse">
  <param name="subapplet.displayname" value="My Java 3D Applet">
  <param name="jnlpNumExtensions" value="1">
  <param name="jnlpExtension1" value="http://download.java.net/media/java3d/webstart/release/java3d-latest.jnlp">
  <param name="progressbar" value="true">
  <param name="noddraw.check" value="true">
</applet>

See more ...

Posted at 19:08 in Java | WriteBacks (1)
WriteBacks

applet+code-でつながるブログリング

applet+codeに関するブログをまとめています。[link]

Posted by at 01/07/2009 10:30:38 PM

Aug 17, 2008

ペンギン in the world

Languagenotation
Japaneseペンギン
Englishpenguin
GermanPinguin
Frenchpingouin
Spanishping?ino
Italianpinguino
Simplified Chinese企鹅
Traditional Chinese企鵝
Russianпингвин
Polishpingwin
CzechTučňák
HungarianPingvin
Dutchpinguïn
Thaiเพนกวิน
Hangul펭귄
Turkishpenguen
Portuguesepinguim
Arabicالبطريق
Persianپنگوئن
Finnishpingviini
Danishpingvin
Swedishpingvin

See more ...

Posted at 13:07 in 雑記 | WriteBacks (0)
WriteBacks

Aug 15, 2008

Maximaでマンデルブロー図形

お盆なので、特にお盆とは関係ないが、ついでに定番のマンデルブロー集合も描いておく。

f(c,z):=c-z^2$
density(x,y):= block([c,z,d,i],
c:x+y*%i,
z:c,
d:0,
for i:0 while (abs(z)<2 and d<40) do (z:f(c,z),d:d+1),
d
)$
plot3d(density,[x,-0.5,1.5],[y,-1.5,1.5],[grid,100,100],
[gnuplot_pm3d,true],[gnuplot_preamble,"set view map;unset surface;set xrange [-0.5:1.5];set yrange [-1.5:1.5];"]);

See more ...

Posted at 19:25 in 数学 | WriteBacks (0)
WriteBacks

Aug 13, 2008

Maximaで微分方程式

せっかくMaxima使い始めたので、慣れるためにさらにMaximaのマニュアルを読み進めて遊んでると、Maximaは微分方程式も解けることがわかった。難しそうな響きを擁して実際難しい、物理学の基本ツールである。一瞬気が進まなかったが、大したものだと思ったので、敬意を表して試してみた。
すると、微分方程式が見えた気分になるようなグラフが出せることがわかり、面白かったので、ここに書き記すことにした。

まず、コマンドラインで微分させてみる。diff()が微分演算だ。

例によってTeXmacsを使って出力を整形している。(%i??)に続く青い部分が入力した部分、(%o??)に続く黒い部分がMaximaの出力だ。diff()の2つ目の引数は何で微分するか、3つ目の引数(省略可)は何回微分するかを意味する。

次に、desolve()を使ってシンプルな微分方程式\frac{d^2}{dx^2}y+\frac{d}{dx}y+y=0を解かせてみる。

desolve()の2つ目の引数は、どの関数について解くかを示す。diffの前の'(アポストロフィー)は、そのLisp関数を評価(実行&展開)せず、そのまま渡すことを指示するものである(desolve()の場合はこれを省略して評価してしまっても良いようだが、後述のode2()では不可)。y(x)の(x)はyがxの関数であることを示し、これを省略するとdy/dxが0になってしまう。
atvalue()を使って、y,y'の初期値を与えてから再度解かせてみる。

入力行の行末に;の代わりに$を付けると、計算結果を出力しない。
念のため、元の方程式に当てはまるかどうか、検算してみる。

%は直前の計算結果を表す。%o7は見るからに0なのだが、なぜか計算してくれなかったので、ratsimp()で整理させた。

常微分方程式(未知関数が1つだけの微分方程式)はode2()を使って解く方法もある。ode2()を使う場合、初期値は一般解の計算結果に対してic1(),ic2(),bc2()などを使って与える。同じ方程式を解かせてみる。

ode2()の2つ目の引数は常微分方程式の未知関数、3つ目の引数はその関数が従属する変数だ。上の場合は、引数からyがxの関数であることがわかるので、yをy(x)と書かなくても良い。

desolve()とode2()の能力を比較するため、解くには少し難しそうな微分方程式xy+\frac{d}{dx}y=0を解かせてみる。

ode2()は難なく解いた。念のため検算しておく。

それに対して、同じものをdesolve()に解かせると、下のような奇妙な式になった。

おそらく、ラプラス変換してx,yについて解いてラプラス逆変換したもの、という意味なのだろうが、これでは使えないだろう。これでも使える解なのかどうかを一応確認してみる。

やはり解けていない、使えない解であることがわかる。常微分方程式を解くのはdesolve()よりode2()の方が優れているようだ。

しかし、desolve()は複数の未知関数を含む連立微分方程式も解ける。
\left\{ \begin{array}{lll} \frac{d}{dx}f(x) & = & \frac{d}{dx}g(x)+\sin{x} \\ \frac{d^2}{dx^2}g(x) & = & \frac{d}{dx}f(x)-\cos{x} \end{array} \right.

e1とe2が方程式、f(x)とg(x)が未知関数だ。desolve()には、上のように、方程式と未知関数をそれぞれリスト形式で与える。

次に、plotdf()を使って微分方程式のグラフを出してみる。plotdf()に与える微分方程式は、dy/dx=F(x,y)または[dx/dt, dy/dt]=[F(x,y), G(x,y)]に限定されている。

load("plotdf");
plotdf(x);
これを実行すると、下のような画面が開く。(FreeBSD 6.xではTcl/Tkをインストールしていても「wishが無い」というエラーが出たが、wishという名のwish?.?へのリンクを作ればOKだった)

たくさんある矢印は、その地点(x,y)におけるdy/dxの向きと大きさを表している。このグラフ上のどこかをクリックすると、その座標を初期値とするグラフが描かれる。下の図は(0,0)辺りをクリックしたもの。

さらに色々な所をクリックする。

次の2つの図は、

plotdf(sin(x));
を実行したものである。

xとyを含む式も試してみる。

plotdf(sin(x)+cos(y));

放物線を横にした、x=y2を出そうとした。両辺をxで微分すると1=2y(dy/dx)なので、dy/dx=1/(2y)である。よりシンプルに、dy/dx=1/yで試そうとして、

plotdf(1/y);
とすると、原因がよくわからないが、Tclのエラーになった。しかし、次のように小さなオフセットを与えると、何か描画された。2つ目以降の引数に色々なオプションを指定しているが、これらが無くてもエラーにはならない。
plotdf(10^-10+1/y,[trajectory_at,0.5,1],[nsteps,200],[xradius,3],[yradius,6],[xcenter,3],[ycenter,0]);

trajectory_atで初期値を(0.5,1)として描画したものだが、左側が何かおかしい。さらに何ヶ所かクリックしてみると、

y=0の近辺でおかしくなるようだ。1/(y+⊿y)-1/yの計算でもしているのだろうか。

x=cos(t), y=sin(t)のグラフは円になるので、[dx/dt, dy/dt]=[-y,x] (=[-sin(t), cos(t)])のグラフも円になる(図は省略)。それよりもっと場を複雑にさせようと考えて、dx/dtとdy/dtを周期関数にしてみた。次のグラフは[dx/dt, dy/dt]=[sin(y),sin(x)]である。

plotdf([sin(y),sin(x)]);

色を変えながら初期値をクリックした。セル状のグラフが無限に続く。
dx/dtとdy/dtを入れ替えると、次のようなグラフになる。
plotdf([sin(x),sin(y)]);

こんなのも作ってみた。

plotdf([
10^-10 - sin(atan(y/x))*cos(%pi/4*sin(8*atan(y/x))) - cos(atan(y/x))*sin(%pi/4*sin(8*atan(y/x))),
10^-10 - sin(atan(y/x))*sin(%pi/4*sin(8*atan(y/x))) + cos(atan(y/x))*cos(%pi/4*sin(8*atan(y/x)))
]);

[dx/dt, dy/dt]=[F(x,y), G(x,y)]で表せない(dx/dtにtが入るとか)方程式も、desolve()で解けるものなら
、いつものplot2d()でグラフにできる。
\left\{\begin{array}{l}\frac{dx}{dt}=3x-7y\\ \frac{dy}{dt}=5x-4y \end{array}\right.
(この方程式は[dx/dt, dy/dt]=[F(x,y), G(x,y)]で表せるが、plotdf()ではおかしなグラフになる)

atvalue(x(t),t=0,1);
atvalue(y(t),t=0,1);
Z:desolve([diff(x(t),t)=3*x(t)-7*y(t), diff(y(t),t)=5*x(t)-4*y(t)],[x(t),y(t)]);
plot2d([parametric, x(t), y(t), [t,0,6]],[nticks,200]),Z;

アステロイド曲線\left\{\begin{array}{l} x''(t)=2y'(t)-3x(t)\\ y''(t)=-2x'(t)-3y(t) \end{array}\right.(実はx=cos(t)^3, y=sin(t)^3から逆算したもの)を描いてみる。

atvalue(x(t),t=0,1);
atvalue(y(t),t=0,0);
atvalue(diff(x(t),t),t=0,0);
atvalue(diff(y(t),t),t=0,0);
Z: desolve([
diff(x(t),t,2)=2*diff(y(t),t)-3*x(t),
diff(y(t),t,2)=-2*diff(x(t),t)-3*y(t)
],[x(t),y(t)]);
plot2d([parametric,x(t),y(t), [t,0,2*%pi]],[nticks,200]),Z;

[parametric, ...]の部分を、それらのリストにすることにより、複数のグラフを1枚に描ける。

eq: [
diff(x(t),t,2)=2*diff(y(t),t)-3*x(t),
diff(y(t),t,2)=-2*diff(x(t),t)-3*y(t)
];
atvalue(x(t), t=0, 1);
atvalue(y(t), t=0, 0);
atvalue(diff(x(t), t), t=0, 0);
atvalue(diff(y(t), t), t=0, 0);
Z1: desolve(eq, [x(t), y(t)]);
atvalue(x(t), t=0,2);
Z2: desolve(eq, [x(t), y(t)]);
atvalue(x(t), t=0,3);
Z3: desolve(eq, [x(t), y(t)]);
L1: [parametric, x(t), y(t), [t, 0, 2*%pi]], Z1;
L2: [parametric, x(t), y(t), [t, 0, 2*%pi]], Z2;
L3: [parametric, x(t), y(t), [t, 0, 2*%pi]], Z3;
plot2d([L1,L2,L3], [nticks,200]);

See more ...

Posted at 20:17 in 数学 | WriteBacks (0)
WriteBacks

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

Jun 22, 2008

ビット数を数えるルーチン

32bitの立っているビット数を高速に数えるアルゴリズムとして、

int bitcount(unsigned long n)
{
    n = ((n&0xaaaaaaaa) >> 1) + (n&0x55555555); 
    n = ((n&0xcccccccc) >> 2) + (n&0x33333333); 
    n = ((n&0xf0f0f0f0) >> 4) + (n&0x0f0f0f0f); 
    n = ((n&0xff00ff00) >> 8) + (n&0x00ff00ff); 
    n = ((n&0xffff0000) >> 16) + (n&0x0000ffff); 
    return n;
}
こういうのが良く知られている。1行目で2ビットずつ足し合わせることで2ビットの中の立っているビット数を求め、2行目でそれらを2つずつ合計して4ビットの中の立っているビット数を求め…というのを32bit使って並列にやっている訳だ。ポイントを理解すれば単純明快、1ビットずつ32回評価するより段違いに速いのは明らかである。

数年前にこれを初めて見た時、独学の計算量節約プログラミングを長年続け、私の頭の中で完成しつつあった独自のプログラム理論は、一気に消し飛んだ。

そんな記憶も彼方の昨日、ネットサーフィンしていて、下のようなコードを見つけた。HAKMEM 169と呼ばれるアルゴリズムらしい。ネット上の色々なC言語実装を参考に勝手に整形した。

int bitcount2(unsigned long n){
    n = n - ((n >> 1) & 033333333333) - ((n >> 2) & 011111111111);
    return ((n + (n >> 3)) & 030707070707) % 077;
}
短い!
1つ目の例と同じ香りが漂うが、演算子の数が少ない。もしかして最強コードか?と思って解読しようとすると、数学的に大変ややこしい。
1行目は、3bitの整数m内の立っているビット数がm - (m>>1) - (m>>2)で計算できる(※1)ことを利用して、3bitずつ並列に、立っているビット数を計算する。2行目は、3bitずつのビット数を2つずつ足し合わせて6bitずつのビット数にしてから(※2)、bが6bitの整数、kが6の倍数の時にb*2^kを6bitの111111で割った余りがbになる(※3)ことを利用して、6bitずつの合計を求めている。

よく読むと、剰余演算を使ってるので、1つ目の例の方が速い。
参考:http://infolab.stanford.edu/~manku/bitcount/bitcount.html
("Parallel"が1つ目の例、"MIT"が2つ目の例に相当する。"Parallel"より速いのはビット数のテーブルを予め用意している。)

元々のMITのHAKMEMのコードは剰余演算のあるアセンブラで書かれており、アセンブラとしては目を疑うほど短い。HAKMEM 169は最速ではないが、コードのコンパクトさと合わせて称えるべきアルゴリズムなのだろう。

See more ...

Posted at 20:11 in PC一般 | WriteBacks (1)
WriteBacks

アセンブラとは

英語表記は、assembler。アセンブリ言語を使って作られたソースを、コンピュータが実行できる形に変換するためのソフトのこと。アセンブリ言語とは、コンピ...[link]

Posted by at 07/07/2008 06:26:22 AM

Jun 21, 2008

C言語の悪習

ある所で、C言語のソースに

if (NULL == pointer){...
という表記があるのを見て、好ましくない書き方だ、と言ったら、一緒にいた2人に、この書き方には理由がある、この書き方もありだ、というようなことを言われた。
このC言語特有の"=="の左側に定数を置く書き方は、知ってると通っぽい、知ってることを自慢したくなる、不思議な匂いを漂わせるようだ。プロはみんなこう書いている、という錯覚を起こさせることもあるようだが、事実は全く逆で、まともな書籍には書かれない、最もわかってる人々は使わない書き方だ。著名なオープンソースではマイナーである。
しかし、これが好ましくない書き方であることを簡潔に説明することは難しく、また100%誤りと言えるものではなく、一応利点もあることから、根強い人気があり、まるで宗教論争のように平行線の議論が繰り返される。

C言語のスタイルでよく議論になってるのを見かけるトピックは他にもある。それらと合わせて、個人的な見解を書いてみる。

(1) 定数==変数
"NULL == pointer"に限らず、

0 == integer
とか
true == boolean
というのも見かける。
肯定派の根拠はもちろん、変数==定数という書き方だと、if文などの条件式で間違って(変数=定数)と書いた場合にコンパイルエラーにならないため、バグの元になりやすいというものだ。
しかし、それが唯一の利点であることが、逆に好ましくない理由を示している。もしその利点が無ければ、標記のような書き方は誰もしないだろう。その理由は、読みにくいからである。
その利点には、可読性を犠牲にするまでの価値があるだろうか。

(計算式=定数)だと定数を左にしなくてもコンパイルエラーになるし、(変数=変数)だとコンパイルエラーにする類似の手段はない。まさかそういう間違いをコンパイルエラーにするために普段から(0+変数==変数)と書く人は居るまい。救われるのは"=="と"="の書き間違い全てではなく、(変数=定数)のみである。
それも、現在広く使われているコンパイラでは、条件式の部分が(変数=定数)になっているとワーニングを出してくれるので、さほど問題ではない。それに、そもそもそんな書き間違いは稀であろう。

変数と定数の比較の時だけ比較対象を左に書き、それ以外の比較では比較対象を右に書く、というのは、コーディングスタイルとして整合性に欠けるのではないか。それとも、定数==変数と書く人は、必ず比較対象を左に書くのだろうか?
C言語以外の、比較対象の定数を比較演算子の右に書ける言語では、定数を左に書く人は居ないだろう。C言語の限られた比較演算の場合だけ、役立つかどうか不明な些細な利点のために、敢えてわかりにくい書き方をするべきだろうか。

(2) (char)NULL

char c = (char)NULL;
というやつである。もちろん正しくは
char c = '\0';
だ。
これの支持者は'\0'より(char)NULLの方が「ヌルキャラクタ」として理解しやすいし、(メジャーな誤用のため)よく見かけるし、(char)NULLが'\0'と異なる値になる処理系はまず存在せず誤動作しないので、間違いだと気付かないのだろう。
しかし、ANSIの言語仕様でNULLは「ヌルポインタ」であり、いかなるオブジェクトも指さないポインタの値と決まっている。C言語的には(int)NULLは0でなくても良いのだから、(char)NULLを'\0'の意味で書くのは誤りである。
それでも、実際に問題が起こらないなら読みやすい方がいいだろう、と言われれば、返す言葉が見つからない。

(3) プログラム終了時でもmalloc()したメモリは必ずfree()すべし
コーディングスタイルやコーディングルールとして、malloc()には必ず対になるfree()を書く、というのがある。習慣的に静的にメモリリークを防ぐためには、基本的には好ましいことであるが、main()を抜ける直前やexit()の直前のfree()は無駄である。争点は、そのようなfree()も書くべきかどうかだ。
free()に限らず、ソケットのclose()やファイルポインタのfclose()など、exit時にOSに返されると規定されているリソース全てについて同じ話が当てはまる。

終了時もfree()すべし派の言い分は、
・その方がプログラムの構造として美しい
・free()して害になることはほとんど無い
・exit()後に処理があると、free()を怠ることにより問題が生じる可能性がある
というのがあったのを記憶している。ソースコードの静的解析ツールを使うと警告が出てしまう、というのもあったかも知れない。3つ目はatexit()やon_exit()のことを言っているが、これは世界にそんなコードが1つあるか無いかというオーダーのレアケースであろう。

そりゃ本当に無害なら書いた方が良い、終了直前のfree()は省略できるというのは知らなくてもいい無駄な知識、だとは思うが、終了直前のfree()も書くとコンパイルされてマシン語コードになってしまうのである。例えばfree()1つくらいなら10~20バイトくらいかも知れないが、それも10個になると100~200バイトにもなる(PowerPCで試した1つの例だと120バイトになった)。それだけ実行ファイルのサイズが無駄に大きくなってしまうのである。

例えばソースコードの読みやすさやメンテナンス性を犠牲にしてまで実行ファイルのサイズが小さくなる書き方を選ぶべきかと問われれば、よほどサイズがクリティカルなプロジェクトで無い限りNoであろうが、終了直前のfree()を省略して害になることはほとんど無かろう。ゴミを作り出すコードを埋め込んでまで、ある種のプログラムの対称性を重視する価値があるだろうか。

自分が書くコード、目の前にあるコードに整っていることを求める気持ちはよくわかる。例えばもし、C言語のソースファイルの最後の閉じ括弧は書かなくても良くて、書かない方がオブジェクトコードが小さくなる、というコンパイラがあったとしたら、個人的には最後の閉じ括弧を書かないのは気持ち悪すぎてできない。書いてコメントアウトするのも同様である。読む用とコンパイル用に閉じ括弧ありのソースとなしのソースを別に作るとしても、コンパイル用のソースの閉じ括弧を消して保存する操作が違和感の極みであり、コードを見ながらの手作業ではかなり後味が悪いだろう。

それでも、私はやはり実利を取るべきであり、ゴミを生むコードは書くべきでないと思う。終了直前のfree()は、どうせまとめて焼却されるゴミは分別するコストが無駄、というアナロジーででも割り切って省くべきではないだろうか。

See more ...

Posted at 23:01 in PC一般 | WriteBacks (3)
WriteBacks

僕自身も、終了直前の free() は呼ばなくてもいいとは思うが、 ゴミを生むコードは書くべきではないというのは言い過ぎだと思う。 以前、組み込み系の OS でファイルディスクリプタを開放せずに、 終了するプログラムを書いたが OS がリソースを開放してくれなかった。 組み込みシステムで使われるような OS は例外的かもしれないが、 OS がプロセスの終了時にリソースの開放を確実にしてくれるのかは疑問。 POSIX あたりで、規定されているのだろうか?

Posted by eternalharvest at 02/02/2010 12:55:29 AM

もちろんexit()時に解放されるリソースはOSや使用するライブラリによって異なります。 私も組み込み系の仕事をしているので、コードはOSによって異なるのが当たり前だと思っています。上の記事もその認識で書いています。 ファイルディスクリプタが指すリソースがプロセス終了時に解放されないことがあるのは普通だと思います。(手元のBSD系OSでもそうでした)

Posted by ynomura at 02/02/2010 08:33:11 PM

http://jbbs.livedoor.jp/bbs/read.cgi/computer/5651/1048816258/150 > if(変数=定数)なんてやっちゃうカスは、ど~せ他んところでバグだらけだろうから所詮カス。 > if(定数==変数)は、私はカスなので↑やっちゃうかも、って宣言してるみたいで、所詮カス。 > そんなコード見たら不安で仕方ない。

Posted by at 02/24/2010 04:06:09 PM

Jun 14, 2008

弾性衝突円盤アプレットの設計メモ

先日公開した弾性衝突円盤アプレットはまた改造して別バージョンを作ろうと思ってるので、自分自身がそのソースコードを読めなくならないよう、設計を書き残す。

●Swing関係
・JAppletを継承するクラスは、盤面クラスのインスタンスを作ってタイマー処理するだけにしている。
 表示画面が複数あっても、盤面のインスタンスを増やすだけで済む。

・タイマーはSwing timerを使っている。汎用のjava.util.Timerを使うのに比べて、タイマーイベントがSwingのイベント発行スレッドにて他のイベントと一緒にシリアライズされるので、タイマーイベントによって処理が割り込まれるということがなく、排他制御がシンプルになる。
 というようなことが、JavaSEのAPI仕様書のjavax.swing.Timerの項に書いてある。
 なお、所々でSwingUtilities.invokeLaterを使ってるのは、同様の効果を期待してのことだが、よく考えると、Swing timerを使ってるので意味が無かった。

・JComponentの盤面の大きさは、JComponent#getPreferredSize()が返すようにしている。レイアウトマネージャはFlowLayoutを使っており、FlowLayoutは余計な調節をしないので、これだけで済むはず。

●描画関係
・円盤の画像は、盤面インスタンスの生成時に全パターン用意するようにしている。(50-58行目)

・盤面の描画は、スプライト的なものを使わず、毎回矩形で塗り潰して全描画、としてるので、せっかくなので、描画の度に色相を少しずつ変えながら塗り潰すことにした。(72-73行目)
 bg_countはそのためだけの変数。
 74行目の青色設定は、状態表示にしか効いていない。

・円盤の座標は中心の座標を保持しているため、描画開始位置は半径を引いた座標になる。

●衝突計算関係
・円盤の移動後に衝突判定と衝突後の速度の計算をしている(127-182行目)が、壁との衝突と他の円盤との衝突が同時に起こった場合、円盤同士の衝突を先に計算すると、壁との衝突の影響が考慮されないので、壁がクッションのような扱いになってしまう。かと言って、壁と他の円盤との衝突を同時に計算しようとすると、複雑な計算になる。
 そもそも、円盤と壁が円盤を挟むように同時に衝突する場合、真ん中の円盤は壁の役割を兼ねるので、円盤同士の衝突の計算は、壁に向かう円盤同士の衝突ではなく、壁に向かう円盤と壁から跳ね返ってきた円盤との衝突とするべきである。
 そのため、壁との衝突を先に処理するようにしている。

・壁との衝突判定は、円盤が壁と深く重なってる場合は衝突状態が続くので、速度ベクトルが壁の方を向いているかどうかを考慮するようにした。(131,138行目)
 円盤同士が深く重なってる場合も衝突状態が続くので、前回の中心間の距離を保存して、距離が縮まっているかどうかを衝突判定の材料にしている。(hit_depth[][]が前回の距離の2乗)

・複数の円盤の衝突が同時に発生した場合の計算は非常に複雑になる。
 1つの円盤が他の2つの円盤に同時に衝突するのはまだいいが、さらに衝突相手がまた別の円盤と同時に衝突する場合は、衝突の影響が全ての繋がる円盤に及ぶので、計算が大変である。
 そのため、1つの円盤は同時に他の1つの円盤としか衝突せず、さらに他の円盤との衝突はその瞬間は無視し、めり込むことを許す、とした。(isHit[]が衝突フラグ)

・衝突による速度変化を一旦別の配列に入れているのは、1つの円盤が同時に複数の円盤と衝突する計算をする場合のことを考えてのものであり、今回の方式では意味が無い。

●円盤の増減関係
・動的に円盤の数を減らせるよう、規定の数の円盤以外は、画面外に出るまで衝突判定の対象外とした。(185-193行目)

・動的に円盤の数を増やす時、既にある円盤と同じ場所から、速度ベクトルを30°ずらしてスタートさせるようにしている。それにより、円盤が分裂する感じになる。
 一気に2つ以上の円盤が増える場合は、分裂する円盤をそれぞれ違う円盤にしている。0番の円盤から分裂させると、いつも同じ円盤から分裂する感じになるので、最後尾の円盤から分裂させるようにしている。(201-207行目)

See more ...

Posted at 20:58 in Java | WriteBacks (0)
WriteBacks

Jun 10, 2008

玄箱HGに移行

10年以上前に購入し、ここ3年連続運転のwebサーバーとして活躍してきたノートパソコン、Mebius MN-340(CPU:Pentium 150MHz)の中から、チュイイイン、ギャアアァァンという、回転刃が金属板を削るような音がしてきたので、先々月から準備し始めていた玄箱に本格的に移行することにした。

初代玄箱の試験用のDebian環境をddでパーティションごとバックアップして、試験用のHDDをまともなHDDに換装して、ddでリストアしようとしたら、Debianが起動しなかった。(ddで取ったHDDのイメージファイルはLinuxでループバックデバイスとして正常にマウントできたのに。なぜだろう?)

仕方なく、Debianの再インストールを始めた。何度かインストールに失敗した後のある夜、寝る前にインストール作業をしていて、操作ミスをして慌てて電源を抜いたら、最悪のタイミングだったようで、ファームウェア異常になって玄箱自体が立ち上がらなくなってしまった。
これを修復するには、JTAGの環境が必要だ。長いこと半田こても握ってないし、そもそも半田こてを用意していない。半田こてを入手しても、その先も険しい道のりだろう。

元々、初代玄箱を試験的に動かしてみて、このMebiusの代わりは荷が重いと感じており、これまでも玄箱HGの購入を少し考えていたので、これまでに調べたことが無駄にならないように、忘れない内に玄箱HGを購入した。

●玄箱HG Debian化メモ
http://www.revulo.com/にあるページの通りの方法で、debian-sarge-2.6.17.3-kuroHG-20060702.tgzを使って、一発で成功した。
具体的な手順は以下の通り。

1. ログイン(ファームウェア1.01の初期状態だと、user=root, pass=kuroadmin)
2. echo -n NGNG > /dev/fl3
3. reboot
 →EMモードで起動する
4. 上記.tgzをtmpimage.tgzという名前にしてzip圧縮、ファームウェア1.01のimage.zipをそれに置き換える
5. ファームウェアダウンロード

●Debianセットアップメモ
- IPアドレス変更
 /etc/network/interfacesを書き換える
- Debian upgrade
 apt-get update
 apt-get upgrade
 を順に実行する
- 時刻設定
 dateコマンドで設定する

- テキストエディタの設定
 デフォルト設定では、何かにつけnanoというエディタが起動するので、
 update-alternatives --all
 を実行して、適切な選択肢でviクローンを選択する。

- IBM JDKのインストール
http://chira-ura.seesaa.net/article/28761629.htmlとそのリンク先に、パッケージの入手方法からalternativesの設定方法まで書かれている。
なお、Tomcatの起動スクリプトに書かれているIBM JDKのインストール先のディレクトリ名は、上記ページの説明とは違って、/usr/lib/j2sdk1.5-ibmとなっている。
上記ページの説明のようにalternativesを適切に設定すれば、複数のJavaをインストールした後のJava VMやJavaコンパイラの切り替えは、それぞれ
update-alternatives --config java
update-alternatives --config javac
というコマンドでできる。

- Tomcat5が起動しない場合の対処
 JavaのセキュリティマネージャをOFFにする。具体的には、/etc/init.d/tomcat5内にて
  TOMCAT5_SECURITY=no
 と設定する。

- Tomcatが使うJava VMについて
javaコマンドに関係なく、/etc/init.d/tomcat5にて優先順が決まっている。
別のものを使う場合は、/etc/default/tomcat5のJAVA_HOMEを設定する。

See more ...

WriteBacks

パソコン バックアップ

パソコン バックアップについての情報です。[link]

Posted by at 07/09/2008 07:59:00 AM

パソコン バックアップ 方法

パソコンのバックアップの方法についての情報です。[link]

Posted by at 07/09/2008 10:21:31 AM

HDD バックアップ

HDD バックアップについての情報です。[link]

Posted by at 07/09/2008 05:06:40 PM

Jun 08, 2008

ImageMagickで半透明画像を作る

前のエントリーに貼っている数式画像は、このサイトの数式画像作成CGIで作った画像を、ImageMagickで白黒反転して半透明にしている。
今回実際に行ったコマンドを書き記す。

convert -negate -resize 75% original.png tmp_nega.png
convert tmp_nega.png -channel alpha -fx "g" tmp_blur.png
convert tmp_blur.png -fx "g>0" result.png

1行目が白黒反転と縦横75%縮小、2行目が各ピクセルのα値(非透過率)の設定、3行目が背景でない部分を白にするコマンドだ。2~3行目では白の度合いを緑のレベルで代用しているが、元画像は白黒なので、赤や青を使っても同じだ。

白さに比例したα値を設定した後に色を白にするのは、そうしないと輝度が下がりすぎてしまうからである。例えば輝度50%のピクセルに50%の透過率を与えると、黒地だと輝度が25%に下がってしまう。黒地のモノクロ画像を半透明にするには、その白さをα値にして、白で塗り潰せば良いのである。
しかし、それだと、暗色の背景で透過合成表示してくれない画像ビューワで見ると真っ白になってしまうため、3行目では白の度合いが0のピクセルだけは黒にしている。"g>0"を"1"にすれば白で塗り潰すことになる。

なお、この例では変換の度に異なる画像ファイルを作るようにしているが、上書きで良いなら同じファイル名を指定しても良い。

参考:
convertの使い方ImageMagickのサイト内)

See more ...

Posted at 00:35 in PC一般 | WriteBacks (0)
WriteBacks

May 31, 2008

円盤の衝突計算

弾性衝突円盤アプレットの衝突計算について、悲しいことに式を求めるのにかなり時間を費やしたので、忘れないように記録する。これが無いと、自分でもソースコードを理解できなくなってしまう。

まず直線上の衝突について考える。物体1,2の質量をm1,m2、衝突前の速度をv1,v2、衝突後の速度をv1',v2'とすると、完全弾性衝突の場合、運動量保存の法則より
m_1v_1+m_2v_2=m_1v'_1+m_2v'_2
-\frac{v'_1-v'_2}{v_1-v_2}=1
が成り立つ。速度変化に注目してこれを整理すると、
v'_1-v_1=\frac{-2m_2(v_1-v_2)}{m_1+m_2}
v'_2-v_2=\frac{-2m_1(v_2-v_1)}{m_1+m_2}
が得られる。これは、
v_{tmp}=\frac{m_1v_1+m_2v_2}{m_1+m_2}---(1)
と置くと、
v'_1-v_1=2(v_{tmp}-v_1)---(2)
v'_2-v_2=2(v_{tmp}-v_2)---(3)
と整理できる。

次に平面上の衝突について考える。
ball_hit.png
図のように円盤と円盤が衝突する場合の速度変化については、円盤の中心同士を結ぶ線上の方向の
速度が衝突によって変化し、衝突方向でない速度成分は変化しないので、衝突方向の速度V1r、V2rについて、上記の式(1)~(3)により速度の変化を求めれば良い。つまり、 Pを円盤の中心間を結ぶベクトル、V1,V2を円盤1,2の速度ベクトル、v1r,v2rを円盤1,2のP方向の速度、V1,v1rの衝突後の速度をそれぞれV1',v1r'とすると、円盤の速度変化は
\vec{V'_1}-\vec{V_1} = \frac{\vec{P}}{|\vec{P}|}(v'_{1r}-v_{1r})
である。

v1rの計算については、θをPとV1の間の角度とすると、
v_{1r} = |\vec{V_1}|\cos{\theta}
であり、PとV1との内積は
\vec{V_1}\cdot\vec{P} = |\vec{V_1}||\vec{P}|\cos{\theta}
だから、
v_{1r} = \frac{\vec{V_1}\cdot\vec{P}}{|\vec{P}|}
で求められる。同様に
v_{2r} = \frac{\vec{V_2}\cdot\vec{P}}{|\vec{P}|}
となる。

直線の場合の式(1)と同様に
v_{tmp}=\frac{m_1v_{1r}+m_2v_{2r}}{m_1+m_2} \left (=\frac{1}{|\vec{P}|}\frac{m_1<br>(\vec{V_1}\cdot\vec{P})+m_2(\vec{V_2}\cdot\vec{P})}{m_1+m_2}\right )
と置くと、(2)より
¥vec{V'_1}-¥vec{V_1} = ¥frac{¥vec{P}}{|¥vec{P}|}(v'_{1r}-v_{1r}) = 2¥frac{¥vec{P}}{|¥vec{P}|}(v_{tmp}-v_{1r})
となる。vtmpとv1rの両方に1/|P|が入っているので、計算上の便宜のため、
v_{1r} = \vec{V_1}\cdot\vec{P}
v_{2r} = \vec{V_2}\cdot\vec{P}
と置き直すと、
¥vec{V'_1}-¥vec{V_1}=2¥frac{¥vec{P}}{|¥vec{P}|^2}(v_{tmp}-v_{1r})
¥vec{V'_2}-¥vec{V_2}=2¥frac{¥vec{P}}{|¥vec{P}|^2}(v_{tmp}-v_{2r})
が得られる。2次元のベクトルをX軸成分とY軸成分に分けて
\vec{P}=(P_x,P_y), \vec{V_n}=(V_{nx},V_{ny})
とスカラー表記にすると、
|\vec{P}|^2 = P_x^2+P_y^2, \vec{V_n}\cdot\vec{P}=V_{nx}P_x+V_{ny}P_y
なので、円盤同士の衝突による速度変化は、
\vec{V'_1}-\vec{V_1} = \left (\frac{2P_x}{P_x^2+P_y^2}(v_{tmp}-v_{1r}),\frac{2P_y}{P_x^2+P_y^2}(v_{tmp}-v_{1r}) \right )
\vec{V'_2}-\vec{V_2} = \left (\frac{2P_x}{P_x^2+P_y^2}(v_{tmp}-v_{2r}),\frac{2P_y}{P_x^2+P_y^2}(v_{tmp}-v_{2r}) \right )
v_{tmp}=\frac{m_1v_{1r}+m_2v_{2r}}{m_1+m_2}
v_{1r} = V_{1x}P_x+V_{1y}P_y
v_{2r} = V_{2x}P_x+V_{2y}P_y
で求められることがわかる。

See more ...

Posted at 22:47 in 数学 | WriteBacks (2)
WriteBacks

突然の質問失礼いたします。 直線の場合の式(1)と同様に・・・と置くとの次の式は左辺がスカラーで右辺がベクトルになっているように思われるのですがいかがでしょうか?

Posted by 谷口 at 10/22/2009 10:40:46 PM

ご指摘ありがとうございます。左辺を書き間違えていました。v1r'-v1rでなく、V1'-V1でした。v1r'-v1rとなってた箇所を修正しました。

Posted by ynomura at 10/25/2009 12:17:14 AM

May 24, 2008

数式画像作成CGI

TeXの数式をPNG画像に変換するCGIを作ってみた。





例:\sum_{k=0}^{n-1}x^k = \frac{1-x^n}{1-x}

See more ...

Posted at 23:31 in PC一般 | WriteBacks (0)
WriteBacks

May 18, 2008

弾性衝突円盤アプレット

何かJavaのGraphics2Dクラスを使ったプログラミングの練習をしようと思って、昔作ったことがある思い出のアプリをJavaで作ってみた。
弾性衝突円盤アプレットのページへ

各円盤は面積に比例した質量を持ち、壁や他の円盤に完全弾性衝突する。
後は説明不要だと思う。
Pentium3だと動作がちょっと重いかも知れない。

大昔、X68000のHuman68k(つまりウィンドウシステムでない)上で動くものを作った時は、アセンブラで書いてスプライトを使って衝突計算アルゴリズムもかなり練って速度を出したのだが、今回、Javaのコードをほとんど工夫せずに書き、描画処理は単純に消して描いての繰り返しで、しかもJavaアプレットとしてブラウザ上で動かすと、遜色ない速度が出たので、驚いた。68000MPUとCeleron 3.2GHzの性能はかくも違うのか。PCの性能の進歩はなんと速いことか。

割とシンプルなコードになったので、ソースコード、説明を添えて公開するつもりだったが、時間が無かったので、とりあえずアプレットだけ先に公開することにした。

円盤の描画は、先にBufferedImageの画像を用意しておいて、背景を塗りつぶしてGraphics#drawImageでペタペタ貼っているだけである。当方の環境ではちらつきが目立たないので、ダブルバッファリングは省略した。
円盤同士の衝突の計算は、円盤の速度の衝突点方向の成分を求めて、高校の物理で習う完全弾性衝突の式を当てはめて、その方向の速度の変化を求めている。

結局Graphics2Dの機能は使わなかった。非矩形の画像の重ね合わせに必要だと思ったのだが、java.awt.imageのBufferedImageクラスのαプレーンを使えば難なくできた。
別の新たなネタを探さねばならない。

See more ...

Posted at 21:24 in Java | WriteBacks (0)
WriteBacks

Apr 27, 2008

MySQL Connector-Jでの再接続

JDBCのDriverManager#getConnectionによるデータベースへの接続は結構時間がかかるので、短時間に何度もDBにアクセスする場合は、Connectionオブジェクトを保持したい。
しかし、MySQLのデフォルト設定の動作では、DBへの接続は8時間使われないと切られてしまう。常駐型のservletなどの連続運転のアプリでは、途中でこのコネクションが切れることがあり得る。この時間を変更することもできるらしいが、いつか切れるので、それで解決にはならない。従って、切れれば再接続する必要がある。

JDBCアプリ側で再接続する場合、DBへの接続が有効か無効かを知る術が問題になる。接続が切れた後のクエリ送信時の例外はSocketクラスのIOExceptionであり、ドライバ(MySQL Connector-J)にてcatchされるので、JDBCアプリではcatchできない。JDKのAPI仕様書を見るとConnection#isClosedというメソッドがあるが、これはJDBCアプリからcloseしたかどうかを返すものであり、MySQL側で接続を切られてもclosedにはならない。

そもそもservletでは、こういう時の定石はconnection poolingらしい。ServletエンジンでDBへの接続を保持して、servletからの要求に応じて接続を貸し出す仕組みで、暇な時に再接続もしてくれる、J2EEに標準装備されている一般的なものである。勿論Tomcatでも使える。
しかし、Tomcatのマニュアルの"JNDI Datasource HOW-TO"のページを見ると、設定方法がかなり複雑である。もっと簡単な解決方法は無いかと思ってConnector-Jのマニュアルを眺めていると、autoReconnectという設定項目が見つかった。DBとの接続が切れていればドライバが自動的に再接続するという設定だ。Obsoleteなproperty(将来削除される設定項目)であり、SQLExceptionを処理できない場合以外は非推奨と書かれているが、今回はとりあえず目的が達成できれば良かったので、

conn = DriverManager.getConnection("jdbc:mysql:/
/localhost/xxxx?user=xxxx&password=xxxx&autoReconnect=true");

という風にautoReconnectの設定を追加した。

実にシンプルかつスマートに解決できたと思っていたが、Tomcatのログを見て、時々エラーが発生していることに気付いた。調べてみると、どうやらautoReconnectの設定をしていても、DBとの接続が切れると、1回はクエリ送信が例外(SQLException)で終わり、それによって再接続され、その次のクエリから接続が有効になるようだ。
そこで、1回はSQLExceptionに対して再試行する、ということも考えたが、そもそもSQLExceptionをcatchして再接続するのならautoReconnectの機能を使う意味が無い、connection poolingの方がスマートだろう、と考えてやめた。

という訳で、connection poolingを試すべく、TomcatのHowToの通りに設定を書き、テストプログラムをコンパイルすると、javax.sql.*が無いというエラーになった。一瞬で確信した通り、プールされたconnectionを取り出すためのDataSourceクラスが1.3のJVMには存在しないのだった(javax.sqlは1.4以降の仕様らしい)。

途方に暮れてさらにConnector-Jのマニュアルを読み進めると、"Common Problems and Solutions"の所に答えが書いてあった。Connection poolingするか、autoReconnectにした上でSQLExceptionをcatchし、MySQLのエラーコードを調べて、回数の上限を決めて再施行するか、のどちらかをする必要があるそうだ。SQLException#getSQLStateでMySQLのエラーコードがわかり、MySQLのinfoによると接続系のエラーが起こるとSQLStateが"08S01"になるらしいので、SQLStateが"08S01"の場合に最大2回クエリを再施行することで解決した。

int retryCount = 3;
boolean transactionCompleted = false;
do{
try{
Statement st = conn.createStatement();
ResultSet rs = st.executeQuery("xxxx");
transactionCompleted = true;
}
catch(SQLException e){
String sqlState = e.getSQLState();
if (sqlState.equals("08S01")){
retryCount--;
}else{
retryCount = 0;
}
}
finally{
...
}
} while(!transactionCompleted && retryCount > 0);

See more ...

Posted at 17:59 in Java | WriteBacks (0)
WriteBacks

Apr 20, 2008

玄箱+Debianで試験運用中

このウェブサイトを支えているPentium1機の動作が不安になってきたことと、ノートPCだが古いので消費電力が気になってきたこともあって、かねてから気になっていた玄箱を(主に玩具としてだが)購入した。
Web上に情報が豊富にあるため、Debian化はすぐにできた。さらに、Debianのパッケージ管理ツール(apt-get/aptitude)がとてもわかりやすく便利なので、webサーバとしてのセットアップも大体楽にできた。玄箱のこともDebianのことも全く何も知らなかったが、結局組み立てから1週間で最低限のwebサーバ化ができた。

MovableType、Tomcatの移行もできたので、試しにしばらく玄箱でこのサイトを動かすことにした。

しかし、早速いくつか問題が起こっている。
CPUの速度は、計算だけのベンチマークテストでは玄箱(PowerPC 200MHz)の方が倍くらい速いようだが、実用上の動作は玄箱の方が遅いケースもある。CGIの動作は常に玄箱の方が遅い。ApacheもPerlも同じバージョンで、設定にも大差ないので、OSが遅いかPerlが遅いくらいしか原因が思いつかない(Pentium 1機のFreeBSDの方はカーネルの再構築もPerlの再コンパイルもしている)。

RAMが少ないのも問題だ。玄箱はRAMが64Mしかないので、Tomcatやmod_perlのMovableTypeを動かすとswapが多発する。例えば、Tomcatを動かしながらcpanを動かすと、時々シェルがキー入力に反応しなくなり、リモートから操作困難になった。

RAM 128Mの玄箱PROを買うべきだったか?

See more ...

WriteBacks

Tomcat5をなるべく実メモリ少なく動かすために、Java VMをいくつかインストールして試してみた。 Java VMの使用RAM領域を完全にスワップアウトしてから小さいservletを動かす方法で、GCJとKaffeとIBM Java VM(for 32bit-ppc ver.142 SR9)(全てJDK1.4.2互換)を比較した所、Tomcat5のservletの動作に必要な実メモリの最小値は、GCJとIBM Java VMが約10M、Kaffeが約6Mだった。 但し、Kaffeは起動時間が半端じゃなく長い。GCJやBM Java VMで2~3分のTomcatの起動に1時間以上かかる。KaffeのJavaコンパイラも強烈に遅いので、残念ながら玄箱では使い物にならない。

Posted by ynomura at 04/26/2008 12:45:34 AM

もう少しMebius MN-340(Pentium 150MHz, RAM 96M)と玄箱の性能を比べてみると、PerlのCGIの動作速度がほぼ互角な以外は、全体的に玄箱の方が1.2~2倍くらい速いようだ。特に、MN-340はHDDが遅いため、swapが比較にならないほど遅いことがわかった。10年以上前に買ったマシンだけのことはある。色々実験していて、64Mのmallocに1分以上かかったのには笑ってしまった(40GのHDDを積んだ玄箱では遅くても4秒くらいで終わる)。RAMを多目に増設していたので、これまでかなり救われていたようだ。

Posted by ynomura at 04/26/2008 08:18:35 PM

Apr 19, 2008

MySQLと日本語文字

MySQLで日本語文字を使うために調べたことを記録する。

●基本事項
サーバーではデータベース毎に文字コードセット(コーディング)を決めることができる(表毎、列毎の指定も可能)。またサーバー~クライアント間のデータ送受信時に使用する文字コードセットは随時変更することができる。
従って、DB上の文字列の文字コードに関わらず、クライアントから所望の文字コードセットで文字列を送受信することが可能である。

MySQLはコンパイル時(configure時)にデフォルトで全ての文字コードセットに対応している訳ではないので、日本語文字コードセットに対応するようにコンパイルする(configure実行時に--with-charsetオプションを付ける)必要がある。

動作中のMySQLで使用できる文字コードセットは、クライアントにて

SHOW CHARSET;
で見ることができる。

●サーバー側の文字コードセット指定
デフォルトの文字コードセットは、my.cnfの[server]エントリーのdefault-character-setで指定できる。
特定のDBに対しては、CREATE DATABASEまたはALTER DATABASE文にてCHARACTER SETオプションで指定できる。

●クライアント側の文字コードセット指定
クライアントで使用する文字コードセットは、

SET NAMES '文字コードセット名';
で指定できる。EUCは'ujis'、UTF-8は'utf8'、Shift-JISは'sjis'である。
実際には、これによって関係する複数の環境変数が同時に更新される。
文字コードセットに関係する環境変数の値は、
SHOW VARIABLES LIKE 'character%';
とすると全て見ることができる。

なお、MySQLのinfoを含め色々な所にmy.cnfの[client]エントリーのdefault-character-setにてデフォルトの文字コードセットを設定できるようなことが書かれているが、筆者の環境(FreeBSD 4.11+MySQL 5.0.27)ではその設定は反映されない。Webを見ると、他でも同じ現象が少なからず発生しているようだ。my.cnfにてこれを設定する場合は、1回mysqlを起動して環境変数が意図通りに変わっているかどうか確認しておくべきだと思う。

●LOAD DATA INFILEで使われる文字コードセット
MySQLのLOAD DATA INFILEでファイル読み出し時に使われる文字コードはcharacter_set_databaseの値であり、上述のSET NAMESでは更新されない。従って、対象ファイルにてデフォルトでない文字コードが含まれる場合は、先に

set character_set_database=文字コードセット名;
のようにして設定する必要がある。
なお、文字コードによっては区切り文字の判定に曖昧さが出る可能性があるため、日本語文字列は""(ダブルクォート)で括ってLOAD DATA INFILE文にFIELDS ENCLOSED BY '"'(シングルクォート2つの間にダブルクォート1つ)を指定するなどした方が良いと思う。
(参考:http://www.hirohama.biz/mysql/2008/01/31-100131.html

See more ...

Posted at 21:27 in PC一般 | WriteBacks (0)
WriteBacks

Apr 17, 2008

Javaと日本語文字

Javaアプレットとサーブレットを作っていて、日本語文字、エンコーディング関連でつまずいたが何とか解決したことを記録する。
ここに書くものは、1つの成功例であって、最善の方法であるとは限らない。

●Servletからappletへの日本語文字列の受け渡し
Applet-servlet間の通信は通常HTTPで行うことになるので、appletからはjava.net.URLConnectionを使うとする。
Servlet側では、javax.servlet.HttpServletResponse#setContentType()でエンコーディングを指定する。
例:

public void doGet(HttpServletRequest req, HttpServletResponse res){
res.setContentType("text/html; charset=utf-8");

Applet側では、InputStreamReaderを使って、送られてきた文字列を読み出す。InputStreamReaderの作成時にエンコーディングを指定する。
例:
URL url = new URL("http://...");
URLConnection uc = url.openConnection();
InputStreamReader isr = new InputStreamReader(uc.getInputStream(), "utf-8");
BufferedReader br = new BufferedReader(isr);
String s = br.readLine();

●ソースコードへの日本語文字の記述
Swingのラベルや標準出力への出力文字列など、プログラムが表示する日本語文字をソースコードに直接埋め込む場合は、javacコマンドでコンパイルする時に、ソースコードで使用しているコーディングの名前を-encodingオプションで指定する。
例:

javac -encoding EUC-JP eucfile.java

コーディング名は、日本語だと"EUC-JP","ISO-2022-JP","UTF-8","Shift_JIS"などがある。
(参考:encoding.doc

●APIからの日本語文字の取得
JDBC等のAPIがJava内部で使用するものでないコーディングで文字コードを返す場合は、文字コードの変換が必要になる。Stringクラスのコンストラクタを使うとJavaのコーディングに変換できる。
例:

ResultSet rs;
String jstr = new String(rs.getBytes(1), "euc-jp");

●LinuxでJavaアプレットを表示すると日本語フォントが出ない(□になる)場合の対処
Mozilla上のJavaアプレットの表示がそうなるのも、appletviewerでそうなるのも共通の原因で、JREに適切にフォントパスを設定すると解決すると思われる。
JREのインストールディレクトリのlib/fonts/の下にfallbackというディレクトリを作成し、日本語を含むフォントファイルへのシンボリックリンクを置く。
例:

cd /usr/java/jdk1.5.0_07/jre/lib/fonts
mkdir fallback
ln -s /usr/share/fonts/ja/TrueType/kochi-gothic.ttf fallback/
ln -s /usr/share/fonts/ja/TrueType/kochi-mincho.ttf fallback/

(参考:JavaSE 5.0のリリースノート (日本語)(English)

See more ...

Posted at 21:56 in Java | WriteBacks (0)
WriteBacks

Apr 13, 2008

JDBC+servlet連携appletを作ってみた

せっかくこの長老マシンに苦労してMySQLとTomcatをインストールしてJDBCとservletが動くようになったのに、これまで簡単なテストプログラムしか動かしておらず、有効活用できていなかったし、あまり勉強にもなっていなかったので、Javaの復習とSwingの勉強を兼ねて、JDBC servletと連携するJavaアプレットを作ってみた。
3択首都当てクイズのページへ
アプリケーションのネタが思いつかなかったので、とりあえず3択クイズにした。

Servletとappletの役割分担は、MVCモデルに従って、クイズの作成、正誤判定、終了判断といったモデルは全てservlet側で行い、applet側は表示だけを行うようにした。すなわち、Mをサーバー側、VとCをクライアント側に分けた。そのため、appletは1問ずつ問題をservletから取得し、解答者の選択をservletに送っている。クイズの終了判定もservletが行うため、servletから問題でなく成績表が送られてくると終了という設計にした。

セッション管理にはjavax.servlet.HttpSessionを使っており、複数のクライアントが同時にクイズを始めても、servletで全てのクライアントの途中経過が管理される。Cookieが無効にされているとサーバー側からセッションIDが送られ、クライアント側でURLのquery部にセッションIDを付けるようにしているので、CookieがOFFでも動くはず。

クイズの問題はMySQLのデータベース上にあり、servletは1問ずつJDBC経由で問題を取得する。1つのセッションの中で問題が重複しないように、出題する問題の番号だけはセッションの新規作成時に全て決めておくようにした。

クライアント側の描画にはSwingを使用した。CardLayoutでペインを切り替えるようにしている。クイズ画面の動く背景は、SwingタイマーとAWTを使って描画している。最初はこの背景描画が非常に重く、試行錯誤して改良している内に、事前に画像ファイルを用意してタイル貼りする方が負荷が軽いことが判明したが、画像ファイルを作成するのが面倒なので、そこまではしなかった。

See more ...

Posted at 20:56 in Java | WriteBacks (2)
WriteBacks

よければjavaのソースを教えてください。

Posted by ななし at 01/11/2012 01:57:35 PM

ここに置いてます。

Posted by ynomura at 01/11/2012 09:33:20 PM

Apr 02, 2008

Tomcat起動時のエラーログ

TomcatでJavaのservletを動かす実験をしていて、ログを見ると、servletは正常に動くのにTomcatが起動時に下のようなもので始まる大量のエラーを出してることに気付いた。

[ERROR] Digester - Parse Error at line 37 column 15: The content of element type
"servlet" must match "(icon?,servlet-name,display-name?,description?,(servlet-c
lass|jsp-file),init-param*,(以下略)

"Parse"とあるので私が追加したweb.xml(XML形式の設定ファイル)がおかしいのだろうとは思ったが、そもそも例えば上で書かれてるような37行目に"servlet"という文字列は無いし、"must match ..."とか言われても、"servlet-name"エレメントと"servlet-class"エレメントしか使ってないので、きちんとmatchしてるように思え、何が悪いのかわからなかった。しばらくしてやっと、上記の右側が正規表現ぽいことから、"servlet"エレメントの中身がicon, servlet-name, display-name, ...がこの順でないといけない、"?"がついてるエレメントはあっても無くても良い、servlet-classまたはjsp-fileのどちらか1つが必須である、という意味だと気付いた。

<servlet>~</servlet>の間に<servlet-name>と<servlet-class>を交互に書いていたのが問題だったようだ。同様に、<web-app>~</web-app>の間では、<servlet>と<servlet-mapping>を交互に置いてはいけないらしい。

See more ...

WriteBacks

Mar 30, 2008

自己出力プログラム

Perlで自分自身を出力するプログラムが作れるという話を聞いたので、考えてみた。

print文1行だけではどう考えても作れないので、コードを

C1; print C2;

(C1,C2は何らかのコード)という形だと仮定する。当然C2は"C1; print C2;"に展開されるコードということになるが、何も無い所から"print"なんて文字列を含む値を生成するモジュールは無いので、これを実現するには、媒介変数を使って"print"を複製するしかない。従って、コードを改めて
C1; $x = "C2printC3"; print "C4$xC5$xC6";

という形に仮定し直す。このコードで表示される文字列は
C4C2printC3C5C2printC3C6

となるが、ここで、C3の部分がダブルクォーテーションで始まり、「$x = "…"」ではなく「$x = q()」または「$x = qq()」を使わないと厳しいことがわかる。従って、コードをさらに
C1; $x = q(C2printC3); print "C4$xC5$xC6";

という形に仮定し直す。このコードで表示される文字列は、やはり
C4C2printC3C5C2printC3C6

となる。

「C4C2」の部分は「C1; $x = q(C2」に一致しないといけないので、C4=「C1; \$x = q(」(C4はダブルクォーテーションの中にあるので"$"はエスケープが必要)となる。
「C3C5C2」の部分は「C3); 」に一致しないといけないので、「C5C2」=「); 」ということになるが、C2はq()の中にあるため")"は使えないので、C2は空、C5=「); 」となる。
「C3C6」の部分は「 "C4\$xC5\$xC6";」に一致しないといけないが、両辺にC6があり、C6の末尾が「C6";」となることから、C6は空ということになる。

仮定したコードからC2,C4,C5,C6を消して書き直してみる。

C1; $x = q(printC3); print "C1; \$x = q($x); $x";

「C3C6」=「 "C4\$xC5\$xC6";」より、C3=「 "C1; \$x = q($x); $x";」となるので、コードは
C1; $x = q(print "C1; \$x = q($x); $x";); print "C1; \$x = q($x); $x";

となる。「C1; 」の部分は""による展開に気をつければ何でも良いということがわかる。

例えばC1を「$a = 0」にすると、コードは

$a = 0; $x = q(print "\$a = 0; \$x = q($x); $x";); print "\$a = 0; \$x = q($x); $x";

となり、「C1; 」を空にすると、コードは
$x = q(print "\$x = q($x); $x";); print "\$x = q($x); $x";

となる。これらのコードを実行すると、コード自身と全く同じ文字列を出力する。

See more ...

Posted at 18:01 in PC一般 | WriteBacks (0)
WriteBacks

Jan 05, 2008

第3正規形とBoyce-Codd正規形

久々にSQLの入門書を読んで復習しているが、これまでRDBの第2正規形や第3正規形をよく理解していなかったので、この機会にきちんと理解しようとしている。
この入門書には簡単な例だけで正規形の正確な意味が書かれていないので、ネットで調べながら勉強していると、ボイス-コッド正規形というものがあることを知った。

第1~3正規形とボイス-コッド正規形の意味を非常に浅く理解した後、ボイス-コッド正規形は必ず第3正規形の条件を満たすというのを読んで、理解するのに結構な時間がかかったので、忘れないように、自分なりの理解で図を書いてみた。

第2正規形とは、「主キーの値が決まれば、他の列の値が決まる」と書いてあることがあるが、私の理解では「主キーまたはいずれかの候補キーの値が全て決まらない限り、他の列の値はどれも決まらない」と書かないと意味が取れないと思う。
第3正規形は、第2正規形でかつ「候補キーでない列の値よって他の列の値が決まることがない」と理解した。
つまり第2正規形は主キーに対して制約がある形、第3正規形は非キー属性にも制約がある形と言えるのではないだろうか。

Boyce-Codd正規形は、「属性間の自明でない依存関係X→Yの全てについて、Xは候補キーかそれを含む集合」ということなので、「候補キーを含まないキーによって定まる属性はない」と理解したが、これだけでは第2~3正規形との関係がよくわからないので、図を書いて考えてみる。

2_3_bcnf.png

上の図において、A~Eは属性(列)で、ABが主キー、BCは(主キー以外の)候補キーであり、第1正規形である(第1正規形の条件を満たす)とする。
ABの値が決まればD,Eの値が決まるので、AB→D、AB→Eの依存関係がある。

B→Dの依存関係があると、候補キーの値が揃わなくても非キー属性の値が決まってしまうので、第2正規形でない。
D→Eの依存関係があると、候補キー以外の属性(非キー属性)によって決まる非キー属性があることになるので、第3正規形でない。(または、AB→D、D→Eで非キー属性が推移的に候補キーに従属するので、第3正規形でない)
C→Dの依存関係は、非キー属性によって決まる非キー属性があることにはならないので、第3正規形の条件を破らないが、候補キーを含まないキーによって定まる属性があることになるので、Boyce-Codd正規形の条件は満たさない。

See more ...

Posted at 11:56 in PC一般 | WriteBacks (0)
WriteBacks

Jan 03, 2008

データベース用語和英対応表

6~7年前に買ったSQLの入門書を、捨てる前に読み返している。この入門書を使って1回SQLを勉強したのだが、実際に使うことが無かったため、全く身に付かず、歳のせいで記憶にも残っていないのだ。実際、MySQLを触りながら復習しようとして、mysqlを起動すると、"select * from (テーブル名);"以外の文はそらでは全く書けなかった。

従って、SQLの入門書とMySQLのマニュアル(MySQL info)とを見ながらSQLを試しているのだが、SQLの入門書が日本語でMySQLのマニュアルが英語であり、筆者にデータベースの基礎知識がないため、日本語と英語の対応が取れない用語がいくつか発生した。

そこで、出くわしたデータベース用語の日本語と英語の対応表を作ることにした。

日本語英語
階層型データベースhierarchical database
ネットワーク型データベースnetwork structure database
リレーショナルデータベースrelational database, RDB
データ定義言語Data Definition Language, DDL
データ操作言語Data Manipulation Language, DML
データ制御言語Data Control Language, DCL
主キーprimary key
候補キーcandidate key
代理キー(代替キー)alternate key
複合キー(連結キー)composite key
外部キーforeign key
非キー属性non-key attribute
正規化normalization
非正規化denormalization
非正規形non-first normal form. NF2
第1正規形first normal form, 1NF
第2正規形second normal form, 2NF
第3正規形third normal form, 3NF
ボイス-コッド正規形Boyce-Codd normal form, BCNF
第4正規形forth normal form, 4NF
第5正規形fifth normal form, 5NF
射影-結合正規形projection-join normal form, PJNF
関係代数relational algebra
集合演算set operation
関係演算relational operation
sum
difference
product
quotient
直積(デカルト積)direct product(Cartesian product)
選択selection
射影projection
結合join
交差結合cross join(Cartesian join)
等結合equi-join
自然結合natural join
内部結合inner join
外部結合outer join
左外部結合left outer join
右外部結合right outer join
自己結合self-join
相関サブクエリーcorrelated subquery
集計関数aggregate function
参照整合性制約referential integrity constraint
一意性制約unique constraint
(2008/1/4 更新)

See more ...

Posted at 21:55 in PC一般 | WriteBacks (0)
WriteBacks