module Cayley_Dickson

!  use fft_mod
  implicit none
  private

  integer, parameter, public     :: wp = selected_real_kind(15)

  public  :: conjg, abs, real, aimag, assignment(=),   &
             operator(+), operator(-), operator(*), operator(/)
             

  type, public :: quater
    complex(kind=wp) :: c(2)
  end type

  type, public :: octonion
    type(quater) :: c(2)
  end type

  type, public :: sedenion
    type(octonion) :: c(2)
  end type

  interface conjg
    module procedure conjg_q, conjg_o, conjg_s
  end interface

  interface abs
    module procedure abs_q, abs_o
  end interface

  interface aimag
    module procedure aimag_q
  end interface

  interface real
    module procedure real_q
  end interface

  interface operator(+)
    module procedure sum_oo, sum_qq
  end interface

  interface operator(-)
    module procedure dif_oo, dif_qq, min_o, min_q
  end interface

  interface operator(*)
    module procedure mul_qq, mul_cq , mul_qc, mul_rq, mul_qr,  &
                     mul_oo, mul_qo, mul_oq, mul_co, mul_oc, mul_ro, mul_or, &
                     mul_ss
  end interface

  interface operator(/)
    module procedure div_qq, div_ss
  end interface

  interface assignment(=)
    module procedure q_equ_c, q_equ_r, o_equ_q, o_equ_c, o_equ_r,  &
                     s_equ_o, s_equ_q, s_equ_c, s_equ_r
  end interface



contains


  elemental subroutine q_equ_c(q, c)
    type(quater), intent(in out) :: q
    complex(wp), intent(in)      :: c
    q = quater([c,(0d0,0d0)])
  end

  elemental subroutine q_equ_r(q, r)
    type(quater), intent(in out)  :: q
    real(wp), intent(in)          :: r
    q%c(1) = r; q%c(2) = 0
  end

  elemental function sum_qq(p, q)  result(r)
    type(quater), intent(in)     :: p, q
    type(quater)                 :: r
    r = quater( p%c+q%c )
  end

  elemental function dif_qq(p, q)  result(r)
    type(quater), intent(in)     :: p, q
    type(quater)                 :: r
    r = quater( p%c-q%c )
  end

  elemental function min_q(q)    result(r)
    type(quater), intent(in)     :: q
    type(quater)                 :: r
    r = quater( -q%c )
  end

  elemental function mul_qq(p, q)  result(r)    !  a c - d' b,  d a + b c`
    type(quater), intent(in)     :: p, q
    type(quater)                 :: r

    r%c(1) = p%c(1) * q%c(1)  -  conjg(q%c(2)) * p%c(2)
    r%c(2) = q%c(2) * p%c(1)  +  p%c(2) * conjg(q%c(1))
  end

  elemental function mul_cq(c, q)  result(r)
    complex(wp), intent(in)      :: c
    type(quater), intent(in)     :: q
    type(quater)                 :: r
    r = quater([c * q%c(1), q%c(2) * c])         !  linear in c
  end

  elemental function mul_qc(q, c)  result(r)
    complex(wp), intent(in)      :: c
    type(quater), intent(in)     :: q
    type(quater)                 :: r
    r = quater([q%c(1) * c, q%c(2) * conjg(c)])  !  linear/anti-linear in c !!
  end

  elemental function mul_rq(x, q)  result(r)
    real(wp), intent(in)         :: x
    type(quater), intent(in)     :: q
    type(quater)                 :: r
    r = quater( x * q%c )
  end

  elemental function mul_qr(q, x)  result(r)
    real(wp), intent(in)         :: x
    type(quater), intent(in)     :: q
    type(quater)                 :: r
    r = quater( q%c * x)
  end

  elemental function div_qq(p, q)  result(r)
    type(quater), intent(in)     :: p, q
    type(quater)                 :: r

    r = p * quater([conjg(q%c(1)), -q%c(2)] / (abs(q%c(1))**2+abs(q%c(2))**2))
  end

  elemental function conjg_q(q)    result(r)
    type(quater), intent(in)     :: q
    type(quater)                 :: r

    r = quater([ conjg(q%c(1)),  -q%c(2) ])
  end

  elemental function abs_q(q)    result(r)
    type(quater), intent(in)     :: q
    real(wp)                     :: r

    r = sqrt(abs(q%c(1))**2 + abs(q%c(2))**2)
  end

  elemental function real_q(q)   result(r)
    type(quater), intent(in)     :: q
    real(wp)                     :: r
    r = real(q%c(1))
  end

  elemental function aimag_q(q)  result(r)   !??? not mapping to base algebra..
    type(quater), intent(in)     :: q
    type(quater)                 :: r

    r = quater([ q%c(1)-real(q%c(1)),  q%c(2) ])    !!! Which pure function??
  end

