Add sectorscalartype (#146)

* Add `Feltype` and `Reltype`

* Simplify some manipulations in terms of `Feltype` and `Reltype`

* Implement review comments

* Change Fscalartype and Rscalartype to sectorscalartype

TensorKitSectors/src/TensorKitSectors.jl

@@ -10,6 +10,7 @@ export Sector, Group, AbstractIrrep
 export Irrep
 
 export Nsymbol, Fsymbol, Rsymbol, Asymbol, Bsymbol
+export sectorscalartype
 export dim, sqrtdim, invsqrtdim, frobeniusschur, twist, fusiontensor, dual
 export otimes, deligneproduct, times
 export FusionStyle, UniqueFusion, MultipleFusion, SimpleFusion, GenericFusion,

TensorKitSectors/src/irreps.jl

@@ -197,22 +197,26 @@ findindex(::SectorValues{SU2Irrep}, s::SU2Irrep) = twice(s.j) + 1
 dim(s::SU2Irrep) = twice(s.j) + 1
 
 FusionStyle(::Type{SU2Irrep}) = SimpleFusion()
+sectorscalartype(::Type{SU2Irrep}) = Float64
 Base.isreal(::Type{SU2Irrep}) = true
 
 Nsymbol(sa::SU2Irrep, sb::SU2Irrep, sc::SU2Irrep) = WignerSymbols.δ(sa.j, sb.j, sc.j)
+
 function Fsymbol(s1::SU2Irrep, s2::SU2Irrep, s3::SU2Irrep,
                  s4::SU2Irrep, s5::SU2Irrep, s6::SU2Irrep)
     if all(==(_su2one), (s1, s2, s3, s4, s5, s6))
         return 1.0
     else
         return sqrtdim(s5) * sqrtdim(s6) *
-               WignerSymbols.racahW(Float64, s1.j, s2.j,
+               WignerSymbols.racahW(sectorscalartype(SU2Irrep), s1.j, s2.j,
                                     s4.j, s3.j, s5.j, s6.j)
     end
 end
+
 function Rsymbol(sa::SU2Irrep, sb::SU2Irrep, sc::SU2Irrep)
-    Nsymbol(sa, sb, sc) || return 0.0
-    return iseven(convert(Int, sa.j + sb.j - sc.j)) ? 1.0 : -1.0
+    Nsymbol(sa, sb, sc) || return zero(sectorscalartype(SU2Irrep))
+    return iseven(convert(Int, sa.j + sb.j - sc.j)) ? one(sectorscalartype(SU2Irrep)) :
+           -one(sectorscalartype(SU2Irrep))
 end
 
 function fusiontensor(a::SU2Irrep, b::SU2Irrep, c::SU2Irrep)
@@ -341,59 +345,62 @@ Base.eltype(::Type{CU1ProdIterator}) = CU1Irrep
 dim(c::CU1Irrep) = ifelse(c.j == zero(HalfInt), 1, 2)
 
 FusionStyle(::Type{CU1Irrep}) = SimpleFusion()
+sectorscalartype(::Type{CU1Irrep}) = Float64
 Base.isreal(::Type{CU1Irrep}) = true
 
 function Nsymbol(a::CU1Irrep, b::CU1Irrep, c::CU1Irrep)
     return ifelse(c.s == 0, (a.j == b.j) & ((a.s == b.s == 2) | (a.s == b.s)),
                   ifelse(c.s == 1, (a.j == b.j) & ((a.s == b.s == 2) | (a.s != b.s)),
                          (c.j == a.j + b.j) | (c.j == abs(a.j - b.j))))
 end
+
 function Fsymbol(a::CU1Irrep, b::CU1Irrep, c::CU1Irrep,
                  d::CU1Irrep, e::CU1Irrep, f::CU1Irrep)
     Nabe = convert(Int, Nsymbol(a, b, e))
     Necd = convert(Int, Nsymbol(e, c, d))
     Nbcf = convert(Int, Nsymbol(b, c, f))
     Nafd = convert(Int, Nsymbol(a, f, d))
 
-    Nabe * Necd * Nbcf * Nafd == 0 && return 0.0
+    T = sectorscalartype(CU1Irrep)
+    Nabe * Necd * Nbcf * Nafd == 0 && return zero(T)
 
     op = CU1Irrep(0, 0)
     om = CU1Irrep(0, 1)
 
     if a == op || b == op || c == op
-        return 1.0
+        return one(T)
     end
     if (a == b == om) || (a == c == om) || (b == c == om)
-        return 1.0
+        return one(T)
     end
     if a == om
         if d.j == zero(HalfInt)
-            return 1.0
+            return one(T)
         else
-            return (d.j == c.j - b.j) ? -1.0 : 1.0
+            return (d.j == c.j - b.j) ? -one(T) : one(T)
         end
     end
     if b == om
-        return (d.j == abs(a.j - c.j)) ? -1.0 : 1.0
+        return (d.j == abs(a.j - c.j)) ? -one(T) : one(T)
     end
     if c == om
-        return (d.j == a.j - b.j) ? -1.0 : 1.0
+        return (d.j == a.j - b.j) ? -one(T) : one(T)
     end
     # from here on, a, b, c are neither 0+ or 0-
-    s = sqrt(2) / 2
+    s = T(sqrt(2) / 2)
     if a == b == c
         if d == a
             if e.j == 0
                 if f.j == 0
-                    return f.s == 1 ? -0.5 : 0.5
+                    return f.s == 1 ? T(-1 // 2) : T(1 // 2)
                 else
                     return e.s == 1 ? -s : s
                 end
             else
-                return f.j == 0 ? s : 0.0
+                return f.j == 0 ? s : zero(T)
             end
         else
-            return 1.0
+            return one(T)
         end
     end
     if a == b # != c
@@ -404,7 +411,7 @@ function Fsymbol(a::CU1Irrep, b::CU1Irrep, c::CU1Irrep,
                 return s
             end
         else
-            return 1.0
+            return one(T)
         end
     end
     if b == c
@@ -415,27 +422,28 @@ function Fsymbol(a::CU1Irrep, b::CU1Irrep, c::CU1Irrep,
                 return f.s == 1 ? -s : s
             end
         else
-            return 1.0
+            return one(T)
         end
     end
     if a == c
         if d == b
             if e.j == f.j
-                return 0.0
+                return zero(T)
             else
-                return 1.0
+                return one(T)
             end
         else
-            return d.s == 1 ? -1.0 : 1.0
+            return d.s == 1 ? -one(T) : one(T)
         end
     end
     if d == om
-        return b.j == a.j + c.j ? -1.0 : 1.0
+        return b.j == a.j + c.j ? -one(T) : one(T)
     end
-    return 1.0
+    return one(T)
 end
+
 function Rsymbol(a::CU1Irrep, b::CU1Irrep, c::CU1Irrep)
-    R = convert(Float64, Nsymbol(a, b, c))
+    R = convert(sectorscalartype(CU1Irrep), Nsymbol(a, b, c))
     return c.s == 1 && a.j > 0 ? -R : R
 end
 

TensorKitSectors/src/sectors.jl

@@ -74,21 +74,28 @@ Return the conjugate label `conj(a)`.
 dual(a::Sector) = conj(a)
 
 """
-    isreal(::Type{<:Sector}) -> Bool
+    sectorscalartype(I::Type{<:Sector}) -> Type
 
-Return whether the topological data (Fsymbol, Rsymbol) of the sector is real or not (in
-which case it is complex).
+Return the scalar type of the topological data (Fsymbol, Rsymbol) of the sector `I`.
 """
-function Base.isreal(I::Type{<:Sector})
-    u = one(I)
-    if BraidingStyle(I) isa HasBraiding
-        return (eltype(Fsymbol(u, u, u, u, u, u)) <: Real) &&
-               (eltype(Rsymbol(u, u, u)) <: Real)
+function sectorscalartype(::Type{I}) where {I<:Sector}
+    if BraidingStyle(I) isa NoBraiding
+        return eltype(Core.Compiler.return_type(Fsymbol, NTuple{6,I}))
     else
-        return (eltype(Fsymbol(u, u, u, u, u, u)) <: Real)
+        Feltype = eltype(Core.Compiler.return_type(Fsymbol, NTuple{6,I}))
+        Reltype = eltype(Core.Compiler.return_type(Rsymbol, NTuple{3,I}))
+        return Base.promote_op(*, Feltype, Reltype)
     end
 end
 
+"""
+    isreal(::Type{<:Sector}) -> Bool
+
+Return whether the topological data (Fsymbol, Rsymbol) of the sector is real or not (in
+which case it is complex).
+"""
+Base.isreal(I::Type{<:Sector}) = sectorscalartype(I) <: Real
+
 # FusionStyle: the most important aspect of Sector
 #---------------------------------------------
 """

src/fusiontrees/manipulations.jl

@@ -17,7 +17,7 @@ function insertat(f₁::FusionTree{I}, i::Int, f₂::FusionTree{I,0}) where {I}
     # this actually removes uncoupled line i, which should be trivial
     (f₁.uncoupled[i] == f₂.coupled && !f₁.isdual[i]) ||
         throw(SectorMismatch("cannot connect $(f₂.uncoupled) to $(f₁.uncoupled[i])"))
-    coeff = Fsymbol(one(I), one(I), one(I), one(I), one(I), one(I))[1, 1, 1, 1]
+    coeff = one(sectorscalartype(I))
 
     uncoupled = TupleTools.deleteat(f₁.uncoupled, i)
     coupled = f₁.coupled
@@ -39,7 +39,7 @@ function insertat(f₁::FusionTree{I}, i, f₂::FusionTree{I,1}) where {I}
     # identity operation
     (f₁.uncoupled[i] == f₂.coupled && !f₁.isdual[i]) ||
         throw(SectorMismatch("cannot connect $(f₂.uncoupled) to $(f₁.uncoupled[i])"))
-    coeff = Fsymbol(one(I), one(I), one(I), one(I), one(I), one(I))[1, 1, 1, 1]
+    coeff = one(sectorscalartype(I))
     isdual′ = TupleTools.setindex(f₁.isdual, f₂.isdual[1], i)
     f = FusionTree{I}(f₁.uncoupled, f₁.coupled, isdual′, f₁.innerlines, f₁.vertices)
     return fusiontreedict(I)(f => coeff)
@@ -59,7 +59,7 @@ function insertat(f₁::FusionTree{I}, i, f₂::FusionTree{I,2}) where {I}
         isdual′ = (isdualb, isdualc, tail(isdual)...)
         inner′ = (uncoupled[1], inner...)
         vertices′ = (f₂.vertices..., f₁.vertices...)
-        coeff = Fsymbol(one(I), one(I), one(I), one(I), one(I), one(I))[1, 1, 1, 1]
+        coeff = one(sectorscalartype(I))
         f′ = FusionTree(uncoupled′, coupled, isdual′, inner′, vertices′)
         return fusiontreedict(I)(f′ => coeff)
     end
@@ -110,8 +110,8 @@ function insertat(f₁::FusionTree{I,N₁}, i, f₂::FusionTree{I,N₂}) where {
     F = fusiontreetype(I, N₁ + N₂ - 1)
     (f₁.uncoupled[i] == f₂.coupled && !f₁.isdual[i]) ||
         throw(SectorMismatch("cannot connect $(f₂.uncoupled) to $(f₁.uncoupled[i])"))
-    coeff = Fsymbol(one(I), one(I), one(I), one(I), one(I), one(I))[1, 1]
-    T = typeof(coeff)
+    T = sectorscalartype(I)
+    coeff = one(T)
     if length(f₁) == 1
         return fusiontreedict(I){F,T}(f₂ => coeff)
     end
@@ -457,8 +457,8 @@ function _recursive_repartition(f₁::FusionTree{I,N₁},
     # precompute the parameters of the return type
     F₁ = fusiontreetype(I, N)
     F₂ = fusiontreetype(I, N₁ + N₂ - N)
-    coeff = @inbounds Fsymbol(one(I), one(I), one(I), one(I), one(I), one(I))[1, 1, 1, 1]
-    T = typeof(coeff)
+    T = sectorscalartype(I)
+    coeff = one(T)
     if N == N₁
         return fusiontreedict(I){Tuple{F₁,F₂},T}((f₁, f₂) => coeff)
     else
@@ -500,8 +500,7 @@ function Base.transpose(f₁::FusionTree{I}, f₂::FusionTree{I},
     p = linearizepermutation(p1, p2, length(f₁), length(f₂))
     @assert iscyclicpermutation(p)
     if usetransposecache[]
-        u = one(I)
-        T = eltype(Fsymbol(u, u, u, u, u, u))
+        T = sectorscalartype(I)
         F₁ = fusiontreetype(I, N₁)
         F₂ = fusiontreetype(I, N₂)
         D = fusiontreedict(I){Tuple{F₁,F₂},T}
@@ -587,8 +586,7 @@ function planar_trace(f₁::FusionTree{I}, f₂::FusionTree{I},
          map(l -> l - count(l .> q′), TupleTools.getindices(linearindex, p2)))
     end
 
-    u = one(I)
-    T = typeof(Fsymbol(u, u, u, u, u, u)[1, 1, 1, 1])
+    T = sectorscalartype(I)
     F₁ = fusiontreetype(I, N₁)
     F₂ = fusiontreetype(I, N₂)
     newtrees = FusionTreeDict{Tuple{F₁,F₂},T}()
@@ -615,8 +613,7 @@ of output trees and corresponding coefficients.
 """
 function planar_trace(f::FusionTree{I,N},
                       q1::IndexTuple{N₃}, q2::IndexTuple{N₃}) where {I<:Sector,N,N₃}
-    u = one(I)
-    T = typeof(Fsymbol(u, u, u, u, u, u)[1, 1, 1, 1])
+    T = sectorscalartype(I)
     F = fusiontreetype(I, N - 2 * N₃)
     newtrees = FusionTreeDict{F,T}()
     N₃ === 0 && return push!(newtrees, f => one(T))
@@ -673,8 +670,7 @@ function elementary_trace(f::FusionTree{I,N}, i) where {I<:Sector,N}
     i < N || f.coupled == one(I) ||
         throw(ArgumentError("Cannot trace outputs i=$N and 1 of fusion tree that couples to non-trivial sector"))
 
-    u = one(I)
-    T = typeof(Fsymbol(u, u, u, u, u, u)[1, 1, 1, 1])
+    T = sectorscalartype(I)
     F = fusiontreetype(I, N - 2)
     newtrees = FusionTreeDict{F,T}()
 
@@ -786,12 +782,7 @@ function artin_braid(f::FusionTree{I,N}, i; inv::Bool=false) where {I<:Sector,N}
     inner_extended = (uncoupled[1], inner..., coupled′)
     vertices = f.vertices
     u = one(I)
-
-    if BraidingStyle(I) isa NoBraiding
-        oneT = Fsymbol(u, u, u, u, u, u)[1, 1, 1, 1]
-    else
-        oneT = Rsymbol(u, u, u)[1, 1] * Fsymbol(u, u, u, u, u, u)[1, 1, 1, 1]
-    end
+    oneT = one(sectorscalartype(I))
 
     if u in (uncoupled[i], uncoupled[i + 1])
         # braiding with trivial sector: simple and always possible
@@ -925,7 +916,7 @@ function braid(f::FusionTree{I,N},
                p::NTuple{N,Int}) where {I<:Sector,N}
     TupleTools.isperm(p) || throw(ArgumentError("not a valid permutation: $p"))
     if FusionStyle(I) isa UniqueFusion && BraidingStyle(I) isa SymmetricBraiding
-        coeff = Rsymbol(one(I), one(I), one(I))
+        coeff = one(sectorscalartype(I))
         for i in 1:N
             for j in 1:(i - 1)
                 if p[j] > p[i]
@@ -940,10 +931,7 @@ function braid(f::FusionTree{I,N},
         f′ = FusionTree{I}(uncoupled′, coupled′, isdual′)
         return fusiontreedict(I)(f′ => coeff)
     else
-        u = one(I)
-        T = BraidingStyle(I) isa NoBraiding ?
-            typeof(Fsymbol(u, u, u, u, u, u)[1, 1, 1, 1]) :
-            typeof(Rsymbol(u, u, u)[1, 1] * Fsymbol(u, u, u, u, u, u)[1, 1, 1, 1])
+        T = sectorscalartype(I)
         coeff = one(T)
         trees = FusionTreeDict(f => coeff)
         newtrees = empty(trees)
@@ -1009,8 +997,7 @@ function braid(f₁::FusionTree{I}, f₂::FusionTree{I},
     if FusionStyle(f₁) isa UniqueFusion &&
        BraidingStyle(f₁) isa SymmetricBraiding
         if usebraidcache_abelian[]
-            u = one(I)
-            T = Int
+            T = Int # do we hardcode this ?
             F₁ = fusiontreetype(I, N₁)
             F₂ = fusiontreetype(I, N₂)
             D = SingletonDict{Tuple{F₁,F₂},T}
@@ -1020,11 +1007,7 @@ function braid(f₁::FusionTree{I}, f₂::FusionTree{I},
         end
     else
         if usebraidcache_nonabelian[]
-            u = one(I)
-            T = BraidingStyle(I) isa NoBraiding ?
-                typeof(Fsymbol(u, u, u, u, u, u)[1, 1, 1, 1]) :
-                typeof(sqrtdim(u) * Fsymbol(u, u, u, u, u, u)[1, 1, 1, 1] *
-                       Rsymbol(u, u, u)[1, 1])
+            T = sectorscalartype(I)
             F₁ = fusiontreetype(I, N₁)
             F₂ = fusiontreetype(I, N₂)
             D = FusionTreeDict{Tuple{F₁,F₂},T}