ソース・プログラム
rot_1~7 それぞれ２変換角パラメータが合計１４自由度を出します。ここでは面倒なので、rot_1..rot_7 の固定順番で、乱数で決めた角度で変換して、左右両辺の結果を比較しています。
八元数の積は、非結合演算なので演算順序に注意が必要です。
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) pure function conj_r(a) result(c)
type(t_r), intent(in) :: a
c%re = a%re
end function conj_r
type(t_r) pure function im_r(a) result(c)
type(t_r), intent(in) :: a
c%re = 0.0_kd
end function im_r
type(t_r) pure function minus_r(a) result(c)
type(t_r), intent(in) :: a
c%re = -a%re
end function minus_r
type(t_r) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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) pure 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
pure 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) pure 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) pure 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
module m_g2
use m_comp
implicit none
contains
pure function so7_1(x, y, a) result(b)
real(kd), intent(in) :: x, y
real(kd), intent(in) :: a(7)
real(kd) :: b(7), r(7, 7)
r = 0.0
r(1, 1) = 1.0
r(2, 2) = cos(x); r(2, 3) = sin(x)
r(3, 2) = -sin(x); r(3, 3) = cos(x)
r(4, 4) = cos(x - y); r(4, 5) = -sin(x - y)
r(5, 4) = sin(x - y); r(5, 5) = cos(x - y)
r(6, 6) = cos(y); r(6, 7) = sin(y)
r(7, 6) = -sin(y); r(7, 7) = cos(y)
b = matmul(r, a)
end function so7_1
pure function so7_2(x, y, a) result(b)
real(kd), intent(in) :: x, y
real(kd), intent(in) :: a(7)
real(kd) :: b(7), r(7, 7)
r = 0.0
r(1, 1) = cos(x); r(1, 3) = -sin(x)
r(2, 2) = 1.0
r(3, 1) = sin(x); r(3, 3) = cos(x)
r(4, 4) = cos(- x + y); r(4, 6) = sin(- x + y)
r(5, 5) = cos(y); r(5, 7) = -sin(y)
r(6, 4) = -sin(- x + y); r(6, 6) = cos(- x + y)
r(7, 5) = sin(y); r(7, 7) = cos(y)
b = matmul(r, a)
end function so7_2
pure function so7_3(x, y, a) result(b)
real(kd), intent(in) :: x, y
real(kd), intent(in) :: a(7)
real(kd) :: b(7), r(7, 7)
r = 0.0
r(1, 1) = cos(x); r(1, 2) = sin(x)
r(2, 1) = -sin(x); r(2, 2) = cos(x)
r(3, 3) = 1.0
r(4, 4) = cos(- x + y); r(4, 7) = sin(- x + y)
r(5, 5) = cos(y); r(5, 6) = sin(y)
r(6, 5) = -sin(y); r(6, 6) = cos(y)
r(7, 4) = -sin(- x + y); r(7, 7) = cos(- x + y)
b = matmul(r, a)
end function so7_3
pure function so7_4(x, y, a) result(b)
real(kd), intent(in) :: x, y
real(kd), intent(in) :: a(7)
real(kd) :: b(7), r(7, 7)
r = 0.0
r(1, 1) = cos(x); r(1, 5) = -sin(x)
r(2, 2) = cos(- x + y); r(2, 6) = -sin(- x + y)
r(3, 3) = cos(y); r(3, 7) = sin(y)
r(4, 4) = 1.0
r(5, 1) = sin(x); r(5, 5) = cos(x)
r(6, 2) = sin(- x + y); r(6, 6) = cos(- x + y)
r(7, 3) = -sin(y); r(7, 7) = cos(y)
b = matmul(r, a)
end function so7_4
pure function so7_5(x, y, a) result(b)
real(kd), intent(in) :: x, y
real(kd), intent(in) :: a(7)
real(kd) :: b(7), r(7, 7)
r = 0.0
r(1, 1) = cos(x); r(1, 4) = sin(x)
r(2, 2) = cos(x + y); r(2, 7) = sin(x + y)
r(3, 3) = cos(y); r(3, 6) = sin(y)
r(4, 1) = -sin(x); r(4, 4) = cos(x)
r(5, 5) = 1.0
r(6, 3) = -sin(y); r(6, 6) = cos(y)
r(7, 2) = -sin(x + y); r(7, 7) = cos(x + y)
b = matmul(r, a)
end function so7_5
pure function so7_6(x, y, a) result(b)
real(kd), intent(in) :: x, y
real(kd), intent(in) :: a(7)
real(kd) :: b(7), r(7, 7)
r = 0.0
r(1, 1) = cos(x - y); r(1, 7) = -sin(x - y)
r(2, 2) = cos(x); r(2, 4) = sin(x)
r(3, 3) = cos(y); r(3, 5) = sin(y)
r(4, 2) = -sin(x); r(4, 4) = cos(x)
r(5, 3) = -sin(y); r(5, 5) = cos(y)
r(6, 6) = 1.0
r(7, 1) = sin(x - y); r(7, 7) = cos(x - y)
b = matmul(r, a)
end function so7_6
pure function so7_7(x, y, a) result(b)
real(kd), intent(in) :: x, y
real(kd), intent(in) :: a(7)
real(kd) :: b(7), r(7, 7)
r = 0.0
r(1, 1) = cos(x); r(1, 6) = -sin(x)
r(2, 2) = cos(x - y); r(2, 5) = -sin(x - y)
r(3, 3) = cos(y); r(3, 4) = -sin(y)
r(4, 3) = sin(y); r(4, 4) = cos(y)
r(5, 2) = sin(x - y); r(5, 5) = cos(x - y)
r(6, 1) = sin(x); r(6, 6) = cos(x)
r(7, 7) = 1.0
b = matmul(r, a)
end function so7_7
type(t_o) pure function rot_1(x, y, a) result(b)
real(kd), intent(in) :: x, y
type(t_o), intent(in) :: a
b = r7_to_o(so7_1(x, y, o_to_r7(a)))
end function rot_1
type(t_o) pure function rot_2(x, y, a) result(b)
real(kd), intent(in) :: x, y
type(t_o), intent(in) :: a
b = r7_to_o(so7_2(x, y, o_to_r7(a)))
end function rot_2
type(t_o) pure function rot_3(x, y, a) result(b)
real(kd), intent(in) :: x, y
type(t_o), intent(in) :: a
b = r7_to_o(so7_3(x, y, o_to_r7(a)))
end function rot_3
type(t_o) pure function rot_4(x, y, a) result(b)
real(kd), intent(in) :: x, y
type(t_o), intent(in) :: a
b = r7_to_o(so7_4(x, y, o_to_r7(a)))
end function rot_4
type(t_o) pure function rot_5(x, y, a) result(b)
real(kd), intent(in) :: x, y
type(t_o), intent(in) :: a
b = r7_to_o(so7_5(x, y, o_to_r7(a)))
end function rot_5
type(t_o) pure function rot_6(x, y, a) result(b)
real(kd), intent(in) :: x, y
type(t_o), intent(in) :: a
b = r7_to_o(so7_6(x, y, o_to_r7(a)))
end function rot_6
type(t_o) pure function rot_7(x, y, a) result(b)
real(kd), intent(in) :: x, y
type(t_o), intent(in) :: a
b = r7_to_o(so7_7(x, y, o_to_r7(a)))
end function rot_7
end module m_g2
program supercomplex
use m_comp
use m_g2
implicit none
real(kd), parameter :: pi = 4 * atan(1.0)
real(kd) :: theta, tx(7), ty(7)
integer :: i
type(t_o) :: ou, ov, uv, u0, v0, uv0
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)
u0 = ou / sqrt(dot(ou, ou))
v0 = ov / sqrt(dot(ov, ov))
uv0 = u0*v0
do i = 1, 100
call random_number(tx)
call random_number(ty)
tx = tx * pi
ty = ty * pi
u0 = rot_1(tx(1), ty(1), u0)
v0 = rot_1(tx(1), ty(1), v0)
u0 = rot_2(tx(2), ty(2), u0)
v0 = rot_2(tx(2), ty(2), v0)
u0 = rot_3(tx(3), ty(3), u0)
v0 = rot_3(tx(3), ty(3), v0)
u0 = rot_4(tx(4), ty(4), u0)
v0 = rot_4(tx(4), ty(4), v0)
u0 = rot_5(tx(5), ty(5), u0)
v0 = rot_5(tx(5), ty(5), v0)
u0 = rot_6(tx(6), ty(6), u0)
v0 = rot_6(tx(6), ty(6), v0)
u0 = rot_7(tx(7), ty(7), u0)
v0 = rot_7(tx(7), ty(7), v0)
uv0 = rot_1(tx(1), ty(1), uv0)
uv0 = rot_2(tx(2), ty(2), uv0)
uv0 = rot_3(tx(3), ty(3), uv0)
uv0 = rot_4(tx(4), ty(4), uv0)
uv0 = rot_5(tx(5), ty(5), uv0)
uv0 = rot_6(tx(6), ty(6), uv0)
uv0 = rot_7(tx(7), ty(7), uv0)
print '(a, 8f8.3)', 'g(u)*g(v)', u0 * v0 + to_o(dot(u0, v0), 0.0_kd, 0.0_kd, 0.0_kd, 0.0_kd, 0.0_kd, 0.0_kd, 0.0_kd)
print '(a, 8f8.3)', ' g(uv) ', uv0
end do
end program supercomplex
