ケーリー＝ディクソン構成二通り

背景

実数・複素数・四元数・八元数

複素数 a+bi を２つの実数（a,b)の組と考えて、その加減乗除の演算規則が適当に与えられている代数系を考えることが出来ます。ここで、a,b に実数じゃなくて複素数を突っ込んだらどうなるかと考えると、うまく適当な演算規則をあたえると、矛盾なく定義が出来て超複素数というべきもの（四元数）が作れることが知られています。この作り方はケーリー=ディクソン構成と言われるものです。

ケーリー=ディクソン構成
積　:(a,b)(c,d):=(ac-~db, da+b~c)
共役:~(a,b):=(~a,-b)

四元数は単なる空想の産物ではなく、自然界では量子力学のスピンのところでパウリのスピン行列として登場してきます。四元数は行列表現があることから推察できるように、交換則が成り立ちません。

調子に乗って進むともう一段階進むと、八元数というものも定義出来て加減乗除の演算規則が作れます。しかし、八元数では結合則も成立せず、だんだんと手足を切り取られてダルマさんになってしまい、この先は加減乗除可能な１６元数というものは作れなくなってしまいます。

ここで、それぞれの数のノルムを考えると、実数ベクトルのノルムを保存するのは回転群、複素数のベクトルのノルムを保存するのはユニタリー群、四元数ベクトルのノルムを保存するのは斜交群（シンプレクティック群）、八元数のノルムを保存するのは例外群（特にG2）とリー群の分類と関係しています。

commutator・associator

交換則からの逸脱具合は、commutator（交換子積） [A,B]=AB-BA、結合則からの逸脱具合は、associator [X,Y,Z]=(XY)Z-X(YZ) で表現されます。 associator は三項演算子になっています。ところで交換子積は非結合積なので、[X,Y,Z]=[[X,Y],Z]-[X,[Y,Z]]となりますが、Jacobiの恒等式によりこの値は[Y,[X,Z]]と分かります。これが交換子積の非結合具合となります。

八元数では、前述のとおり一般には結合則は満たさないのですが、部分的に二つの元の演算では結合則が成り立ちます。すなわち[x,x,y]=[x,y,x]=[y,x,x]=0 となります。

プログラム

ケーリー=ディクソン構成は、実数からはじまって同じアルゴリズムが三回反復されて、複素数・四元数・八元数が作られます。Fortran90/95 の場合、派生型を使って下から組み上げていく形で順々に定義してゆき、最後に演算子オーバーロードで二項演算子などを定義する形が素直だと思いますが、ほぼ同じ内容の繰り返しが必要となります。これに対し、Fortran2003/08では無制限ポインタを使って、上側から下側の実数に向かってゆく再帰的な定義で反復なしに１回の記述で定義することが出来ます。（ただし無制限多態ポインタ使用のため見苦しいプログラムになります。）

以下に両方式のプログラムを示します。ここで簡単のため、加減乗除のうち加算と乗算のみを定義しいます。減算と除算は、加算と乗算から定義されるのでここでは省きました。

ソースプログラム・出力結果

Fortran90/95

 (533.0000,0.0000000E+00)
   533.0000      0.0000000E+00 cc
   533.0000      0.0000000E+00  0.0000000E+00  0.0000000E+00
  -6.000000       13.00000       32.00000       1.000000       50.00000
  -1.000000      -12.00000       67.00000

   204.0000      0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00
  0.0000000E+00  0.0000000E+00  0.0000000E+00
