fortran66のブログ

fortran について書きます。

Fortran から X-Window に点を打つ

X-window グラフィック

これまで Win32 用のグラフィックルーチンは作ってきましたが、X-Window 用も出来ないものかと前々から考えておりました。Fortran だけで実現しようとすると、沢山のインターフェースを用意する必要が出てきますが、Linux 系では C コンパイラは必ず入っているので、C 言語も使うことにしてハードルを下げることにします。

Bash on Windows で実行して、とりあえず、点を打てれば良しとします。


参考ページ:
linux - X11/Xlib.h not found in Ubuntu - Stack Overflow

  • X-Window Graphics の基礎

http://www.cc.kyoto-su.ac.jp/~isida/Pdfs/x.pdf

実行結果

モンテカルロ法で円周率計算
f:id:fortran66:20180707233919p:plain

gfortran / intel fortran どちらでも実行できました。

  • gfortran 実行例
gfortran plot.c plot.f90 -lX11

c コンパイラgnu のものを使います。

gfortran -c -o plot.lib plot.c
ifort plot.f90 plot.lib -lX11

ソース・プログラム

C 言語プログラム plot.c

X window の open/close と 点打ちの3つのルーチンを用意します。

#include <X11/Xlib.h>
#include <X11/Xutil.h>

static Display *d;
static Window  w;
static GC      gc;
static unsigned long white, black;

void X_open(int nx, int ny)
{
  // Open a display.
  d = XOpenDisplay(0);
  if ( !d ) return;
  white = WhitePixel(d, DefaultScreen(d));
  black = BlackPixel(d, DefaultScreen(d));

  // Create the window
  w = XCreateSimpleWindow(d, DefaultRootWindow(d), 0, 0, nx, ny, 0, black, white);
  XMapWindow(d, w);
  gc = XCreateGC(d, w, 0, 0);
}

void X_point(int ix, int iy)
{
  XDrawPoint(d, w, gc, ix, iy);
  XFlush(d);
}

void X_close(void)
{
  XFreeGC(d, gc);
  XDestroyWindow(d, w);
  XFlush(d);
  XCloseDisplay(d);
}
Fortran プログラム plot.f90

C 言語の3つのルーチンに対するインターフェースを用意します。そこから先は Fortran だけで話が済みます。

module m_plot
    implicit none
    interface
         subroutine Xopen(nx, ny) bind(c, name = 'X_open')
             integer, value :: nx, ny
         end subroutine Xopen

         subroutine Xpoint(ix, iy) bind(c, name = 'X_point')
             integer, value :: ix, iy
         end subroutine Xpoint

         subroutine Xclose() bind(c, name = 'X_close')
         end subroutine Xclose
     end interface
end module m_plot


program plot
    use m_plot
    implicit none
    integer, parameter :: n = 10000, nx = 500, ny = nx
    real :: x(n), y(n)
    integer :: i, ix, iy

    call random_seed()
    call random_number(x)
    call random_number(y)

    call Xopen(nx, ny)
    do ix = 0, nx
        iy = sqrt(real(ny**2 - ix**2))
        call Xpoint(ix, iy)
    end do
    do iy = 0, ny
        ix = sqrt(real(nx**2 - iy**2))
        call Xpoint(ix, iy)
    end do

    do i = 1, n
        ix = nx * x(i)
        iy = ny * y(i)
        call Xpoint(ix, iy)
    end do
    print *, 'Monte-Carlo Pi=', 4 * count(x**2 + y**2 < 1.0) / real(n)

    print *, 'press enter'
    read *
    call Xclose()
end program plot

補足[H30.7.10]

環境によっては、窓の立ち上がりが遅くて?待たねばならないようです。XFlush も乱用できません。

#include <X11/Xlib.h>
#include <X11/Xutil.h>
static Display* d;
static Window   w;
static GC       gc;
unsigned long white, black;

