fortran66のブログ

fortran について書きます。

【メモ帳】オイラーの定数 γ

オイラーの定数 γ 速めの収束

irresistible integrals をパラパラ眺めていたら、§9.4.1 に \sum_{k=1}^n{1\over k}-\log(n+0.5) がやはりオイラー定数 γ に収束して、しかも収束が速いとありました。

Irresistible Integrals: Symbolics, Analysis and Experiments in the Evaluation of Integrals

Irresistible Integrals: Symbolics, Analysis and Experiments in the Evaluation of Integrals

本来の定義の \sum_{k=1}^n{1\over k}-\log(n) は収束が遅いので有名で、収束誤差が \mathcal{O}(1/2n) なのですが、\log に 0.5 足したものは \mathcal{O}(1/24n^2) の誤差で収束するとのことです。証明もそれほど面倒ではない模様です。

こんな、積分の中点公式みたいなもので加速されて出てくるなんて。お前はまだガンマを知らない、という気持ちです。

お前はまだグンマを知らない DVD

お前はまだグンマを知らない DVD

以下で、計算して確かめてみます。

実行結果

1000 まで足して、小数点以下 7 桁位まで求まっています。比較として本来の定義で計算した値も求めています。こちらは小数点以下 3 桁ほどまで求まっていて、理論通りの結果となっています。

文献値 0.57721 56649 (小数点以下 10 桁)

         100  0.577219790140490       0.582207331651529
         200  0.577216701374822       0.579713581573410
         300  0.577216126324240       0.578881405643301
         400  0.577215924668091       0.578465144068524
         500  0.577215831235245       0.578215331568328
         600  0.577215780449559       0.578048766753451
         700  0.577215749814171       0.577929780547829
         800  0.577215729924379       0.577840534693221
         900  0.577215716284744       0.577771117576444
        1000  0.577215706526555       0.577715581568206
続行するには何かキーを押してください . . .

プログラム

    program euler_const
        implicit none
        integer, parameter :: kd = kind(1.0d0)
        real(kd) :: g, s
        integer :: i
        s = 0.0_kd
        do i = 1, 10**3
            s = s + 1.0_kd / real(i, kind = kd)
            if (mod(i, 10**2) == 0) print *, i, s - log(i + 0.5_kd), s - log(i + 0.0_kd)
        end do
    end program euler_const

【ニュース】Fortran2018 対応 Modern Fortran Explained 11月に出る!

Modern Fortran Explained ~Incorporating Fortran 2018~ 5th edition

(ごっすエディション)

New to this Edition:

Based on the text of the first edition, all the older language levels have been incorporated into a uniform exposition that has then been extensively revised

Four new chapters describe in detail the features of the latest standard, Fortran 2018

表紙が赤くなって、新たに4章が追加される模様です。Fortran 95 Explained 青黄版のような感じでしょうか?

global.oup.com

Modern Fortran Explained (Numerical Mathematics and Scientific Computation)

Modern Fortran Explained (Numerical Mathematics and Scientific Computation)

  • 作者: Michael Metcalf,John Reid,Malcolm Cohen
  • 出版社/メーカー: Oxford University Press, U.S.A.
  • 発売日: 2011/05/08
  • メディア: ペーパーバック
  • この商品を含むブログを見る

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

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

Table of Contents
1. Whence Fortran?
2. Language elements
3. Expressions and assignments
4. Control constructs
5. Program units and procedures
6. Allocation of data
7. Array features
8. Specification statements
9. Intrinsic procedures and modules
10. Data transfer
11. Edit descriptors
12. Operations on external files
13. Advanced type parameter features
14. Procedure pointers
15. Object-oriented programming
16. Submodules
17. Coarrays
18. Floating-point exception handling
19. Interoperability with C
20. Fortran 2018 coarray enhancements
21. Fortran 2018 enhancements to interoperability with C
22. Fortran 2018 conformance with ISO/IEC/IEEE 60559:2011
23. Minor Fortran 2018 features

【ニュース】Arjen Markus 氏の modern Fortran 講義 その他

2018 Workshop on Fortran Modernization for Scientific Application (ForMoSA)

