fortran66のブログ

fortran について書きます。

【メモ帳】ニュースなど

多倍長計算ルーチン GNU MPFR の Fortran binding

GNU MPFR Library

github.com

Fortran 次期規格の Generics について

The State of Fortran Genericseverythingfunctional.wordpress.com

考えているのは単なる型引数でしょうかね。

ジョセフ・ナイ   ウクライナ戦争に見る8つの教訓

Fifth, information warfare makes a difference. As RAND’s John Arquilla pointed out two decades ago, the outcomes of modern warfare depend not only on whose army wins, but also on “whose story wins.”

宣伝戦の重要性は第一次大戦後に欧州視察に出かけた近衛文麿も述べていたが、その後をみると嘘のスケールが小さすぎる。

www.project-syndicate.org

mobile.twitter.com

【メモ帳】メビウス関数 μ(n) その2

計算法

以前メビウス関数の記事を書きましたが、Fortran での計算法を別のやり方でやってみます。

fortran66.hatenablog.com

ここでは素数の因数ごとに符号反転、素数の二乗以上のべきが混じるものは 0 にするという、愚直な方法でまず一覧表を作ってそれを利用することにします。

プログラム

算法:エラトステネスの篩で n 以下の素数の表を作り、これに基づき素数ごとにその倍数のメビウス関数の表の符号を反転し、さらに二乗三乗・・・の倍数のメビウス関数の表を 0 にクリアしてゆきます。

求まったメビウス関数の表を元に、メビウスの反転公式を使って素朴に計算したゼータ関数の逆数を 2,3 の場合に計算してみます。

    module moebius_m
        implicit none
    contains
        function iprime(n) result(ires) ! primes under n
            integer, intent(in) :: n
            integer, allocatable :: ires(:)
            integer :: itab(n), i
            itab(1) = 0
            forall(i = 2:n) itab(i) = i
            do i = 2, int(n**0.5)    
                if (itab(i) /= 0) itab(i**2::i) = 0
            end do
            ires = pack(itab, itab /= 0)
        end function iprime 
        
        subroutine moebius_table(moebius)  
            integer, intent(out) :: moebius(:)
            integer, allocatable :: iprime_table(:)
            integer :: i, j, k, n
            n = size(moebius)
            iprime_table = iprime(n) 
            moebius = 1
            do i = 1, size(iprime_table)
                k = iprime_table(i)
                moebius(k::k) = -moebius(k::k)
                j = k**2
                do while(j <= n)
                    moebius(j::j) = 0 
                    j = j * k
                end do
            end do
        end subroutine moebius_table
    end module moebius_m
    

    program MoebiusP
        use moebius_m
        implicit none
        real, parameter :: pi = 4 * atan(1.0)
        integer, allocatable :: m(:)
        integer :: i, n
        real :: x, s
        integer :: mu(10**5)
        call moebius_table(mu)

        
   ! zeta(2) = pi^2/6
        x = 0.0
        s = 0.0
        do i = 1, 10000
            x = x + mu(i) / real(i)**2
            s = s + 1.0   / real(i)**2
        end do
        print *, x, 6.0 / pi**2, s, pi**2 / 6
        print *, 'zeta(2)^-1 * zeta(2) =', x * s
        print *
        
    ! zeta(3) = Apery's constant
        x = 0.0
        s = 0.0
        do i = 1, 1000
            x = x + mu(i) / real(i)**3
            s = s + 1.0   / real(i)**3
        end do
        print *, x, 1.0 / 1.20205690, s, 1.20205690
        print *, 'zeta(3)^-1 * zeta(3) =', x * s
        
    end program MoebiusP

計算結果

以前の結果を再現しています。一気にテーブルを作るので前より計算結果が早く出てきます。

  0.6079294      0.6079271       1.644725       1.644934
 zeta(2)^-1 * zeta(2) =  0.9998769

  0.8319076      0.8319074       1.202051       1.202057
 zeta(3)^-1 * zeta(3) =  0.9999951

オンライン整数列辞典 oeis.org

【メモ帳】exp(-x^4) の積分

(1/4)!

Youtube のおすすめで (1/4)! を計算する動画が出てきたのですが、 (1/4)! は  \int_0^\infty \exp(-x^4)dx に等しいので見てみたところレムニスケート周率で表わされることが学べたので、メモっておくことにします。


www.youtube.com

exp(-x4) の積分

 y=x^4 で置換すると  dy = 4x^3dx なので、
\int_0^\infty \exp(-x^4)dx=\int_0^\infty \exp(-y)y^{-{3\over4}}dy/4={1\over4}\int_0^\infty \exp(-y)y^{{1\over4}-1}dy={1\over4}\Gamma({1\over4})
が示される。一般に  ({1\over n})!=\int_0^\infty \exp(-x^n)dx となる。

fortran66.hatenablog.com

