fortran66のブログ

fortran について書きます。

非結合代数　八元数 (octonion) 計算

八元数 (Octonion)

Octonion とは八個のタマネギのことではなく、八元数を意味します。

f:id:fortran66:20120313142501j:plain
タマネギ部隊

八元数は非結合代数の中でも簡単なものされていて、加減乗除の出来る数となっています。

加減乗除のできる数は、実数・複素数・四元数・八元数ときて打ち止めとなっています。ケーリー・ディクソン構成によって拡張して作ってゆけます。

ここで複素数は共役が一致せず、四元数はさらに交換則が成り立たず、八元数ではそのうえ結合則も成り立たなくなっています。

fortran66.hatenablog.com

自然界でいえば、四元数はスピンと関係しており、八元数はリー群の例外群に関係して、その中の G2 群は希土類やウラン系元素のｆ 電子軌道に現れています。

計算

ここでは、以前作ったルーチンを拡張して、チェックの計算をしてみます。

八元数公式のチェック

まずは八元数の基本的な公式が成り立つことを確かめます。その後ネットで検索して出てきた計算を再現してみます。

四元数の交換則からのズレの平均www.johndcook.com

SF 作家のグレッグ・イーガンがコメント欄に現れて、解析解を求めています。（グレッグ・イーガンの小説は読んだことありませんｗ　同じ数学系 SF 作家のルディ・ラッカーのホワイトライトは読んだことありますが、内容の記憶はありませんｗ）

順列都市〔上〕

順列都市〔上〕

順列都市〔下〕

順列都市〔下〕

ホワイト・ライト (ハヤカワ文庫SF)

ホワイト・ライト (ハヤカワ文庫SF)

10^7 回平均

f:id:fortran66:20190215031444p:plain
四元数交換則からのズレ平均 （横軸x100)

八元数の結合則からのズレの平均

www.johndcook.com

f:id:fortran66:20190215031549p:plain
八元数結合則からのズレ平均 （横軸x100）

上記平均計算では、正規分布に基づく乱数を基に規格化して単位ベクトルの生成を行っています。正規分布を用いると多次元球の表面付近にベクトルが生成されるからいいのかもしれません。

一様分布で計算すると平均も分布も微妙にずれます。一様分布を多次元化すると、対角方向に分布が集中していると思われるので、まぁもっともかもしれません。

今後の予定としては、G2 群の表現行列を見てみることにします。

群と位相 (基礎数学選書 (5))

群と位相 (基礎数学選書 (5))

群と表現 (基礎数学選書 (10))

群と表現 (基礎数学選書 (10))

例外型単純リー群

例外型単純リー群

四元数と八元数―幾何、算術、そして対称性

四元数と八元数―幾何、算術、そして対称性

実行結果

 |a*b|=|a||b|   62016.0000000000        62016.0000000000
 (au,av)=(a,a)(u,v)=(ua,va)  -432.000000000000       -432.000000000000
  -432.000000000000
 (au,bv)+(bu,av)=2(a,b)(u,v)  -4032.00000000000       -4032.00000000000
 ~~u = u
 ~~u      :   0.000   2.000  -1.000   1.000  -2.000   2.000  -1.000  -2.000
   u      :   0.000   2.000  -1.000   1.000  -2.000   2.000  -1.000  -2.000

 ~(u+v)=~u+~v
 ~(u+v)   :   0.000  -1.000   0.000   1.000   1.000  -1.000  -1.000   3.000
 ~u + ~v  :   0.000  -1.000   0.000   1.000   1.000  -1.000  -1.000   3.000

 (au,v)=(u,~av)
 (au,  v) : -29.000
 (u ,~av) : -29.000

 (u,v)=(v,u)=1/2(~uv+~vu)=1/2(u~v+v~u)
 (u,v)    :  -9.000
 (v,u)    :  -9.000
 1/2(~uv+~vu):  -9.000
 1/2(u~v+v~u):  -9.000

 a(~au) = (a~a)u
 a(~au)   :  -48.000 -96.000-144.000-192.000 240.000 288.000 336.000 384.000
 (a~a)u   :  -48.000 -96.000-144.000-192.000 240.000 288.000 336.000 384.000

 associator: (a, a, u)=(a, u, a)=(u, a, a)=0
 (a, a, u) :   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000
 (a, u, a) :   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000
 (u, a, a) :   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000

 associator: (a, ~a, u)=(a, u, ~a)=(u, a, ~a)=0
 (a, ~a, u):   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000
 (a, u, ~a):   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000
 (u, a, ~a):   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000

 associator: (a, b, c)=(b, c, a)=(c, a, b)=-(b, a, c)
 (a, b, c) :   0.000-440.000 496.000 -72.000  64.000-640.000 256.000 256.000
 (b, c, a) :   0.000-440.000 496.000 -72.000  64.000-640.000 256.000 256.000
 (c, a, b) :   0.000-440.000 496.000 -72.000  64.000-640.000 256.000 256.000
 (b, a, c) :   0.000 440.000-496.000  72.000 -64.000 640.000-256.000-256.000

 ~b(au)+~a(bu)=2(a,b)u=(ua)~b+(ub)~a
 ~b(au)+~a(bu):   0.000 896.000-448.000 448.000-896.000 896.000-448.000-896.000
    2(a,b)u   :   0.000 896.000-448.000 448.000-896.000 896.000-448.000-896.000
 (ua)~b+(ub)~a:   0.000 896.000-448.000 448.000-896.000 896.000-448.000-896.000

 Maufang: (au)(va)=a(uv)a
 (au)(va)  :-374.000 -76.000 -18.000 328.000 462.000 -28.000-230.000-220.000
  a(uv)a   :-374.000 -76.000 -18.000 328.000 462.000 -28.000-230.000-220.000

 cross product: Re(u)=Re(v)=0; u X v = uv + (u,v)
 u X v     :   0.000  -4.000  -4.000   2.000   6.000  -3.000  -6.000  -7.000
 uv + (u,v):   0.000  -4.000  -4.000   2.000   6.000  -3.000  -6.000  -7.000

 commutator   1.13197636058180
 associator   1.09479514480162
