2019-02-07

【ニュース】intel fortran v.19.0.2 出る

intel fortran v.19.0.2

software.intel.com

Changes in Update 2 (Intel® Fortran Compiler 19.0.2)
Intel® Fortran Compiler 19.0 Update 2 includes functional and security updates. Users should update to the latest version.

BBC Micro Bit を FORTRAN 77 で

youtu.be

Fortran ハンドブック

作者: 田口俊弘
出版社/メーカー: 技術評論社
発売日: 2015/07/10
メディア: 単行本（ソフトカバー）
この商品を含むブログを見る

Modern Fortran Explained: Incorporating Fortran 2018 (Numerical Mathematics and Scientific Computation)

作者: Michael Metcalf,John Reid,Malcolm Cohen
出版社/メーカー: Oxford Univ Pr
発売日: 2018/11/06
メディア: ペーパーバック
この商品を含むブログを見る

Fortran90/95による実践プログラミング

作者: 安田清和,水野正隆,小野英樹
出版社/メーカー: 大阪大学出版会
発売日: 2014/03/28
メディア: 単行本（ソフトカバー）
この商品を含むブログを見る

2019-02-04

【ニュース】メモ

Fortran 規格目安箱（サーベイ）：最終報告

英 BCS サイトにて
www.fortran.bcs.org

[01/19] Fortran Standardisation Benefits Survey - Final Report
The Final Report on the Fortran Standardisation Benefits Survey is now available. The Survey was developed by members of the Fortran specialist group committee and was open from 30 June 2018 to 3 January 2019.

http://www.fortran.bcs.org/2019/FortranBenefitsSurvey_final.pdf

ハッカーニュースに Fortran ネタ

news.ycombinator.com

ベンチマークで Rust 厨が湧くパターンｗいつも通りインライン・アセンブラで草

謎の音楽　Fortran 77

www.youtube.com

2019-01-30

【ニュース】ACM SIGPLAN Fortran Forum 2018年12月号出る！　その他

Compiler Support for the Fortran 2003/08/18 のみ！

ようやく出た ACM SIGPLAN Fortran Forum。

しかし記事は、Ian D. Chivers & Jane Sleightholme のコンパイラ対応状況のみ。

山里は　冬ぞ寂しさ　まさりける
　　　人目も草も　かれぬと思へば

とはいうものの、Fortran 2018 対応状況も加わって進歩を感じます。ARM, Fujitsu, NEC も加わって頑張っています。

dl.acm.org

サーバーレス Fortran

よく分かりませんが Fortran が動いたｗ

その後、それをAWSにアップロードし、CSVデータをアップロードした段階でトリガーが実行され、AWSのLambdaのファンクションとして、このFortranのアプリケーションが実行されたことを確認した。

thinkit.co.jp
f:id:fortran66:20190130011935j:plain

2019-01-18

【新刊】fortran本【予定】

CRC Press Fortran 新刊予定

逃げ水のように、出版予定年月日が伸びてゆく CRC press の Fotran 本。蜃気楼のようにサイトのリストに現れたり削られたりします。

ブレグジット失敗でのポンド暴落を待っているのですが、中々しぶとい。ドイツ銀行株空売りで EU も潰しておきたいところ。

Damian Rouson 本

Modern Fortran: Tricks of the Software Engineering Trade - CRC Press Book

トルコ行進曲

Computational Methods in Science and Engineering: Models, Algorithms, Coding, and Analysis with GNU Fortran - CRC Press Book

トルコと言えば、なぜかギュレン師系の学校が日本にもそこそこあるのが笑える。

その他　HPC 本リスト

CRC Press Online - Series: Chapman & Hall/CRC Computational Science

2019-01-17

【乞食速報】Springer 数学電子本 3,799円引き

新年乞食

数学・統計学のイーブックに有効な3,799円のギフトクーポン！ | Springer

数学・統計学のイーブックに有効な3,799円のギフトクーポン！
キャンペーン期間は2019年1月31日までです*
12,000以上ものタイトルをご用意
クーポンコードMATH19E を入力後ショッピングカートにて適用

ご注文金額は、他のディスカウントも適用された状態でご購入手続き時に表示される額が7,599円以上となる必要があります。

見逃していたｗ

2018-12-31

【メモ帳】反切の起源

狩野直喜「支那 文学史」

支那文学史―上古より六朝まで (1970年)

作者: 狩野直喜
出版社/メーカー: みすず書房
発売日: 1970
メディア: ?
この商品を含むブログを見る

