From a27ac9855482ee56688bb2e5d69f4abd1080d29d Mon Sep 17 00:00:00 2001 From: Alberto Ramos Date: Sun, 3 Oct 2021 09:35:25 +0200 Subject: [PATCH] SU2 routines re-organized in different files --- src/Groups/AlgebraSU2.jl | 113 +++++++++++++++++++++++++++++++++++++ src/Groups/GroupSU2.jl | 119 --------------------------------------- src/Groups/Groups.jl | 13 ++++- src/Groups/SU2Types.jl | 26 +++++++++ 4 files changed, 149 insertions(+), 122 deletions(-) create mode 100644 src/Groups/AlgebraSU2.jl create mode 100644 src/Groups/SU2Types.jl diff --git a/src/Groups/AlgebraSU2.jl b/src/Groups/AlgebraSU2.jl new file mode 100644 index 0000000..dea936f --- /dev/null +++ b/src/Groups/AlgebraSU2.jl @@ -0,0 +1,113 @@ +### +### "THE BEER-WARE LICENSE": +### Alberto Ramos wrote this file. As long as you retain this +### notice you can do whatever you want with this stuff. If we meet some +### day, and you think this stuff is worth it, you can buy me a beer in +### return. +### +### file: AlgebraSU2.jl +### created: Sun Oct 3 09:24:25 2021 +### + +SU2alg(x::T) where T <: AbstractFloat = SU2alg{T}(x,0.0,0.0) +SU2alg(v::Vector{T}) where T <: AbstractFloat = SU2alg{T}(v[1],v[2],v[3]) +projalg(g::SU2{T}) where T <: AbstractFloat = SU2alg{T}(imag(g.t1), real(g.t2), imag(g.t2)) +dot(a::SU2alg{T}, b::SU2alg{T}) where T <: AbstractFloat = a.t1*b.t1 + a.t2*b.t2 + a.t3*b.t3 +norm(a::SU2alg{T}) where T <: AbstractFloat = sqrt(a.t1^2 + a.t2^2 + a.t3^2) +norm2(a::SU2alg{T}) where T <: AbstractFloat = a.t1^2 + a.t2^2 + a.t3^2 + +Base.:+(a::SU2alg{T}) where T <: AbstractFloat = SU2alg{T}(a.t1,a.t2,a.t3) +Base.:-(a::SU2alg{T}) where T <: AbstractFloat = SU2alg{T}(-a.t1,-a.t2,-a.t3) +Base.:+(a::SU2alg{T},b::SU2alg{T}) where T <: AbstractFloat = SU2alg{T}(a.t1+b.t1,a.t2+b.t2,a.t3+b.t3) +Base.:-(a::SU2alg{T},b::SU2alg{T}) where T <: AbstractFloat = SU2alg{T}(a.t1-b.t1,a.t2-b.t2,a.t3-b.t3) + +Base.:*(a::SU2alg{T},b::Number) where T <: AbstractFloat = SU2alg{T}(a.t1*b,a.t2*b,a.t3*b) +Base.:*(b::Number,a::SU2alg{T}) where T <: AbstractFloat = SU2alg{T}(a.t1*b,a.t2*b,a.t3*b) +Base.:/(a::SU2alg{T},b::Number) where T <: AbstractFloat = SU2alg{T}(a.t1/b,a.t2/b,a.t3/b) + +""" + function Base.exp(a::T, t::Number=1) where {T <: Algebra} + +Computes `exp(a)` +""" +function Base.exp(a::SU2alg{T}) where T <: AbstractFloat + + rm = sqrt( a.t1^2+a.t2^2+a.t3^2 )/2.0 + if (abs(rm) < 0.05) + rms = rm^2/2.0 + ca = 1.0 - rms *(1.0 - (rms/6.0 )*(1.0 - rms/15.0)) + sa = 0.5 - rms/6.0*(1.0 - (rms/10.0)*(1.0 - rms/21.0)) + else + ca = CUDA.cos(rm) + sa = CUDA.sin(rm)/(2.0*rm) + end + + t1 = complex(ca,sa*a.t1) + t2 = complex(sa*a.t2,sa*a.t3) + return SU2{T}(t1,t2) +end + +function Base.exp(a::SU2alg{T}, t::T) where T <: AbstractFloat + + rm = t*sqrt( a.t1^2+a.t2^2+a.t3^2 )/2.0 + if (abs(rm) < 0.05) + rms = rm^2/2.0 + ca = 1.0 - rms *(1.0 - (rms/6.0 )*(1.0 - rms/15.0)) + sa = t*(0.5 - rms/6.0*(1.0 - (rms/10.0)*(1.0 - rms/21.0))) + else + ca = CUDA.cos(rm) + sa = t*CUDA.sin(rm)/(2.0*rm) + end + + t1 = complex(ca,sa*a.t1) + t2 = complex(sa*a.t2,sa*a.t3) + return SU2{T}(t1,t2) +end + + +""" + function expm(g::G, a::A) where {G <: Algebra, A <: Algebra} + +Computes `exp(a)*g` + +""" +function expm(g::SU2{T}, a::SU2alg{T}) where T <: AbstractFloat + + rm = sqrt( a.t1^2+a.t2^2+a.t3^2 )/2.0 + if (abs(rm) < 0.05) + rms = rm^2/2.0 + ca = 1.0 - rms *(1.0 - (rms/6.0 )*(1.0 - rms/15.0)) + sa = 0.5 - rms/6.0*(1.0 - (rms/10.0)*(1.0 - rms/21.0)) + else + ca = CUDA.cos(rm) + sa = CUDA.sin(rm)/(2.0*rm) + end + + t1 = complex(ca,sa*a.t1)*g.t1-complex(sa*a.t2,sa*a.t3)*conj(g.t2) + t2 = complex(ca,sa*a.t1)*g.t2+complex(sa*a.t2,sa*a.t3)*conj(g.t1) + return SU2{T}(t1,t2) +end + +""" + function expm(g::SU2, a::SU2alg, t::Float64) + +Computes `exp(t*a)*g` + +""" +function expm(g::SU2{T}, a::SU2alg{T}, t::T) where T <: AbstractFloat + + rm = t*sqrt( a.t1^2+a.t2^2+a.t3^2 )/2.0 + if (abs(rm) < 0.05) + rms = rm^2/2.0 + ca = 1.0 - rms *(1.0 - (rms/6.0 )*(1.0 - rms/15.0)) + sa = t*(0.5 - rms/6.0*(1.0 - (rms/10.0)*(1.0 - rms/21.0))) + else + ca = CUDA.cos(rm) + sa = t*CUDA.sin(rm)/(2.0*rm) + end + + t1 = complex(ca,sa*a.t1)*g.t1-complex(sa*a.t2,sa*a.t3)*conj(g.t2) + t2 = complex(ca,sa*a.t1)*g.t2+complex(sa*a.t2,sa*a.t3)*conj(g.t1) + return SU2{T}(t1,t2) + +end diff --git a/src/Groups/GroupSU2.jl b/src/Groups/GroupSU2.jl index 737c872..7174153 100644 --- a/src/Groups/GroupSU2.jl +++ b/src/Groups/GroupSU2.jl @@ -15,20 +15,12 @@ # https://en.wikipedia.org/wiki/Cayley%E2%80%93Dickson_construction using CUDA, Random -import Base.:*, Base.:+, Base.:-,Base.:/,Base.:\,Base.exp,Base.zero,Base.one -import Random.rand -struct SU2{T} <: Group - t1::Complex{T} - t2::Complex{T} -end SU2(a::T, b::T) where T <: AbstractFloat = SU2{T}(complex(a), complex(b)) inverse(b::SU2{T}) where T <: AbstractFloat = SU2{T}(conj(b.t1), -b.t2) dag(a::SU2{T}) where T <: AbstractFloat = inverse(a) norm(a::SU2{T}) where T <: AbstractFloat = sqrt(abs2(a.t1) + abs2(a.t2)) norm2(a::SU2{T}) where T <: AbstractFloat = abs2(a.t1) + abs2(a.t2) tr(g::SU2{T}) where T <: AbstractFloat = complex(2.0*real(g.t1), 0.0) -Base.one(::Type{SU2{T}}) where T <: AbstractFloat = SU2{T}(one(T),zero(T)) -Random.rand(rng::AbstractRNG, ::Random.SamplerType{SU2{T}}) where T <: AbstractFloat = exp(SU2alg{T}(randn(rng,T),randn(rng,T),randn(rng,T))) """ function normalize(a::T) where {T <: Group} @@ -47,29 +39,6 @@ Base.:*(a::SU2{T},b::SU2{T}) where T <: AbstractFloat = SU2{T}(a.t1*b.t1-a.t2*co Base.:/(a::SU2{T},b::SU2{T}) where T <: AbstractFloat = SU2{T}(a.t1*conj(b.t1)+a.t2*conj(b.t2),-a.t1*b.t2+a.t2*b.t1) Base.:\(a::SU2{T},b::SU2{T}) where T <: AbstractFloat = SU2{T}(conj(a.t1)*b.t1+a.t2*conj(b.t2),conj(a.t1)*b.t2-a.t2*conj(b.t1)) -struct SU2alg{T} <: Algebra - t1::T - t2::T - t3::T -end -SU2alg(x::T) where T <: AbstractFloat = SU2alg{T}(x,0.0,0.0) -SU2alg(v::Vector{T}) where T <: AbstractFloat = SU2alg{T}(v[1],v[2],v[3]) -projalg(g::SU2{T}) where T <: AbstractFloat = SU2alg{T}(imag(g.t1), real(g.t2), imag(g.t2)) -dot(a::SU2alg{T}, b::SU2alg{T}) where T <: AbstractFloat = a.t1*b.t1 + a.t2*b.t2 + a.t3*b.t3 -norm(a::SU2alg{T}) where T <: AbstractFloat = sqrt(a.t1^2 + a.t2^2 + a.t3^2) -norm2(a::SU2alg{T}) where T <: AbstractFloat = a.t1^2 + a.t2^2 + a.t3^2 -Base.zero(::Type{SU2alg{T}}) where T <: AbstractFloat = SU2alg{T}(zero(T),zero(T),zero(T)) -Random.rand(rng::AbstractRNG, ::Random.SamplerType{SU2alg{T}}) where T <: AbstractFloat = SU2alg{T}(randn(rng,T),randn(rng,T),randn(rng,T)) - -Base.:+(a::SU2alg{T}) where T <: AbstractFloat = SU2alg{T}(a.t1,a.t2,a.t3) -Base.:-(a::SU2alg{T}) where T <: AbstractFloat = SU2alg{T}(-a.t1,-a.t2,-a.t3) -Base.:+(a::SU2alg{T},b::SU2alg{T}) where T <: AbstractFloat = SU2alg{T}(a.t1+b.t1,a.t2+b.t2,a.t3+b.t3) -Base.:-(a::SU2alg{T},b::SU2alg{T}) where T <: AbstractFloat = SU2alg{T}(a.t1-b.t1,a.t2-b.t2,a.t3-b.t3) - -Base.:*(a::SU2alg{T},b::Number) where T <: AbstractFloat = SU2alg{T}(a.t1*b,a.t2*b,a.t3*b) -Base.:*(b::Number,a::SU2alg{T}) where T <: AbstractFloat = SU2alg{T}(a.t1*b,a.t2*b,a.t3*b) -Base.:/(a::SU2alg{T},b::Number) where T <: AbstractFloat = SU2alg{T}(a.t1/b,a.t2/b,a.t3/b) - function isgroup(a::SU2{T}) where T <: AbstractFloat tol = 1.0E-10 if (abs2(a.t1) + abs2(a.t2) - 1.0 < 1.0E-10) @@ -79,91 +48,3 @@ function isgroup(a::SU2{T}) where T <: AbstractFloat end end -""" - function Base.exp(a::T, t::Number=1) where {T <: Algebra} - -Computes `exp(a)` -""" -function Base.exp(a::SU2alg{T}) where T <: AbstractFloat - - rm = sqrt( a.t1^2+a.t2^2+a.t3^2 )/2.0 - if (abs(rm) < 0.05) - rms = rm^2/2.0 - ca = 1.0 - rms *(1.0 - (rms/6.0 )*(1.0 - rms/15.0)) - sa = 0.5 - rms/6.0*(1.0 - (rms/10.0)*(1.0 - rms/21.0)) - else - ca = CUDA.cos(rm) - sa = CUDA.sin(rm)/(2.0*rm) - end - - t1 = complex(ca,sa*a.t1) - t2 = complex(sa*a.t2,sa*a.t3) - return SU2{T}(t1,t2) -end - -function Base.exp(a::SU2alg{T}, t::T) where T <: AbstractFloat - - rm = t*sqrt( a.t1^2+a.t2^2+a.t3^2 )/2.0 - if (abs(rm) < 0.05) - rms = rm^2/2.0 - ca = 1.0 - rms *(1.0 - (rms/6.0 )*(1.0 - rms/15.0)) - sa = t*(0.5 - rms/6.0*(1.0 - (rms/10.0)*(1.0 - rms/21.0))) - else - ca = CUDA.cos(rm) - sa = t*CUDA.sin(rm)/(2.0*rm) - end - - t1 = complex(ca,sa*a.t1) - t2 = complex(sa*a.t2,sa*a.t3) - return SU2{T}(t1,t2) -end - - -""" - function expm(g::G, a::A) where {G <: Algebra, A <: Algebra} - -Computes `exp(a)*g` - -""" -function expm(g::SU2{T}, a::SU2alg{T}) where T <: AbstractFloat - - rm = sqrt( a.t1^2+a.t2^2+a.t3^2 )/2.0 - if (abs(rm) < 0.05) - rms = rm^2/2.0 - ca = 1.0 - rms *(1.0 - (rms/6.0 )*(1.0 - rms/15.0)) - sa = 0.5 - rms/6.0*(1.0 - (rms/10.0)*(1.0 - rms/21.0)) - else - ca = CUDA.cos(rm) - sa = CUDA.sin(rm)/(2.0*rm) - end - - t1 = complex(ca,sa*a.t1)*g.t1-complex(sa*a.t2,sa*a.t3)*conj(g.t2) - t2 = complex(ca,sa*a.t1)*g.t2+complex(sa*a.t2,sa*a.t3)*conj(g.t1) - return SU2{T}(t1,t2) -end - -""" - function expm(g::SU2, a::SU2alg, t::Float64) - -Computes `exp(t*a)*g` - -""" -function expm(g::SU2{T}, a::SU2alg{T}, t::T) where T <: AbstractFloat - - rm = t*sqrt( a.t1^2+a.t2^2+a.t3^2 )/2.0 - if (abs(rm) < 0.05) - rms = rm^2/2.0 - ca = 1.0 - rms *(1.0 - (rms/6.0 )*(1.0 - rms/15.0)) - sa = t*(0.5 - rms/6.0*(1.0 - (rms/10.0)*(1.0 - rms/21.0))) - else - ca = CUDA.cos(rm) - sa = t*CUDA.sin(rm)/(2.0*rm) - end - - t1 = complex(ca,sa*a.t1)*g.t1-complex(sa*a.t2,sa*a.t3)*conj(g.t2) - t2 = complex(ca,sa*a.t1)*g.t2+complex(sa*a.t2,sa*a.t3)*conj(g.t1) - return SU2{T}(t1,t2) - -end - -export SU2, SU2alg, inverse, dag, tr, projalg, expm, exp, norm, norm2, isgroup diff --git a/src/Groups/Groups.jl b/src/Groups/Groups.jl index e6dfa1b..138f535 100644 --- a/src/Groups/Groups.jl +++ b/src/Groups/Groups.jl @@ -12,8 +12,8 @@ module Groups -using Random -import Base.:*, Base.:+, Base.:-,Base.:/,Base.:\,Base.one,Base.zero +using CUDA, Random +import Base.:*, Base.:+, Base.:-,Base.:/,Base.:\,Base.exp,Base.one,Base.zero import Random.rand abstract type Group end @@ -21,9 +21,16 @@ abstract type Algebra end export Group, Algebra -include("GroupSU2.jl") +## +# SU(2) and 2x2 matrix operations +## +include("SU2Types.jl") export SU2, SU2alg +include("GroupSU2.jl") +include("AlgebraSU2.jl") +## END SU(2) + ## # SU(3) and 3x3 matrix operations ## diff --git a/src/Groups/SU2Types.jl b/src/Groups/SU2Types.jl new file mode 100644 index 0000000..e6e5340 --- /dev/null +++ b/src/Groups/SU2Types.jl @@ -0,0 +1,26 @@ +### +### "THE BEER-WARE LICENSE": +### Alberto Ramos wrote this file. As long as you retain this +### notice you can do whatever you want with this stuff. If we meet some +### day, and you think this stuff is worth it, you can buy me a beer in +### return. +### +### file: SU2Types.jl +### created: Sun Oct 3 09:22:48 2021 +### + +struct SU2{T} <: Group + t1::Complex{T} + t2::Complex{T} +end + +struct SU2alg{T} <: Algebra + t1::T + t2::T + t3::T +end + +Base.zero(::Type{SU2alg{T}}) where T <: AbstractFloat = SU2alg{T}(zero(T),zero(T),zero(T)) +Random.rand(rng::AbstractRNG, ::Random.SamplerType{SU2alg{T}}) where T <: AbstractFloat = SU2alg{T}(randn(rng,T),randn(rng,T),randn(rng,T)) +Base.one(::Type{SU2{T}}) where T <: AbstractFloat = SU2{T}(one(T),zero(T)) +Random.rand(rng::AbstractRNG, ::Random.SamplerType{SU2{T}}) where T <: AbstractFloat = exp(SU2alg{T}(randn(rng,T),randn(rng,T),randn(rng,T)))