fortran66のブログ

fortran について書きます。

ホーアの「Hints on programming language design」に於ける FORTRAN 評

C. A. R. Hoare の FORTRAN IV (FORTRAN66) 時代の FORTRAN 評です。

f:id:fortran66:20170829014010p:plain

Hints on programming language design.

皮肉が効きすぎて、イマイチ最後のオチが分かりませんw


ASA Standard FORTRAN
前書き
W. P. Heising, History and summary of FORTRAN standardization development for the ASA
History and summary of FORTRAN standardization development for the ASA
本文
S. Gorn, FORTRAN vs. Basic FORTRAN: a programming language for informational processing on automatic data processing systems
FORTRAN vs. Basic FORTRAN







バッカスの考えた FORTRAN に、AHO の書いたコンパイラ、ホーアがいちゃもん・・・



f:id:fortran66:20170829014441p:plain
AMAZON 本 ☟


[追記:H29/9/3]
ホーア曰く
Charles Anthony Richard Hoare

I don't know what the programming language of the year 2000 will look like, but I know it will be called FORTRAN.

You can't teach beginning programmers top-down design because they don't know which way is up.

これ好きかもw

Interesting Quotes より

oberon と golang

Google 社の go 言語は、oberon の C 表記版だということが言われているようです。

開発者の一人で N. Wirth の下で学位を取得した Robert Griesemer 氏がそれに関連する講演を行っています。

GopherCon 2015: Robert Griesemer - The Evolution of Go


GopherCon 2015: Robert Griesemer - The Evolution of Go

講演スライド
The Evolution of Go

Modula-2, Oberon
f:id:fortran66:20170828220345p:plain

OO や Generics について
f:id:fortran66:20170828220119p:plain

type の継承メカニズムがなく、メソッドがポインタのようなもので、コルーチン的コンカレント処理がある点に鑑みると、Oberon よりも Modula/Modula-2 の延長のような気もします。

記事 Link

Multiprocessing in Python with Fortran and OpenMP

マルチスレッド・マルチコアの利用が難しい Python から、Fortran 経由で OpenMP を利用しようという記事です。

www.admin-magazine.com

OpenMP は元々は、Fortran の並列化コンパイラ指示行が各ベンダーごとにばらばらだったのを統一しようとして生まれた規格です。後に C 言語が乗ってきて、最近では母屋を取られて C 言語規格の方が主導するようになりました。元々共有メモリー用なのに意外にオーバーヘッドが重かった気がします。

N. Wirth "From modula to Oberon"

Oberon で落とした機能についてその根拠を述べています。

onlinelibrary.wiley.com

https://www.research-collection.ethz.ch/bitstream/handle/20.500.11850/68917/eth-3189-01.pdf

ガンマ関数

Paul J. Nahin "Inside Interesting Integrals"


Paul J. Nahin "Inside Interesting Integrals" という本を眺めていると、\int_0^\infty e^{-y^2}dy, \int_0^\infty e^{-y^3}dy, \int_0^\infty e^{-y^4}dy, ・・・, \int_0^\infty e^{-y^n}dy という積分がでてきて、\int_0^\infty e^{-y^2}dy=\sqrt{\pi}/2は Gaussian の積分で常識として、2 乗以外は、3 乗も 4 乗も使われるのを見たこともないし、どうでもいいですわ~と思ったのですが、これを y^n=x の変数変換をすると、おなじみのガンマ関数が出てきてぎゃふんとなりました。

ny^{n-1}dy=dx よりdy={1\over n}x^{{1\over n}-1}dxなので、\int_0^\infty x^{{1\over n}-1}e^{-x}{1\over n}dx={1\over n}\Gamma({1\over n})=\Gamma({1\over n}+1) となります。

これはガンマ関数を提示するのに、なかなか面白い導入部ではないかと思いました。さて、e^{-y^n}は、n が 1 より大きい正の整数の時は急激に減少する関数なので、積分範囲の上限が +∞ でも、実際には小さな値まで取れば十分になっています。

数値計算

ここでは数値積分を計算して、それを Fortran2008 から組み込み関数になったガンマ関数と比較してみます。

ソース・プログラム

台形積分をリチャードソン補外することで、シンプソン法相当の計算をしています。はじめからシンプソン法をやればいいかもしれませんが、もう一段やればロンバーグ法相当になるので、まぁいいでしょうw