tcevents.chem.uzh.ch

上記タイトルの集まりがあったようで、 Modern Fortran in Practice の著者 Arjen Markus 氏による、Modern Fortran 紹介の他、プリプロセッサを使ったメタプログラミングや、ユニット・テスト、ソースコード文書化など最近の流行りのテーマの実技もあったようです。

Modern Fortran in Practice

Modern Fortran in Practice

Arjen Markus 氏のスライドが以下のページで見られます。
tcevents.chem.uzh.ch

ユニット・テスト等

Fortran におけるユニット・テストや文書化については、2016年の英国計算機協会 Fortran 年会での講演にも出てきているので、英欧方面で流行り出しているのかもしれません。
BCS Fortran Specialist Group Annual General Meeting 2016 Agenda

悪質な言論テロに対してなぜ報道がないのか?

www.sankei.com

 5月30日の各紙夕刊(産経は同31日朝刊)には、犯人が殺人未遂容疑で再逮捕され、ここで韓国籍の男の実名が李宰弦であることが明らかになる。

 犯行の動機は産経が比較的詳しく報道している。同30日の記事では「男は出頭時、『NHKの報道内容に腹が立ってやった』などと話し、その後『日本のメディアに腹が立った』とも話したとされる」とある。

 また、同31日の記事では「李容疑者は19日に出頭した際、『無責任な報道をする日本のメディアへのメッセージだ』『日本のメディアのトップはNHKだと思った』などと話していたが、現在は黙秘しているという」と報じた。

富士通 ARM チップ試作品完成

ソフトバンクがからんだだけで胡散臭く見える ARM・・・ むしろ DMM がからみ始めた PEZY に期待w
cloud.watch.impress.co.jp

NEC SX 賞を貰う

中国で昔でいう姑娘のことを MM と呼んでいたような...ポコペン
it.impressbm.co.jp

【太陽黒点】宇宙天気予報センターサイト・リニューアル

太陽黒点情報

先ごろ、宇宙天気予報センターのサイトがリニューアルされて、太陽黒点のページに行くのが難しくなりました。

以下にリンクをメモっておきます。
黒点情報 -- 宇宙天気情報センター(NICT)


ここ2周期ほど、太陽黒点数が減少して、かつての小氷期を思わせる様相です。これから、極小に向かうようなのでどの位0記録が続くか楽しみです。

これとは別に、太陽の極地方がどうなっているのか、気になっているのですが、だれかロケット飛ばして見てきてください。

【寝言】支那思想史

狩野直喜中国哲学史」

ここしばらく気分転換に少しづつ読んでいました。数年前に、古本屋で千円ばかりで買いました。

武内義雄の「支那思想史」と合わせて読むと、師匠の狩野直喜が提起した問などに、武内義雄が研究を進めて答えている感じで面白いと感じました。支那思想史の簡潔さも、このようなまとめの上に可能になったものかと思います。

以前に読んだときは、宋学以降がつまらなくて、清朝考証学者の列挙羅列にもうんざりしましたが、今回読んでみると東漢後漢)から西漢前漢)へ逆行してゆく清朝考証学の流れも分かって、康有為の孔子教の出所などもはっきりし、興味深く読めました。

共産中国もいずれは儒教へ回帰してゆくと思いますが、戻ろうにもどこへ戻っていいのやら?w 私としては孔子教復活で馬鹿をやる日を楽しみにしています。

本邦の江戸時代の宋学批判者は、漢代を飛ばして戦国期の儒学ないし孔子の考えそのものに遡ろうとしましたが、支那人は遡るのも、もっぱら漢代の注疏に留まったのが面白いところです。

宋学のいいところは、量が多く意味不明で難解な五経は読まなくてよくて、四書だけを読めば十分という隠された真の動機。ハウツー本・ダイジェスト本に通ずるいい感じの思想。ダメなのは、仏教からアイデア丸パクリなのに、それを姑息に隠そうとして、無駄に理屈が複雑なところ。

中国哲学史 (1953年)

中国哲学史 (1953年)

中国哲学史 (岩波オンデマンドブックス)|岩波オンデマンドブックス