続行するには何かキーを押してください . . .

    module m_comp
      implicit none
      interface operator(+) 
        procedure :: plus_r, plus_c, plus_h, plus_o
      end interface
      interface operator(*) 
        procedure :: prod_r, prod_c, prod_h, prod_o
      end interface
      interface operator(-) 
        procedure :: minus_r, minus_c, minus_h, minus_o
      end interface
      interface operator(.cc.) 
        procedure :: conj_r, conj_c, conj_h, conj_o
      end interface
      interface conj 
        procedure :: conj_r, conj_c, conj_h, conj_o
      end interface
      type :: t_r
        real :: r
      end type t_r    

      type :: t_c
        type(t_r) :: re, im
      end type t_c    

      type :: t_h
        type(t_c) :: re, im
      end type t_h    

      type :: t_o
        type(t_h) :: re, im
      end type t_o    

    contains  
    
      type(t_r) function conj_r(a) result(c)
        type(t_r), intent(in) :: a
        c%r = a%r 
      end function conj_r
      type(t_r) function minus_r(a) result(c)
        type(t_r), intent(in) :: a
        c%r = -a%r 
      end function minus_r
      type(t_r) function plus_r(a, b) result(c)
        type(t_r), intent(in) :: a, b
        c%r = a%r + b%r
      end function plus_r
      type(t_r) function prod_r(a, b) result(c)
        type (t_r), intent(in) :: a, b
        c%r = a%r * b%r
      end function prod_r

      type(t_c) function conj_c(a) result(c)
        type(t_c), intent(in) :: a
        c%re =  a%re 
        c%im = -a%im 
      end function conj_c
      type(t_c) function minus_c(a) result(c)
        type(t_c), intent(in) :: a
        c%re = -a%re 
        c%im = -a%im 
      end function minus_c
      type(t_c) function plus_c(a, b) result(c)
        type (t_c), intent(in) :: a, b
        c%re = a%re + b%re
        c%im = a%im + b%im
      end function plus_c
      type(t_c) function prod_c(a, b) result(c)
        type (t_c), intent(in) :: a, b
        c%re = a%re * b%re + (-b%im * a%im)
        c%im = b%im * a%re +   a%im * b%re
      end function prod_c
      type(t_c) function to_c(r0, r1)
        real, intent(in) :: r0, r1
        to_c = t_c(t_r(r0), t_r(r1))
      end function to_c  
      
      type(t_h) function conj_h(a) result(c)
        type(t_h), intent(in) :: a
        c%re = .cc.a%re 
        c%im =    -a%im 
      end function conj_h
      type(t_h) function minus_h(a) result(c)
        type(t_h), intent(in) :: a
        c%re = -a%re 
        c%im = -a%im 
      end function minus_h
      type(t_h) function plus_h(a, b) result(c)
        type (t_h), intent(in) :: a, b
        c%re = a%re + b%re
        c%im = a%im + b%im
      end function plus_h
      type(t_h) function prod_h(a, b) result(c)
        type (t_h), intent(in) :: a, b
        c%re = a%re * b%re + (-(.cc.b%im) * a%im)
        c%im = b%im * a%re + a%im * (.cc.b%re)
      end function prod_h
      type(t_h) function to_h(r0, r1, r2, r3)
        real, intent(in) :: r0, r1, r2, r3
        to_h = t_h(to_c(r0, r1), to_c(r2, r3))
      end function to_h 
      
      type(t_o) function conj_o(a) result(c)
        type(t_o), intent(in) :: a
        c%re = .cc.a%re 
        c%im =    -a%im 
      end function conj_o
      type(t_o) function minus_o(a) result(c)
        type(t_o), intent(in) :: a
        c%re = -a%re 
        c%im = -a%im 
      end function minus_o
      type(t_o) function plus_o(a, b) result(c)
        type (t_o), intent(in) :: a, b
        c%re = a%re + b%re
        c%im = a%im + b%im
      end function plus_o
      type(t_o) function prod_o(a, b) result(c)
        type (t_o), intent(in) :: a, b
        c%re = a%re * b%re + (-(.cc.b%im) * a%im)
        c%im = b%im * a%re + a%im * (.cc.b%re)
      end function prod_o
      type(t_o) function to_o(r0, r1, r2, r3, r4, r5, r6, r7)
        real, intent(in) :: r0, r1, r2, r3, r4, r5, r6, r7
        to_o = t_o(to_h(r0, r1, r2, r3), to_h(r4, r5, r6, r7))
      end function to_o  
    end module m_comp
    
    program supercomplex
      use m_comp
      implicit none
      complex :: a, b, c
      type(t_r) :: ra, rb, rc
      type(t_c) :: ca, cb, cc 
      type(t_h) :: ha, hb, hc 
      type(t_o) :: oa, ob, oc 
      ra = t_r(2.0)
      rb = t_r(3.0)
      
      a = (2.0, 3.0)
      b = (4.0, 5.0)
      c = a * b
      ca = to_c(2.0, 3.0)
      cb = to_c(4.0, 5.0)
      cc = ca * cb
      print *, conjg(c) * c
      print *, (.cc.cc) * cc, 'cc'
      
      ha = t_h(ca, to_c(0.0, 0.0))
      hb = t_h(cb, to_c(0.0, 0.0))
      hc = ha * hb
      print *, (.cc. hc) * hc
      hc = to_h(0.0, 0.0, 0.0, 0.0)

      oa = to_o(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0)
      ob = t_o(hb, hc) 
      oc = oa * ob
      print *, oc
      print *
      print *, (.cc. oa) * oa
    end program supercomplex

