2007-01-01から1年間の記事一覧

関数の速度(7)

並列処理を使うと速くなることがある。 具体的にはSSEを使うと速くなることがある。 http://www.intel.co.jp/jp/download/index.htm SSEというのは、floatを4つ並べて一つの命令でまとめて計算する技術である。例えば、 x[0] /= y[0]; x[1] /= y[1]; x[2] /=…

関数の速度(6)

atan2 とにかくこれを速くしたい。 場合分けはあとまわしにして、まずは、0 前回のatanがそのまま使えて、 y /= x; double y2 = y * y; return y * (((256 * y2 + 5943) * y2 + 19250) * y2 + 15015) / (((1225 * y2 + 11025) * y2 + 24255) * y2 + 15015);y…

関数の速度(5)

前回求めた有理関数を展開すると、 h1とh2が近そうだが、実際に[1,2]で最大誤差を見てみると、 0 : error = 0.00520393 1 : error = 0.000249046 2 : error = 0.000186153 3 : error = 0.00129726 4 : error = 0.0901862 と、h2がわずかによい。n = 5だと、 …

関数の速度(4)

log(1+x)/xを有理関数で近似するとき、テーラー展開で4次まであわせると、 となったが、4次まであわせるには、 という形などが考えられる。 その中から最適なものを選ぶべきではないだろうか。 一般に、近似したい関数を、 として、n次までの近似を求めたい…

何回コインを投げれば連続10回表が出るんだろう(4)