!------    Now octonions:     ----------------------------------------!

  elemental subroutine o_equ_q(o, q)
    type(octonion), intent(in out) :: o
    type(quater), intent(in)       :: q
    o = octonion([ q, quater(0) ])             ! quater(0) ???
  end

  elemental subroutine o_equ_c(o, c)
    type(octonion), intent(in out) :: o
    complex(wp), intent(in)        :: c
    o%c(1) = c;  o%c(2) = 0d0 
  end

  elemental subroutine o_equ_r(o, r)
    type(octonion), intent(in out) :: o
    real(wp), intent(in)           :: r
    o%c(1) = r;  o%c(2) = 0d0 
  end

  elemental function conjg_o(o)   result(r)
    type(octonion), intent(in)    :: o
    type(octonion)                :: r

    r = octonion([ conjg(o%c(1)),  -o%c(2) ])
  end

  elemental function abs_o(o)    result(r)
    type(octonion), intent(in)   :: o
    real(wp)                     :: r
    r = sqrt(abs(o%c(1))**2 + abs(o%c(2))**2)  !??? inefficient?
  end

  elemental function sum_oo(p, q)  result(r)
    type(octonion), intent(in)     :: p, q
    type(octonion)                 :: r
    r = octonion( p%c+q%c )
  end

  elemental function dif_oo(p, q)  result(r)
    type(octonion), intent(in)     :: p, q
    type(octonion)                 :: r
    r = octonion( p%c-q%c )
  end

  elemental function min_o(p)      result(r)
    type(octonion), intent(in)     :: p
    type(octonion)                 :: r
    r = octonion( -p%c )
  end

  elemental function mul_oo(p, q)  result(r)    !  a c - d' b,  d a + b c`
    type(octonion), intent(in)     :: p, q
    type(octonion)                 :: r

    r%c(1) = p%c(1) * q%c(1)  -  conjg(q%c(2)) * p%c(2)
    r%c(2) = q%c(2) * p%c(1)  +  p%c(2) * conjg(q%c(1))
  end

  elemental function mul_qo(p, q)  result(r)
    type(quater), intent(in)       :: p
    type(octonion), intent(in)     :: q
    type(octonion)                 :: r
  end

  elemental function mul_oq(q, p)  result(r)
    type(quater), intent(in)       :: p
    type(octonion), intent(in)     :: q
    type(octonion)                 :: r
    r = octonion([q%c(1) * p, q%c(2) * conjg(p)]) 
  end

  elemental function mul_co(c, q)  result(r)
    complex(wp), intent(in)        :: c
    type(octonion), intent(in)     :: q
    type(octonion)                 :: r
    r = octonion([c * q%c(1), q%c(2) * c])    ! NB: C and Q do not commute !!!
  end

  elemental function mul_oc(q, c)  result(r)
    complex(wp), intent(in)        :: c
    type(octonion), intent(in)     :: q
    type(octonion)                 :: r
    r = octonion([q%c(1) * c, q%c(2) * conjg(c)]) 
  end

  elemental function mul_ro(x, q)  result(r)
    real(wp), intent(in)           :: x
    type(octonion), intent(in)     :: q
    type(octonion)                 :: r
    r = octonion( x * q%c )
  end

  elemental function mul_or(q, x)  result(r)
    real(wp), intent(in)           :: x
    type(octonion), intent(in)     :: q
    type(octonion)                 :: r
    r = octonion( q%c * x)
  end

!------    Now sedenions:     ----------------------------------------!

  elemental subroutine s_equ_o(s, o)
    type(sedenion), intent(in out) :: s
    type(octonion), intent(in)     :: o