Fortan2003/08

ケーリー=ディクソン構成を再帰的に実現しますが、再帰の出発値に当たるのは、組を作る数の型として intrinsic な real 型が現れた時になります。 intrinsic な real 型は継承できないなので、real 型を他の派生型と同等に扱うには、無制限多態ポインタ class(*) を使う必要があります。しかし class(*) は代入するにも出力するにも何するにも select type で型を確定する必要があるので厄介です。このため出力用のユーザー定義ルーチンも用意しておきました。（楽な抜け道があるのかもしれませんが、私の今の知識ではこれで精一杯）

   -3.000000
   -5.000000       10.00000
 (-5.000000,10.00000)
   -60.00000       12.00000       30.00000       24.00000
   -60.00000      -12.00000      -30.00000      -24.00000
   -1000.000      -1200.000      -1400.000      -1600.000       608.0000
   400.0000       600.0000       800.0000
   -1000.000      -1200.000      -1400.000      -1600.000       608.0000
   400.0000       600.0000       800.0000

   -1000.000       1264.000       1488.000       1712.000      -104.0000
  -288.0000      -472.0000      -656.0000
   -1000.000       1264.000       1488.000       1712.000      -104.0000
  -288.0000      -472.0000      -656.0000

   -1000.000      -1184.000      -1368.000      -1552.000      -104.0000
   528.0000       752.0000       976.0000
   -1000.000      -1184.000      -1368.000      -1552.000      -104.0000
   528.0000       752.0000       976.0000

    8489664.      0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00
  0.0000000E+00  0.0000000E+00  0.0000000E+00

    350.0000      -390.0000      -198.0000      -186.0000       218.0000
  -254.0000      -222.0000      -270.0000
    350.0000      -42.00000      -158.0000      -470.0000       218.0000
  -206.0000      -190.0000      -294.0000
続行するには何かキーを押してください . . .

ここではノルムが実部のみになること、また八元数の非結合性や部分結合性等を確認しています。

module m_sc
  implicit none
  
  type :: t_sc
    class(*), allocatable :: re, im
  contains
    generic :: assignment(=) => assign_sc
    generic :: operator(-) => minus 
    generic :: operator(.cc.) => conj
    generic :: operator(+) => plus
    generic :: operator(*) => prod    
    procedure :: minus, plus, conj, prod
    generic :: write(formatted) => pr_sc
    procedure :: assign_sc, pr_sc
  end type t_sc
  
