fortran66のブログ

fortran について書きます。

【メモ帳】Stirling の公式

ガンマ関数

Stirling の公式  n! = \sqrt{2\pi n}({n\over e})^n をまた考えてみます。 \Gamma(n+1)=n! で、また \Gamma(x)=\int_0^\infty x^{n-1}\exp(-x) dx であるので、  n! =\int_0^\infty x^{n}\exp(-x) dx を考えることにします。

この積分被積分関数x が大きくなると増大する x^n とそれを減衰させて抑え込む \exp(-x) の積になっていて、x の増加につれてまず増えた後に尻すぼみに減る山のような形をしています。山の頂上が来るのは微係数が 0 になる点なので  n x^{n-1}\exp(-x)-x^n\exp(-x)=0 より  x=n の所と分かります。これより山の頂上の高さは n^n\exp(-n) と求まります。これの対数を取ると  n! \simeq n\log(n)-n となり、この時点で階乗の対数の近似によく使われる近似が得られています。

ここで山の形を Gaussian で近似することを考えます。山の高さは分かっているので、山の幅を決めればいいことになります。なぜ Gaussian なのか、Lorentzian じゃダメなのか!と言いたくなるでしょうが、図を書いてみると裾がストンと落ちているので、Gaussian かなという気もします。しかしそれよりも、Stirling の公式の  \sqrt(2\pi n) という形が Gaussian の積分を思わせます。もう少し真面目に考えると、被積分関数を鞍点法のノリで  \exp の形にまとめた後 \exp の中身を x=n で Taylor 展開してみると、定数項は山の高さに、1回微分は山頂だから 0 で消えて、二次の項が残って Gaussian の形になるので、この項で切ればいい気がします。

jupyter wolfram script にて

式の形からして、左右対称な釣り鐘型にはならんだろ Gaussian じゃダメなんでは?という猜疑心の深い方のために、3次の項まで Taylor 展開して左右対称からの崩れも見てみることにします。

f[x_] := x^n Exp[-x]
g[x_] := Log[f[x]]
Series[g[x], {x, n, 3}]

Exp[Normal[%]]

n = 50
g2[x_]:= -n -(x-n)^2 / (2 n)
Plot[{x^n * Exp[-x], n^n * Exp[ g2[x] ]}, {x, 0, 2 * n}]

g3[x_]:= -n -(x-n)^2 / (2 n) + (x-n)^3 / (3 n^2)
Plot[{x^n * Exp[-x], n^n * Exp[ g3[x] ]}, {x, 0, 2 * n}]

青が被積分関数、橙が近似によるもの。

以上の結果から、山の形を Gaussian で近似して面積を求めたのが Stirling の公式ということが分かります。Gaussian の幅は \sigma^2=2n となります。また 3次の項までとるとGaussian の左右対称からのズレがおおよそ修正されることも分かります。

【メモ帳】バッカス講演

Programming in America in the 1950s – Some Personal Impressions

同名の記事について、以前記事にしていました。 fortran66.hatenablog.com

This lecture’s transcript was included in the edited volume from the conference, viz. Backus, J., “Programming in America in the 1950s – Some Personal Impressions,” in Metropolis, N., and Howlett, J., Rota, Gian-Carlo, A History of Computing in the Twentieth Century, New York: Academic Press, 1980, pp. 125 – 135.

www.semanticscholar.org


www.youtube.com

thenewstack.io

Software という単語の起源

質疑(31:48~)の所で、software と firmware という単語の発明者が、Ascher Opler だと聴衆中の人が言っています。 1960年頃の Rand Corporation での集まりで出たようです。確かではないが Software という言葉は 1960 年に出来た言葉だろうとのことです。

ただ wikipedia をみますと、firmware は Ascher Opler で正しいようですが、software の方は 1953 年の Rand Corporation Research Memorandum に出てくるのが工学的な文脈で確認できる最古の出版例だそうです。こうしてみると Rand Corporation 界隈で出てきたもののような気がします。

en.wikipedia.org

en.wikipedia.org

Paul Niquette が起源を主張しているようです。

www.niquette.com

www.nytimes.com

訃報に訂正が出ていて Turkey 説は正しくないようです。但し、M.J. Lorenzo 著の以下の本によると、単語そのものの初出ではないが 1958 年に書かれた記事の中で今日的な意味で software という単語を使っているようです。

【寝言】月齢1.4

新しい月

