乱数の質その2

モンテカルロ法だとよくわからなかったglibcとMelsenne Twisterの違い。ならば単純にヒストグラムを取ったらどうなるという訳でテストです。
一様乱数だと処理が面倒なので、0から9999までの乱数を発生させて、各数値の発生回数をカウントします。

全体で1010回の試行の結果は次の通り。

分散の値はMTの方が大きくなっているけど、そういうものなのかな??

実際的にはホワイトノイズ(例のシャーってノイズ)にA/Dコンを接続したものや量子力学的現象を利用するとかいろいろな方法が提案されていますが、何が一番なんだろうな。というか、乱数の質ってなんだか哲学的な話になってきたので簡単には結論は出せないんでしょうけど。

乱数の質?

乱数を使用してあれこれという場面は結構ありますが、質がどうなのかのチェックはしていなかったと思い、ちょっとテストした結果です。

比較にはLinux標準のglibcに含まれるrand()と Mersenne Twisterと呼ばれる手法によるgenrand_real1()を利用して一様乱数を発生させてどんな感じなのかを調べます。ちなみにMersenne Twisterは非常に質が良いとされている最近出てきた乱数生成手法。
ただヒストグラムを取っても面白くないので、モンテカルロ法で円周率を計算させてみます。

プログラム自体は単純なので省略。各10億(109)点(20億組)から円周率を計算することを1000回繰り返すので、1兆(1012)組の乱数を発生させることになります。
MT.hの中身を見ていないのですが、けっこう面倒なことをしているようで、rand()が360分かかったのに対して、genrand()の方は415分が所要時間。

で、結果はどうなのかというと、
円周率計算結果
と大差無しの結果が得られました。なんか肩透かしだな。というよりrand()もなかなかやるじゃないかという訳で安心して使用できそうだという結論。それともあと4桁くらい余計に生成させると差が出るんだろうか?  ただ4桁増やすと、結果出るのが4150000分後って288日もかかるんでやる気は無いけど。

追記:
もしかして、109個ってのが多すぎだったかも。3桁くらい減らしてみたら、違いが出るか?

多角形からの円周率

円周率の算出方法として内接する多角形の一辺の長さの合計から割り出す方法があります。

基本から考えると、こんな感じかな。
まず円に内接する多角形を考える。この多角形は中心を頂点とするn個の二等辺三角形に分割でき、底辺の長さの合計が円周に近似される。
さて角数をnとすると、各三角形の頂角θは

であり、各三角形の底角θ1

頂点から底辺の中央に向けて垂線を引くと、底辺の1/2の点で交わるので、円の半径をrとすれば、底辺の1/2の長さr1

となるので、n角形の辺の長さの合計Rは

と表される。
元の円の半径rを0.5とすれば、
となる。

例えば、4角形ならば4cos(90-180/4)=2.828…、5角形なら5cos(90-180/5)=2.938…
120角形ならば、120cos(90-180/120)=3.141233…
といった具合。

つまり、

とまぁ、Texを使えば数式をそれらしく書けるのが嬉しくて下らないことをだらだら書きました。

現代では三角関数が使えるので、それこそ電卓でも持ち出せば1000角形だろうが一瞬で計算OKと。それでは意味はないので、幾何的に処理することを考えないと。
幾何的に処理となると、三平方の定理を使ってやればいいのかな? ちょっと図を書いて考えてみよう。




KiCAD日本語対応その後

メインシステムであるFedoraが30になったので、そちらの対応に追われていました。
Fedora30ではKiCadが5.1.2となっています。例によって.srpmを調達してnewstroke_font.cppを入替え。無事にLinux版はビルドに成功しました。

それにしても、バージョンアップの度にこれだと面倒だな… Windows版だとWinbuilderにまた変化が有ったようなのでこちらは手強そう。

ところで、この日本語対応版KiCADのバイナリって需要があるのかな? 本家kicad.jpも活動が停まっているように見えるしな。

Fedora30

新しい版の発表が先日ありました。予定では7日とのことで、連休明けだなぁと思って居ましたが、運良く休み期間中に発表になったので早速入替え。半年に一度の煩わしい楽しい時間の始まりです。

とりあえず最小限インストールして、何時のように必要なものをインストール。
さすがに半年に一度インストールしている関係で環境構築はすっかり慣れし親しんだルーチンワークで、約2時間ほどで8割方元の環境が戻りました。あとは開発系のパッケージだけど、なにが必要化はよく掴んでないので、必要と言われるままに入れていくしか無いな。

例によって再ビルドするパッケージが山積みとなりますが、楽しくしばらく遊べそう。

演算精度その2

単純に整数の逆数を足し込んでいくだけのライプニッツ級数は収束が遅いことがわかりました。なので、べき乗項の有るオイラー級数だったらどうなのかをさっそくテスト。大して計算式に変化はありませんのでちょいちょいと変更してコンパイル。

#include <stdio.h>
#include <quadmath.h>
#include <math.h>

void main(void)
{
    char    c1[80];
    long    loop;
    __float128  res=0.0Q,r1=1.0Q;
    
    for (loop=1;loop<100000000;loop++) {
        res=res+1.0Q/(r1*r1);
        r1=r1+2.0Q;
    };
    res=sqrtq(res*8.0Q);
    quadmath_snprintf(c1,sizeof(c1),"%.40Qf",res);
    printf("loop=%ld\tresult=%s\n",loop,c1);
}