contains
  
  recursive function minus(a) result(c)
    class(t_sc), intent(in) :: a
    class(t_sc), allocatable:: c
    allocate(t_sc::c)
    select type (re => a%re)
    type is (real)
      c%re = -re
    class is (t_sc)
      allocate(c%re, source = -re)  
      select type (im => a%im)
      class is (t_sc)
      allocate(c%im, source = -im)  
      end select    
    end select
  end function minus 
  recursive function conj(a) result(c)
    class(t_sc), intent(in) :: a
    class(t_sc), allocatable:: c
    allocate(t_sc::c)
    select type (re => a%re)
    type is (real)
      c%re = re
    class is (t_sc)
      allocate(c%re, source = conj(re))  
      select type (im => a%im)
      class is (t_sc)
      allocate(c%im, source = -im)  
      end select    
    end select
  end function conj

  recursive function plus(a, b) result(c)
    class(t_sc), intent(in) :: a, b
    class(t_sc), allocatable:: c
    allocate(t_sc::c)
    select type (ar => a%re)
    type is (real)
      select type (br => b%re)
      type is (real)
      c%re = ar + br
      end select
    class is (t_sc)
      select type (br => b%re)
      class is (t_sc)
      allocate(c%re, source = ar + br)
      end select
      select type (ai => a%im)
      class is (t_sc)
      select type (bi => b%im)
      class is (t_sc)
      allocate(c%im, source = ai + bi)
      end select
      end select
    end select
  end function plus
  
  recursive function prod(a, b) result(c)
    class(t_sc), intent(in) :: a, b
    class(t_sc), allocatable:: c
    allocate(t_sc::c)
    select type (ar => a%re)
    type is (real)
      select type (br => b%re)
      type is (real)
      c%re = ar * br
      end select
    class is (t_sc)
      select type (br => b%re)
      class is (t_sc)
      select type (ai => a%im)
      class is (t_sc)
      select type (bi => b%im)
      class is (t_sc)
      allocate(c%re, source = ar * br + conj(-bi) * ai)
      allocate(c%im, source = bi * ar + ai * conj(br) )
      end select
      end select
      end select
    end select
  end function prod
  
  subroutine assign_sc(c, a)
    class(t_sc), intent(out) :: c
    real, intent(in) :: a(:)
    allocate(t_sc::c%re)
    allocate(t_sc::c%im)
    select case (size(a))
    case (1) ! R
      allocate(c%re, source = a(1))  
    case (2) ! C    
      select type (cr => c%re)
      type is (t_sc)
        cr = a(1:1)
      end select  
      select type (ci => c%im)
        type is (t_sc)  
        ci = a(2:2)
      end select
    case (4) ! H    
      select type (cr => c%re)
      type is (t_sc)
        cr = a(1:2)
      end select  
      select type (ci => c%im)
        type is (t_sc)  
        ci = a(3:4)
      end select
    case (8) ! O    
      select type (cr => c%re)
      type is (t_sc)
        cr = a(1:4)
      end select  
      select type (ci => c%im)
        type is (t_sc)  
        ci = a(5:8)
      end select
    case default
      stop 'assign_c: error!'
    end select    
  end subroutine assign_sc
  
  recursive subroutine pr_sc(dtv, unit, iotype, vlist, iostat, iomsg)
    class(t_sc), intent(in) :: dtv
    integer, intent(in) :: unit
    character (len = *), intent(in) :: iotype
    integer, intent(in) :: vlist(:)
    integer, intent(out) :: iostat
    character (len = *), intent(in out) :: iomsg
    select type (re => dtv%re)
    type is (real)
      write(unit, *, iostat = iostat) re
    type is (t_sc)
      call re%pr_sc(unit, iotype, vlist, iostat, iomsg)    
      select type (im => dtv%im)
      type is (t_sc)    
      call im%pr_sc(unit, iotype, vlist, iostat, iomsg)    
      end select
    end select      
  end subroutine pr_sc
end module m_sc


program hello
  use m_sc
  implicit none
  class(t_sc), allocatable :: ra, rb, rc, ca, cb, cc 
  class(t_sc), allocatable :: ha, hb, hc, oa, ob, oc, od 
  allocate(ra, rb)
  
  ra = [1.0]
  rb = [2.0]
  rc = -ra + (-rb)
  print *, rc
  
  allocate(ca, cb)
  ca = [1.0, 2.0]
  cb = [3.0, 4.0]
  cc = ca * cb
  print *, cc
  print *, (1.0, 2.0) * (3.0, 4.0)
!  print *, conj(cc)
  
  allocate(ha, hb, oa, ob)
  ha = [1.0, 2.0, 3.0, 4.0]
  hb = [5.0, 6.0, 7.0, 8.0]
  hc = ha * hb
  print *, hc
  print *, conj(hc)
  
  oa = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]
  ob = [5.0, 6.0, 7.0, 8.0, -1.0, -2.0, -3.0, -4.0]
  oc = (oa * oa) * ob
  od = oa * (oa * ob)
  print *, oc
  print *, od
  print *
  oc = (oa * ob) * oa
  od = oa * (ob * oa)
  print *, oc
  print *, od
  print *
  oc = (ob * oa) * oa
  od = ob * (oa * oa)
  print *, oc
  print *, od
  print *
  print *, conj(od) * od
  print *

  oa = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]
  ob = [5.0, 6.0, 7.0, 8.0, -1.0, -2.0, -3.0, -4.0]
  oc = [1.0, 1.0, 1.0, 1.0, -2.0, -2.0, -1.0, -1.0]

  print *, (oa * ob) * oc
  print *, oa * (ob * oc)    