中国哲学史 (岩波オンデマンドブックス)|岩波オンデマンドブックス

中国思想史 (岩波全書 73)

中国思想史 (岩波全書 73)

孔子の考えについての寝言:

孔子が周公旦に見た新機軸は自我の発見(内省の発明)だと思われます。当時としては画期的な能力であったのでしょう。(Rubyラーが reflection を自慢するくらいw)

「忠」は己の内面の省察、「恕」は対称性による他人の内面の推察にあたると思われます。

ここで、他人の内心の省察は友愛w方向へ向かうと前提しこれを「仁」としたと考えます。

この前提の根拠は二つあって、親子の愛情の「孝」と兄弟の愛情の「悌」の存在でしょう。(素朴な周初には、アメリカの西部劇に見られるように、○○一家とか○○兄弟とかが当然の協力単位だったのだろうと思います。)

また、孔子はこれらは人間に先験的に与えられているものとして、証明不用なほど自明と見做していたのでしょう。それゆえ「性」を語らなかった事が説明されます。

但し、孔子は親兄弟が骨肉相争う春秋時代を遺憾に思って活動しているので、前提が成立していないという矛盾が生じますw

結局、自我の内省(及びその外延)である忠(恕)の発見・発明が最も大事で、仁は忠恕した結果なすもの。何をなすかは孝悌を前提とした善行。これを親兄弟を離れて外挿してゆけば他人との友、そしてさらにその先に及んで広がってゆくもの。と、私には思われます。

【メモ帳】

Exercise 1.4.3

入力するのが大変w

t1:=33422213193503283445319840060700101890113888695441601636800
t2:=47865783109185901637528055700883208512000182209549068756000
t3:=63273506018045330510555274728827082768779925144537753208000
t4:=77218653725969794800710549093404104300057079699419429079500

FactorInteger[t1]
{{2,6},{5,2},{7,2},{11,1},{13,2},{19,1},{23,1},{31,1},{37,1},{47,1},{59,1},{61,1},{97,1},{113,1},{127,1},{131,1},{137,1},{139,1},{149,1},{151,1},{157,1},{163,1},{167,1},{173,1},{179,1},{181,1},{191,1},{193,1},{197,1},{199,1}}

FactorInteger[200!/110!/90!/t1]
{{1,1}}

200!/110!90! 以後2ズレ

Irresistible Integrals: Symbolics, Analysis and Experiments in the Evaluation of Integrals

Irresistible Integrals: Symbolics, Analysis and Experiments in the Evaluation of Integrals

暇つぶしに読むにはメンドイ感じ
冒頭一番最初の Fibonacci の漸化式で、出だしから本文とインデックスが1ずれていて草。しかもあとから母関数を求めるところでもこのズレが効いてて困る。

【メモ帳】Wolfram 言語で分割数

Wolfram Development Platform にて

さっぱりわからんw
でも一応出来たから確認クイズ付きのチュートリアルは優秀w

Hardy & Ramanujan と Rademacher の式による計算。微分を計算機がやってくれる。GCDも用意してある。

f:id:fortran66:20180620040628p:plain

D[Sinh[Pi/k*Sqrt[2/3*(z-1/24)]]/Sqrt[z-1/24],z]  

(ここで画面に出るメニューの Simplify を押して、出てくる式をコピペして定義しないとうまくいかない。)
f[z_,k_]:=Simplify[ 云々

f0[h_,k_, n_]:=Exp[-2*Pi*I*n*h/k+Pi*I*Sum[(j/k)*(h*j/k-Floor[h*j/k]-1/2),{j,1,k-1}]]  

ff[h_, k_,n_]:=If [GCD[h,k]==1,f0[h,k,n], 0]  

F[k_,n_]:=Sum[ff[h,k,n], {h,1,k}]  

P[n_]:=Sum[Sqrt[k]/(Pi*Sqrt[2])*F[k,n]*f[n,k] ,{k,1,10}] 

Table[Round[N[P[i]]],{i,1,10}]

Round[N[P[1000],1000]]

分割数 p(1000)=24061467864032622473692149727991


fortran66.hatenablog.com