void X_open(int nx, int ny){
    // Open a display.
    d = XOpenDisplay(0);
    if ( !d ) return;
    //
    white = WhitePixel(d, DefaultScreen(d));
    black = BlackPixel(d, DefaultScreen(d));
    // Create the window
    w = XCreateSimpleWindow(d, DefaultRootWindow(d), 0, 0, nx, ny, 0, black, white);
    XMapWindow(d, w);
    gc = XCreateGC(d, w, 0, 0);
    XFlush(d);
}

void X_point(int ix, int iy){
     XDrawPoint(d, w, gc, ix, iy);
//     XFlush(d);
}

void X_flush(){
     XFlush(d);
}

void X_close(void){
    XFreeGC(d, gc);
    XDestroyWindow(d, w);
    XFlush(d);
    XCloseDisplay(d);
}
module m_plot
    implicit none
    interface
         subroutine Xopen(nx, ny) bind(c, name = 'X_open')
             integer, value :: nx, ny
         end subroutine Xopen

         subroutine Xpoint(ix, iy) bind(c, name = 'X_point')
             integer, value :: ix, iy
         end subroutine Xpoint

         subroutine Xclose() bind(c, name = 'X_close')
         end subroutine Xclose

         subroutine Xflush() bind(c, name = 'X_flush')
         end subroutine Xflush
     end interface
end module m_plot


program plot
    use m_plot
    implicit none
    integer, parameter :: n = 10000, nx = 500, ny = nx
    real :: x(n), y(n)
    integer :: i, ix, iy

    call random_seed()
    call random_number(x)
    call random_number(y)

    call Xopen(nx, ny)
    call sleep(1)                              ! non-standard
    do ix = 0, nx
        iy = sqrt(real(ny**2 - ix**2))
        call Xpoint(ix, iy)
    end do
    do iy = 0, ny
        ix = sqrt(real(nx**2 - iy**2))
        call Xpoint(ix, iy)
    end do
    call Xflush()

    do i = 1, n
        ix = nx * x(i)
        iy = ny * y(i)
        call Xpoint(ix, iy)
    end do
    call Xflush()
    print *, 'Monte-Carlo Pi=', 4 * count(x**2 + y**2 < 1.0) / real(n)

    print *, 'press enter'
    read *
    call Xclose()
end program plot

【朗報】R. W. Numrich の coarray 本、出版予定復活

Parallel Programming with Co-arrays

10年くらい前から延び延びになって、出版社(Chapman and Hall)のサイトからも消えていた Robert Numrich の coarray 本が復活していました。Fortran 2018 規格も固まったところで、満を持しての出版でしょうか。

Numrich は coarray の発案者なので、内容は確かなのではないかと思われます。

www.crcpress.com

September 12, 2018 の出版予定ですが、何度も延期されてきたので信じないで待ちましょう。

www.youtube.com


Features

  • Provides a tutorial on the coarray model
  • Presents technical specification of the coarray model in enough detail for programmers to write standard-conforming code
  • Includes examples from linear algebra and examples for irregular problems based on graphs
  • Covers object-oriented design using coarrays
  • Contains exercises with solutions
  • Offers Fortran code samples on a supplementary website

Parallel Programming with Co-arrays (Chapman & Hall/CRC Computational Science)

Parallel Programming with Co-arrays (Chapman & Hall/CRC Computational Science)

おまけ

なお同時期に誘導ミサイル本の改訂版が出るようです。Fortran プログラムも付属するようなので、興味ある人は買いましょう。

MATLAB® and Fortran programs are included to illustrate computational approaches for missile guidance

www.crcpress.com


Modern Missile Guidance, Second Edition

Modern Missile Guidance, Second Edition

【メモ帳】ディンキン図の証明で出てくる妙な式に

ディンキン図

能書き

物理のための Lie 群とかいうような本を読んでいますと、三次元の回転群や二次元ユニタリー群などの Lie 代数を見た後、より一般的なルート図を出してきて、次に回転群やユニタリー群の他にどのようなルート図が可能なのかをソ連の数学者ディンキンによるディンキン図によって証明して分類を尽くします。