end program Hello

a(:) は必ずしも a ならず

配列の全配列操作を記述するとき、スカラーと配列を明示的に区別するために、配列名のあとに (:) をつけることが推奨されることがありますが、文脈によって微妙に意味が違う場合もあり、また副プログラム呼び出しの引数にしたときの挙動など、コンパイラの実装にも依存するので、時には注意が必要かと思います。

特に Fortran 2003 で、allocatable 属性の配列に代入するときの挙動がそれまでのものと変わって、自動再割り付けが行われたりするので危険です。

Cray が 2005 年に出した、Fortran 2003 and Beyond という Fortran 2003 規格と Fortran 2008 規格案を紹介した文書に、いい実例が出ていたので少し変えて紹介します。

ソース・プログラム

(:) があるとたとえ全部分を指定していても部分配列となり、自動再割り付けは行われません。（文字列の場合も同様です）

    program Console6
      implicit none
      block 
        character(len = :), allocatable :: string
        allocate(character(16)::string)
        string = '0123456789abcdef'
        !
        print *, len(string), ':', string, ':' ! length = 16
        string(:) = 'pad'
        print *, len(string), ':', string, ':' ! length = 16
        string = 'short'
        print *, len(string), ':', string, ':' ! lenght =  5
      end block
      !
      print *
      !
      block
        real, allocatable :: a(:), b(:), c(:)
        allocate(a(10), b(20))
        a = 1.10
        b = 1.20
        c = a
        print *, size(c)  ! size(c) = 10
        c = b
        print *, size(c)  ! size(c) = 20
        c(:) = a(:)       ! error size-mismatch
      end block
    end program Console6

実行結果

          16 :0123456789abcdef:
          16 :pad             :
           5 :short:

          10
          20
forrtl: severe (408): fort: (2): Subscript #1 of the array A has value 11 which is greater than the upper bound of 10

Image              PC        Routine            Line        Source
libifcoremdd.dll   5A1B4138  Unknown               Unknown  Unknown
Console6.exe       00C323AF  _MAIN__                    26  Console6.f90
Console6.exe       00C327DF  Unknown               Unknown  Unknown
Console6.exe       00C3521E  Unknown               Unknown  Unknown
Console6.exe       00C35100  Unknown               Unknown  Unknown
Console6.exe       00C34FAD  Unknown               Unknown  Unknown
Console6.exe       00C35238  Unknown               Unknown  Unknown
KERNEL32.DLL       73F68744  Unknown               Unknown  Unknown
ntdll.dll          774C587D  Unknown               Unknown  Unknown
ntdll.dll          774C584D  Unknown               Unknown  Unknown
続行するには何かキーを押してください . . .

４・２３

こまった時はシェイクスピア

今日はシェークスピアの誕生日にして命日ではありませんか。

四月馬鹿 「How Thou Canst Maketh a Fine Program in Fortran」

今更ですが、四月馬鹿のジョーク記事で、18世紀風？の古雅な英語で書かれた Fortran インストールおよび、”Hello World!" ならぬ "Good morrow, and well met, O world!" 表示までの紹介記事がありました。

www.digitalocean.com

ちょっと面白かったです。作者はクラウドコンピューティングサービス会社のテクニカルライターの方のようです。

Fortran 2008 の pointer function で作った連想配列

Fortran 2008 の新機能 pointer function を用いて、簡単な連想配列 (associative array) を作ってみます。

ここでは古典的な手法で、文字列をインデックスとした整数型連想配列を作ってみます。

まず文字列を適当な整数に変換し、あらかじめ用意した配列のサイズの剰余をとって、それを実際の引数とします。このとき、単射性が無いため数値が重なることがあります（いわゆる衝突）、これが起きた場合は、適当な数を足して空いている場所が見つかるまでずらします。

文字列を数値に変換するために、まず文字列を文字配列に直し、次にそれをASCIIコードに直し、総和を取って、適当な素数で剰余を取っています。また衝突時にずらす数も素数に取れば、全ての空きを探すことになります（全射）。

