fortran66のブログ

fortran について書きます。

【メモ帳】Backus の質疑応答から

1978年 "The History of Fortran I, II, and III" 講演の質疑応答より

この質疑応答が中々いい所を突いていて面白いのです。

J.A.N. Lee. Transcript of question and answer session: "The History of Fortran I, II, and III". www.softwarepreservation.org

I-N を整数型にした理由

fortran66.hatenablog.com

Destination = Source という代入の向き

今のコンピュータ言語のほとんどは、右辺から左辺への代入になっていますが、FORTRAN の影響によるものだと思います。なぜその順序に決定したか、バッカスが答えています。

COBOL で有名なグレース・ホッパーの同時期の言語では左辺から右辺への代入になっていたようです。

f:id:fortran66:20190820013008p:plain
代入について

アセンブラでは、ソース → デスティネーション と デスティネーション ← ソース の記法の両方があって未だ統一されていません。昔は 80 系と 68 系で逆になっているという話だったと思うですが、80 系が支配的になった後でも統一が成りませんでした。だから、右辺から左辺への代入の向きは自明ではないと思うのです。

FORMAT が実行時にインタープリタ方式で処理される理由

FORTRAN は(書式付きの)I/O が遅いと昔から腐されてきましたが、その理由として FORMAT が実行時にインタープリタ方式で解釈されていることが考えられます。実際 FORTRAN66 以降では、実行時の動的 FORMAT 書き換えが可能です。

コンパイル時に確定させていいはずなのに、なぜ実行時に解釈するのかを、バッカスは担当者に答えさせていまます。理由は、そうするつもりだったけど手が回らなかったためのようです。

f:id:fortran66:20190820014348p:plain
なぜ FORMAT がインタープリタ処理

内部ファイルで、文字列と数値の相互変換をするとき柔軟に行けるんですが、書式付きI/Oでフォーマットを使うので他言語とくらべてどうしても遅くなってしまいます。

f:id:fortran66:20190824022523p:plain
メチルメタフィジー

【メモ帳】雑記

Backus 記憶違い

www.semanticscholar.org

Backus が 1950年代のコンピュータ事情を回顧していますが、手紙等を調べた結果、色々記憶違いがあって自分がそれまで公言していたことに誤りがあったと正直に述べています。

ワシントン (子どもの伝記全集 8)

ワシントン (子どもの伝記全集 8)

Fortran wiki ページ NECFujitsu 抜け

fortranwiki.org

NECFujitsu が表から抜けているようなので、けしからんです。

gfortran-9.2

mag.osdn.jp

9.1 の bugfix の模様です。 gcc.gnu.org

gcc.gnu.org

gcc.gnu.org

9.1 と文書は変わらないようです。

Fortran
Asynchronous I/O is now fully supported. The program needs to be linked against the pthreads library to use it, otherwise the I/O is done synchronously. For systems which do not support POSIX condition variables, such as AIX, all I/O is still done synchronously.

The BACK argument for MINLOC and MAXLOC has been implemented.

The FINDLOC intrinsic function has been implemented.

The IS_CONTIGUOUS intrinsic function has been implemented.

Direct access to the real and imaginary parts of a complex variable via c%re and c%im has been implemented.

Type parameter inquiry via str%len and a%kind has been implemented.

C descriptors and the ISO_Fortran_binding.h source file have been implemented.

The MAX and MIN intrinsics are no longer guaranteed to return any particular value in case one of the arguments is a NaN. Note that this conforms to the Fortran standard and to what other Fortran compilers do. If there is a need to handle that case in some specific way, one needs to explicitly check for NaN's before calling MAX or MIN, e.g. by using the IEEE_IS_NAN function from the intrinsic module IEEE_ARITHMETIC.

A new command-line option -fdec-include, set also by the -fdec option, has been added to increase compatibility with legacy code. With this option, an INCLUDE directive is also parsed as a statement, which allows the directive to be spread across multiple source lines with line continuations.

A new BUILTIN directive, has been added. The purpose of the directive is to provide an API between the GCC compiler and the GNU C Library which would define vector implementations of math routines.

よくよく見ると、ISO_Fortran_binding.h が利用できるようになっているようです。

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

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

Fortran は省エネ!

http://cyber.dabamos.de/blog/#article-2019-06-30

thenewstack.io

Energy consumed Run-time
Imperative 125J 5585ms
Object-Oriented 879J 32965ms
Functional 1367J 42740ms
Scripting 2320J 88322 ms

ベンチマーク競争サイトのプログラムは、ムキ勢が高度抽象記法やランタイム安全チェック機能を捨てて、インライン=アセンブラやら使うのであまりあてにならない気がw もっと素直にその言語が宣伝している機能で書かれている方がよろしいかと。

【メモ帳】Julia 言語で BMP write

Julia 言語よく分からんw