休みなので読みさしだった狩野直喜の「支那文学史」の残りを斜め読みしていましたが、六朝文学の所はあんまり興味もないうえに、内容も人物・作品の列挙羅列だし、疲れました。

ただ当時の政治情勢が内容に反映してデカダン趣味に走ったという背景などが述べられていて、この様な因果付けは大変勉強になります。

第四篇「六朝文學」第七章　梁の文學　の冒頭で反切について書いてあったのが興味を引きました。

　それから前にも一寸述べた通り、声韻と云ふ事が此の時代に矢釜敷なつて来た。すなわち詩の意味よりも寧ろ之れを表はす文字の響きに関して、面倒な規則が設けらるゝ様になった。尤此等は決して偶然に起つた事にあらず、已に魏の孫炎反切を云へり。一体漢までは学者が或る字の音を示すには「讀若某」と云ふのみ。しかるに炎に至り初めて反切を創め、或字の音をあらはすに○○反とか○○切と云ひ、二字の音をして以つて一の音を出すの方法を起したり（顔子家訓音辞篇）。

（めんどくさいので漢字は現代文字でｗ）

漢字は表音文字ではないので、音をあらわすのに反切では、漢字二文字で前のは子音、後ろのは母音でローマ字みたいに読んでいますが、その起源が三国時代ににあるようです。仏教伝来で、天竺の文字を読むために発達したという説が上記引用のあとに紹介されています。

丁度、その反切が出来た魏の頃の魏志倭人伝に我が国の記録とともに『邪馬台国』が出てくるわけですが、「邪馬台」の子音をとれば「ヤマト」になるので、これは大和朝廷の音に過ぎないと思うのですが、邪馬台国には大して興味が無いので、どうでもいいですｗ

ただとりあえず、反切が魏の時代にはあったということで、邪馬台＝ヤマト音説が少なくとも矛盾しないことが分かってよろしかったです。

[追記：Ｒ元10.12]
武内義雄「支那学研究法」（武内義雄全集第九巻100頁）に、邪馬台の台（臺）は古音では韻書で「ト」と発音する部に入っているというようなことが書いてあって、普通にそのものずばり「ヤマト」読みしていました。

漢文研究法: 中国学入門講義 (東洋文庫)

作者: 狩野直禎,狩野直喜
出版社/メーカー: 平凡社
発売日: 2018/07/12
メディア: 単行本
この商品を含むブログを見る

中国哲学史 (1953年)

作者: 狩野直喜
出版社/メーカー: 岩波書店
発売日: 1953/12/05
メディア: 単行本
この商品を含むブログを見る

春秋研究

作者: 狩野直喜
出版社/メーカー: みすず書房
発売日: 1994/11/01
メディア: 単行本
この商品を含むブログを見る

2018-12-24

【メモ帳】ガンマ関数とスターリングの公式　

ガンマ関数

正の整数値に対してのみ定義されている階乗を、定義域を実数に広げて、微分可能な連続関数として拡張した物がガンマ関数 $\Gamma(n) = \int_0^\infty x^{n-1}e^{-x}dx$ ですが、階乗とは $n!=\Gamma(n+1)$ で結びついています。*1

ぱっと見、なんでこれが階乗になるのか意味不明な式ですが、 $n!$ の場合の被積分関数を少し見みてみることにして、これを $f(x)=x^ne^{-x}$ と置いておきます。

この被積分関数 $f(x)$ は、無限大に発散してゆく $x$ のべき $x^n$ と尻すぼみに減衰させる指数関数 $\exp(-x)$ の積になっています。べき関数よりも減衰関数の方が強いので、 $x$ を０から増やしてゆくと、被積分関数 $f(x)$ は初めは勢いよく増えますが、どこかで天井を迎えてその後０に向かって減衰してゆくことになります。値は常に正になっています。

その天井の x 座標は、 $f(x)$ の微係数が 0 になる所ですが、それは $f'(x) = e^{-x}x^{n-1}(n-x)$ なので $x = n$ と分かります。ここで、整数 $n$ がちゃんと形をとどめて出てくるので少し納得です。

fortran66.hatenablog.com

これが半無間区間で積分されるわけですが、どうせ被積分関数のピーク位置以外の寄与はカスみたいに小さいというのがパターンなので、今回もそんなところだろうと予想して、被積分関数の値 $f(n)=e^{-n}n^n$ を数値的に計算してガンマ関数と比較して見ることにします。

プログラム１