(これに引き続いて、ウェイトを定義して、ルートとウェイトで議論を続けるのかと思いきや、突如どこからともなくヤング図とワイルの指標公式が降ってきて、ルート図を全く使わずユニタリー群を分類し始めたりして、そんならはじめっからこっちをやれよと言いたくなる展開でずっこけるのですが、すぐに素粒子の分類に入って誤魔化されます。そもそも、有限群のところで群は指標に尽きると言っておきながら、連続群に入るや指標のことをすっかり忘れてルート・ウェイトで話は尽きるようなそぶりを見せつつ、結局天下りのワイルの指標公式で指標に戻る所が、ドリフや植木等なみのずっこけ劇です。ヤング図とワイル公式から行きたければ↓)

The Theory of Group Characters and Matrix Representations of Groups (AMS Chelsea Publishing)

The Theory of Group Characters and Matrix Representations of Groups (AMS Chelsea Publishing)

The Theory of Group Representations (Phoenix Edition)

The Theory of Group Representations (Phoenix Edition)

そもそも、実ベクトル、複素ベクトル、四元数ベクトルのノルムを保存する線形変換として、回転群、ユニタリー群、シンプレクティック群は素直に出てきます。これらの群は比較的古くから知られていて、ワイルの本の題名にちなんで古典群と呼ばれるようです。一方、八元数の場合は結合則が成り立たないため素直に線形変換を作れず、例外群と呼ばれる少数の場合だけが可能になります。例外群は、ディンキン図によってすべてが出尽くしていることが証明されます。その証明はそれほど難しくない不等式によるのですが、その中で妙な式が出てきます。(八元数から先は、体をなさないので考えなくて良い。)

古典群 -不変式と表現- (シュプリンガー数学クラシックス 第)

古典群 -不変式と表現- (シュプリンガー数学クラシックス 第)


丁度、以下のブログにその証明がありますが、その中でも「重み2の辺を1本だけ持つ場合」のところで『ここでなぜか  u=\sum_{i=1}^p i \cdot u_i を考える。』と言ってます。ホントになぜこう置くのか昔から妙に思っておりました。

arxiv.hatenablog.com

はじめに勉強した時は、鉄のカーテンの向こう側の妙なテクニックだなと思っていたのですが、よくよく考えるとこう置くことで、一番きつい条件の不等式を解いていることを示せなければ、分類を尽くしていないように思います。そこで、この妙な式(ベクトルの置き方)について考えてみます。

眠くなってきたので式だけw

上記ブログ 『重み2の辺を1本だけ持つ場合』参照のこと。


解くべきは、不等式 \langle u, v \rangle^2 \leq |u|^2\,|v|^2 を用いての条件出し。

ここではベクトルの係数c_i, d_iを任意実係数としてより一般的に考える。

u=\sum_{i=1}^p c_i\, u_i

v=\sum_{i=1}^q d_i\, v_i
すると、 \langle
 u_i,u_{i+1} \rangle =-{1\over 2}
の条件の下で、

\begin{align}
\left|u\right|^2&=\langle u,u\rangle\\
 &=\sum_{i=1}^{p}c_i^2\langle u_i,u_i\rangle +\sum_{i=1}^{p-1}c_i\,c_{i+1}2\langle u_i,u_{i+1}\rangle \\
&=\sum_{i=1}^{p}c_i^2-\sum_{i=1}^{p-1}c_i\,c_{i+1}\\
&={1\over2}(c_1^2+\sum_{i=1}^{p-1}(c_i-c_{i+1})^2+c_p^2)\\
\end{align}
 v の方も同様。

内積\langle u,v\rangle \langle  u_p, v_q\rangle =-{1\over \sqrt{2}}より、
 
\begin{align}
\langle u,v\rangle&=\langle\sum_{i=1}^pc_i u_i, \sum_{i=1}^q d_i v_i\rangle\\
&=c_p\, d_q \langle u_p, v_q\rangle\\
&=-{c_p\,d_q\over\sqrt{2}}\\
\end{align}

不等式

不等式の条件を出来るだけ厳しくするには、ベクトル長(ノルム)u, vをできる限り短く小さい値に、また内積は出来るだけ大きい値にすればよい(積c_p\,d_qを大に)。