型変換規則がよく分からない。構造体をまとめて書きだす方法も分からない。構造体に初期値を与えられるのかも分からない。ASCII 文字型も分からない。

頭部に定義部が無いと軽い感じになるのだなと思いました。xxx...end の block 構造が重いわけではないのかと。

struct Bmp_file_header
    bfType::AbstractString # 'BM'  B=42 M=4d 
    bfSize::Int32 # file size in bytes
    bfReserved1::Int16
    bfReserved2::Int16
    bfOffBits::Int32
end

struct Bmp_info_header
    biSize::Int32          # = Z'28' ! size of bmp_info_header ; 40bytes 
    biWidth::Int32
    biHeight::Int32
    biPlanes::Int16        # = 1 ! always 1
    biBitCount::Int16
    biCompression::Int32   # = 0 ! 0:nocompression, 1:8bitRLE, 2:4bitRLE, 3:bitfield
    biSizeImage::Int32
    biXPelsPerMeter::Int32 # = 3780 ! 96dpi
    biYPelsPerMeter::Int32 # = 3780 ! 96dpi 
    biClrUsed::Int32       # = 0
    biClrImportant::Int32  # = 0 
end

struct Rgb
    ib::UInt8
    ig::UInt8
    ir::UInt8
end

function wr_bmp(bmp, fn)
    ny, nx = size(bmp)
    Array{Rgb}(undef, ny, nx)
    bmp_file_header = Bmp_file_header("BM", 14 + 40 + 0 + (3 * nx + mod(ny, 4)) * ny, 0, 0, 14 + 40) 
    bmp_info_header = Bmp_info_header(0x28, nx, ny, 1, 24, 0, (3 * nx + mod(ny, 4)) * ny, 3780, 3780, 0, 0)
    open(fn * ".bmp", "w") do file
        write(file, bmp_file_header.bfType)
        write(file, bmp_file_header.bfSize)
        write(file, bmp_file_header.bfReserved1)
        write(file, bmp_file_header.bfReserved2)
        write(file, bmp_file_header.bfOffBits)
        
        write(file, bmp_info_header.biSize)
        write(file, bmp_info_header.biWidth)
        write(file, bmp_info_header.biHeight)
        write(file, bmp_info_header.biPlanes)
        write(file, bmp_info_header.biBitCount)
        write(file, bmp_info_header.biCompression)
        write(file, bmp_info_header.biSizeImage)
        write(file, bmp_info_header.biXPelsPerMeter)
        write(file, bmp_info_header.biYPelsPerMeter)
        write(file, bmp_info_header.biClrUsed)
        write(file, bmp_info_header.biClrImportant)
        
        for iy = ny:-1:1
            for ix = 1:nx
                write(file, bmp[iy, ix].ib, bmp[iy, ix].ig, bmp[iy, ix].ir)
            end
            for ix =1:nx % 4
                write(file, convert(UInt8, 0))
            end    
        end    
    end    
end
nx = 250; ny = 250
B = Array{Rgb}(undef, ny, nx)
for ix = 1:nx
    for iy = 1:ny
        B[ix, iy] = Rgb(convert(UInt8, 0), convert(UInt8, 0), convert(UInt8, 255))
    end
end

wr_bmp(B, "test")

f:id:fortran66:20190816002555p:plain
bmp -> png 手動変換

Julia 1.0 Programming Complete Reference Guide: Discover Julia, a high-performance language for technical computing

Julia 1.0 Programming Complete Reference Guide: Discover Julia, a high-performance language for technical computing

【メモ帳】mybinder 用 gfortran-9 Dockerfile

binder で gfortran-9

gfortran-9 は、まだ add-apt-repository ppa:ubuntu-toolchain-r/test しなければならないので、binder 備え付けの apt.txt では対応できず、Dockerfile を書く必要があります。しかし mybinder サイトの説明を読んでも幼女には難しいので、Opencoarray の Dockerfile を参考にさせていただきました。

Binder Jupyter lab

Binder Jupyter Notebook

f66blog.github.io

f:id:fortran66:20190813204130p:plain
Mybinder Jupyter Gfortran-9 Docker

FROM jupyter/scipy-notebook

USER root