amazon AA

    program p_Gamma
      implicit none
      real, parameter :: pi = 4 * atan(1.0)
      integer, parameter :: n = 2**8  
      integer :: i
      real :: a, b, x(n), y(n), s, t, h, f
      !
      a = 0.0
      b = 15.0
      h = (b - a) / (n - 1)
      forall (i = 1:n) x(i) = a + h * (i - 1)
      !  
      do i = 1, 10
        f = real(i)
        y = exp(-x**f)
        s = richard(h, y) ! ~ simpson's method
        t = gamma(1.0 / f + 1) !gamma(1.0/i) / i
        print *, i, t, s, abs(t - s) / t  
      end do    
    contains
      real function trap(y) ! Trapezoidal rule
        real, intent(in) :: y(:)
        integer :: n
        n = size(y)
        trap = 0.5 *(y(1) + y(n)) + sum(y(2:n-1))
      end function trap
  
      real function richard(h, y) ! Richardson extrapolation
        real, intent(in) :: h, y(:)
        real :: s1, s2
        s1 =     h * trap(y)
        s2 = 2 * h * trap(y(::2))
        richard = (4 * s1 - s2) / 3 
      end function richard
    end program p_Gamma

実行結果

Γ(1/n+1) を\int_0^\infty e^{-y^n}dy積分を数値的に計算して、組み込み関数の値と比較しています。積分範囲の∞を 15 にしたとき、相対誤差は単精度の限界のあたりに来ています。

           1   1.000000      0.9999995      4.7683716E-07
           2  0.8862270      0.8862269      6.7256636E-08
           3  0.8929795      0.8929799      4.0048832E-07
           4  0.9064025      0.9064024      6.5759579E-08
           5  0.9181687      0.9181687      0.0000000E+00
           6  0.9277194      0.9277192      1.9274572E-07
           7  0.9354376      0.9354377      1.2743693E-07
           8  0.9417427      0.9417428      1.2658371E-07
           9  0.9469653      0.9469649      4.4059956E-07
          10  0.9513507      0.9513475      3.3832430E-06
続行するには何かキーを押してください . . .

無印 MODULA

無印 Modula は Oberon に近いミニマリズム

ここで無印 Modula とは、Modula-2 ではなく、初代の Modula のことを指します。

link.springer.com
この記事に、pascal と modula の相違点がまとめられています。これによりますと、この時点ですでに variant record, FOR 文, 集合型(set type), 部分範囲型(subrange type) などが省かれており、Oberon での切り捨てがこの時点ですでに萌芽していたことが分かります。

f:id:fortran66:20170822023624p:plain

FOR 文を捨てて、WHILE 一本さらしに巻いて突き進むのが 1976 年の時点というのが興味深いです。

また E.W.Dijkstra の『構造化プログラミング』の中の C.A.R.Hoare の「データ構造化序論」は Wirth の PASCAL と密接な関係があるわけですが*1、理由はともあれ早くもその中の直和(variant 型)、べき集合(集合型)、数え上げ(部分範囲型)が要らない子になっているのも面白いです。

Structured programming



Module の use

Fortran では Module を引用して利用するときに、use モジュール名 の形で宣言しますが、use は modula に現れていることが分かりました(ただしモジュール名ではなく、モジュール内の個別要素名)。Modula-2 では import を用いて、今の Python と全く同じ形式で引用を行っています。

無印 Modula 仕様書

無印 Modula の仕様書が ETH のデジタルアーカイブにあって読むことが出来ます。

MODULA
a language for modular multiprogramming
MODULA

*1:Recollections about the development of Pascal p.103 "Whereas the basic framework of Pascal stems from ALGOL W, many of the new features emerged from suggestions made by C. A. R. Hoare, including enumeration, subrange, set, and file types. The form of COBOL-like record types was due to Hoare, as well as the idea to represent a computer's "logical words" by a well-known abstraction--namely, sets (of small integers). These "bits and pieces" were typically presented and discussed during meetings of the IFIP Working Group 2.1 (ALGOL), and thereafter appeared as communications in the ALGOL Bulletin. They were collected in Hoare's Notes on Data Structuring [Hoare 1972]."

メモ帳 oberon & javascript

Oberon で分割数計算

Oberonは何故かロシアで人気があったようです。
[追記 H29/8/22]
link.springer.com


ロシア製 Oberon → javascript 変換機
Oberon online compiler
出力がイマイチ? JS.alert しかない?


Fortran 90 出始めの頃や、ELF90 や F で命令を大文字で書く習慣があったのは、 MODULA-2 や Oberon の影響だったのではという気がしてきました。もっとも、Fortran 2003 までは小文字は使ってもいいけど、処理系が必ずしも用意する必要もないということだったらしいので、それが理由だという説もありますが。

MODULE test;
IMPORT JS;
VAR i: INTEGER;