ところでこのベクトル長の式は、数学的に解くのはめんどくさいが、物理的には両端が壁に固定された、ばねで結ばれたpないしq個の質点の作るポテンシャルエネルギー(ばね定数・質量は1、係数c_i, d_iは壁からの距離と見なす)に対応している(?w)。したがって、ある両端の長さに対する最低エネルギー状態は、質点がそれぞれ等間隔に並んだ場合に実現される。これは、c_i, d_iiに比例することを意味する。比例間隔は不等式の両辺でキャンセルされるので、結局u,vは、u=\sum_{i=1}^p i\,u_iv=\sum_{i=1}^q i\,v_iとおいて構わないことになる。こうして妙な式が出てくる。

間違ってたらゴメンw

群と表現 (理工系の基礎数学 9)

群と表現 (理工系の基礎数学 9)

群と表現 (岩波オンデマンドブックス)

群と表現 (岩波オンデマンドブックス)

【寝言】最近の盗掘竹簡と論語の成立事情

盗掘竹簡

この頃、佐藤信弥「中国古代史研究の最前線」を電車の中で読んでいましたが、とても面白かったです。落合淳思の殷代甲骨文に対して、西周金文が専門のようで、とても興味深く読めました。

一番興味深かったのは、漢代の武帝の少し後に即位し、悪宦官霍光によって廃された廃帝の墓からでた論語に、現行論語には無くて斉論にしかないとされる知道篇があったという所で、続報が是非とも知りたいところです。漢の武帝はオカルト好きで、オカルトと親和性の高い斉方面の儒者を贔屓にした点に鑑みてありうるかなと思えます。

近年、支那大陸では盗掘が盛んで得体のしれない竹簡がゾロゾロと闇市に出回っているようで、それを公的な博物館や国立大学などが買い集めているようです。マルクス=レーニン主義唯物論に基づけば、祟りも不敬も関係ないようです。レーニン廟からレーニンのミイラが盗まれる日も遠くないでしょう。

支那人は文字を崇拝して、文盲の苦力ですら道端に文字の書いてあるたばこの包み紙が落ちていると、それを拾って大事にすると、戦前の支那通の人が書いていますが、そのせいか支那人には考古学的な資料が出ても司馬遷史記を優先するような傾向がありましたが、昨今は高い金を出して盗掘品を買いあさって元を取る気になったか、その辺も変わってきているようです。


以前は、浅野祐一、湯浅邦弘が熱心に新出土竹簡の真正性から論じて、内容を紹介していましたが、新しい世代は真正性は良しとして、より内容に踏み込んで論じるようになっているようです。

新しい世代は、長きにわたったマルクス主義の軛から解放されて、のびのび明快に論じているのが気持ちよいです。

浅野・湯浅本などによると、盗掘された竹簡は水に浸ってブヨブヨになっているとかで、支那そばの上にのってる支那竹みたいなものかと空想していましたが、「中国古代史研究の最前線」には、泥に埋もれたチャーシューの塊のような写真が載っていました。

古代中国の宇宙論

古代中国の宇宙論

諸子百家「再発見」―掘り起こされる古代中国思想

諸子百家「再発見」―掘り起こされる古代中国思想

諸子百家 (講談社学術文庫)

諸子百家 (講談社学術文庫)

湯浅邦弘の論語本も、冒頭に論語の成立事情に関して、最近の盗掘竹簡の知見を交えつつ論じていますが、「中国古代史研究の最前線」はこれよりも竹簡に関しては新情報が載っているようでした。

論語 - 真意を読む (中公新書)

論語 - 真意を読む (中公新書)

【乞食速報】Springer 電子本 14.99 ユーロ

Springer セール

Springer のセールも最近頻度が落ちているうえに、値段も以前の $10 から 15EUR にかわってお得度が落ちています。

1842年まで遡及する膨大な知の遺産にアクセス!
✔ 50,000点以上もの優れたタイトルが14.99€*
✔ 割引適用にはクーポンコードBOOKARCHIVES18を入力してください
✔ 本キャンペーンは2018年7月18日まで有効です

www.springer.com

Fortran の古いのや、PASCAL, Modula-2, Oberon 等のN. Wirth 言語本も買えます。