レムニスケート周率 \varpi

定義により、 \varpi=2 \int_0^1{dy\over\sqrt{1-y^4}}。その値は、

2.622057554292119810464839589891119413682754951431623162816821703...

oeis.org

Youtube での導出

  • Gamma 関数の相反公式より

 \Gamma({1\over4}) \Gamma({1-1\over4})=\Gamma({1\over4}) \Gamma({3\over4})={\pi\over\sin({\pi\over4})}

  • Beta 関数の性質から

B(1/4,1/2)={\Gamma({1\over4})\Gamma({1\over2})\over\Gamma({1\over4}+{1\over2})}={\Gamma({1\over4})\sqrt{2}\over\Gamma({3\over4})}

  • Beta 関数の定義から

ここでも  y=x^4 で置換する。

B(1/4,1/2)=\int_0^1y^{{1\over4}-1}(1-y)^{{1\over2}-1}dy=\int_0^1x^{-3}(1-x^4)^{-{1\over2}}4x^3dx=4\int_0^1{1\over\sqrt{1-x^4}}dx=2\varpi

以上のことから、

B(1/4,1/2)={\Gamma({1\over4})^2\sqrt{\pi}\over\sqrt{2}\pi}={\Gamma({1\over4})^2\over\sqrt{2\pi}}

つまり


\Gamma({1\over4})^2=2\varpi\sqrt{2\pi}

したがって


({1\over4})!={1\over4}\Gamma({1\over4})={1\over4}\sqrt{2\varpi\sqrt{}2\pi}

結論


\int_0^\infty \exp(-x^4)dx=({1\over4})!={1\over4}\sqrt{2\varpi\sqrt{}2\pi}

寝言

\Gamma({1\over4})^2 が出てくるので、\int_0^\infty \exp(-x^2)dx の時のように、積分を二乗して球座標に行って解けないかとも思いましたが、今日は方角が悪くてやる気が出ません!Bata 関数の出方を見る限り難しい気もします。

Fortran で検算

[tex: \exp(-x4)] の積分は台形積分で、急速に減少する関数なので上限を 2 までにします。分割数は 32 。

    program Lemniscate
        implicit none
        real, parameter :: pi = 4 * atan(1.0)
        real, parameter :: aLemniscate = 2.622057554292119810
        integer, parameter :: n = 32
        integer :: i
        real :: x(0:n), y(0:n), s, xmax = 2.0, h 
        h = xmax / n
        
        print *, '1/4gamma(1/4)=              ', gamma(5.0/4.0), gamma(1.0/4.0)/4.0
        
        forall(i = 0:n) x(i) = h * i 
        y = exp(-x**4)
        
        s = (y(0) + y(n)) / 2 + sum(y(1:n-1)) ! Trapezoidal rule
        
        print *, '\int_0^\infty \exp(-x^4)dx= ', h * s

        print *, '1/4sqrt(2Lsqrt(2pi))=       ', sqrt(2 * aLemniscate * sqrt(2 * pi)) / 4
    end program Lemniscate

実行結果

 1/4gamma(1/4)=                0.9064025      0.9064025
 \int_0^\infty \exp(-x^4)dx=   0.9064025
 1/4sqrt(2Lsqrt(2pi))=         0.9064025

【メモ帳】実数フォーマット戦国時代

48 bit への回帰

昔の CDC の実数は 48 bit 長だったと思います。

www.hpcwire.com

IBM も 360 シリーズに統一される前は、現在の 8 bit = 1byte ではなく、1 文字 6 bit で科学技術用は単精度浮動小数点数が 1 語 36 bit 2 進型で精度が今より良くなっていました。その後 360 シリーズで 32 bit 16 進 単精度になって精度がた落ちで倍精度必須的な感じになり、IEEE754 で 2 進 32 bit に統一されて現在に至っています。

モリー移動の律速化と、AI 用途での必要精度の低下で、歴史が循環しているようです。   

Microsoft仮数部を束に

そろそろ初代 Thinking Machine ばりに 1 bit 演算機で任意長計算するのが復活してもいい頃です。

悲しみ

安倍ちゃんのテロによる暗殺がショックでやる気が起こりません。

日本のサヨクやメディアはこの期に及んで安倍元首相への誹謗中傷捏造歪曲に励んでいてがっかりします。

https://sa.kapamilya.com/absnews/abscbnnews/media/2019/reuters/06/28/20190628-putin-brings-own-mug-rtr.jpg

news.abs-cbn.com

【黙祷】安倍元首相テロリストに暗殺さる!

黙祷

これほどの人物の命を、大切に守ってゆこうとしない社会にショックです。残念でなりません。

妬み僻み嫉みで足の引っ張り合いしかしないサヨク・メディアが憎しみを煽ったことが巡り巡った結果だと思います。

