fortran66のブログ

fortran について書きます。

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

背景

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

複素数 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