ソース・プログラム

    module m_hash
      implicit none
     ! private
      integer, parameter :: nhash = 17, nd = 13
      type :: t_key
        character(len = :), allocatable :: key  
      end type t_key
      type (t_key) :: keyc(nhash)  
      integer, target :: vals(nhash) = 0  
    contains
      integer function ihash(text)
        character(len = *), intent(in) :: text
        ihash = mod(sum(iachar(transfer(text, ' ', size = len_trim(text)))), nhash) + 1
      end function ihash  
       
      function ia(text) result(ires)
        character(len = *), intent(in) :: text
        integer, pointer :: ires 
        integer :: key, loc, i
        key = ihash(text)
        do i = 1, nhash
          if (.not. allocated(keyc(key)%key)) then
            keyc(key)%key = trim(text)
            ires => vals(key)
            exit
          else if (keyc(key)%key == trim(text)) then 
            ires => vals(key)
            exit
          else ! collision
            key = mod(key + nd, nhash) + 1  
          end if
        end do  
        if (i > nhash) stop 'associative array exhausted!' 
      end function ia      
    end module m_hash
    
    program Hash
      use m_hash
      implicit none
      
      ia('a') = 41
      print *, ia('a')
      ia('a') = 100
      print *, ia('a')
      ia('b') = 200
      print *, ia('a') + ia('b')
      
      
      ia('FORTRAN77') = 1978
      ia('Fortran90') = 1991
      ia('Donald Trump') = 1
      ia('Steve Bannon') = 4
      ia('Milo Yiannopoulos') = 10
      
      block
        integer :: i
        do i = 1, nhash
          print '(i3,a,a20,a,i10)', i, ':', keyc(i)%key, '=>', ia(keyc(i)%key)
        end do    
      end block
    end program Hash

実行結果

          41
         100
         300
  1:                    =>         0
  2:           Fortran90=>      1991
  3:                    =>         0
  4:                    =>         0
  5:           FORTRAN77=>      1978
  6:                    =>         0
  7:        Donald Trump=>         1
  8:                    =>         0
  9:                    =>         0
 10:   Milo Yiannopoulos=>        10
 11:                    =>         0
 12:                    =>         0
 13:                   a=>       100
 14:                   b=>       200
 15:                    =>         0
 16:                    =>         0
 17:        Steve Bannon=>         4
続行するには何かキーを押してください . . .

変数っぽく扱える pointer function

Fortran 2008 の新機能に、pointer を返り値としてもつ関数が、あたかも変数のように扱えるというものがあります。Modern Fortran Explained の 20.5.2 pointer functions denoting variables や The new features of Fortran 2008 の 6.2 Pointer functions に記述があります。

連想配列のようなものに適しているのではないかと思います。

実例１

    module m_sub
      implicit none
    contains
      function storage(key) result(loc)
        integer, intent(in) :: key
        integer, pointer :: loc
        integer, target :: m(100) = 0
        loc => m(key)
      end function storage
    end module m_sub

    program Console2
      use m_sub
      implicit none
      integer :: i
      do i = 1, 10
        storage(i) = i**2
      end do
      do i = 1, 10
        print *, storage(i)
      end do        
    end program Console2

実行例１

           1
           4
           9
          16
          25
          36
          49
          64
          81
         100
続行するには何かキーを押してください . . .

実例２

フィボナッチ数列を再帰で求めますが、一度求めた値はテーブルに記録して再利用します。フィボナッチ数の計算では意味がないのですが、テーブルのインデックスはとびとびの値でよくなっています。

二重再帰を使っています。

    module m_fib
      implicit none
      integer, private, parameter :: nmax = 1000
      integer, private :: nx = 0
      integer, private :: keys(nmax) = -huge(0)  
      integer, private, target :: vals(nmax) = 0  
    contains
      recursive function fib(n) result(ires)
        integer, intent(in) :: n
        integer :: ires
        select case(n)
          case (1:2)
            ires = 1
          case (3:)
            ires = fib_table(n)
          case default
            ires = 0
        end select     
      end function fib
    
      recursive function fib_table(key) result(ires)
        integer, intent(in) :: key
        integer, pointer :: ires 
        integer :: loc
        loc = findloc(keys(:nx), key, dim = 1) 
        if (loc <= 0) then 
          nx = nx + 1  
          keys (nx) = key
          vals(nx) = fib(key - 1) + fib(key - 2)
          loc = nx
        end if
        ires => vals(loc)
      end function fib_table
    end module m_fib
    
    program fibonacci
      use m_fib
      implicit none
      integer :: i
      do i = 1, 40
        print *, i, fib(i)
      end do   
    end program fibonacci