六月三十日の夕方、夕焼け雲が綺麗なので写真を撮ったのですが、後から拡大して見てみると月齢 1 台の月が写っていました。ネットの月齢サイトで調べてみると月齢 1.3~1.4 程度と思われます。

こんなに細い月を見るのはしばらくぶりです。うれしいです。1年の折り返しに当たり、天皇陛下のお恵みか!w

【寝言】Moonglow 

Moonglow / thema from picnic

picnic の映画は見てないんですが曲が好きで。メドレーと書かれたりもしますが、今風に言うとマッシュアップのようになっていて、同時に二つの曲の旋律が流れるところが良かったりします。

この前、月が綺麗に輝いていたし。

夏の夜のようなけだるい感じがいいんですが、イメージを反転させたチャチャチャ版もいい。


www.youtube.com


www.youtube.com


www.youtube.com


www.youtube.com


www.youtube.com


www.youtube.com


www.youtube.com

【メモ帳】arcsin

導出メモ

すぐ忘れるw

\displaystyle \sin(x)={\exp(ix)-\exp(-ix)\over 2i}

\displaystyle 2i\exp(ix)\sin(x)=\exp(2ix)-1

\displaystyle \exp(2ix)-2i\exp(ix)\sin(x)-1=0

これを  \exp(ix)二次方程式と見て解くと

\displaystyle \exp(ix)=i\sin(x)+\sqrt{-\sin^2(x)+1}

但しプラスの解だけを取った。ここで両辺の対数を取ると、

\displaystyle ix=\log\left(i\sin(x)+\sqrt{1-\sin^2(x)}\right)

となる。

さて  \sin(x)=z と置き換えることを考えると  x=\arcsin(z) なので、

\displaystyle \arcsin(z) =-i\log(iz+\sqrt{1-z^2})

が求まる。

check

    program arcsine
        implicit none
        real :: x
        complex :: c, ci = (0.0, 1.0)
        integer :: i
        do i = 1, 10
            call random_number(x)
            print *, x, ':', asin(x), -ci * log(ci*x+sqrt(1.0-x*x))
        end do
        
        print *
        
        do i = 1, 10
            c = cmplx(i*i*i, i*i) * 1.0e-2
            print *, c, ':', asin(c), -ci * log(ci*c+sqrt(1.0-c*c))
        end do    
    end program arcsine
 3.9208680E-07 :  3.9208680E-07 (3.9208680E-07,-7.6866031E-14)
 2.5480442E-02 :  2.5483200E-02 (2.5483200E-02,-1.2659941E-08)
 0.3525161     :  0.3602585     (0.3602585,-5.6824399E-09)
 0.6669145     :  0.7300602     (0.7300602,-8.5801748E-09)
 0.9630555     :   1.298127     (1.298127,-7.1473925E-09)
 0.8382882     :  0.9941360     (0.9941360,-6.4942043E-09)
 0.3353550     :  0.3419821     (0.3419821,-3.6983820E-09)
 0.9153272     :   1.156320     (1.156320,2.0078712E-09)
 0.7958636     :  0.9204326     (0.9204326,-2.7873030E-08)
 0.8326931     :  0.9839536     (0.9839535,-2.3087523E-09)

(9.9999998E-03,9.9999998E-03) : (9.9996664E-03,1.0000333E-02) (9.9996664E-03,1.0000329E-02)
(7.9999998E-02,3.9999999E-02) : (8.0021039E-02,4.0117647E-02) (8.0021039E-02,4.0117655E-02)
(0.2700000,8.9999996E-02) : (0.2721770,9.3304269E-02) (0.2721770,9.3304262E-02)
(0.6400000,0.1600000) : (0.6775884,0.2039517) (0.6775884,0.2039518)
(1.250000,0.2500000) : (1.270154,0.7668216) (1.270154,0.7668217)
(2.160000,0.3600000) : (1.385768,1.424145) (1.385768,1.424145)
(3.430000,0.4900000) : (1.422694,1.915248) (1.422694,1.915249)
(5.120000,0.6400000) : (1.444065,2.324834) (1.444065,2.324834)
(7.290000,0.8100000) : (1.459105,2.681223) (1.459105,2.681221)
(10.00000,1.000000) : (1.470634,2.998273) (1.470634,2.998266)

【寝言】夏至

至点 


www.youtube.com

Euler の γ

 \displaystyle \gamma := \lim_{n\rightarrow\infty}\left(\sum_{i=1}^n{1\over i} - \log(n)\right)