被積分関数のピーク値と積分した値に対応するガンマ関数の比を計算します。ガンマ関数は $(n-1)!$ と $n!$ 相当の二つの場合を採っておきます。 $e^{-n}n^n/\Gamma(n)$ と $e^{-n}n^n/\Gamma(n+1)$

    program gam1
      use, intrinsic :: iso_fortran_env 
      implicit none
      integer, parameter :: kq = real128
      real(kq) :: x
      integer  :: i
      do i = 1, 1000, 100
         x = i 
         print '(i3,2(x,g0))', i, exp(-x) * x**x / gamma(x), exp(-x) * x**x / gamma(x + 1)  
      end do    
    end program gam1

計算結果１

  1 .367879441171442321595523770161461 .367879441171442321595523770161461
101 4.00601365311257804891572563285707 .396635015159661192961953032956148E-001
201 5.65363854114232006542228612995131 .281275549310563187334442096017479E-001
301 6.91947392501882452722897920740193 .229882854651788190273388013534928E-001
401 7.98715292661699106820555066172287 .199180870987954889481435178596563E-001
501 8.92805158365303910433045178190687 .178204622428204373339929177283547E-001
601 9.77883421755885529035262693536473 .162709387979348673716349865813030E-001
701 10.5613016450838296424350939106298 .150660508489070323001927159923414E-001
801 11.2896672383911749737148370845054 .140944659655320536500809451741626E-001
901 11.9738080975208380196990875864777 .132894651470819511872353913279397E-001

計算結果を見ると、被積分関数のピーク値は $(n-1)!$ よりは大きく、積分後の値 $(n)!$ よりも（当たり前ですが）小さくなっています。そこで、戯れにその中間の $(n-1/2)!$ 相当の値と比較してみることにします。

プログラム２

    program gam2
      use, intrinsic :: iso_fortran_env 
      implicit none
      integer, parameter :: kq = real128
      real(kq) :: x
      integer  :: i
      do i = 1, 1000, 100
         x = i 
         print '(i3,x,g0)', i, exp(-x) * x**x / gamma(x + 0.5_kq)  
      end do    
    end program gam2

計算結果２

  1 .415107497420594703340268249441337
101 .399106893561339198388218895228938
201 .399024988331960555536648899830374
301 .398997508756508333925078382047327
401 .398983735395490147411070127245623
501 .398975460605820122182191678946810
601 .398969939583726224779254994017756
701 .398965993792399175599028757808330
801 .398963033242696722186803273419373
901 .398960729877773834988078478722955

すると、うまい具合に一定比　0.39896... に収束しているように見えます。0.39896 という値を電卓で適当に叩いて、逆数を取って二乗すると 6.2826... という値になったので多分 $2\pi$ であろうと思われます。

すなわち数値的には、ガンマ関数の階乗相当地点での被積分関数のピーク値は $f(n) / \Gamma(n+{1\over2}) \sim {1\over\sqrt{2\pi}}$ と求まりました。少しずれてはいますが、確かに被積分関数のピーク値が階乗と関係深いことが見て取れます。

スターリングの公式

スターリングの公式は、俗には階乗の近似式 $n!\sim n(\ln(n)-1)$ もしくはこれにもう少し係数などが付いた物を指します。

今、求めた式を少し移項した上で両辺の対数を取ってみると、たしかに $\ln\Gamma(n+{1\over2})=(n-{1\over2})!\sim n\ln(n)-n+\ln(\sqrt{2\pi})$ となってよく似た式になります。

あるいは $1/2$ だけずらすことにすると $n!=\Gamma(n+1)$ なので、 $n\rightarrow n+{1\over2}$ と置き換えて、 $n!=\sqrt{2\pi}f(n+1/2)\sim\sqrt{2\pi}e^{-(n+1/2)}(n+1/2)^{(n+1/2)}$ となります。

数値的に $\ln\Gamma(n+{1\over2})\sim n\ln(n)-n+\ln(\sqrt{2\pi})$ の差を取って比較してみます。

プログラム３

差 $\ln\Gamma(n+{1\over2})-\left({1\over2}\ln{2\pi}-n+n\ln(n)\right)$ とその差に $n$ を掛けたものを出力します。 $n$ を掛けるのはスターリングの公式の近似の高次項からの類比による戯れです。

    program gam3
      use, intrinsic :: iso_fortran_env 
      implicit none
      integer, parameter :: kq = real128
      real(kq), parameter :: pi = 4 * atan(1.0_kq)
      real(kq) :: x
      integer  :: i
      do i = 1, 1000, 100
         x = i 
         print '(i3,2(x,g0))', i, log(gamma(x + 0.5_kq)) - ( log(2 * pi) / 2 - x + x * log(x) ), & 
                                ( log(gamma(x + 0.5_kq)) - ( log(2 * pi) / 2 - x + x * log(x) ) ) * i 
      end do    
    end program gam3