RUN apt-get update -y && \
    apt-get install -y --no-install-recommends software-properties-common && \
    add-apt-repository ppa:ubuntu-toolchain-r/test -y && \
    apt-get update -y && \
    apt-get install -y --no-install-recommends build-essential \
      gcc-9>=9.1.0 \
      gfortran-9>=9.1.0 \
      g++-9>=9.1.0 \
      ${transientBuildDeps} && \
    update-alternatives --install /usr/bin/gcc gcc /usr/bin/gcc-9 60 \
      --slave /usr/bin/gfortran gfortran /usr/bin/gfortran-9 && \
    update-alternatives --set gcc "/usr/bin/gcc-9" && \
    gcc --version && \
    gfortran --version && \
    apt-get clean && \
    apt-get purge -y --auto-remove ${transientBuildDeps} && \
    rm -rf /var/lib/apt/lists/* /var/log/* /tmp/*        

ARG NB_USER
ARG NB_UID
ENV USER ${NB_USER}
ENV HOME /home/${NB_USER}

USER ${NB_USER}
RUN cd ${HOME} && \
    git clone https://github.com/f66blog/binder_gf9.git && \
    cd binder_gf9  && \
    pip install --user ./jupyter-gfort-kernel  && \
    jupyter kernelspec install ./jupyter-gfort-kernel/gfort_spec --user 

WORKDIR ${HOME}/binder_gf9

【メモ帳】分割数 p(721) チェック

分割数 p(721)

ハーディーの「ラマヌジャン その生涯と業績に想起された主題による十二の講義 」に 721 の分割数計算の級数和が出ていたので(講義VIII, p.193)、検算のため比較してみました。

f:id:fortran66:20190810183240p:plain
分割数 p(721)

実行結果

多分、本の方は総和の小数点以下の0の数が1個少ないのではないかと思います。

     1       161061755750279601828302117.84821                                -0.00000
     2                     -124192062781.96844                                 0.00000
     3                           -706763.61926                                 0.00000
     4                              2169.16830                                -0.00000
     5                                 0.00000                                 0.00000
     6                                14.20704                                -0.00000
     7                                 6.07837                                -0.00000
     8                                 0.18926                                 0.00000
     9                                 0.04914                                -0.00000
    10                                 0.00000                                 0.00000
    11                                 0.08814                                 0.00000
    12                                -0.03525                                -0.00000
    13                                 0.03248                                 0.00000
    14                                -0.00687                                 0.00000
    15                                 0.00000                                -0.00000
    16                                -0.01133                                -0.00000
    17                                 0.00000                                -0.00000
    18                                -0.00554                                 0.00000
    19                                -0.00860                                -0.00000
    20                                -0.00000                                -0.00000
    21                                -0.00526                                 0.00000
 
p(721)       161061755750279477635534762.00040                                -0.00000

Hardy-Ramanujan-Rademacher の公式

無限和を 21 で止めています。

\displaystyle{
p(n)=\sum_{k=1}^\infty{\sqrt{k}\over\pi\sqrt{2}}\sum_{0\lt h\leq k\\ (h,k)=1}e^{{-2\pi i n h\over k}+\pi i\sum_{j=1}^{k-1}{j\over k}({h j\over k}-\lfloor{h j\over k}\rfloor-{1\over2})}
\left. {d\over dz}{\sinh\left({\pi\over k}\sqrt{{2\over3}(z-{1\over24})}\right)\over
\sqrt{z-{1\over24}} }\right|_{z=n}
}

プログラム

プログラムは、以前 Fortran で書いたものを用いることにします。倍精度では p(250) 以降で桁が足りなくなるので、4倍精度で計算することにします。フォーマット文や、無限級数の打ち切りなどを調節しました。

fortran66.hatenablog.com

    module m_part
      implicit none
      integer , parameter :: kd = kind(1.0q0)
      real(kd), parameter :: pi = 4 * atan(1.0_kd)
    contains  
      complex(kd) function cpartition(n)
        integer, intent(in) :: n
        integer, parameter :: kmax = 21
        real(kd) :: f, s, t, u
        complex(kd) :: c
        integer :: j, m, k
        cpartition = (0.0_kd, 0.0_kd)
        do k = 1, kmax
          f = sqrt(k / 2.0_kd) / pi 
          c = (0.0_kd, 0.0_kd)
          do m = 1, k
            if (igcd(m, k) == 1) then
              s = -real(2 * n * m, kd) / k + s_hk(m, k)
              c = c + exp(cmplx(0.0_kd, pi * s, kind = kd)) 
            end if
          end do
          cpartition = cpartition + f * c * da(k, n) 
          print '(i6, 2f40.5)', k, f * c * da(k, n) 
        end do  
      contains
        real(kd) function s_hk(h, k)
          integer, intent(in) :: h, k
          real(kd) :: t, u
          integer :: j
          s_hk = 0.0_kd
          do j = 1, k - 1
             t = real(m * j, kd) / k
             u = real(    j, kd) / k
             s_hk = s_hk + (u - floor(u) - 0.5_kd) * (t - floor(t) - 0.5_kd)                
          end do
        end function s_hk
      end function cpartition
  
      integer pure function igcd(ia, ib)
        integer, intent(in) :: ia, ib
        integer :: m(2)
        m = [ia, ib]
        do 
          m = [m(2), mod(m(1), m(2))]
          if (m(2) == 0) exit
        end do
        igcd = m(1)
      end function igcd  
      
      real(kd) pure function da(k, n) ! d/dx sinh( f * xs ) / xs ; xs = sqrt(x - 1/24)
        integer, intent(in) :: k, n 
        real(kd) :: xs, f
        xs = sqrt( real(n, kd) - 1.0_kd / 24.0_kd )
        f  = pi / k * sqrt(2.0_kd / 3.0_kd)
        da = ( f * cosh(f * xs) - sinh(f * xs) / xs ) / (2 * xs * xs) 
      end function da
    end module m_part
  
    program Pn
      use m_part
      implicit none
      integer :: i
      complex(kind = kd) :: z
      do i = 721, 721 
          z = cpartition(i)
        print *
        print '(a, i3, a, 2f40.5)', 'p(', i, ')', z
      end do  
    end program Pn

【寝言】支那人大暴れw

海外で中共派対香港派で争いw

最近ネットは朝鮮ニュースばかりで面白くないですが、一段落したのか、ようやく支那ネタも復活してきました。

阿片窟から這い出てきたような支那人がオースラリア・ニュージーランドで大暴れのようです。

youtu.be

傑作なのは、80年代よろしくラジカセみたいなのを担いだおっさんが能天気なドリフの音楽みたいなのを流しているところです。

youtu.be

www.youtube.com

朝鮮

シェークスピアの「嵐(テンペスト)」に、キャリバンという半獣人みたいしょうもないのが出てきます。魔女(みたいなの)に虐げられていたのを、メインキャラのプロスペロに救われたのですが、自由の身になったはずのキャリバンはプロスペロに感謝するどころか娘を襲おうとします。そして懲罰されると逆恨みして憎しみを連ねるけど、自由にどこへでも行けばいいのに下男としてプロスペロの下で働きつつ裏切るとか、近代的合理主義や理性では訳の分からない存在ですが、韓国を連想させます。

ところが西洋カルチュラル・スタディーズに、プロスペロとキャリバンを欧州植民地主義とそこから解放された土人として読み解くというのがあって草生えました。

テンペスト―シェイクスピア全集〈8〉 (ちくま文庫)

テンペスト―シェイクスピア全集〈8〉 (ちくま文庫)

The tempest - Signet Classic Edition

The tempest - Signet Classic Edition

嵐

シャーロック・ホームズの「プライオリ学校」では、公爵様が若い頃に売春婦の妾に産ませた子供を、素性を隠して秘書として手元に置いていましたが、何度も不品行を繰り返し、そのたびに恥を忍び大金を払ってもみ消してたものの、妾の子は反省も感謝もせず、ついに逆恨みを募らせて公爵の正規の嫡男を誘拐し、爵位と財産を奪おうとする話でしたが、なんとなく韓国を連想させます。

公爵のとっつあんは、妾の子にかつての妾の面影を見出して若気の至りと、甘やかしてしまうのは性格の弱さだとして、何の関係もない嫡男や周囲の人間が命の危険にまでさらされてしまっては堪ったもんじゃないと思います。

シャーロック・ホームズの帰還 (新潮文庫)

シャーロック・ホームズの帰還 (新潮文庫)

The Return of Sherlock Holmes

The Return of Sherlock Holmes

H.G.ウェルズの「モロー博士の島」では、モロー博士に手術されて獣たちが人間のようにされますが、獣人たちはモロー博士を恐れ敬うとともに、文明化のストレスから憎しみを募らせ、モロー博士の死と共に獣性が勝って、野蛮へ戻ることに恐怖しつつも結局狂気に支配されてゆきます。ビクトリア朝末期からエドワード朝にかけて、大英帝国繁栄の終わりの予感から反進化論・退化ブームが来たらしいのですが、その一つと思われます。

モロー博士の島 他九篇 (岩波文庫)

モロー博士の島 他九篇 (岩波文庫)

The Island Of Dr. Moreau: A First Unabridged Edition (Annotated) By H.G. Wells.

The Island Of Dr. Moreau: A First Unabridged Edition (Annotated) By H.G. Wells.

ノーベル賞作家 W.ゴールディングの「蠅の王」でも、遭難したイギリスの学童が、陰気臭くも野獣化してゆく姿を描いていました。

ロビンソン・クルーソー」や「スイスのロビンソン一家」、「冒険ダン吉」みたいに平静日常を失わず、創意工夫で土人をてなづけて面白おかしく愉快にやってもらいたいものです。

蠅の王 (新潮文庫)

蠅の王 (新潮文庫)

Lord of the Flies

Lord of the Flies

完訳ロビンソン・クルーソー (中公文庫)

完訳ロビンソン・クルーソー (中公文庫)

スイスのロビンソン〈上〉 (1950年) (岩波文庫)

スイスのロビンソン〈上〉 (1950年) (岩波文庫)