fortran66のブログ

fortran について書きます。

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

新年乞食

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

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

【メモ帳】反切の起源

狩野直喜支那文学史

[asin:B000J942NI:detail]

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

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

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

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

(めんどくさいので漢字は現代文字でw)

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

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

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

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

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

中国哲学史 (1953年)

中国哲学史 (1953年)

春秋研究

春秋研究

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

ガンマ関数

正の整数値に対してのみ定義されている階乗を、定義域を実数に広げて、微分可能な連続関数として拡張した物がガンマ関数 \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 を0から増やしてゆくと、被積分関数 f(x) は初めは勢いよく増えますが、どこかで天井を迎えてその後0に向かって減衰してゆくことになります。値は常に正になっています。

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

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

プログラム1

被積分関数のピーク値と積分した値に対応するガンマ関数の比を計算します。ガンマ関数は (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

  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)! よりは大きく、積分後の値 (n+1)! よりも(当たり前ですが)小さくなっています。そこで、戯れにその中間の (n+1/2)! 相当の値と比較してみることにします。

プログラム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

計算結果2

  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}) の差を取って比較してみます。

プログラム3

\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 であることも分かります。

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

プログラム4

\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

計算結果4

  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}\right) となっていると考えられます。

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

結論

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

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

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

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

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

【寝言】永谷園東海道五十三次カード

永谷園 東海道五十三次カード

永谷園東海道五十三次カードを集めているのですが、今日お茶漬けを買いに行ったら、なんと欅坂46カードとかいう胡乱なカードに替わっていて憤死。ダブりも入れてまだ三十枚に達していないというのに、この仕打ち。やんぬるかな!やんぬるかな!

梅茶漬けとわさび茶漬けは東海道五十三次でしたが買う気しないし。

www.nagatanien.co.jp

この怒り。
Hey!Hey! カード詐欺100億円中国にあげちゃう!

【乞食速報】Springer 含紙本 3,799円引き! & iSUS 冬のプレゼント その他

Springer 年末セール

ブレグジット、仏黄巾の賊、独ババァ乱心で ポンド・EUR 安来ないかな?もっと安くなれ。

お好きな1冊が最大3,799円割引でご購入いただけます*
www.springer.com

iSUS 冬のプレゼント

www.xlsoft.com

久しぶりに Modern Fortran MEAP Update

OO の章
livebook.manning.com

ACM SIGPLAN Fortran Forum 12月号まだ出ず!

ブレグジットに鑑みて?
ACM SIGPLAN Fortran Forum

來來支那人倶楽部

先月のブルームバーグ主催の昭南島での米支対話。次期大統領選に民主党から出馬を目指しているという噂のブルームバーグ王岐山キッシンジャー、ポールソンを集めた会議だから、支那べったりの話が出るのかと思いきやw 金に狡い元GSのポールソンが中共No.2 王岐山の目の前で米国内の空気を読んで支那戦不可避原稿を棒読みw

Hank Paulson Opening Remarks at New Economy Forum

NHK 漢詩の時間 蘇東坡の巻

蘇東坡は左遷されても悲憤慷慨せず、飄々としているのがいい。そのせいで懲りてないとさらに左遷されたりしてるがw 

前・後赤壁の賦もエピソードも主題もいい内容で、天下の名文と言われるだけのことはある。支那文明もここで終わりにして十分。蛇足が長すぎる。

www4.nhk.or.jp

最近のNHK漢詩の時間は、ほんとつまらん先生のつまらん講義で、時折現代中国語の動物の鳴き声みたいな読みを聞かされたり迷惑。何の意味があるのか。江戸屋猫八のウグイスの鳴きまねでも放送した方がよっぽどいい。

今年の前期は長恨歌をテーマに半年も講義をしながら、史実にも碌に触れず、楊貴妃可哀想、玄宗可哀想のありきたりの感傷で、せっかくのいいテーマが台無しだった。

後期の宋代の詩は予想外につまらなくないが、なんで宋詞を無視するのかよく分からない。

【メモ帳】IBM FORTRAN IV H EXTENDED は四倍精度、非同期 I/O 等あり

IBM 370 用 FORTRAN IV H (1974)

http://web.eah-jena.de/~kleine/history/languages/GC28-6515-10-FORTRAN-IV-Language.pdf

APPENDIX に拡張機能の記述あり。
ASYNCHRONOUS I/O は WAIT など Fortran 2003 で導入された命令がすでに入っています。最近の処理系には大抵入っている四倍精度も、すでに70年代初頭には導入されていたようです。

但しこれは 370 のみのようです。
また Hercules のエミュレータには 370 もありますが、FORTRAN IV H コンパイラは 360 用のものと同じようで、EXTENDED ではなく四倍精度などは使えないようです。


Testing the new VM/370 Sixpack Beta 1.3 distribution - M99

OPTIMIZATION TECHNIQUES FOR FORTRAN IV (G AND H) PROGRAMS WRITTEN FOR, THE IBM 360 UNDER OS

(1971)
https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19730023354.pdf

【遅ニュース】新刊 Chivers & Sleightholme

Introduction to Programming with Fortran

I. Chivers と J. Sleightholme のIntroduction to Programming with Fortran の第4版が夏に出ていたようです。

表紙が前と似通っているので新版だと気づいていませんでした。

www.springer.com


古い部分を残したまま、色々と新しい情報を付け足した模様で、ページが 300 頁ばかり増えているようです。

Introduction to Programming with Fortran

Introduction to Programming with Fortran

Introduction to Programming with Fortran: With Coverage of Fortran 90, 95, 2003, 2008 and 77

Introduction to Programming with Fortran: With Coverage of Fortran 90, 95, 2003, 2008 and 77

Introduction to Programming with Fortran: With Coverage of Fortran 90, 95, 2003, 2008 and 77

Introduction to Programming with Fortran: With Coverage of Fortran 90, 95, 2003, 2008 and 77

Introduction to Programming with Fortran: with coverage of Fortran 90, 95, 2003 and 77

Introduction to Programming with Fortran: with coverage of Fortran 90, 95, 2003 and 77