続行するには何かキーを押してください . . .

ソース・プログラム

    module m_comp
      implicit none
      integer, parameter :: kd = kind(0.0d0)
      
      interface operator(+) 
        procedure :: add_r, add_c, add_h, add_o
      end interface
      interface operator(*) 
        procedure :: mul_r, mul_c, mul_h, mul_o
      end interface
      interface operator(/) 
        procedure :: scale_r, scale_c, scale_h, scale_o
        procedure :: div_r, div_c, div_h, div_o
      end interface
      interface operator(-) 
        procedure :: minus_r, minus_c, minus_h, minus_o
        procedure :: sub_r, sub_c, sub_h, sub_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
      interface dot 
        procedure :: dot_r, dot_c, dot_h, dot_o
      end interface
      interface cross 
        procedure :: cross_r, cross_c, cross_h, cross_o
      end interface
      interface im 
        procedure :: im_r, im_c, im_h, im_o
      end interface

      type :: t_r
        real(kd) :: re
      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%re = a%re 
      end function conj_r
      type(t_r) function im_r(a) result(c)
        type(t_r), intent(in) :: a
        c%re = 0.0_kd
      end function im_r
      type(t_r) function minus_r(a) result(c)
        type(t_r), intent(in) :: a
        c%re = -a%re 
      end function minus_r
      type(t_r) function add_r(a, b) result(c)
        type(t_r), intent(in) :: a, b
        c%re = a%re + b%re
      end function add_r
      type(t_r) function sub_r(a, b) result(c)
        type(t_r), intent(in) :: a, b
        c%re = a%re - b%re
      end function sub_r
      type(t_r) function div_r(a, b) result(c)
        type (t_r), intent(in) :: a, b
        c = a * conj(b) / dot(b, b)
      end function div_r
      type(t_r) function mul_r(a, b) result(c)
        type (t_r), intent(in) :: a, b
        c%re = a%re * b%re
      end function mul_r
      type(t_r) function scale_r(a, b) result(c)
        type (t_r), intent(in) :: a
        real(kd), intent(in) :: b
        c%re = a%re / b
      end function scale_r
      real(kd) function dot_r(a, b) result(c)
        type (t_r), intent(in) :: a, b
        c = a%re * b%re
      end function dot_r
      type(t_r) function cross_r(a, b) result(c)
        type (t_r), intent(in) :: a, b
        c = (conj(b) * a - conj(a) * b) / 2.0_kd
      end function cross_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 im_c(a) result(c)
        type(t_c), intent(in) :: a
        c%re = im(a%re)
        c%im = a%im
      end function im_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 add_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 add_c
      type(t_c) function sub_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 sub_c
      type(t_c) function mul_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 mul_c
      type(t_c) function div_c(a, b) result(c)
        type (t_c), intent(in) :: a, b
        c = a * conj(b) / dot(b, b)
      end function div_c
      type(t_c) function to_c(r0, r1)
        real(kd), intent(in) :: r0, r1
        to_c = t_c(t_r(r0), t_r(r1))
      end function to_c  
      type(t_c) function scale_c(a, b) result(c)
        type (t_c), intent(in) :: a
        real(kd), intent(in) :: b
        c%re = a%re / b
        c%im = a%im / b
      end function scale_c
      real(kd) function dot_c(a, b) result(c)
        type (t_c), intent(in) :: a, b
        c = dot_r(a%re, b%re) + dot_r(a%im, b%im)
      end function dot_c
      type(t_c) function cross_c(a, b) result(c)
        type (t_c), intent(in) :: a, b
        c = (conj(b) * a - conj(a) * b) / 2.0_kd
      end function cross_c
      
      type(t_h) function conj_h(a) result(c)
        type(t_h), intent(in) :: a
        c%re = conj(a%re) 
        c%im =     -a%im 
      end function conj_h
      type(t_h) function im_h(a) result(c)
        type(t_h), intent(in) :: a
        c%re = im(a%re)
        c%im = a%im
      end function im_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 add_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 add_h
      type(t_h) function sub_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 sub_h
      type(t_h) function mul_h(a, b) result(c)
        type (t_h), intent(in) :: a, b
        c%re = a%re * b%re - conj(b%im) * a%im
        c%im = b%im * a%re + a%im * conj(b%re)
      end function mul_h
      type(t_h) function div_h(a, b) result(c)
        type (t_h), intent(in) :: a, b
        c = a * conj(b) / dot(b, b)
      end function div_h
      type(t_h) function to_h(r0, r1, r2, r3)
        real(kd), intent(in) :: r0, r1, r2, r3
        to_h = t_h(to_c(r0, r1), to_c(r2, r3))
      end function to_h
      type(t_h) function scale_h(a, b) result(c)
        type (t_h), intent(in) :: a
        real(kd), intent(in) :: b
        c%re = a%re / b
        c%im = a%im / b
      end function scale_h
      real(kd) function dot_h(a, b) result(c)
        type (t_h), intent(in) :: a, b
        c = dot_c(a%re, b%re) + dot_c(a%im, b%im)
      end function dot_h
      type(t_h) function cross_h(a, b) result(c)
        type (t_h), intent(in) :: a, b
        c = (conj(b) * a - conj(a) * b) / 2.0_kd
      end function cross_h
      type(t_h) function to_h_arr(a) result(c)
        real(kd), intent(in) :: a(4)
        c = to_h(a(1), a(2), a(3), a(4))
      end function to_h_arr
      
      type(t_o) function conj_o(a) result(c)
        type(t_o), intent(in) :: a
        c%re = conj(a%re) 
        c%im =     -a%im 
      end function conj_o
      type(t_o) function im_o(a) result(c)
        type(t_o), intent(in) :: a
        c%re = im(a%re)
        c%im = a%im
      end function im_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 add_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 add_o
      type(t_o) function sub_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 sub_o
      type(t_o) function mul_o(a, b) result(c)
        type (t_o), intent(in) :: a, b
        c%re = a%re * b%re - conj(b%im) * a%im
        c%im = b%im * a%re + a%im * conj(b%re)
      end function mul_o
      type(t_o) function div_o(a, b) result(c)
        type (t_o), intent(in) :: a, b
        c = a * conj(b) / dot(b, b)
      end function div_o
      type(t_o) function scale_o(a, b) result(c)
        type (t_o), intent(in) :: a
        real(kd), intent(in) :: b
        c%re = a%re / b
        c%im = a%im / b
      end function scale_o
      type(t_o) function cross_o(a, b) result(c)
        type (t_o), intent(in) :: a, b
        c = (conj(b) * a - conj(a) * b) / 2.0_kd
      end function cross_o
          
      type(t_o) function arr8_to_o(r)
        real(kd), intent(in) :: r(8)
        arr8_to_o = t_o(to_h(r(1), r(2), r(3), r(4)), to_h(r(5), r(6), r(7), r(8)))
      end function arr8_to_o

      type(t_o) function to_o(r0, r1, r2, r3, r4, r5, r6, r7)
        real(kd), 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  
      type(t_o) function r7_to_o(a) result(b)
        real(kd), intent(in) :: a(7)
        b = t_o(to_h(0.0_kd, a(1), a(2), a(3)), to_h(a(4), a(5), a(6), a(7)))
      end function r7_to_o  
      function o_to_r7(a) result(b)
        type(t_o), intent(in) :: a
        real(kd) :: b(7)
        associate(hrcr=>a%re%re, hrci=>a%re%im, hicr=>a%im%re, hici=>a%im%im)
          b = [hrcr%im%re, hrci%re%re, hrci%im%re, hicr%re%re, hicr%im%re, hici%re%re, hici%im%re]
        end associate
      end function o_to_r7  
      real(kd) function dot_o(a, b) result(c)
        type (t_o), intent(in) :: a, b
        c = dot_h(a%re, b%re) + dot_h(a%im, b%im)
      end function dot_o
      
      type (t_o) function associator_o(a, b, c) result(d)
        type (t_o), intent(in) :: a, b, c
        d = (a * b) * c - a * (b * c)  
      end function associator_o
      
    end module m_comp
    
    program supercomplex
      use m_comp
      implicit none
      type(t_r) :: ra, rb
      type(t_o) :: oa, ob, oc, ou, ov, ow, oz1, oz2 

      oa = to_o(1.0_kd,  2.0_kd, 3.0_kd,  4.0_kd,  3.0_kd,  2.0_kd,  1.0_kd,  2.0_kd)
      ob = to_o(9.0_kd, 10.0_kd,11.0_kd, 12.0_kd, 13.0_kd, 14.0_kd, 15.0_kd, 16.0_kd)
      ou = to_o(0.0_kd,  2.0_kd,-1.0_kd,  1.0_kd, -2.0_kd,  2.0_kd, -1.0_kd, -2.0_kd)
      ov = to_o(0.0_kd, -1.0_kd, 1.0_kd, -2.0_kd,  1.0_kd, -1.0_kd,  2.0_kd, -1.0_kd)
      ow = to_o(0.0_kd, -4.0_kd, 2.0_kd,  1.0_kd,  3.0_kd, -5.0_kd,  9.0_kd,  6.0_kd)
     
      oc = oa * ob
      print *, '|a*b|=|a||b|', dot(oc, oc), dot(oa, oa) * dot(ob,ob)
      
      oc = to_o(-1.0_kd,  -2.0_kd,  -3.0_kd,  -4.0_kd,  5.0_kd,  6.0_kd,  7.0_kd,  8.0_kd)
      print *, '(au,av)=(a,a)(u,v)=(ua,va)', dot(oa*ou, oa*ov), dot(oa, oa) * dot(ou, ov), dot(ou*oa, ov*oa)
      print *, '(au,bv)+(bu,av)=2(a,b)(u,v)', dot(oa*ou, ob*ov) + dot(ob*ou, oa*ov), 2 * dot(oa, ob)*dot(ou,ov) 

      print *, '~~u = u'
      print '(a, 8f8.3)', ' ~~u      :' , conj(conj(ou))
      print '(a, 8f8.3)', '   u      :' , ou
      print *

      print *, '~(u+v)=~u+~v'
      print '(a, 8f8.3)', ' ~(u+v)   :', conj(ou + ov)
      print '(a, 8f8.3)', ' ~u + ~v  :', conj(ou)+conj(ov)
      print *

      print *, '(au,v)=(u,~av)'
      print '(a, 8f8.3)', ' (au,  v) :', dot(oa*ou, ov)
      print '(a, 8f8.3)', ' (u ,~av) :', dot(ou, conj(oa)*ov)
      print *
      
      print *, '(u,v)=(v,u)=1/2(~uv+~vu)=1/2(u~v+v~u)'
      print '(a, 8f8.3)', ' (u,v)    :', dot(ou, ov)
      print '(a, 8f8.3)', ' (v,u)    :', dot(ov, ou)
      oz1 = conj(ou)*ov+conj(ov)*ou
      oz2 = ou*conj(ov)+ov*conj(ou)
      ra = oz1%re%re%re
      rb = oz2%re%re%re
      print '(a, 8f8.3)', ' 1/2(~uv+~vu):', ra%re/2
      print '(a, 8f8.3)', ' 1/2(u~v+v~u):', rb%re/2
      print *

      print *, 'a(~au) = (a~a)u '
      print '(a, 8f8.3)', ' a(~au)   : ', oa * (conj(oa) * oc)
      print '(a, 8f8.3)', ' (a~a)u   : ', (oa * conj(oa)) * oc
      print *
      
      print *, 'associator: (a, a, u)=(a, u, a)=(u, a, a)=0'
      print '(a, 8f8.3)', ' (a, a, u) :', associator_o(oa, oa, ou) 
      print '(a, 8f8.3)', ' (a, u, a) :', associator_o(oa, ou, oa) 
      print '(a, 8f8.3)', ' (u, a, a) :', associator_o(ou, oa, oa) 
      print *
      
      print *, 'associator: (a, ~a, u)=(a, u, ~a)=(u, a, ~a)=0'
      print '(a, 8f8.3)', ' (a, ~a, u):', associator_o(oa, conj(oa), ou) 
      print '(a, 8f8.3)', ' (a, u, ~a):', associator_o(oa, ou, conj(oa)) 
      print '(a, 8f8.3)', ' (u, a, ~a):', associator_o(ou, oa, conj(oa)) 
      print *
      
      print *, 'associator: (a, b, c)=(b, c, a)=(c, a, b)=-(b, a, c)'
      print '(a, 8f8.3)', ' (a, b, c) :', associator_o(oa, ob, oc) 
      print '(a, 8f8.3)', ' (b, c, a) :', associator_o(ob, oc, oa) 
      print '(a, 8f8.3)', ' (c, a, b) :', associator_o(oc, oa, ob) 
      print '(a, 8f8.3)', ' (b, a, c) :', associator_o(ob, oa, oc) 
      print *
      
      print *, '~b(au)+~a(bu)=2(a,b)u=(ua)~b+(ub)~a' 
      print '(a, 8f8.3)', ' ~b(au)+~a(bu):', conj(ob)*(oa*ou)+conj(oa)*(ob*ou) 
      print '(a, 8f8.3)', '    2(a,b)u   :', to_o(2*dot(oa,ob),0.0_kd,0.0_kd,0.0_kd,0.0_kd,0.0_kd,0.0_kd,0.0_kd)*ou 
      print '(a, 8f8.3)', ' (ua)~b+(ub)~a:', (ou*oa)*conj(ob)+(ou*ob)*conj(oa)
      print *
      
      print *, 'Maufang: (au)(va)=a(uv)a'
      print '(a, 8f8.3)', ' (au)(va)  :', (oa*ou)*(ov*oa)
      print '(a, 8f8.3)', '  a(uv)a   :', oa*(ou*ov)*oa
      print *
      
      print *, 'cross product: Re(u)=Re(v)=0; u X v = uv + (u,v)  '
      print '(a, 8f8.3)', ' u X v     :', cross(ou, ov)
      print '(a, 8f8.3)', ' uv + (u,v):', ou*ov + to_o(dot(ou,ov), 0.0_kd,0.0_kd,0.0_kd,0.0_kd,0.0_kd,0.0_kd,0.0_kd) 
      print *
      
      open(9, file='quartcomm.csv')
      call random_seed()
      block
         real, parameter :: pi = 4 * atan(1.0_kd) 
         type (t_h) :: h(2), com         
         integer :: i, n = 10**7, m(201)
         real(kd) :: x(4), y(4), p(4), q(4), r1(4), r2(4), s, t
         m = 0
         do i = 1, n
           call random_number(x)
           call random_number(y)
           p = sqrt(-2.0_kd*log(x))
           q = 2 * pi * y
           r1 = p * cos(q)
           r2 = p * sin(q)
           h(1) = to_h_arr(r1)
           h(2) = to_h_arr(r2)
           h(1) = h(1) / sqrt(dot(h(1), h(1)))
           h(2) = h(2) / sqrt(dot(h(2), h(2)))
           
           com = h(1)*h(2) - h(2)*h(1)
           t = sqrt(dot(com, com))
           s = s + t
           m(100*t+1) = m(100*t+1) + 1
         end do    
         print *, 'commutator', s / n
         
        do i = 1, size(m)
            write(9, *), i / 100.0, ',', m(i)
        end do    
      end block    
      close(9)

      open(9, file='octassoc.csv')
      block
        real, parameter :: pi = 4 * atan(1.0_kd) 
        type (t_o) :: o(3), assc
        integer :: i, j, n = 10**7, m(200)
        real(kd) :: x(24), y(24), r(24), t(24), r8(8, 3), s, w
        call random_seed()
        m = 0
        s = 0.0_kd
        do i = 1, n
          call random_number(x)
          call random_number(y)
          r = sqrt(-2*log(x))
          t = 2*pi*y
          r8 = reshape(r*cos(t), [8, 3])
          do j = 1, 3
            o(j) = arr8_to_o(r8(:, j))
            o(j) = o(j) / sqrt( dot(o(j), o(j)) )
          end do 
          assc = associator_o(o(1), o(2), o(3))
          w = sqrt(dot(assc, assc))
          s = s + w
          m(100*w+1) = m(100*w+1) + 1
        end do    
        print *, 'associator', s / n  
        do i = 1, size(m)
          write(9, *), i, ',', m(i)
        end do    
      end block
      close(9)
    end program supercomplex