【寝言】津田左右吉「論語と孔子の思想」

論語孔子の思想

近頃、暑いので斜め読みしておりました。

前半の論語の成立事情に関する考察は、先行研究の結果に基づいているはずですが、さっぱりその言及も参照文献もなく、どこまでが事実として認められていることで、どこからが著者の推測・意見か判然としないので、あまり読んでいて面白くなかったです。

前漢までの論語の成立を論じていますが、事実関係としては、武内義雄の「論語の研究」などから引いてきたと思われるようなものが多いです。

しかし、後半の第五編「論語儒学の学」当たりになりますと、支那人の思考様式そのものを批評して、否定にかかって面白くなってきます。支那文明そのもののに起因して儒教という虚飾がまかり通っていると明快に述べております。最後は儒家・法家・道家すべてを近代人的な合理主義から否定して、小気味いいです。

これは漢学者などの著作には見られないもので、支那人自己批判は望めないにしても、誰かが言わなくてはならないものであるのだから、自由世界の文明国民たる日本人がどうしても言わなくてはならないことでしょう。

この点については、結語「論語の研究の方法と態度」にも明示されており、はなはだもっともだと思われるところです。

しかし、儒教を捨てても共産主義マルクスレーニン主義)に嵌る支那人の愚昧さには草生えるw しかもその虚飾性といい、やってることは儒教時代と何ら変わりがない。

論語之研究 (1940年)

論語之研究 (1940年)

論語之研究 (1972年)

論語之研究 (1972年)



支那人とインデアン

支那人は、王を神格化するのではなく、禹や鯀のように神話の神の方を引きずり落として王としたという考察は興味深いです。支那人の宗教心のなさは、自然の崇高さへの敬意の欠如の産物であろうというのは、さもありなんという感じです。

ところで、仰韶文化の焼き物の装飾をみていると、現代アメリカ・インデアンの砂絵に見られる装飾とよく似たものがあるので、興味深いものだなと思います。

f:id:fortran66:20180702024927j:plain
https://kknews.cc/zh-mo/culture/lpr2zye.html

f:id:fortran66:20180702025044j:plain
◆ネイティブアメリカン インディアン砂絵サンドペイント4方向*ホピMESAアパッチ - ANTIQUE BITTE

【メモ帳】Carrige-Control について その他

FORTRAN carriage-control

ラインプリンタ時代の名残りで、昔の FORTRAN の FORMAT 文の先頭にいつもついていた、FORMAT(1H , とか (1H1, とか (1H+, とか、あの出力制御文字について、Clive Page 氏などがまとめをしてくれていたようです。

FORTRAN carriage-control

規格からの廃止は Fortranメインフレーム)の行志向が、近年の Unix のバイト志向のために追いやられてしまったことが背景にあるようです。また出力装置が carriage-control に対応するかどうかは規格上曖昧だったため、コンソールでの対応もまちまちになっているようです。

処理系ごとにコンソールへの WRITE 文出力で、行頭が1文字空いたり、空かなかったりするのはそのためなのでしょう。

Intel Fortran では vms 互換オプションで有効になります。

補足

Clive Page 氏は、Professional Programmer's Guide to Fortran77の著者だと思います。77 向けの高度な利用法が書かれた、実戦的な内容のよい本だった記憶があります。かなり以前から著者によって ps ファイルも無料公開されており PDF 版もネットに落ちています。

https://www.star.le.ac.uk/~cgp/prof77.html

Professional Programmer's Guide to Fortran 77 (Professional Programmers Guides)

Professional Programmer's Guide to Fortran 77 (Professional Programmers Guides)


amazon(日)には、ぺ氏の別の著作 Fortran 77 (Pocket Programming Guide) がありますが、著者名が間違っていて、Ulive Page になっています。先頭文字が 90 度回転してますw

Pocket Guide: Fortran 77 (Pocket Programming Guide)

Pocket Guide: Fortran 77 (Pocket Programming Guide)

こちらはポケット・リファレンス本で、内容的にはまだ悟りが開けていない普通な感じだったような気がします。