PROCEDURE part(n, m:INTEGER): INTEGER;
VAR res: INTEGER;
BEGIN  
  IF n = 1 THEN 
   res := 1 
  ELSIF m = 1 THEN 
   res := 1
  ELSIF n < 0 THEN 
   res := 0 
  ELSIF m < 0 THEN 
   res := 0 
  ELSE
   res := part(n - m, m) + part(n, m - 1)
  END;
  RETURN res
END part;

PROCEDURE partition(n :INTEGER): INTEGER;
BEGIN  
  RETURN part(n, n)
END partition;
 
BEGIN
 i := 0;
 WHILE i < 10 DO
    i := i + 1;
   JS.alert(partition (i));
 END
END test.

javascript で分割数計算

var pp = function (n, m) {
  if (n == 0) {return 1};
  if (m == 1) {return 1};
  if (n <= 0) {return 0};
  if (m <= 0) {return 0};
  return pp(n - m, m) + pp(n, m - 1);
}

var partition = function (n) { return pp(n, n) };

for(i = 1; i < 11; i++) {
  console.log(partition(i));
};
1 
2 
3 
5 
7 
11 
15 
22 
30 
42 
undefined

メモ帳 S-function の正規化 2

Littlewood の S-function の standardization 2

夏休みの勉強のつづき
fortran66.hatenablog.com


分割{λ}= [\lambda_1,\lambda_2,..., \lambda_p] を S-function とみて、{λ}が弱い減少列になっていない時に弱い減少列に直す。

規則1 尻の 0 は取り去る。
規則2 尻に負数が出たら 0。(ここでは係数項 0 とする)
規則3 隣り合う数の入れ替えは\{...\lambda_i,\lambda_{i+1},...\}=−\{...\lambda_{i+1}−1,\lambda_i+1,...\}


規則3 より  \lambda_i = \lambda_{i+1} - 1 になっているところがあると、\{\lambda\}=-\{\lambda\}となって S-function は 0 になる。

ソース・プログラム

書き直したプログラム

main = print (standard (1, [-5,1,1,1,1,1])) 

type Partition = [Int]
type PartitionF = (Int, Partition) -- Partition with a (phase)factor


standard :: PartitionF -> PartitionF
standard (f, xs) = if is_std xs then (f, xs)
                                else (standard . std) (f, xs)

is_std :: Partition -> Bool 
is_std []  = True
is_std xs  = if last xs <= 0 then False 
                             else and [x >= y|(x,y) <- zip xs (tail xs)]

std :: PartitionF -> PartitionF
std (f, [])= (f, [])
std (f, xs) = if (last xs) <  0  then (0, [])
         else if (last xs) == 0  then standard (f, init xs)
         else std' (f, xs)

std' :: PartitionF -> PartitionF                           
std' (f, []) = (f, [])
std' (f, [x]) = (f, [x])
std' (f, x:y:xs) = if x == (y-1) then (0, [])
              else if x < y      then (-f, (y-1):(x+1):xs) 
              else (f', x:xs')
    where 
        (f', xs') = std' (f, y:xs)

実行結果

[-5,1,1,1,1,1] = -[0]= -[] = -1

sh-4.3$ main                                                                                                                                                                 
(-1,[])    

対称群の指標を目指して (マーナハン・中山方式)

Murnaghan-Nakayama の式の部品

チェック用
\chi^{\{31^3\}}_{(6)}=-1, \chi^{\{21^4\}}_{(6)}=+1

生意気に Maybe モナドを使ってみましたが、くるんだ数値を取り出すのに一苦労w

ソース・プログラム

import Data.Maybe

xs3111_6  = [[-3,1,1,1], [3,-5,1,1], [3,1,-5,1], [3,1,1,-5]]
xs21111_6 = [[-4,1,1,1,1], [2,-5,1,1,1], [2,1,-5,1,1], [2,1,1,-5,1], [2,1,1,1,-5]]

main = print ((cmn . mn') xs3111_6) --xs21111_6) 

type Partition = [Int]
type PartitionF = (Int, Partition) -- Partition with a (phase)factor

mn' :: [Partition] -> [PartitionF]
mn' xss = map (\xs -> standard (1, xs)) xss 

cmn :: [PartitionF] -> Int
cmn []   = 0
cmn (x:xs) =  fromJust (ipf0 x)  +  (cmn xs)

ipf0 :: PartitionF -> Maybe Int
ipf0 (f, []) = Just f
ipf0 (f, _)  = Nothing 

実行結果

sh-4.3$ main                                                                                                                                                                 
-1     

amazon AA この Hutton の本は薄くて例題の話題は豊富で面白い。みなまで言わない典型的英国風の教科書か。