収束が遅い(~O(1/n) ) Euler の γ ですが、引くlog(n) を定義からずらして log(n+1/2) にすると 定義より早く O(1/24n2) で収束するのが草。

計算による比較

小さい項からの和

    program euler_gamma
        implicit none
        integer, parameter :: kd = kind(0.0d0)
        real(kd) :: s
        integer :: i, n = 10**7
        s = 0.0_kd
        do i = n, 1, -1
            s = s + 1.0_kd / i
        end do
        print *, s - log(real(n)), s - log(n + 0.5_kd), 1.0_kd / real(n)**2 / 24
   end  program euler_gamma    

カスを集めるのに小さい方から足しました。

定義通り、log(n+1/2)、誤差見積もり

  0.577215967910643       0.577215664901544       4.166666650965334E-016

組み手和(Kahan の和)

   block      
        real(kd) :: s, t, u, v 
        integer :: i, n = 10**7
        s = 0.0_kd
        v = 0.0_kd
        do i = 1, n
            u = 1.0_kd / i - v
            t = s + u
            v = (t - s) - u
            s = t
        end do
        print *, s - log(real(n)), s - log(n + 0.5_kd), 1.0_kd / real(n)**2 / 24
    end block    

カス項が効いているようなので、もう少し真面目な総和計算するとより良い値が得られました。

  0.577215967910632       0.577215664901534       4.166666650965334E-016

Wolfram language

N[EulerGamma, 20]
0.57721566490153286061

もう少し収束の速いもの

 \displaystyle \gamma = {1\over2}+\sum_{n=1}^{\infty}\sum_{m=2^{n-1}}^{2^n-1}{1\over2m(2m+1)(2m+2)}

    program euler_gamma
        implicit none
        integer, parameter :: kd = kind(0.0d0)
        integer :: i, j, n, m, m0, m1, nmax
        real(kd) :: s, t, x2
        nmax = 27
        s = 0.5_kd
        m1 = 1
        do n = 1, nmax
            t = 0.0_kd
            m0 = m1
            m1 = 2**n
            do m = m0, m1 - 1 
                x2 = 2 * m
                t = t + n / (x2 * (x2 + 1) * (x2 + 2))   
            end do 
            s = s + t
            print *, n, m1, s, t
            if (t < 1.0e-16_kd) exit
        end do    
        print *, 'gamma =', s
    end program euler_gamma
           1           2  0.541666666666667       4.166666666666666E-002
           2           4  0.564285714285714       2.261904761904762E-002
           3           8  0.572991591741592       8.705877455877456E-003
           4          16  0.575914184398886       2.922592657293993E-003
           5          32  0.576829154090275       9.149696913897526E-004
           6          64  0.577103770405349       2.746163150734348E-004
           7         128  0.577183875992416       8.010558706751096E-005
           8         256  0.577206763957738       2.288796532201553E-005
           9         512  0.577213201244027       6.437286288142040E-006
          10        1024  0.577214989382304       1.788138277449749E-006
          11        2048  0.577215481120550       4.917382461222898E-007
          12        4096  0.577215615230996       1.341104457486260E-007
          13        8192  0.577215651552576       3.632158007173262E-008
          14       16384  0.577215661331463       9.778887010694266E-009
          15       32768  0.577215663950808       2.619344739581377E-009
          16       65536  0.577215664649300       6.984919308599662E-010
          17      131072  0.577215664834837       1.855369191549276E-010
          18      262144  0.577215664883949       4.911271389529140E-011
          19      524288  0.577215664896910       1.296029950023478E-011
          20     1048576  0.577215664900320       3.410605131646578E-012
          21     2097152  0.577215664901215       8.952838470575934E-013
          22     4194304  0.577215664901450       2.344791028008255E-013
          23     8388608  0.577215664901511       6.128431095930439E-014
          24    16777216  0.577215664901527       1.598721155460254E-014
          25    33554432  0.577215664901531       4.163336342344230E-015
          26    67108864  0.577215664901532       1.082467449009531E-015
          27   134217728  0.577215664901533       2.810252031082286E-016
 gamma =  0.577215664901533

思ったより収束が早くないです。

【メモ帳】ロシア語で考えるんだ!

Fortran 講座 メモリー確保位置

モリーの割り付けについて詳しく説明されているようですが、ロシア語なのでよく分かりません。

よくみるとウクライナの方のようです。


www.youtube.com