実行結果

           1           1
           2           1
           3           2
           4           3
           5           5
           6           8
           7          13
           8          21
           9          34
          10          55
          11          89
          12         144
          13         233
          14         377
          15         610
          16         987
          17        1597
          18        2584
          19        4181
          20        6765
          21       10946
          22       17711
          23       28657
          24       46368
          25       75025
          26      121393
          27      196418
          28      317811
          29      514229
          30      832040
          31     1346269
          32     2178309
          33     3524578
          34     5702887
          35     9227465
          36    14930352
          37    24157817
          38    39088169
続行するには何かキーを押してください . . .

参考　配列を利用したフィボナッチ計算結果再利用型

    module m_fib2
      implicit none
    contains
      recursive function fib(n) result(ires)
        integer, intent(in) :: n
        integer :: ires  
        select case(n)
        case (1:2)
          ires = 1
        case (3:)  
          ires = fib(n - 1) + fib(n - 2)
        case default
          ires = 0
        end select  
      end function fib

      recursive function fib4(n) result(ires)
        integer, intent(in) :: n
        integer :: ires
        integer, save :: stor(100) = 0
        select case(n)
          case(1:2)  
            ires = 1
          case(3:)  
            if (stor(n) == 0) stor(n) = fib4(n - 1) + fib4(n - 2)
            ires = stor(n)
          case default
            ires = 0
        end select    
      end function fib4
    end module m_fib2
    
    program fibonacci
      use m_fib2
      implicit none
      integer :: i
      do i = 1, 38
        print *, i, fib(i)
      end do  
      pause
      do i = 1, 38
        print *, i, fib4(i)
      end do   
    end program fibonacci

派生型の再帰的割り付け成分

Fortran 2008 では、派生型の成分に自分自身を allocatable 属性で持てます。つまり再帰的に派生型を定義できます。再帰的定義は古典的なリスト構造によく使われます。

Fortran 2003 までは、派生型の再帰的な成分は pointer 型に限られていました。pointer 型の成分に記憶領域の割り付けをする場合、解放時の処理を利用者側が自分でやらなければなりませんでした。記憶領域の解放漏れが起きないよう慎重に解放処理を用意する必要がありました。

ところが allocatable 属性の場合は、記憶領域解放時の処理を、処理系側が引き受けてくれるので、利用者側の負担が大幅に減ります。

根元を解放すれば、そこから連なる子・孫・末代まで皆、処理系側が解放してくれます。根元を解放するとそこから連なる記憶領域がアクセス不能の死に領域になってしまう pointer 型との大きな違いです。

ソースプログラム

古典的な一次元単方向リストを定義して、整数値を入れてゆきます。つぎに根元からリストをたぐって代入した整数値を出力します。最後に、リストの根元を解放します。この時、ファイナライザ（いわゆるデストラクタ）を用意して、いまわの際に成分を出力させています。

    module m_final
      type :: t_list
        integer :: ival
        type (t_list), allocatable :: next
      contains
        final :: destruct
      end type t_list  
    contains
      subroutine destruct(this)
        type (t_list), intent(in) :: this
        print *, this%ival
      end subroutine destruct
    end module m_final
    
    program final
      use m_final
      implicit none
      type (t_list), allocatable, target :: root
      type (t_list), pointer :: last
      integer :: i
      allocate(root, source = t_list(0))
      last => root
      do i = 1, 10
        allocate(last%next, source = t_list(i))
        last => last%next
      end do    
      last => root
      
      print *, 'print list'
      do 
        if (.not. associated(last)) exit
        print *, last%ival
        last => last%next
      end do  
      
      print *, 'finalizer'
      deallocate(root)
    end program final

実行結果

ファイナライズする順番は枝葉の末端から根元に向かってではなく、根元から枝葉の末端に向かって解放されてゆくようです。

 print list
           0
           1
           2
           3
           4
           5
           6
           7
           8
           9
          10
 finalizer
           0
           1
           2
           3
           4
           5
           6
           7
           8
           9
          10
続行するには何かキーを押してください . . .

参考
・M.Metcalf, J.Reid and M.Cohen: Modern Fortran Explained　20.3.1
・J.Reid: The new features of Fortran 2008