計算結果

  1 -.397207708399179641258481821872649E-001 -.397207708399179641258481821872649E-001
101 -.412538895125284184356880268682438E-003 -.416664284076537026200449071369262E-001
201 -.207296549782888548831458473423111E-003 -.416666065063605983151231531580453E-001
301 -.138427374882831965831612061863506E-003 -.416666398397324217153152306209154E-001
401 -.103906861724177446060175670243060E-003 -.416666515513951558701304437674669E-001
501 -.831669800064454703178056020480413E-004 -.416666569832291806292206066260687E-001
601 -.693288850874847850447310043197465E-004 -.416666599375783558118833335961676E-001
701 -.594388897581945996700061322962287E-004 -.416666617204944143686742987396564E-001
801 -.520183057158577580123800189897250E-004 -.416666628784020641679163952107697E-001
901 -.462449097365583306083968409790626E-004 -.416666636726390558781655537221354E-001

計算結果をみると、差は結構小さくよく近似されていることが分かります。また高次の補正項が $-0.04166666... / n = -1/24n$ であることも分かります。

そこで調子に乗って、もう一段同じ手続きを進めます。

プログラム４

差 $\ln\Gamma(n+{1\over2})-\left({1\over2}\ln{2\pi}-n+n\ln(n)-{1\over24n}\right)$ とその差に $n^3$ を掛けたものを出力します。 $n^3$ は、計算結果からの逆算です。

    program gam4
      use, intrinsic :: iso_fortran_env 
      implicit none
      integer, parameter :: kq = real128
      real(kq), parameter :: pi = 4 * atan(1.0_kq)
      real(kq) :: x
      integer  :: i
      do i = 1, 1000, 100
         x = i 
         print '(i3,2(x,g0))', i, log(gamma(x + 0.5_kq)) - ( log(2 * pi) / 2 - x + x * log(x) - 1.0_kq / (24 * x) ), & 
                                ( log(gamma(x + 0.5_kq)) - ( log(2 * pi) / 2 - x + x * log(x) - 1.0_kq / (24 * x) ) ) * i**3 
      end do    
    end program gam4

計算結果４

  1 .194589582674870254081848447940181E-002 .194589582674870254081848447940181E-002
101 .235900012835689724514385849088463E-008 .243048019124623958856896254701693E-002
201 .299305005315181808524918841287283E-009 .243053652546747070948926446747635E-002
301 .891260273918649549369722899779121E-010 .243054706952683739145563255973093E-002
401 .376939439171983446954764391899663E-010 .243055077420759382117630006617249E-002
501 .193282185350048823554048332146143E-010 .243055249243288499852057323939240E-002
601 .111964863741344148635793726845698E-010 .243055342696905859466697952881662E-002
701 .705588053109556685876677566004617E-011 .243055399095148574762536148061177E-002
801 .472941897939918885290992188350704E-011 .243055435722780067273990289109438E-002
901 .332300511741232087996684894263910E-011 .243055460846571599191284911257351E-002

計算結果をみると、差は $10^{-11}$ 程度とかなり小さくなっていることが分かります。また高次の補正項は電卓をいじっていると $-0.00230554... / n^3 = -7/2880n^3$ であることが分かります。

すなわち $\ln{(n-{1\over2})!}=\ln\Gamma(n+{1\over2})\sim\left({1\over2}\ln{2\pi}-n+n\ln(n)-{1\over24n}+{7\over2880n^3}\right)$ となっていると考えられます。

ここで近似は 10^-17 程度と倍精度を超えるところまで来ており、これ以上は同じやり方を繰り返しても数値誤差のためにうまくゆきませんでした。これを普通のスターリングの公式の補正項と比較してみると、べきも二次の（逆数）項がないなど、こちらの方が補正項が小さいことが分かります。

結論

以上、ガンマ関数が被積分関数のピーク値の寄与しかないと考えることで、おおむね直観的に理解でき、しかもそこから出てくるスターリングの公式に似た式は、~~元のスターリングの公式よりも補正項が小さくなっており、より良い近似になっていることも分かりました~~。

ここでは数値的に補正項を類推しましたが、被積分関数をピーク値での１の幅の短冊積分になっていると思ってやれば、オイラー＝マクローリンの公式で導出できると思われます。

補足[H31.1.26]　休みの徒然に逆さまに迎えに行くと

スターリングの公式より、ベルヌーイ数を使った展開から $\log\Gamma(z)\simeq{1\over2}\log2\pi-z+(z-{1\over2})\log z +{1\over12z}-{1\over360z^3}+\cdots$ なので、 $z$ を $z+{1\over2}$ にまとめ直すとここで求めた係数が出る。

