ソース・プログラム
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