!    s = sedenion([ o, octonion(0d0) ])       !  Why not???
    s%c(1) = o;  s%c(2) = 0d0 
  end

  elemental subroutine s_equ_q(s, q)
    type(sedenion), intent(in out) :: s
    type(quater), intent(in)       :: q
    s%c(1) = q;  s%c(2) = 0d0 
  end

  elemental subroutine s_equ_c(s, c)
    type(sedenion), intent(in out) :: s
    complex(wp), intent(in)        :: c
    s%c(1) = c;  s%c(2) = 0d0 
  end

  elemental subroutine s_equ_r(s, r)
    type(sedenion), intent(in out) :: s
    real(wp), intent(in)           :: r
!    s = sedenion([ r, 0d0 ])                 !  Why not???
    s%c(1) = r;  s%c(2) = 0d0 
  end

  elemental function conjg_s(s)   result(r)
    type(sedenion), intent(in)    :: s
    type(sedenion)                :: r

    r = sedenion([ conjg(s%c(1)),  -s%c(2) ])
  end

  elemental function mul_ss(p, q)  result(r)    !  a c - d' b,  d a + b c`
    type(sedenion), intent(in)     :: p, q
    type(sedenion)                 :: r

    r%c(1) = p%c(1) * q%c(1)  -  conjg(q%c(2)) * p%c(2)
    r%c(2) = q%c(2) * p%c(1)  +  p%c(2) * conjg(q%c(1))
  end

  elemental function div_ss(p, q)  result(r)    !  a c - d' b,  d a + b c`
    type(sedenion), intent(in)     :: p, q
    type(sedenion)                 :: r

    r = p * sedenion(   [conjg(q%c(1)), -q%c(2)] *   &
           ( 1/(abs(q%c(1))**2+abs(q%c(2))**2) )   )
  end


end module


program test_CD
  use Cayley_Dickson
  implicit none 

  complex(wp), parameter  :: Ima = (0d0,1d0)
 
  type(quater)   :: q, i, j, k
  type(octonion) :: a, b
  type(sedenion) :: s, t, u

  q = quater([Ima, Ima])

  print *, "Hello", q
  print *, "Hello", q*q, abs(q)
  print *, conjg(q)*q

  q = Ima
  print *, "Ima=", q
  print *, "0=", quater(0)

  i = Ima; j = quater([0,1]); k = i*j

  print *, "i=", i
  print *, "j=", j
  print *, "k=", k


  a = k; b = j                ! (a,b) = (e3 + e10)
  s = sedenion([a,b])

  print "(16(1x,i5))", [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]
  print "(16(1x,f5.2))", a, b
  print "(16(1x,f5.2))", s

  a = 0d0; a%c(2) = j
  b = 0d0; b%c(2) = -k        ! (a,b) = (e6 - e15)
  t = sedenion([a,b])
  print "(16(1x,f5.2))", t

  print *, "s*t, t*s:"
  print "(16(1x,f5.2))", s*t     ! (e3 + e10)×(e6 - e15)
  print "(16(1x,f5.2))", t*s     ! (e6 - e15)×(e3 + e10)

  u = 1d0
  print *, "1, 1/s, s*(1/s), (1/s)*s:"
  print "(16(1x,f5.2))", u
  print "(16(1x,f5.2))", u/s
  print "(16(1x,f5.2))", s*(u/s)
  print "(16(1x,f5.2))", (u/s)*s

  print *, "s*t*1/t:"
  print "(16(1x,f5.2))", s*t*u/t
  print *, "s*(t*1/t):"
  print "(16(1x,f5.2))", s*(t*u/t)
  print *, "s*(t*(1/t)):"
  print "(16(1x,f5.2))", s*(t*u/t)

  s = sedenion((i+j-k)*s%c-t%c); s = s/sedenion(u%c+s%c+t%c)
  t = t/s*sedenion([t%c(1)*(j+2d0*k+3d0*i), t%c(2)])
  print *, "s, t:"
  print "(16(1x,f5.2))", s
  print "(16(1x,f5.2))", t

  print *, "s*(1/s), (1/s)*s, t*(1/t), (1/t)*t:"
  print "(16(1x,f5.2))", s*(u/s), (u/s)*s, t*(u/t), (u/t)*t

end