$(z-{1\over2})\log(z-{1\over2})=(z-{1\over2})\log z(1-{1\over2z})=(z-{1\over2})(\log z + \log(1-{1\over2z}))$

対数関数をテイラー展開してやれば
$(z-{1\over2})\log(1-{1\over2z})=(z-{1\over2})(\log z - ( {1\over2z}+{1\over2(2z)^2}+{1\over3(2z)^3}+\cdots))$

開いて移項すれば
$(z-{1\over2})\log z=(z-{1\over2})\log(z-{1\over2})+{1\over2}+({1\over8z}-{1\over4z})+({1\over24z^2}-{1\over16z^2})+({1\over64z^3}-{1\over48z^3})+\cdots$

これを元の式に代入すると、
$\log\Gamma(z)\simeq{1\over2}\log2\pi-(z-{1\over2})+(z-{1\over2})\log(z-{1\over2})+({1\over8}-{1\over4}+{1\over12}){1\over z}+({1\over24}-{1\over16}){1\over z^2}+({1\over64}-{1\over48}-{1\over360}){1\over z^3}+\cdots$

係数を整理して、
$\log\Gamma(z)\simeq{1\over2}\log2\pi-(z-{1\over2})+(z-{1\over2})\log(z-{1\over2})-{1\over24z}-{1\over48z^2}-{23\over2880z^3}+\cdots$

ここで ${1\over24z}$ の項の $z$ を $z-{1\over2}$ の形にすると、
${1\over24(z-{1\over2})}={1\over24z}(1+{1\over2z}+{1\over4z^2})+\cdots={1\over24z}+{1\over48z^2}+{1\over96z^3}+\cdots$
なので、これによって置き換えると二乗の項は消え、また ${1\over96}={30\over2880}$ なので、
$\log\Gamma(z)\simeq{1\over2}\log2\pi-(z-{1\over2})+(z-{1\over2})\log(z-{1\over2})-{1\over24(z-{1\over2})}+{7\over2880z^3}+\cdots$
となる。 $z$ の三次以上の項を無視するなら ${1\over z^3}$ を ${1\over (z-{1\over2})^3}$ に置き換えてもいい。

$z\rightarrow z+{1\over2}$ の置き換えをすると、数値的に求めた式が出る。

補足：
アルティンの本の第三章にガンマ関数の漸近形の議論があります。

*1:なおこの形に一意に決まる条件があるようです。

ガンマ関数入門 (はじめよう数学)

作者:エミールアルティン
日本評論社

Amazon

fortran66のブログ

fortran について書きます。

【ニュース】intel fortran v.19.0.2 出る

intel fortran v.19.0.2

BBC Micro Bit を FORTRAN 77 で

【ニュース】メモ

Fortran 規格目安箱（サーベイ）：最終報告

ハッカーニュースに Fortran ネタ

謎の音楽　Fortran 77

【ニュース】ACM SIGPLAN Fortran Forum 2018年12月号出る！　その他

Compiler Support for the Fortran 2003/08/18 のみ！

サーバーレス Fortran

【新刊】fortran本【予定】

CRC Press Fortran 新刊予定

Damian Rouson 本

トルコ行進曲

その他　HPC 本リスト

【乞食速報】Springer 数学電子本 3,799円引き

新年乞食

【メモ帳】反切の起源

狩野直喜「支那 文学史」

【メモ帳】ガンマ関数とスターリングの公式

ガンマ関数

プログラム１

計算結果１

プログラム２

計算結果２

スターリングの公式

プログラム３

計算結果

プログラム４

計算結果４

結論

補足[H31.1.26]　休みの徒然に逆さまに迎えに行くと

intel fortran v.19.0.2

BBC Micro Bit を FORTRAN 77 で

Fortran 規格目安箱（サーベイ）：最終報告

ハッカーニュースに Fortran ネタ

謎の音楽 Fortran 77

Compiler Support for the Fortran 2003/08/18 のみ！

サーバーレス Fortran

CRC Press Fortran 新刊予定

Damian Rouson 本

その他 HPC 本リスト

新年乞食

狩野直喜「支那文学史」

ガンマ関数

プログラム１

計算結果１

プログラム２

計算結果２

スターリングの公式

プログラム３

計算結果

プログラム４

計算結果４

結論

補足[H31.1.26] 休みの徒然に逆さまに迎えに行くと

謎の音楽　Fortran 77

その他　HPC 本リスト

補足[H31.1.26]　休みの徒然に逆さまに迎えに行くと