分布はだいたい分かったので、これで擬似乱数を調べてみよう。 100万回行って、得られた回数を100ごとに区切って、度数をみる。 VC2005EEで、rand()で得られた値の特定のビットを拾う。 if(rand() & bit) { ...bitが1の結果は、 〜100 56000 〜200 24000 〜3…

何回コインを投げれば連続10回表が出るんだろう(3)

もう一度数学的に考えてみよう。 p(m,n) = 2-n-1(1 - p(n,n) - ... - p(m-n-1,n)) がm > 2nで成り立つ。 確率の総和が1であるから、 p(m,n) = 2-n-1(p(m-n,n) + ...) これをすべて並べて、 p(2n+1,n) = 2-n-1(p(n+1,n) + p(n+2,n) + ...) p(2n+2,n) = 2-n-1(…

何回コインを投げれば連続10回表が出るんだろう(2)

この計算は、浮動小数点数だと誤差がたまってくるが、幸い2進数では有限小数になるので、ある程度の桁までは正確に計算できる。 Arrayで、2048進数で考えるとプログラミングがやさしくなる。 例えば、3 + 11/2048 + 13/2048^2 → [ 3, 11, 13 ]とする。 10万…

何回コインを投げれば連続10回表が出るんだろう(1)

http://blogs.wankuma.com/episteme/archive/2007/12/13/112757.aspx この分布を擬似乱数で出すのは難しい、擬似乱数の質を調べるのによい題材、ということだと思う。 しかし、それにはまず正しい分布を知ることが必要。 というか、私にはそちらを出すほうが…

関数の速度(3)

前々回、expがsin/cosより速いのが気になった。 また、logがatanより速いのも気になった。 exp(x) 35ns sin(x) 45ns cos(x) 45ns log(x) 43ns atan(x) 87ns 数学的には、expとsin/cosは同類である。 logとatanもテーラー展開すれば似たようなものである。 そ…

関数の速度(2)

前回、 double x = 0.5; for(int i = 0; i sum += exp(x); x += dx; // dx = 1e-8 }とした、すなわち、 exp(0.5) + exp(0.50000001) + exp(0.50000002) + ... を計算しているのは、確実にexpをN回計算するため。 これを、 double x = 0.5; for(int i = 0; i …

関数の速度(1)

私はいつも実行速度と戦っている。 遅さの原因がどこにあるかと探ると、 そこにはたいてい atan2 がある。 テーブルを作って速くしたりしているが、劇的には速くならない。 やはり肝心なのはatan2を使わないこと。 atan2の使用を避ける方法をいつも考えてい…

JavaScriptのStringは値?

http://blog.livedoor.jp/dankogai/archives/50963275.htmlえー! http://d.hatena.ne.jp/odz/20071206/1196927108ここにあるようにReadOnlyなだけだと思ってた。 var str = "a"; var str2 = str; // アドレスをコピー var str3 = str.bold(); // 実体をコピ…

フィボナッチ数列の計算法(3)

前に、http://blog.livedoor.jp/dankogai/archives/50958771.htmlを見て、フィボナッチ数の計算がO(1)なのはちょっと気持ち悪いがいいんじゃない?と思ったが、この続きのような記事を見て、ハッとした。http://blog.livedoor.jp/dankogai/archives/50962361…

フィボナッチ数列の計算法(2)

前回の宿題。 べき乗計算に指数の計算を使う(powを使う)場合、80ビット浮動小数点数を使うと、フィボナッチ数列を何項目まで正しく求められるだろうか。 D言語でIntelの場合、realは80ビット浮動小数点数となる。これを用いた計算と、64ビット整数で前回の…

フィボナッチ数列の計算法(1)

http://blog.livedoor.jp/dankogai/archives/50958771.html O(1)の計算方法が紹介されているが、roundをつかうのはちょっと気持ち悪い(64ビット整数を求めるのに、80ビットの浮動小数点数の計算では追いつかない?)という人のために、O(logn)で計算する方…

zzの逆関数(5)

一般の c に対して、 を解く。 ここで、 とおく。 ただし、aは虚部が実部の絶対値より大きく、 aの絶対値が1より大きくなるように虚部を調整しておく。 そうすると、 以前と同じようなグラフが描けて、同じように解ける。 以下、ソース。 ただし、あまり検証…

zzの逆関数(4)

を解く。 同じように、 を解けばよい。 すなわち、の、となる点を取ればよい。 は、同じグラフで、となる点を取ればよい。

zzの逆関数(3)

前回、グラフで解の存在だけ示したが、 それだけではさみしいので、具体的数値を求めてみよう。 グラフによると、wが0に近いところに解があるので、 それを求めてみる。 下より、wは小さいので、 これを上に代入して、 Xはそれなりに大きいので、log(2πn)と…

zzの逆関数(2)

まず、zz = 1 の解を考える。 とする(x, y, s, tは実数)。 より、 nは整数。 nが0なら、uが0で、zは1となる。 nは正とする。負なら共役複素数を考えればよい。 再び、 とおくと(v, wは実数)、 最後は多価だが、これだけ考えればよい。 これから、 これの…

zzの逆関数(1)

前に読んだこの本の巻末で、 zzの逆関数について言及されていた。xのx乗のはなし (はじめよう数学)作者: 土基善文,上野健爾,高橋陽一郎,浪川幸彦出版社/メーカー: 日本評論社発売日: 2002/07メディア: 単行本この商品を含むブログ (3件) を見る考察はされて…

複素数(11)

あとのためにまとめておこう。 D言語での複素数 あれこれ http://d.hatena.ne.jp/inamori/20071005/p1exp http://d.hatena.ne.jp/inamori/20071008/p1cos/sin/tan http://d.hatena.ne.jp/inamori/20071009/p1log http://d.hatena.ne.jp/inamori/20071011/p1s…

複素数(10)

cdouble exp(cdouble z)などを実装してきたが、 cdouble以外のcfloat、crealも同じコードで済む。 何度も同じコードを書きたくないとき、 D言語ではmixinを使うらしい。 import std.cstream;void main(char[][] args) { dout.writefln(pow(2.0f, 5)); dout.w…

複素数(9)

純虚数の場合 引数が純虚数だと速くなる可能性があるので、 それ用にも書いてみる。 ここまでは一価だからやさしい。 log 虚部は[-π, π]だから、 sqrt 符号はyの正負による。 acos 実部が[0, π]だから、-を取る。 asin 実部が[-π/2, π/2]だから、 -を取ると…

複素数(8)

べき乗 powの値域は非常に複雑。 だが、log zは2πiの周期を持つから、alog zは2πaiの周期になる。 aによって3種類に分かれる。 1)aが整数の場合 一価関数となる。 2)aが整数でない有理数の場合 a = m / nと表すと(m, nは互いに素、n > 0)、 n価関数となる…

複素数(7)

atan より、 ただし、logが2π周期だから、 値は、実数がπ周期。C99では、http://docs.hp.com/ja/B2355-60104-06/catan.3M.html値域は、[-π/2, π/2]。 素直に計算すればよい。 import std.cstream; import std.math;void main(char[][] args) { dout.writefln…

複素数(6)

acos をzとおくと、eiZの2次方程式となり、 これを解くと、 より、 ただし、logが2π周期で、頭に±がついている。 主値は、C++では定義されていない? C99ではあるらしい。http://docs.hp.com/ja/B2355-60104-06/cacos.3M.html値域は、実部が[0, π]だというこ…

複素数(5)

atan2 前回のlogで、 log(-1+0i); // 0+3.14159i log(-1-0.0i); // 0+3.14159iとなる。これは、atan2が次のようになるためである。 atan2(0.0, -1); // 3.14159 atan2(-0.0, -1); // -3.14159 sqrt z = reiθ とすると、 sqrt(z) = sqrt(r)eiθ/2 となる。 主…

複素数(4)

対数関数 指数関数の逆関数だから、 とおいて、 atan2は便宜的に使った。 だから、 ところが、atan2は周期2πの多価関数だから、logも同様となる。 例えば、を考えると、 となるが、同時に、 となるため、は多価関数である。 プログラムとして、 返り値が複数…

複素数(3)

三角関数 三角関数は指数関数とは姉妹のようなものである。 オイラーの公式より、 これらから、 z = x + iy を代入して、 sin, tanは、 これらより、 import std.cstream; import std.math;void main(char[][] args) { cdouble c = 1 + 0.1i; dout.writefln(…

複素数(2)

指数関数 これはやさしい。 オイラーの公式、 より、 import std.cstream; import std.math;void main(char[][] args) { cdouble c = 1 + 0.1i; dout.writefln(exp(c)); // 2.7047+0.271375i }cdouble exp(cdouble c) { double im = c.im; return std.math.e…