diff --git a/src/WignerSymbols.jl b/src/WignerSymbols.jl
index 22fff82..5b6839e 100644
--- a/src/WignerSymbols.jl
+++ b/src/WignerSymbols.jl
@@ -10,10 +10,10 @@ clebschgordan(j₁, m₁, j₂, m₂, j₃, m₃ = m₁+m₂) = clebschgordan(Fl
 function clebschgordan(T::Type{<:AbstractFloat}, j₁, m₁, j₂, m₂, j₃, m₃ = m₁+m₂)
     s = wigner3j(T, j₁, j₂, j₃, m₁, m₂, -m₃)
     iszero(s) && return s
-    if isodd(j₁ - j₂ + m₃)
+    if isodd(convert(Int,j₁ - j₂ + m₃))
         s *= -one(s)
     end
-    s /= sqrt(convert(T, j₃+j₃+1))
+    s *= sqrt(convert(T, j₃+j₃+one(j₃)))
     return s
 end
 
@@ -34,8 +34,8 @@ function wigner3j(T::Type{<:AbstractFloat}, j₁, j₂, j₃, m₁, m₂, m₃ =
     β₂ = convert(Int, j₁ - m₁ )
     β₃ = convert(Int, j₂ + m₂ )
 
-    # extra sign in definition
-    sgn = isodd(j₁ - j₃ - m₃) ? -sgn : sgn
+    # extra sign in definition: α₁ - α₂ = j₁ + m₂ - j₂ + m₁ = j₁ - j₂ + m₃
+    sgn = isodd(α₁ - α₂) ? -sgn : sgn
 
     # dictionary lookup or compute
     if haskey(Wigner3j, (β₁, β₂, β₃, α₁, α₂))
@@ -44,7 +44,8 @@ function wigner3j(T::Type{<:AbstractFloat}, j₁, j₂, j₃, m₁, m₂, m₃ =
         s1 = Δ²(j₁, j₂, j₃)
         s2 = prod(map(primefactorial, (β₂, β₁ - α₁, β₁ - α₂, β₃, β₃ - α₁, β₂ - α₂)))
 
-        snum, rnum = splitsquare(s1.num * s2)
+
+        snum, rnum = splitsquare(s1.num*s2)
         sden, rden = splitsquare(s1.den)
         s = convert(BigInt, snum) // convert(BigInt, sden)
         r = convert(BigInt, rnum) // convert(BigInt, rden)
@@ -118,6 +119,7 @@ function δ(j₁, j₂, j₃)
 end
 
 # triangle coefficient
+Δ(j₁, j₂, j₃) = Δ(Float64, j₁, j₂, j₃)
 function Δ(T::Type{<:AbstractFloat}, j₁, j₂, j₃)
     if !δ(j₁, j₂, j₃)
         throw(DomainError())
@@ -139,7 +141,7 @@ end
 
 # reorder parameters determining the 3j symbol to canonical order:
 # j₁ >= j₂ >= j₃ and m₁ >= 0 or m₁ == 0 && m₂ >= 0
-function reorder3j(j₁, j₂, j₃, m₁, m₂, m₃, sign = one(UInt8))
+function reorder3j(j₁, j₂, j₃, m₁, m₂, m₃, sign = one(Int8))
     if j₁ < j₂
         return reorder3j(j₂, j₁, j₃, m₂, m₁, m₃, -sign)
     elseif j₂ < j₃
@@ -150,7 +152,7 @@ function reorder3j(j₁, j₂, j₃, m₁, m₂, m₃, sign = one(UInt8))
         return reorder3j(j₁, j₂, j₃, -m₁, -m₂, -m₃, -sign)
     else
         # sign doesn't matter if total J=j₁ + j₂ + j₃ is even
-        if iseven(j₁ + j₂ + j₃)
+        if iseven(convert(UInt,j₁ + j₂ + j₃))
             sign = one(sign)
         end
         return (j₁, j₂, j₃, m₁, m₂, m₃, sign)
@@ -183,14 +185,14 @@ function compute3jseries(β₁, β₂, β₃, α₁, α₂)
     nums = Vector{T}(length(krange))
     dens = Vector{T}(length(krange))
     for (i, k) in enumerate(krange)
-        num = iseven(k) ? primefactorial(1) : -primefactorial(1)
+        num = iseven(k) ? one(T) : -one(T)
         den = primefactorial(k)*primefactorial(k-α₁)*primefactorial(k-α₂)*
             primefactorial(β₁-k)*primefactorial(β₂-k)*primefactorial(β₃-k)
         nums[i], dens[i] = divgcd!(num, den)
     end
     den = commondenominator!(nums, dens)
     totalnum = sumlist!(nums)
-    totalden = convert(BigInt, PrimeFactorization(den))
+    totalden = convert(BigInt, den)
     return totalnum//totalden
 end
 
@@ -209,7 +211,7 @@ function compute6jseries(β₁, β₂, β₃, α₁, α₂, α₃, α₄)
     end
     den = commondenominator!(nums, dens)
     totalnum = sumlist!(nums)
-    totalden = convert(BigInt, PrimeFactorization(den))
+    totalden = convert(BigInt, den)
     return totalnum//totalden
 end
 
diff --git a/src/primefactorization.jl b/src/primefactorization.jl
index 19a4693..25d81a8 100644
--- a/src/primefactorization.jl
+++ b/src/primefactorization.jl
@@ -56,12 +56,9 @@ struct PrimeFactorization{T<:Unsigned} <: Integer
 end
 PrimeFactorization(powers::Vector{T}) where {T<:Unsigned} = PrimeFactorization{T}(powers, one(Int8))
 
-Base.copy(a::PrimeFactorization) = PrimeFactorization(copy(a.powers), a.sign)
-Base.:(==)(a::P, b::P) where {P<:PrimeFactorization} = a.powers == b.powers && a.sign == b.sign
-
 # define our own factor function, returning an instance of PrimeFactorization
 function primefactor(n::Integer)
-    iszero(n) && return PrimeFactorization(UInt8[], zero(UInt8))
+    iszero(n) && return PrimeFactorization(UInt8[], zero(Int8))
     sn = n < 0 ? -one(Int8) : one(Int8)
     n = abs(n)
     m = length(factortable)
@@ -104,7 +101,16 @@ function primefactorial(n::Integer)
     @inbounds return PrimeFactorization(copy(factorialtable[n+1]))
 end
 
-# Conversion from PrimeFactorization to `BigInt` using `bigprime`
+# Methods for PrimeFactorization:
+Base.copy(a::PrimeFactorization) = PrimeFactorization(copy(a.powers), a.sign)
+
+Base.one(::Type{PrimeFactorization{T}}) where {T<:Unsigned} = PrimeFactorization(Vector{T}())
+Base.zero(::Type{PrimeFactorization{T}}) where {T<:Unsigned} = PrimeFactorization(Vector{T}(), zero(Int8))
+
+Base.promote_rule(P::Type{<:PrimeFactorization},::Type{<:Integer}) = P
+Base.promote_rule(P::Type{<:PrimeFactorization},::Type{BigInt}) = BigInt
+
+Base.convert(P::Type{<:PrimeFactorization}, n::Integer) = convert(P, primefactor(n))
 function Base.convert(::Type{BigInt}, a::PrimeFactorization)
     A = big(a.sign)
     for (n, e) in enumerate(a.powers)
@@ -114,6 +120,24 @@ function Base.convert(::Type{BigInt}, a::PrimeFactorization)
     end
     return A
 end
+Base.convert(::Type{PrimeFactorization{T}}, a::PrimeFactorization{T}) where {T<:Unsigned} = a
+Base.convert(::Type{PrimeFactorization{T1}}, a::PrimeFactorization{T2}) where {T1<:Unsigned, T2<:Unsigned} = PrimeFactorization(map(T1, a.powers), a.sign)
+
+Base.:(==)(a::PrimeFactorization, b::PrimeFactorization) = a.powers == b.powers && a.sign == b.sign
+function Base.:<(a::PrimeFactorization, b::PrimeFactorization)
+    if a.sign != b.sign
+        return a.sign < b.sign
+    elseif a.sign < 0
+        return <(-b, -a)
+    else
+        ag, bg = divgcd(a, b)
+        if length(ag.powers) <= length(bg.powers) && all(k->ag.powers[k]<bg.powers[k], 1:length(ag.powers))
+            return true
+        else
+            return convert(BigInt, ag) < convert(BigInt, bg)
+        end
+    end
+end
 
 # Methods for PrimeFactorization: only fast multiplication, and lcm and gcd.
 # Addition and subtraction will require conversion to BigInt
@@ -181,15 +205,15 @@ end
 # given a list of numerators and denominators, compute the common denominator and
 # the rescaled numerator after putting all fractions at the same common denominator
 function commondenominator!(nums::Vector{P}, dens::Vector{P}) where {P<:PrimeFactorization}
+    isempty(nums) && return one(P)
     # accumulate lcm of denominator
-    den = copy(dens[1])
+    den = PrimeFactorization(copy(dens[1].powers))
     for i = 2:length(dens)
         _vmax!(den.powers, dens[i].powers)
     end
     # rescale numerators
     for i = 1:length(nums)
-        powers = _vsub!(_vadd!(nums[i].powers, den.powers), dens[i].powers)
-        nums[i] = PrimeFactorization(powers, nums[i].sign)
+        _vsub!(_vadd!(nums[i].powers, den.powers), dens[i].powers)
     end
     return den
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index e051425..665dbb5 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,96 @@
 using WignerSymbols
 using Base.Test
 
-# write your own tests here
-@test 1 == 2
+smalljlist = 0:1//2:10
+largejlist = 0:1//2:1000
+# test triangle coefficient
+for j1 in smalljlist, j2 in smalljlist
+    for j3 = abs(j1-j2):(j1+j2)
+        @test Δ(j1,j2,j3) ≈ sqrt(factorial(float(j1+j2-j3))*factorial(float(j1-j2+j3))*
+                                    factorial(float(j2+j3-j1))/factorial(float(j1+j2+j3+1)))
+    end
+end
+
+# test 3j:
+#--------
+# test cg orthogonality
+for j1 in smalljlist, j2 in smalljlist
+    d1 = 2*j1+1
+    d2 = 2*j2+1
+    M = zeros(d1*d2, d1*d2)
+    ind1 = 1
+    for m1 in -j1:j1, m2 in -j2:j2
+        ind2 = 1
+        for j3 in abs(j1-j2):(j1+j2)
+            for m3 in -j3:j3
+                M[ind1,ind2] = clebschgordan(j1,m1,j2,m2,j3,m3)
+                ind2 += 1
+            end
+        end
+        ind1 += 1
+    end
+    @test M'*M ≈ one(M)
+end
+
+# test recurrence relations: Phys Rev E 57, 7274 (1998)
+for k = 1:10
+    j2 = convert(BigFloat,rand(0:1//2:1000))
+    j3 = convert(BigFloat,rand(0:1//2:1000))
+    m2 = convert(BigFloat,rand(-j2:j2))
+    m3 = convert(BigFloat,rand(-j3:j3))
+
+    for j in max(abs(j2-j3),abs(m2+m3))+1:(j2+j3)-1
+        X = j*sqrt(((j+1)^2-(j2-j3)^2)*((j2+j3+1)^2-(j+1)^2)*((j+1)^2-(m2+m3)^2))
+        Y = (2*j+1)*((m2+m3)*(j2*(j2+1)-j3*(j3+1)) - (m2-m3)*j*(j+1))
+        Z = (j+1)*sqrt((j^2-(j2-j3)^2)*((j2+j3+1)^2-j^2)*(j^2-(m2+m3)^2))
+        tol = 10*max(abs(X),abs(Y),abs(Z))*eps(BigFloat)
+        @test (X*wigner3j(BigFloat,j+1,j2,j3,-m2-m3,m2,m3) + Z*wigner3j(BigFloat,j-1,j2,j3,-m2-m3,m2,m3))≈(-Y*wigner3j(BigFloat,j,j2,j3,-m2-m3,m2,m3)) atol=tol
+    end
+end
+
+# test 6j
+#----------
+# test orthogonality
+for j1 in smalljlist, j2 in smalljlist, j4 in smalljlist
+    for j5 in max(abs(j1-j2-j4),abs(j1-j2+j4),abs(j1+j2-j4)):(j1+j2+j4)
+        j6range = max(abs(j2-j4),abs(j1-j5)):min((j2+j4),(j1+j5))
+        j3range = max(abs(j1-j2),abs(j4-j5)):min((j1+j2),(j4+j5))
+        @test length(j6range) == length(j3range)
+        M = zeros(length(j3range),length(j6range))
+        for (k2,j6) in enumerate(j6range)
+            for (k1,j3) in enumerate(j3range)
+                M[k1,k2] = sqrt(2*j3+1)*sqrt(2*j6+1)*wigner6j(j1,j2,j3,j4,j5,j6)
+            end
+        end
+        @test M'*M ≈ one(M)
+    end
+end
+
+# test special case
+for j1 in smalljlist, j2 in smalljlist
+    j6 = 0
+    j4 = j2
+    j5 = j1
+    for j3 in abs(j1-j2):(j1+j2)
+        @test wigner6j(j1,j2,j3,j4,j5,j6) ≈ (-1)^(j1+j2+j3)/sqrt((2*j1+1)*(2*j2+1))
+    end
+end
+
+# test recurrence relations: Phys Rev E 57, 7274 (1998)
+for k = 1:10
+    j2 = convert(BigFloat,rand(largejlist))
+    j3 = convert(BigFloat,rand(largejlist))
+    l1 = convert(BigFloat,rand(largejlist))
+    l2 = convert(BigFloat,rand(abs(l1-j3):(l1+j3)))
+    l3 = convert(BigFloat,rand(abs(l1-j2):min(l1+j2)))
+
+    for j in intersect(abs(j2-j3):(j2+j3), abs(l2-l3):(l2+l3))
+        X = j*sqrt(((j+1)^2-(j2-j3)^2)*((j2+j3+1)^2-(j+1)^2)*((j+1)^2-(l2-l3)^2)*((l2+l3+1)^2 - (j+1)^2))
+        Y = (2*j+1)*( j*(j+1)*( -j*(j+1) + j2*(j2+1) + j3*(j3+1) - 2*l1*(l1+1)) +
+                        l2*(l2+1)*( j*(j+1) + j2*(j2+1) - j3*(j3+1) ) +
+                        l3*(l3+1)*( j*(j+1) - j2*(j2+1) + j3*(j3+1) ) )
+        Z = (j+1)*sqrt((j^2-(j2-j3)^2)*((j2+j3+1)^2-j^2)*(j^2-(l2-l3)^2)*((l2+l3+1)^2 - j^2))
+        tol = 10*max(abs(X),abs(Y),abs(Z))*eps(BigFloat)
+        @test (X*wigner6j(BigFloat,j+1,j2,j3,l1,l2,l3) + Z*wigner6j(BigFloat,j-1,j2,j3,l1,l2,l3))≈(-Y*wigner6j(BigFloat,j,j2,j3,l1,l2,l3)) atol=tol
+    end
+end