結果は、
$ gcc -o paicalc2.o paicalc2.c -lquadmath
$ ./paicalc2.o
 loop=10000000 result=3.1415926217588012757267165560912027187844
$ ./paicalc2.o
 loop=100000000 result=3.1415926504066943431811710021935384220221

と言う訳でライプニッツ級数の1億項の計算に比べて1000万項の計算で1桁精度が向上しました。項数を同じにした場合は2桁の向上ですね。でもそんなものか…

実際にはその後いろんな方法が考案されていて、そちらを使えばってのが身も蓋もない言い方。実際ガウス・ルジャンドル法だと3回目には18桁正しい値が得られるし、昨今の多数桁計算にはマーチンの式が使用されている。
とは言え、今から300年前に手計算だけでこれにたどり着いたのだからすごいことだと思うな。

演算精度の話

大先輩のブログでライプニッツ級数で円周率の計算って話が出ていました。もともとはプログラム電卓でどこまでやれるのかということが主眼なので、でかいPCを持ち出すのは野暮の極みといったところですが、ついつい悪乗り。

gccでもVer.4.8以降だと4倍精度小数点つまり128bit floatが使用できるとのことでしたので、さっそくテスト。
ちなみにソース(paicalc.c)はこんな感じ。面倒なので定義通りの記述。なんかねぇ…

#include <stdio.h>
#include <quadmath.h>
#include <math.h>

void main(void)
{
    char    c1[80];
    long    loop;
    __float128  res=0.0Q,r1=1.0Q,r2=3.0Q;
    
    for (loop=1;loop<100000000;loop++) {
        res=res+1.0Q/r1-1.0Q/r2;
        r1=r1+4.0Q;
        r2=r2+4.0Q;
    };
    res=res*4.0Q;
    quadmath_snprintf(c1,sizeof(c1),"%.40Qf",res);
    printf("loop=%ld\tresult=%s\n",loop,c1);
}
コンパイルして実行した結果は、
$ gcc -o paicalc.o paicalc.c -lquadmath
$ ./paicalc.o
 loop=100000000 result=3.1415926485897931884626429145303967199477

円周率の頭の方の分かっている数値は3.141592653589793…なので4倍精度だからといって7桁程度しか出てませんが、式から言えば当たり前ですね。
それとこのquadmathライブラリは当たり前ですが、ソフト処理なのでしょうから通常の64bit演算に比べると大幅に時間がかかります。倍精度の場合は0.76秒で済みましたが、4倍精度の場合は19.4秒と25倍以上も時間がかかっています。

オイラー式だとべき乗級数なのでもうちょっと収束が早いかな。ただ結果が円周率の2乗なので、開平演算の誤差が乗ってきそう… まぁ明日試してみますか。

Windows10初期化その4

結局分かったのは、CADソフトであるDraftsightの2019SP0という最新版を入れるとおかしくなると言うこと。

レジストリの例のCacheを消しても改善しないので、一旦アンインストール。そしてShellExtensionのCacheの中身を消すことで問題は解決するのですが、なんか気分が良くない。一方旧版であるDraftsight2018を入れる分には問題が起きないので、2019版はSP1が出るまでお預けです。
Readme.txtを読んでないもので、2018->2019の変化がなんなのかよくわかりません。なので2018のままでも特に支障は有りませんが。

まずはメーカーに文句を言うべきなんでしょうけど、どうも英文しか受け付けない雰囲気なのでちょっと面倒そうです。このあたりは片手間で日本でも売っている製品を利用する点でのリスクですね。

PDF関係のツールと音響関係のツールを戻したところで、クローン作成で一段落。
あとはVisiualStudioだけど最新の2019を入れるべきか2017にしておくかで悩み中。

Windows10 初期化その3

またまたずっこけたWindows10。今まで普通に使えていたアプリケーションをインストールしただけでおかしくなったのが納得のいかない点なのです。

で、バックアップから書き戻すつもりが最新版、つまりずっこけた状態をバックアップしてしまうというドジを踏んでしまったので、やむなく本腰を入れて原因追及を開始。

ファイルを右クリックして「プログラムから開く」で無反応となるトラブルなので、PathかShell系の設定に異常を疑い調査。
Pathには特に不審点は無いので次にShell系の調査。

Regeditで\HKEY_CURRENT_USER\Software\Microsoft\Windows\CurrentVersion\
の中から怪しそうなShell Extensionをチェック。すると中にCachedという項があってずらずらと値が並んでいます。残念ながら、全部Hash化された値なので何に該当するのかはさっぱり。なのでエイヤっと全部削除して再起動。その結果トラブルが解決したのでした。

さて、怪しいアプリケーションを再インストールしたら現象が再現するんだろうか? そして再現した場合、同じ処置で解決するのだろうか? なんかもやもやします。

連休はこんなことに時間を使うつもりはなかったんだがなぁ…

Windows初期化その2

使用頻度の高い順に戻していって、8割方終わったところで、異常がないことを確認。この先また異常が出ると面倒なのでここで一旦クローンを作成。
これで万一この後こけても復旧の手間がだいぶ省ける。まぁ何も起きないに越したことはないのですが。

CADソフトのライセンスをいったん戻しておくのを忘れたので、どういう扱いになるのかがすごく心配。ダメと言われるとまた$99の支払いが発生するので、どうにかして回避しなきゃ。

4/27追記
心配していたことが再発。なんだろうなぁ… とりあえず怪しいやつは特定できた気がするので、そいつを除いて全体を再構築しなきゃ。