Juliaで任意精度演算

Juliaで階乗を計算する関数は、その名もfactorial()です。たとえば20!を計算するには以下のようにすれば一発です。

julia> factorial(20)
2432902008176640000

ただ、デフォルトでは64ビット符合付き整数での演算なので、その範囲を超える場合にはエラーになってしまいます。

julia> factorial(21)
ERROR: OverflowError()
Stacktrace:
 [1] factorial_lookup(::Int64, ::Array{Int64,1}, ::Int64) at ./combinatorics.jl:31
 [2] factorial(::Int64) at ./combinatorics.jl:39

big()を使ってあげると最高の精度で演算ができるようにBigIntやBigFloatに変換してくれます。

julia> typeof(21)
Int64

julia> typeof(big(21))
BigInt

そうすることでfactorial()もInt64ではなくBigIntでの計算に切り替わりますので、正しく階乗を計算できます。

julia> factorial(big(21))
51090942171709440000

BigIntにしてしまえば、もっと大きな数についても階乗が計算できます。（しかもJuliaは高速！）

julia> factorial(big(1000))
402387260077093773543702433923003985719374864210714632543799910429938512398629020592044208486969404800479988610197196058631666872994808558901323829669944590997424504087073759918823627727188732519779505950995276120874975462497043601418278094646496291056393887437886487337119181045825783647849977012476632889835955735432513185323958463075557409114262417474349347553428646576611667797396668820291207379143853719588249808126867838374559731746136085379534524221586593201928090878297308431392844403281231558611036976801357304216168747609675871348312025478589320767169132448426236131412508780208000261683151027341827977704784635868170164365024153691398281264810213092761244896359928705114964975419909342221566832572080821333186116811553615836546984046708975602900950537616475847728421889679646244945160765353408198901385442487984959953319101723355556602139450399736280750137837615307127761926849034352625200015888535147331611702103968175921510907788019393178114194545257223865541461062892187960223838971476088506276862967146674697562911234082439208160153780889893964518263243671616762179168909779911903754031274622289988005195444414282012187361745992642956581746628302955570299024324153181617210465832036786906117260158783520751516284225540265170483304226143974286933061690897968482590125458327168226458066526769958652682272807075781391858178889652208164348344825993266043367660176999612831860788386150279465955131156552036093988180612138558600301435694527224206344631797460594682573103790084024432438465657245014402821885252470935190620929023136493273497565513958720559654228749774011413346962715422845862377387538230483865688976461927383814900140767310446640259899490222221765904339901886018566526485061799702356193897017860040811889729918311021171229845901641921068884387121855646124960798722908519296819372388642614839657382291123125024186649353143970137428531926649875337218940694281434118520158014123344828015051399694290153483077644569099073152433278288269864602789864321139083506217095002597389863554277196742822248757586765752344220207573630569498825087968928162753848863396909959826280956121450994871701244516461260379029309120889086942028510640182154399457156805941872748998094254742173582401063677404595741785160829230135358081840096996372524230560855903700624271243416909004153690105933983835777939410970027753472000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

ピンクノイズを作るPython関数

Voss法でピンクノイズを作るプログラムをMatlabからPythonに移植しました。ベクトルではなくループで書いているところなんて、Python的・numpy的にも汚い書き方だと思うので、適当に修正して使って下さい。ピンクノイズについては「DSP Generation of Pink Noise」が詳しく、下記プログラムもこのページを全面的に参考にさせていただいています。

Audacityで10秒間のピンクノイズのスペクトルを見てみたところ、周波数特性もほぼ-10 dB/decade（≈ -3 dB/octave）になってるみたいですし、まぁまぁなんじゃないでしょうか。

Kindleの進化

僕が持っているKindleは、アメリカ国外にも初めて販売開始したKindle 2 Internationalというもの。100カ国ぐらいで3G回線を使って電子書籍のダウンロードができるモデル。当時はiPadが出る前だったので「おおタブレットですか、かっこいい」などと言われたものです。英語しか読めない*1ですし、ウェブブラウザがSSLに対応していないようで見られないサイトも増えてきましたが、The Pragmatic Bookshelf (PragProg)やManning PublicationsやApressがいいかんじのコンピュータ関連の電子書籍を出してくれるので、いまだに現役で使っています。

で、この間のAmazonセール期間に、家族のために新しいKindleを買いました。なんと第8世代だそうです。隣り合わせで写真を撮ったら「ひいおじいちゃんと孫」のようだったので載せておきます。それでも画面サイズは一緒のままなのですね。でもコントラストは上がっていたし、反応速度も上がっていて、e-inkのにじみも少なく、再安価モデルでも僕には十分な品質に見えました。（USB接続してPCからファイル転送できるのかはいずれ試さないと！）

*1:ファームウェアをいじると日本語も表示できますがフォントが気に入らないので英語版として使っています。

昔のウェブブラウザで

昔のウェブブラウザでURL欄に「about:version」と入力すると、そのブラウザのバージョン情報が表示されるとかいうのがあった気がするのだけど、Safariでそれをやっても検索結果が出てしまう。「about:blank」は動いたが、「about:plugins」とかもあった気がする…。気のせいだろうか、それともNCSA MosaicとかNetscapeとかの時代のものなのか？

Swiftを試してみた（簡単な分布関数）

昨日の続きで（チラシの裏のように）書いてみています。今日は二項分布と正規分布。あとは、累積分布関数の逆関数（norminvとbinoinv）に加えてt分布とF分布があると日常的な統計処理には使えそう。（たぶん仕事はSwiftでなくRでやると思うけど）

関数名はMatlabと合わせてあります。

Swiftを試してみた（簡単な記述統計関数）

新しく使う言語を試すときには音響プログラムを書いてみたり、統計プログラムを書いてみたりします。Pythonが簡単で面白かったので、Swiftはどんなだろうと試してみます。今回はSwiftで簡単な統計関数を作ってみました。と言っても平均や分散程度ですし、エラーチェックなども入っていないお遊びコードです。これを参考にはしないで下さいね。（実行環境：Xcode 8.2.1、Swift 3.0.2）

Twitterのような近況報告

1時間だけの数学の講義資料を作るのに10時間くらいかかった。効率が低いままだから、いつまでたっても週末も仕事漬けなのだ。とは言っても昨日は学会会場までバイク乗ったし、今日は10kmほどランニングしたし、仕事だけではないが。これからOblivionを30分ほどやる。（132字）