疫病、戦争、暗殺、次は飢餓でしょうか?グローバリゼーションや中国の膨張とともに、理性が衰退し時代が逆行している気がします。歴史を振り返るに、文明の崩壊は気候変動を背景とした民族大移動で引き起こされているので、嫌な感じしかしません。

日本の周囲には、人権抑圧的独裁国家や古代奴隷制のような国家、テロリスト賛美の犯罪結社のような国ともいえないような醜い集団などが多数存在しており、国内にもこれらに同調する人物がいます。

その様な環境下でも明るく広い視野を失わない安倍元首相のような政治家は貴重な存在でした。安倍元首相の命が失われたのは本当に大きな損失だと思います。とても残念で悲しいです。

【寝言】ニュースなど

COBOL/Fortran が Latin 語扱いw

www.nytimes.com

ウクライナ侵攻後の Yandex に関する記事

Nginx はどうなのか知りたい。

www.nytimes.com

BL エロ同人界隈

イギリス政界の中心地、ロンドン・ウェストミンスターにある会員制クラブ「カールトン」での出来事が、辞任理由だった。

しかし実情は、酒を飲みすぎただけではなかった。議員は、この会員制クラブで男性2人につかみかかり、少なくとも1人の局部に触れたとされている。

ピンチャー氏に対しては同じような性的問題行為の指摘が、過去にも繰り返されていた。それでもジョンソン首相は、副幹事長の職を与えた。副幹事長とは、下院議員の間の規律維持を担うほか、議員たちの私生活の相談にも乗る立場だ。

この後、首相官邸は報道陣に、ジョンソン氏は確かに知っていたが「忘れていた」のだと説明した。

www.bbc.com

ピンチャーだけに、摘まみたくなったのでしょう。名は体を表す。

学士会報 第九五五号

サヨク臭漂う記事の多い笑える号でした。理系記事では乳がん検診で散乱問題を解くのにおっぱいを円錐で近似しているのが印象的でした。

マルクスエコロジーに関心があって、持続可能なエコ社会主義を構想していたことを発見したという記事か冒頭に来ていて、茶を吹き出しそうになりました。10年くらい前に耶蘇の聖書の単語を十個置きに読むと未来の予言が書いてあるとか、そういう本が全米大ヒットしていた気がしますが、同レベルの滑稽さです。エルビス・プレスリーもまだ生きていそうです。

ただ SDGs は大衆の阿片であるとか、サヨちゃんの脳内汁が垣間見えて面白い所もあります。SDGs は「グリーン資本主義」の免罪符と罵っています。で、結論が「脱成長」なのですが、それは日本がここしばらく成し遂げている事なのですが、なぜか日本は二周遅れだ!と高説を垂れて呉れます。ギャグマンガよりも面白いです。

  • 日本の目指すべき新しい資本主義

読むだけ無駄な SDGs 系の寝言です。

グローバル資本主義プーチン大統領が冷戦構造を復活させてくれることでようやく終わりそうです。仕上げに中国を完全に封じ込めた上で、特別蠱毒作戦で中露で戦争をさせるのが良かろうと思います。支那人の性格からしてロシアが弱体化すれば必ず弱みにつけ込んで攻め込んで返り討ちに会って負けると思われます。

  • ロシアは何をめぐってウクライナ・米欧と対立しているのか

最近の米メディアなどではNATOの東方拡大問題はブッシュ政権が悪い説が多く、親子のどっちなのか判然としませんが、ブッシュ(父)政権のベーカー国務長官が拡大しないと口約束をした説と、ブッシュ(子)政権のときジョージア(旧グルジア)のNATO加盟による拡大を提起したのがロシアを刺激した説があって、親子でリバーシブルのようです。

この記事ではブッシュ(子)が悪いことになっています。一方、ウクライナ侵攻前の下斗米伸夫の記事「プーチン・ロシア考」には、クリントンが選挙の東欧移民票目当てで東方不拡大合意を反故にしたと書いてあります。リベラル・メディアは女性部下へのセクハラ+パワハラおじさんビル・クリントンが大好きなのでブッシュが悪いと言いそうですが、ソ連ロシアの専門家の言うことも少しも信用ならんので、何がどうなのか全く見当がつきません。

ともあれ日ソ中立条約を破ったロシアが約束にこだわるのもはなはだ滑稽で笑えます。

読者のお便り欄では、米軍による空襲の思い出話がでていて、戦闘機から無差別機銃掃射受けたことなどが書かれています。そろそろ昭和百年ですねw

【寝言】Wendy Carlos の動画 (旧 Walter Carlos)

Wendy Carlos demonstrates her Moog Synthesizer in 1970

昔、Walter Carlos がいつの間にか Wendy Carlos に性転換していてびっくりしましたが、1970 年の時点で十分カマっぽいです。


www.youtube.com

89 年には限界突破していたようです。


www.youtube.com

Switched on Bach


www.youtube.com