diff --git a/Project.toml b/Project.toml
index 359b24e..cdc34fc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -53,7 +53,7 @@ GraphPlot = "0.6.0"
 Graphs = "1.12.0"
 Makie = "0.19.11"
 NetworkLayout = "0.4.9"
-Oscar = "1.2.2"
+Oscar = "1.5.0"
 StatsBase = "0.34.4"
 WGLMakie = "0.8.15"
 julia = "≥ 1.10"
diff --git a/ext/MakieExt/LDPC/cycles.jl b/ext/MakieExt/LDPC/cycles.jl
index f74df9f..4f441ce 100644
--- a/ext/MakieExt/LDPC/cycles.jl
+++ b/ext/MakieExt/LDPC/cycles.jl
@@ -21,7 +21,7 @@ enumerated. An empty figure and dictionary are returned when there are no cycles
   already cached.
 - Run `using Makie` to activate this extension.
 """
-function CodingTheory.simple_cycle_length_distribution_plot(L::AbstractLDPCCode; len::Int = -1)
+function CodingTheory.simple_cycle_length_distribution_plot(L::AbstractLDPCCode; len::Int = 16)
     dist = simple_cycle_length_distribution(L, len = len)
     x_data = collect(keys(dist))
     y_data = collect(values(dist))
diff --git a/src/Classical/cyclic_code.jl b/src/Classical/cyclic_code.jl
index 7c93198..b71307d 100644
--- a/src/Classical/cyclic_code.jl
+++ b/src/Classical/cyclic_code.jl
@@ -550,6 +550,42 @@ function _generator_matrix(F::FqField, n::Int, k::Int, g::FqPolyRingElem)
     return G
 end
 
+function _classify_factors(poly::fpPolyRingElem)
+    # the problem with reverse(poly) == poly is that coefficents(poly) is not a fixed size
+    n = degree(poly)
+    F = base_ring(parent(poly))
+    F0 = F(0)
+    facs = [x[1] for x in collect(factor(poly))]
+    self_reciprocal_facs = Vector{fpPolyRingElem}()
+    pairs = Vector{Tuple{fpPolyRingElem, fpPolyRingElem}}()
+    temp = zeros(F, 1, n + 1)
+    temp2 = zeros(F, 1, n + 1)
+    for f in facs
+        f_coeffs = collect(coefficients(f))
+        temp[1, end - length(f_coeffs) + 1:end] .= reverse(f_coeffs)
+        temp2[1, end - length(f_coeffs) + 1:end] .= f_coeffs
+        if temp == temp2
+            push!(self_reciprocal_facs, f)
+        else
+            for f2 in facs
+                if f != f2
+                    f_coeffs = collect(coefficients(f2))
+                    temp2[1, 1:end - length(f_coeffs)] .= F0
+                    temp2[1, end - length(f_coeffs) + 1:end] .= f_coeffs
+                    if temp == temp2
+                        push!(pairs, (f, f2))
+                        break
+                    end
+                end
+            end
+        end
+        temp[1, :] .= F0
+        temp2[1, :] .= F0
+    end
+
+    return self_reciprocal_facs, pairs
+end
+
 # TODO: make flat optional throughout
  """
     defining_set(nums::Vector{Int}, q::Int, n::Int, flat::Bool = true)
diff --git a/src/Classical/linear_code.jl b/src/Classical/linear_code.jl
index 0f02922..2c734e9 100644
--- a/src/Classical/linear_code.jl
+++ b/src/Classical/linear_code.jl
@@ -852,7 +852,7 @@ function _are_perm_equivalent_exhaustive_search(C1::AbstractLinearCode, C2::Abst
     sym = symmetric_group(nc)
     for e in collect(sym) 
         P = permutation_matrix(C1.F, e)
-        if G1 * P == G2 
+        if G1 * P == G2
             return (true, P)
         end
     end
diff --git a/src/Classical/weight_dist.jl b/src/Classical/weight_dist.jl
index 5fc6093..1e40c67 100644
--- a/src/Classical/weight_dist.jl
+++ b/src/Classical/weight_dist.jl
@@ -419,6 +419,7 @@ end
 
 """
     minimum_distance_Gray(C::AbstractLinearCode; alg::Symbol = :Zimmermann, v::Bool = false)
+    
 Return the minimum distance of `C` using a deterministic algorithm based on enumerating
 constant weight codewords of the binary reflected Gray code. If a word of minimum weight
 is found before the lower and upper bounds cross, it is returned; otherwise, the zero vector 
@@ -427,7 +428,7 @@ is returned.
 show_progress will display a progress meter for each iteration of a weight that takes longer than 
     a second
 """
-function minimum_distance_Gray(C::AbstractLinearCode; alg::Symbol = :auto, v::Bool = false, 
+function minimum_distance_Gray(C::AbstractLinearCode; alg::Symbol = :auto, verbose::Bool = false, 
     show_progress = true)
 
     ord_F = Int(order(C.F))
@@ -443,39 +444,50 @@ function minimum_distance_Gray(C::AbstractLinearCode; alg::Symbol = :auto, v::Bo
 
     if alg == :auto
         if typeof(C) <: AbstractCyclicCode
-            v && println("Detected a cyclic code, using Chen's adaption.")
+            verbose && println("Detected a cyclic code, using Chen's adaption.")
             alg = :Chen
             # TODO: fix this case
         elseif typeof(C) <: AbstractQuasiCyclicCode
-            v && println("Detected a quasi-cyclic code, using White's adaption.")
+            verbose && println("Detected a quasi-cyclic code, using White's adaption.")
             alg = :White
         else
-            v && println("Using Zimmermann's algorithm.")
+            verbose && println("Using Zimmermann's algorithm.")
             alg = :Zimmermann
         end
     end
+    # TODO should never hit this case with the else in there?
     alg == :auto && throw(ErrorException("Could not determine minimum distance algo automatically"))
 
-    if alg in (:Brouwer, :Zimmermann) 
-        return _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg = alg, verbose = v, show_progress = show_progress)
+    # really this should work for all just the same?
+    if alg in (:Brouwer, :Zimmermann)
+        return _minimum_distance_BZ(C::AbstractLinearCode, alg, verbose, show_progress)
     end
     println("Warning: old enumeration algorithm selected. Performance will be slow") # TODO remove when all code updated
     return _minimum_distance_enumeration_with_matrix_multiply(C::AbstractLinearCode; info_set_alg = alg)
 end
 
-function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zimmermann,
-    verbose::Bool = false, dbg = Dict(), show_progress=false)
+# this is a private function, so there's no point in having keyword arguments
+function _minimum_distance_BZ(C::AbstractLinearCode, info_set_alg::Symbol, verbose::Bool,
+    show_progress::Bool)
+    
+    # you never pass this in anywhere so moving out
+    dbg::Dict{String, Int} = Dict()
     dbg_key_exit_r = "exit_r"
 
-    ord_F = Int(order(C.F))
-    ord_F == 2 || throw(ArgumentError("Currently only implemented for binary codes."))
+    # removing since checked in main function
+    # ord_F = Int(order(C.F))
+    # ord_F == 2 || throw(ArgumentError("Currently only implemented for binary codes."))
+    # put this at user entry point, although if 2^k is large, 2^(n - k) maybe small enough to do
+    # so we should put this in an if statement above if we are doing G or H
     C.k < 2^16 || throw(DomainError("The given linear code has length k >= 2^16 which is not supported"))
+
     info_set_alg ∈ (:Brouwer, :Zimmermann) || throw(ArgumentError("Unknown information set algorithm. Expected `:Brouwer`, `:Zimmermann`"))
 
     generator_matrix(C, true) # ensure G_stand exists
     if _has_empty_vec(C.G, :cols) 
         #TODO err string can instruct the user to construct a new code without 0 cols and tell them the function for that
-        throw(ArgumentError("Codes with standard form of generator matrix having 0 columns not supported")) 
+        throw(ArgumentError("Codes with standard form of generator matrix having 0 columns not supported"))
+        # is that actually possible? I thought we remove this in the linear code constructor? If we don't we should and then remove this here
     end
     # generate if not pre-stored
     parity_check_matrix(C)
@@ -494,20 +506,22 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim
     end
 
     if haskey(dbg, dbg_key_exit_r)
-        verbose && println("dbg Dict: largest message weight searched stored @key=$dbg_key_exit_r")
+        verbose && println("dbg Dict: largest message weight searched stored @key = $dbg_key_exit_r")
         dbg[dbg_key_exit_r] = -1
     end
 
     k, n = size(C.G)
-    A_mats_trunc = [Matrix{UInt16}(undef, k, n-k) for _ in 1:length(A_mats)]
-    for i in 1:size(A_mats, 1) 
-       A_mats_trunc[i] = deepcopy(A_mats[i][k+1 : n, :])
+    A_mats_trunc = [Matrix{UInt16}(undef, k, n - k) for _ in 1:length(A_mats)]
+    for i in 1:size(A_mats, 1)
+       A_mats_trunc[i] = deepcopy(A_mats[i][k + 1 : n, :])
     end
 
+    # I don't remember this case in the paper
     if info_set_alg == :Brouwer && rnks[h] != k
         println("Rank of last matrix too small")
         return
     end
+    # I think we have to remove this from the verbose statement to work
     if verbose
         print("Generated $h information sets with ranks: ")
         for i in 1:h
@@ -532,10 +546,10 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim
         (!triply_even_flag && !doubly_even_flag && even_flag) && println("Detected an even code.")
     end
 
-    # initial_perm_ind will match the permutation we use for the 'found' vector if the found vector is nonzero. To simplify the code below we're going to choose an initial permutation arbitrarily.  
-    initial_perm_ind = 1 
-    # following loop is the r=1 case of the enumeration. We do this case here because we want to make a good guess at the terminating r before we start multiple threads
-    for (j, g) in enumerate(A_mats) # loop over the A_mats rather than the original G because it would add another case to deal with later 
+    # initial_perm_ind will match the permutation we use for the 'found' vector if the found vector is nonzero. To simplify the code below we're going to choose an initial permutation arbitrarily.
+    initial_perm_ind = 1
+    # following loop is the r = 1 case of the enumeration. We do this case here because we want to make a good guess at the terminating r before we start multiple threads
+    for (j, g) in enumerate(A_mats) # loop over the A_mats rather than the original G because it would add another case to deal with later
         # can make this faster with dots and views
         w, i = _min_wt_col(g)
         if w <= C.u_bound
@@ -546,13 +560,15 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim
     end
 
     verbose && println("Current upper bound: $(C.u_bound)")
+    # I'm actually surpised this works since Oscar now requires [:, 1:1] for most things
     found = A_mats[1][:, 1]
 
     l = 0
-    if verbose 
+    # not entirely sure what's going on here since it's never used later, I guess just for testing, but if we want it printed it never ends up printed with the current comments
+    if verbose
         _, _, b_rnks = information_sets(C.G, :Brouwer, permute = true, only_A = false)
         b_h = length(b_rnks)
-        b_lower_bounds = [information_set_lower_bound(r+1, n, k, l, [k - 0 for i in 1:b_h], :Brouwer, even = even_flag, doubly_even = doubly_even_flag, triply_even = triply_even_flag) for r in 1:k-1]
+        b_lower_bounds = [information_set_lower_bound(r + 1, n, k, l, [k - 0 for i in 1:b_h], :Brouwer, even = even_flag, doubly_even = doubly_even_flag, triply_even = triply_even_flag) for r in 1:k - 1]
         b_r_term = findfirst(x -> x ≥ C.u_bound, b_lower_bounds)
 
         # _, _, z_rnks = information_sets(G, :Zimmermann, permute = true, only_A = false)
@@ -563,38 +579,41 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim
         # verbose && println("Predicted termination weight based on current upper bound: Brouwer $b_r_term Zimm $z_r_term")
     end
 
-    #Note the r+1 here. 
-    lower_bounds_for_prediction = [information_set_lower_bound(r+1, n, k, l, rank_defs, info_set_alg, even = even_flag, doubly_even = doubly_even_flag, triply_even = triply_even_flag) for r in 1:k-1]
+    # Note the r + 1 here
+    lower_bounds_for_prediction = [information_set_lower_bound(r + 1, n, k, l, rank_defs, info_set_alg, even = even_flag, doubly_even = doubly_even_flag, triply_even = triply_even_flag) for r in 1:k - 1]
     r_term = findfirst(x -> x ≥ C.u_bound, lower_bounds_for_prediction)
     if isnothing(r_term)
-        raise(DomainError("invalid termination r")) 
+        raise(DomainError("Invalid termination r")) 
     end
     verbose && println("Predicted termination weight based on current upper bound: $r_term")
 
-    #In the main loop we check if lower bound > upper bound before we enumerate and so the lower bounds for the loop use r not r+1
-    lower_bounds = [information_set_lower_bound(r, n, k, l, rank_defs, info_set_alg, even = even_flag, doubly_even = doubly_even_flag, triply_even = triply_even_flag) for r in 1:k-1]
+    # In the main loop we check if lower bound > upper bound before we enumerate and so the lower bounds for the loop use r not r + 1
+    lower_bounds = [information_set_lower_bound(r, n, k, l, rank_defs, info_set_alg, even = even_flag, doubly_even = doubly_even_flag, triply_even = triply_even_flag) for r in 1:k - 1]
 
     predicted_work_factor = fld(n, k) * sum([binomial(k, i) for i in 1:r_term])
     verbose && println("Predicted work factor: $predicted_work_factor")
     if show_progress 
-        prog_bar = Progress(predicted_work_factor, dt=1.0, showspeed=true) # updates no faster than once every 1s
+        prog_bar = Progress(predicted_work_factor, dt = 1.0, showspeed = true) # updates no faster than once every 1s
     end
-    weight_sum_bound = min(2 * C.u_bound + 5, n-k)
+    # don't understand this but okay
+    weight_sum_bound = min(2 * C.u_bound + 5, n - k)
     verbose && println("Codeword weights initially checked on first $weight_sum_bound entries")
 
     num_thrds = Threads.nthreads()
     verbose && println("Number of threads ", num_thrds)
     for r in 2:k
+        # TODO this case should never happen given the error above?
         if r > 2^16
             verbose && println("Warning: Reached an r larger than 2^16") 
         end
         C.l_bound < lower_bounds[r] && (C.l_bound = lower_bounds[r];)
+        # I assume we want to uncomment these?
         # an even code can't have have an odd minimum weight
         # (!triply_even_flag && !doubly_even_flag && even_flag) && (C.l_bound += C.l_bound % 2;)
         # (!triply_even_flag && doubly_even_flag) && (C.l_bound += 4 - C.l_bound % 4;)
         # triply_even_flag && (C.l_bound += 8 - C.l_bound % 8;)
         if C.l_bound >= C.u_bound
-            dbg[dbg_key_exit_r] = r-1
+            dbg[dbg_key_exit_r] = r - 1
             break
         end
         verbose && println("r: $r")
@@ -608,6 +627,7 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim
             end
             i_count > 0 && println("$i_count of the original $h information sets no longer contribute to the lower bound")
         end
+        # wait... does this work for nonbinary? cause we do throw an error above if it's not binary
         p = Int(characteristic(C.F))
 
         uppers = [C.u_bound for _ in 1:num_thrds]
@@ -630,17 +650,19 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim
                     first_vec[i] = 1
                 end
             else
+                # we shouldn't need the CT. here anymore
                 CodingTheory._subset_unrank_to_vec!(init_rank, UInt64(r), first_vec)
             end
  
-            # as in White Algo 7.1 we loop over matrices first 
+            # I couldn't figure out why we want to do this but sure
+            # as in White Algo 7.1 we loop over matrices first
             for i in 1:h
                 if keep_going[] == false
                     verbose && println(thrd_stop_msg)
                     break
                 end
 
-                c_itr = zeros(UInt16, C.n - C.k) 
+                c_itr = zeros(UInt16, C.n - C.k)
                 is_first = true
                 curr_mat = A_mats_trunc[i]
                 count = UInt128(0)
@@ -650,11 +672,12 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim
                         println(thrd_stop_msg)
                         break
                     end
-                    show_progress && ProgressMeter.next!(prog_bar) 
+                    show_progress && ProgressMeter.next!(prog_bar)
+
                     if r - rank_defs[i] > 0
                         if is_first 
                             LinearAlgebra.mul!(c_itr, curr_mat, first_vec)
-                            @inbounds @simd for j in eachindex(c_itr) 
+                            @inbounds @simd for j in eachindex(c_itr)
                                 c_itr[j] %= p
                             end
                             is_first = false
@@ -675,21 +698,24 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim
                             verbose && @assert w != 0
                             if uppers[ind] > w 
                                 subset_vec_full = zeros(Int, k)
+                                # CT.
                                 CodingTheory._subset_unrank_to_vec!(UInt128(init_rank + count), UInt64(r), subset_vec_full)
 
                                 uppers[ind] = w 
                                 founds[ind] = vcat(subset_vec_full, c_itr)
-                                verbose && @assert size(founds[ind], 1) == C.n "found vector has length $(size(founds[ind], 1)) but should be n=$(C.n)"
-                                exit_thread_indicator_vec[ind] = i 
+                                # is it possible to hit this assert if programmed correctly and checked on one or two runs?
+                                verbose && @assert size(founds[ind], 1) == C.n "found vector has length $(size(founds[ind], 1)) but should be n = $(C.n)"
+                                exit_thread_indicator_vec[ind] = i
 
-                                println("Adjusting (local) upper bound: $w for c_itr=$(Int.(c_itr))")
+                                println("Adjusting (local) upper bound: $w for c_itr = $(Int.(c_itr))")
                                 if C.l_bound == uppers[ind]
-                                    println("early exit")
+                                    verbose && println("Early exit")
                                     Threads.atomic_cas!(keep_going, true, false)
                                 else
                                     r_term = findfirst(x -> x ≥ C.u_bound, lower_bounds)
                                     isnothing(r_term) && (r_term = k;)
                                     verbose && println("Updated termination weight: $r_term")
+                                    # can we update the progress meter size here or just maybe add an amount equal to the amount now skipped?
                                 end
                             end
                         end
@@ -697,7 +723,7 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim
                     count = add!(count, count, 1)
                 end
             end
-            loc = argmin(uppers) 
+            loc = argmin(uppers)
             C.u_bound = uppers[loc]
             found = founds[loc]
             initial_perm_ind = exit_thread_indicator_vec[loc]
@@ -705,11 +731,12 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim
     end
 
     C.d = C.u_bound
+    # why would this not always be the proper wt?
     y = matrix(C.F, 1, n, perms_mats[initial_perm_ind] * found) # weight(y) >= C.d, with equality not being the typical case
     verbose && @assert iszero(C.H * transpose(y))
     if dbg[dbg_key_exit_r] == -1
         dbg[dbg_key_exit_r] = r
-    end 
+    end
     show_progress && ProgressMeter.finish!(prog_bar)
     verbose && println("Computation complete")
     return C.u_bound, y
diff --git a/src/CodingTheory.jl b/src/CodingTheory.jl
index 746e530..63fb7b8 100644
--- a/src/CodingTheory.jl
+++ b/src/CodingTheory.jl
@@ -54,6 +54,8 @@ const CTPolyRingElem = PolyRingElem{<:CTFieldElem}
 const CTGroupAlgebra = GroupAlgebraElem{fpFieldElem, GroupAlgebra{fpFieldElem, FinGenAbGroup, FinGenAbGroupElem}}
 const CTChainComplex = Union{ComplexOfMorphisms{AbstractAlgebra.FPModule{fpFieldElem}}} # residue and group algebras later
 const CTPolyMatrix = Union{AbstractAlgebra.Generic.MatSpaceElem{fpPolyRingElem}, AbstractAlgebra.Generic.MatSpaceElem{FqPolyRingElem}}
+const CTLRPolyElem = AbstractAlgebra.Generic.LaurentMPolyWrap{fpFieldElem, fpMPolyRingElem,
+AbstractAlgebra.Generic.LaurentMPolyWrapRing{fpFieldElem, fpMPolyRing}}
 
 include("Classical/types.jl")
 export AbstractCode, AbstractNonadditiveCode, AbstractNonlinearCode, AbstractAdditiveCode,
@@ -72,7 +74,8 @@ include("Quantum/types.jl")
 export AbstractSubsystemCode, AbstractSubsystemCodeCSS, AbstractStabilizerCode, AbstractStabilizerCodeCSS,
     AbstractGraphStateSubsystem, AbstractGraphStateSubsystemCSS, AbstractGraphStateStabilizer,
     AbstractGraphStateStabilizerCSS, AbstractHypergraphProductCode, AbstractEASubsystemCode,
-    AbstractEASubsystemCodeCSS, AbstractEAStabilizerCode, AbstractEAStabilizerCodeCSS 
+    AbstractEASubsystemCodeCSS, AbstractEAStabilizerCode, AbstractEAStabilizerCodeCSS, AbstractGeneralizedToricCode
+
 # misc
 export LogicalTrait, GaugeTrait, CSSTrait, HasLogicals, HasNoLogicals, HasGauges, HasNoGauges,
     IsCSS, IsNotCSS, copy, ChainComplex
@@ -95,7 +98,7 @@ export kronecker_product, Hamming_weight, weight, wt, Hamming_distance, distance
     is_valid_bipartition, extract_bipartition, is_Hermitian_self_orthogonal,
     row_supports, row_supports_symplectic, strongly_lower_triangular_reduction,
     residue_polynomial_to_circulant_matrix, group_algebra_element_to_circulant_matrix,
-    load_alist, extended_binomial
+    load_alist, extended_binomial, mult_order, prepare_mult_order
     # load_alist, _rref_non_pivot_cols
     # , _min_wt_row
     # , circ_shift
@@ -473,4 +476,12 @@ export copying, gauging, thickening_and_choose_heights, coning, quantum_weight_r
 include("Quantum/homological_measurements.jl")
 export homological_measurement, Cheeger_constant
 
+#############################
+# Quantum/GeneralizedToricCode.jl
+#############################
+
+include("Quantum/GeneralizedToricCode.jl")
+export GeneralizedToricCode, FiniteGeneralizedToricCode, maximum_dimension,
+    Laurent_polynomial_ring, defining_polynomials, twist_vectors
+
 end
diff --git a/src/LDPC/MP_decoders.jl b/src/LDPC/MP_decoders.jl
index b79679c..2276311 100644
--- a/src/LDPC/MP_decoders.jl
+++ b/src/LDPC/MP_decoders.jl
@@ -1139,12 +1139,12 @@ function _message_passing_fast(H_Int::Matrix{UInt8}, v::Matrix{UInt8}, syndrome_
             current_bits[i] = var_to_check_messages[i, curr_iter] >= 0 ? 0 : 1
         end
 
-        LinearAlgebra.mul!(syn, H_Int, current_bits)
-        if syndrome_based
-            all(syn[i] % 2 == v[i, 1] for i in 1:num_check) && return true, current_bits, iter, var_to_check_messages[:, curr_iter]
-        else
-            all(iszero(syn[i] % 2) for i in 1:num_check) && return true, current_bits, iter, var_to_check_messages[:, curr_iter]
-        end
+        # LinearAlgebra.mul!(syn, H_Int, current_bits)
+        # if syndrome_based
+        #     all(syn[i] % 2 == v[i, 1] for i in 1:num_check) && return true, current_bits, iter, var_to_check_messages[:, curr_iter]
+        # else
+        #     all(iszero(syn[i] % 2) for i in 1:num_check) && return true, current_bits, iter, var_to_check_messages[:, curr_iter]
+        # end
 
         if schedule == :parallel
             temp_iter = curr_iter
diff --git a/src/LDPC/cycles.jl b/src/LDPC/cycles.jl
index c200ff0..ab780ab 100644
--- a/src/LDPC/cycles.jl
+++ b/src/LDPC/cycles.jl
@@ -333,7 +333,7 @@ is returned when there is no cycles.
   columns of `parity_check_matrix(L)` then top-to-bottom by rows.
 """
 function enumerate_simple_cycles(L::AbstractLDPCCode; len::Int = 16)
-    ispositive(len) || throw(DomainError("Cycle length parameter must be a positive integer"))
+    is_positive(len) || throw(DomainError("Cycle length parameter must be a positive integer"))
 
     # TODO add this variable to struct
     if len > L.max_cyc_len
diff --git a/src/Lattices/Voronoi.jl b/src/Lattices/Voronoi.jl
new file mode 100644
index 0000000..ecec9b4
--- /dev/null
+++ b/src/Lattices/Voronoi.jl
@@ -0,0 +1,9 @@
+# Copyright (c) 2025 David Marquis
+# All rights reserved.
+#
+# This source code is licensed under the BSD-style license found in the
+# LICENSE file in the root directory of this source tree.
+
+#############################
+     # general functions
+#############################
diff --git a/src/Lattices/lattices.jl b/src/Lattices/lattices.jl
new file mode 100644
index 0000000..737060b
--- /dev/null
+++ b/src/Lattices/lattices.jl
@@ -0,0 +1,232 @@
+# Copyright (c) 2025 Eric Sabo, David Marquis
+# All rights reserved.
+#
+# This source code is licensed under the BSD-style license found in the
+# LICENSE file in the root directory of this source tree.
+
+#############################
+        # constructors
+#############################
+
+# TODO do we even want to accept a Gram matrix here?
+"""
+    Lattice(B::MatrixElem{<:RingElem}; type::Symbol = :basis)
+
+Return a lattice using the columns of `B` as a basis.
+"""
+function Lattice(B::MatrixElem{<:RingElem}; type::Symbol = :basis)
+    if type == :basis
+        rank(B) == ncols(B) || throw(ArgumentError("Basis must have full column rank."))
+        return Lattice(B)
+    elseif type == :Gram
+        # which decomp to use here?
+
+        return Lattice(B)
+    else
+        throw(ArgumentError("Parameter `type` must be either `:basis` or `:Gram`."))
+    end
+end
+
+"""
+    Lattice(V::Vector{MatrixElem{<:RingElem}})
+
+Return a lattice using the elements of the vectors as a basis.
+"""
+function Lattice(V::Vector{MatrixElem{<:RingElem}})
+    # check equal sizes and base rings
+    # stack rows into columns of matrix
+
+    rank(B) == ncols(B) || throw(ArgumentError("Basis must have full column rank."))
+    return Lattice(B)
+end
+
+"""
+
+
+"""
+function Lattice(G::MatrixElem{<:RingElem})
+
+
+    return Lattice(B)
+end
+
+"""
+    construction_A(C::AbstractLinearCode)
+
+Return the lattice given by applying Construction A to the code `C`.
+"""
+function construction_A(C::AbstractLinearCode)
+
+
+end
+
+"""
+    construction_B(C::AbstractLinearCode)
+
+Return the lattice given by applying Construction B to the code `C`.
+"""
+function construction_B(C::AbstractLinearCode)
+
+
+end
+
+"""
+    construction_C(C::AbstractLinearCode)
+
+Return the lattice given by applying Construction C to the code `C`.
+"""
+function construction_C(C::AbstractLinearCode)
+
+
+end
+
+#############################
+      # getter functions
+#############################
+
+"""
+    basis(L::Lattice)
+
+Return the basis of the lattice.
+"""
+basis(L::Lattice) = L.B
+
+"""
+    rank(L::Lattice)
+
+Return the rank of the lattice.
+"""
+rank(L::Lattice) = size(L.B, 1)
+
+"""
+    dimension(L::Lattice)
+    dim(L::Lattice)
+
+Return the dimension of the lattice.
+"""
+dimension(L::Lattice) = size(L.B, 2)
+dim(L::Lattice) = dimension(L)
+
+"""
+    is_full_rank(L::Lattice)
+
+Return `true` if the lattice is full rank; otherwise return `false`.
+"""
+is_full_rank(L::Lattice) = dimension(L) == rank(L)
+
+"""
+    discriminant(L::Lattice)
+    disc(L::Lattice)
+
+Return the discriminant of the lattice.
+"""
+discriminant(L::Lattice) = det(Gram_matrix(L))
+disc(L::Lattice) = discriminant(L)
+
+"""
+    Gram_matrix(L::Lattice)
+
+Return the Gram matrix of the basis of the lattice.
+"""
+Gram_matrix(L::Lattice) = gram(L.B)
+
+"""
+    volume(L::Lattice)
+    vol(L::Lattice)
+    determinate(L::Lattice)
+    det(L::Lattice)
+
+Return the volumne of the lattice.
+"""
+volume(L::Lattice) = sqrt(disc(L))
+vol(L::Lattice) = volume(L)
+determinate(L::Lattice) = volume(L)
+det(L::Lattice) = volume(L)
+
+"""
+    check_matrix(L::Lattice)
+
+Return the check matrix of the lattice.
+"""
+check_matrix(L::Lattice) = inverse(L.B)
+
+"""
+    shortest_vector(L::Lattice)
+
+"""
+shortest_vector(L::Lattice)
+
+"""
+    successive_minima(L::Lattice)
+
+"""
+successive_minima(L::Lattice)
+
+"""
+    dual_basis(L::Lattice)
+
+Return a dual basis for the lattice.
+"""
+function dual_basis(L::Lattice)
+    Ginv = inverse(Gram_matrix(L))
+    new_basis = Ginv * B
+    # b^*_i = \sum_j (G_inv)_ij b_j
+    return new_basis
+end
+
+"""
+    dual(L::Lattice)
+
+Return the dual lattice.
+"""
+dual(L::Lattice) = Lattice(dual_basis(L))
+
+"""
+    is_unimodular(L::Lattice)
+
+Return `true` if the lattie is unimodular; otherwise return `false`.
+"""
+function is_unimodular(L::Lattice)
+    # Gram matrix has integer entries and disc(L) == 1
+end
+
+#############################
+     # general functions
+#############################
+
+"""
+
+
+"""
+function in(v::MatrixElem{<:RingElem}, L::Lattice)
+    # check sizes and rings
+    # setup Bx = v
+    # true if solution exists and is in ZZ
+end
+
+function is_sublattice(L1::Lattice, L2::Lattice)
+
+end
+
+function is_equal(L1::Lattice, L2::Lattice)
+    # L1 ⊆ L2 && L2 ⊆ L1
+
+end
+
+"""
+
+
+"""
+function coding_gain(L::Lattice)
+    
+end
+
+"""
+
+
+"""
+function covering_radius(L::Lattice)
+
+end
+
+# SVP, CVP
diff --git a/src/Lattices/misc_known_lattices.jl b/src/Lattices/misc_known_lattices.jl
new file mode 100644
index 0000000..be4da91
--- /dev/null
+++ b/src/Lattices/misc_known_lattices.jl
@@ -0,0 +1,5 @@
+# Copyright (c) 2025 David Marquis
+# All rights reserved.
+#
+# This source code is licensed under the BSD-style license found in the
+# LICENSE file in the root directory of this source tree.
diff --git a/src/Lattices/types.jl b/src/Lattices/types.jl
new file mode 100644
index 0000000..bdd559a
--- /dev/null
+++ b/src/Lattices/types.jl
@@ -0,0 +1,19 @@
+# Copyright (c) 2025 Eric Sabo, David Marquis
+# All rights reserved.
+#
+# This source code is licensed under the BSD-style license found in the
+# LICENSE file in the root directory of this source tree.
+
+#############################
+      # abstract types
+#############################
+
+abstract type AbstractLattice end
+
+#############################
+      # concrete types
+#############################
+
+struct Lattice <: AbstractLattice
+    B::MatrixElem{<:RingElem}
+end
diff --git a/src/Quantum/GeneralizedToricCode.jl b/src/Quantum/GeneralizedToricCode.jl
new file mode 100644
index 0000000..be5312a
--- /dev/null
+++ b/src/Quantum/GeneralizedToricCode.jl
@@ -0,0 +1,520 @@
+# Copyright (c) 2025 Eric Sabo
+# All rights reserved.
+#
+# This source code is licensed under the BSD-style license found in the
+# LICENSE file in the root directory of this source tree.
+
+#############################
+        # constructors
+#############################
+
+function _BB_ansatz_check(f::T) where T <: Union{MPolyQuoRingElem{FqMPolyRingElem}, MPolyQuoRingElem{fpMPolyRingElem}}
+
+    count_x = 0
+    count_y = 0
+    for e in exponents(f.f)
+        if e[1] != 0 && e[2] != 0
+            return false
+        end
+        if e[1] != 0
+            count_x += 1
+        end
+        if e[2] != 0
+            count_y += 1
+        end
+    end
+    if count_x > 1 && count_y > 1
+        return false
+    end
+    return true
+end
+
+"""
+    BivariateBicycleCode(a::MPolyQuoRingElem{FqMPolyRingElem}, b::MPolyQuoRingElem{FqMPolyRingElem})
+
+Return the bivariate bicycle code defined by the residue ring elements `a` and `b`.
+"""
+function BivariateBicycleCode(a::T, b::T) where T <: Union{MPolyQuoRingElem{FqMPolyRingElem},
+    MPolyQuoRingElem{fpMPolyRingElem}}
+
+    R = parent(a)
+    R == parent(b) || throw(DomainError("Polynomials must have the same parent."))
+    F = base_ring(base_ring(a))
+    order(F) == 2 || throw(DomainError("This code family is currently only defined over binary fields."))
+    length(symbols(parent(a))) == 2 || throw(DomainError("Polynomials must be over two variables."))
+    g = gens(modulus(R))
+    length(g) == 2 || throw(DomainError("Residue rings must have only two generators."))
+
+    m = -1
+    l = -1
+    for g1 in g
+        exps = collect(exponents(g1))
+        length(exps) == 2 || throw(ArgumentError("Moduli of the incorrect form."))
+        iszero(exps[2]) || throw(ArgumentError("Moduli of the incorrect form."))
+        !iszero(exps[1][1]) && !iszero(exps[1][2]) && throw(ArgumentError("Moduli of the incorrect form."))
+        if iszero(exps[1][1])
+            m = exps[1][2]
+        else
+            l = exps[1][1]
+        end
+    end
+
+    if !_BB_ansatz_check(a) || !_BB_ansatz_check(b)
+        throw(ArgumentError("Polynomials do not satisfy the bivariate bicycle code ansatz."))
+    end
+
+    # already has the modulus built into R
+    I = ideal(R, [a, b])
+    Q, _ = quo(R, I)
+    k_dim = 2 * vector_space_dimension(Q)
+
+    return BivariateBicycleCode(R, F, 2 * l * m, k_dim, a, b, l, m)
+end
+
+"""
+    CoprimeBivariateBicycleCode(a::MPolyQuoRingElem{FqMPolyRingElem}, b::MPolyQuoRingElem{FqMPolyRingElem})
+
+Return the coprime bivariate bicycle code defined by the residue ring elements `a` and `b`.
+
+# Note
+
+- This is defined in https://arxiv.org/pdf/2408.10001.
+"""
+function CoprimeBivariateBicycleCode(a::ResElem, b::ResElem)
+    R = parent(a)
+    S = base_ring(a)
+    R == parent(b) || throw(DomainError("Polynomials must have the same parent."))
+    F = base_ring(S)
+    order(F) == 2 || throw(DomainError("This code family is currently only defined over binary fields."))
+    length(gens(S)) == 1 || throw(DomainError("Polynomials must be over one variable."))
+    f = modulus(R)
+    deg_P = degree(f)
+    f == gen(S)^deg_P - 1 || throw(ArgumentError("Residue ring not of the form π^(l * m) - 1."))
+
+    # BUG this is not particularly true since l or m could be factored itself and yet still be coprime
+    facs = Nemo.factor(deg_P)
+    length(facs) == 2 || throw(ArgumentError("Residue ring not of the form π^(l * m) - 1."))
+    k = collect(keys(facs.fac))
+    v = collect(values(facs.fac))
+    l = k[1]^v[1]
+    m = k[2]^v[2]
+
+    k_dim = 2 * degree(gcd(a, b, f))
+    return CoprimeBivariateBicycleCode(R, F, 2 * l * m, k_dim, a, b, l, m)
+end
+
+"""
+    GeneralizedToricCode(f::CTLRPolyElem, g::CTLRPolyElem)
+
+Return the generalized toric code defined by the polynomials `f` and `g`.
+"""
+function GeneralizedToricCode(f::CTLRPolyElem, g::CTLRPolyElem)
+    LR = parent(f)
+    LR == parent(g) || throw(ArgumentError("The polynomials must be over the same ring."))
+    length(symbols(LR)) == 2 || throw(ArgumentError("The polynomials must be over a Laurent polynomial ring in two variables."))
+    # TODO check polynomial ansatz here
+
+    R2 = Oscar._polyringquo(LR)
+    R = codomain(R2)
+    (x, y) = gens(LR)
+    I = ideal(LR, [f, g])
+    II = ideal(R, R2.(gens(I)))
+    Q, _ = quo(R, II)
+    n = 2 * length(monomial_basis(Q))
+
+    return GeneralizedToricCode(LR, base_ring(LR), n, f, g)
+end
+
+"""
+    FiniteGeneralizedToricCode(f::CTLRPolyElem, g::CTLRPolyElem, a1::Tuple{Int, Int}, a2::Tuple{Int, Int})
+
+Return the generalized toric code defined by the polynomials `f` and `g` and twist vectors `a1`, `a2`.
+"""
+function FiniteGeneralizedToricCode(f::CTLRPolyElem, g::CTLRPolyElem, a1::Tuple{Int, Int},
+    a2::Tuple{Int, Int})
+
+    LR = parent(f)
+    LR == parent(g) || throw(ArgumentError("The polynomials must be over the same ring."))
+    length(symbols(LR)) == 2 || throw(ArgumentError("The polynomials must be over a Laurent polynomial ring in two variables."))
+    # TODO check polynomial ansatz here
+
+    R2 = Oscar._polyringquo(LR)
+    R = codomain(R2)
+    (x, y) = gens(LR)
+    I = ideal(LR, [f, g, x^a1[1] * y^a1[2] - 1, x^a2[1] * y^a2[2] - 1])
+    II = ideal(R, R2.(gens(I)))
+    Q, _ = quo(R, II)
+    k_dim = 2 * vector_space_dimension(Q)
+    n = 2 * length(monomial_basis(Q))
+
+    return FiniteGeneralizedToricCode(LR, base_ring(LR), n, k_dim, f, g, a1, a2)
+end
+
+"""
+    MonomialCode(f::CTLRPolyElem, g::CTLRPolyElem)
+
+Return the monomial code defined by the polynomials `f` and `g`.
+"""
+function MonomialCode(f::CTLRPolyElem, g::CTLRPolyElem)
+    LR = parent(f)
+    LR == parent(g) || throw(ArgumentError("The polynomials must be over the same ring."))
+    length(symbols(LR)) == 2 || throw(ArgumentError("The polynomials must be over a Laurent polynomial ring in two variables."))
+
+    R2 = Oscar._polyringquo(LR)
+    R = codomain(R2)
+    (x, y) = gens(LR)
+    I = ideal(LR, [f, g])
+    II = ideal(R, R2.(gens(I)))
+    Q, _ = quo(R, II)
+    n = 2 * length(monomial_basis(Q))
+
+    return MonomialCode(LR, base_ring(LR), n, f, g)
+end
+
+"""
+    FiniteMonomialCode(f::CTLRPolyElem, g::CTLRPolyElem, a1::Tuple{Int, Int}, a2::Tuple{Int, Int})
+
+Return the monomial code defined by the polynomials `f` and `g` and twist vectors `a1`, `a2`.
+"""
+function FiniteMonomialCode(f::CTLRPolyElem, g::CTLRPolyElem, a1::Tuple{Int, Int},
+    a2::Tuple{Int, Int})
+
+    LR = parent(f)
+    LR == parent(g) || throw(ArgumentError("The polynomials must be over the same ring."))
+    length(symbols(LR)) == 2 || throw(ArgumentError("The polynomials must be over a Laurent polynomial ring in two variables."))
+
+    R2 = Oscar._polyringquo(LR)
+    R = codomain(R2)
+    (x, y) = gens(LR)
+    I = ideal(LR, [f, g, x^a1[1] * y^a1[2] - 1, x^a2[1] * y^a2[2] - 1])
+    II = ideal(R, R2.(gens(I)))
+    Q, _ = quo(R, II)
+    k_dim = 2 * vector_space_dimension(Q)
+    n = 2 * length(monomial_basis(Q))
+
+    return FiniteMonomialCode(LR, base_ring(LR), n, k_dim, f, g, a1, a2)
+end
+
+#############################
+      # getter functions
+#############################
+
+polynomial_ring(S::AbstractBivariateBicycleCode) = S.R
+
+Laurent_polynomial_ring(S::Union{AbstractMonomialCode, AbstractGeneralizedToricCode}) = S.LR
+
+field(S::AbstractMonomialCode) = S.F
+
+"""
+    defining_polynomials(S::AbstractMonomialCode) -> Tuple{CTLRPolyElem, CTLRPolyElem}
+
+Return the polynomials defining the monomial code `S`.
+"""
+defining_polynomials(S::AbstractMonomialCode) = S.f, S.g
+
+"""
+    twist_vectors(S::AbstractMonomialCode) -> Tuple{Tuple{Int, Int}, Tuple{Int, Int}}
+
+Return the twist vectors of the monomial code `S` if they are defined.
+"""
+function twist_vectors(S::AbstractMonomialCode)
+    if hasproperty(S, :a1) && hasproperty(S, :a2)
+        return S.a1, S.a2
+    elseif hasproperty(S, :l) && hasproperty(S, :m)
+        return (S.l, 0), (0, S.m)
+    else
+        throw(ArgumentError("The twist vectors are not defined for this code.."))
+    end
+end
+
+length(S::AbstractMonomialCode) = S.n
+
+"""
+    dimension(S::AbstractMonomialCode) -> Int
+
+Return the dimension of the monomial code `S` if it is defined.
+"""
+function dimension(S::AbstractMonomialCode)
+    if hasproperty(S, :k)
+        return S.k
+    else
+        throw(ArgumentError("The dimension is not defined for this code. Use `maximum_dimension` instead."))
+    end
+end
+
+
+#############################
+      # setter functions
+#############################
+
+#############################
+     # general functions
+#############################
+
+"""
+    maximum_dimension(S::Union{MonomialCode, GeneralizedToricCode}) -> Int
+
+Return the maximum dimension of the monomial code `S`.
+"""
+function maximum_dimension(S::Union{MonomialCode, GeneralizedToricCode})
+    R2 = Oscar._polyringquo(S.LR)
+    R = codomain(R2)
+    (x, y) = gens(S.LR)
+    I = ideal(S.LR, [S.f, S.g])
+    II = ideal(R, R2.(gens(I)))
+    Q, ϕ = quo(R, II)
+    return 2 * vector_space_dimension(Q)
+end
+
+"""
+    CSSCode(S::BivariateBicycleCode)
+
+Return the CSS code defined by the bivariate bicycle code `S`.
+"""
+function CSSCode(S::BivariateBicycleCode)
+    x = matrix(S.F, [mod1(i + 1, S.l) == j ? 1 : 0 for i in 1:S.l, j in 1:S.l]) ⊗ identity_matrix(S.F, S.m)
+    y = identity_matrix(S.F, S.l) ⊗ matrix(S.F, [mod1(i + 1, S.m) == j ? 1 : 0 for i in 1:S.m, j in 1:S.m])
+
+    A = zero_matrix(S.F, S.l * S.m, S.l * S.m)
+    for ex in exponents(lift(S.a))
+        # iszero(ex[1]) || iszero(ex[2]) || throw(ArgumentError("Polynomial `a` must not have any `xy` terms"))
+        power, which = findmax(ex)
+        if which == 1
+            A += x^power
+        elseif which == 2
+            A += y^power
+        end
+    end
+
+    B = zero_matrix(S.F, S.l * S.m, S.l * S.m)
+    for ex in exponents(lift(S.b))
+        # iszero(ex[1]) || iszero(ex[2]) || throw(ArgumentError("Polynomial `b` must not have any `xy` terms"))
+        power, which = findmax(ex)
+        if which == 1
+            B += x^power
+        elseif which == 2
+            B += y^power
+        end
+    end
+
+
+    return CSSCode(hcat(A, B), hcat(transpose(B), transpose(A)))
+end
+
+"""
+    CSSCode(S::CoprimeBivariateBicycleCode)
+
+Return the coprime bivariate bicycle code defined by the residue ring elements `a` and `b`.
+
+# Note
+
+- This is defined in https://arxiv.org/pdf/2408.10001.
+"""
+function CSSCode(S::CoprimeBivariateBicycleCode)
+    deg_P = degree(modulus(S.R))
+    x = matrix(S.F, [mod1(i + 1, S.l) == j ? 1 : 0 for i in 1:S.l, j in 1:S.l]) ⊗ identity_matrix(S.F, S.m)
+    y = identity_matrix(S.F, S.l) ⊗ matrix(S.F, [mod1(i + 1, S.m) == j ? 1 : 0 for i in 1:S.m, j in 1:S.m])
+
+    P = x * y
+    A = zero_matrix(S.F, deg_P, deg_P)
+    exps = findall(i -> !is_zero(i), collect(coefficients(lift(S.a)))) .- 1
+    for ex in exps
+        A += P^ex
+    end
+
+    B = zero_matrix(S.F, deg_P, deg_P)
+    exps = findall(i -> !is_zero(i), collect(coefficients(lift(S.b)))) .- 1
+    for ex in exps
+        B += P^ex
+    end
+
+    return CSSCode(hcat(A, B), hcat(transpose(B), transpose(A)))
+end
+
+"""
+    CSSCode(S::FiniteGeneralizedToricCode)
+
+Return the CSS code defined by the finite generalized toric code `S`.
+"""
+function CSSCode(S::FiniteGeneralizedToricCode)
+    R2 = Oscar._polyringquo(S.LR)
+    R = codomain(R2)
+    (x, y) = gens(S.LR)
+    swap = hom(S.LR, S.LR, [x^-1, y^-1])
+    f_anti = swap(S.f)
+    g_anti = swap(S.g)
+    f_R2 = R2(S.f)
+    g_R2 = R2(S.g)
+    f_anti_R2 = R2(f_anti)
+    g_anti_R2 = R2(g_anti)
+    I = ideal(S.LR, [x^S.a1[1] * y^S.a1[2] - 1, x^S.a2[1] * y^S.a2[2] - 1])
+    II = ideal(R, R2.(gens(I)))
+    Q, ϕ = quo(R, II)
+
+    mono = monomial_basis(Q)
+    len_mon = length(mono)
+    n = 2 * len_mon
+    LR_edge_index = Dict{fpMPolyRingElem, Int}(mono[i] => i for i in 1:len_mon)
+    TB_edge_index = Dict{fpMPolyRingElem, Int}(mono[i] => i + len_mon for i in 1:len_mon)
+
+    Fone = S.F(1)
+    row = 1
+    X_stabs = zero_matrix(S.F, len_mon, n)
+    Z_stabs = zero_matrix(S.F, len_mon, n)
+    for edge in mono
+        f_shift = simplify(ϕ(edge * f_R2))
+        for term in terms(f_shift.f)
+            # X_12
+            X_stabs[row, TB_edge_index[term]] = Fone
+        end
+        g_shift = simplify(ϕ(edge * g_R2))
+        for term in terms(g_shift.f)
+            # X_14
+            X_stabs[row, LR_edge_index[term]] = Fone
+        end
+
+        g_shift = simplify(ϕ(edge * g_anti_R2))
+        for term in terms(g_shift.f)
+            # Z_12
+            Z_stabs[row, TB_edge_index[term]] = Fone
+        end
+        f_shift = simplify(ϕ(edge * f_anti_R2))
+        for term in terms(f_shift.f)
+            # Z_14
+            Z_stabs[row, LR_edge_index[term]] = Fone
+        end
+        row += 1
+    end
+
+    return CSSCode(X_stabs, Z_stabs)
+end
+
+# TODO add show methods for other codes
+function show(io::IO, S::AbstractGeneralizedToricCode)
+    if isa(S, FiniteGeneralizedToricCode)
+        println(io, "Finite Generalized Toric Code:")
+        println(io, "\tf: $(S.f)")
+        println(io, "\tg: $(S.g)")
+        println(io, "\ta1: $(S.a1)")
+        println(io, "\ta2: $(S.a2)")
+    else
+        println(io, "Generalized Toric Code:")
+        println(io, "\tf: $(S.f)")
+        println(io, "\tg: $(S.g)")
+    end
+end
+
+"""
+    is_pure(S::AbstractMonomialCode) -> Bool
+
+Return `true` if the logicals of the monomial code `S` are pure.
+"""
+function is_pure(S::AbstractMonomialCode)
+    # TODO return true if l and m are odd?
+    I_f = ideal(S.R, [S.f])
+    I_g = ideal(S.R, [S.g])
+    return I_f ∩ I_g == I_f * I_g
+end
+
+# """
+#     is_principal(S::AbstractMonomialCode) -> Bool
+
+# Return `true` if the logicals of the monomial code `S` are principal.
+# """
+# function is_principal(S::AbstractMonomialCode)
+#     TODO return true if l and m are odd?
+#     TODO: check if pure
+#     I_f = ideal(S.R, [S.f])
+#     I_g = ideal(S.R, [S.g])
+#     # TODO bug here with types and annihilator
+#     return is_principal(annihilator(I_f)) && is_principal(annihilator(I_g))
+# end
+
+#############################
+         # 3D Trial
+#############################
+
+# this one is really the same as above
+function Generalized3DToricCode(f::CTLRPolyElem, g::CTLRPolyElem)
+    LR = parent(f)
+    LR == parent(g) || throw(ArgumentError("The polynomials must be over the same ring."))
+    length(symbols(LR)) == 3 || throw(ArgumentError("The polynomials must be over a Laurent polynomial ring in three variables."))
+    return Generalized3DToricCode(LR, base_ring(LR), f, g)
+end
+
+function FiniteGeneralized3DToricCode(f::CTLRPolyElem, g::CTLRPolyElem, a1::Tuple{Int, Int},
+    a2::Tuple{Int, Int}, l::Int)
+
+    l > 0 || throw(DomainError(l, "Parameter l must be positive"))
+    LR = parent(f)
+    LR == parent(g) || throw(ArgumentError("The polynomials must be over the same ring."))
+    length(symbols(LR)) == 3 || throw(ArgumentError("The polynomials must be over a Laurent polynomial ring in three variables."))
+    return FiniteGeneralized3DToricCode(LR, base_ring(LR), f, g, a1, a2, l)
+end
+
+function maximum_dimension(S::AbstractGeneralized3DToricCode)
+    R2 = Oscar._polyringquo(S.LR)
+    R = codomain(R2)
+    (x, y, z) = gens(S.LR)
+    if isa(S, FiniteGeneralized3DToricCode)
+        I = ideal(S.LR, [S.f, S.g, x^S.a1[1] * y^S.a1[2] - 1, x^S.a2[1] * y^S.a2[2] - 1, z^S.l - 1])
+    else
+        I = ideal(S.LR, [S.f, S.g])
+    end
+    II = ideal(R, R2.(gens(I)))
+    Q, ϕ = quo(R, II)
+    return 2 * vector_space_dimension(Q)
+end
+
+function CSSCode(S::FiniteGeneralized3DToricCode)
+    R2 = Oscar._polyringquo(S.LR)
+    R = codomain(R2)
+    (x, y, z) = gens(S.LR)
+    swap = hom(S.LR, S.LR, [x^-1, y^-1, z^-1])
+    f_anti = swap(S.f)
+    g_anti = swap(S.g)
+    f_R2 = R2(S.f)
+    g_R2 = R2(S.g)
+    f_anti_R2 = R2(f_anti)
+    g_anti_R2 = R2(g_anti)
+    I = ideal(S.LR, [x^S.a1[1] * y^S.a1[2] - 1, x^S.a2[1] * y^S.a2[2] - 1, z^S.l - 1])
+    II = ideal(R, R2.(gens(I)))
+    Q, ϕ = quo(R, II)
+
+    mono = monomial_basis(Q)
+    len_mon = length(mono)
+    n = 2 * len_mon
+    LR_edge_index = Dict{fpMPolyRingElem, Int}(mono[i] => i for i in 1:len_mon)
+    TB_edge_index = Dict{fpMPolyRingElem, Int}(mono[i] => i + len_mon for i in 1:len_mon)
+
+    Fone = S.F(1)
+    row = 1
+    X_stabs = zero_matrix(S.F, len_mon, n)
+    Z_stabs = zero_matrix(S.F, len_mon, n)
+    for edge in mono
+        f_shift = simplify(ϕ(edge * f_R2))
+        for term in terms(f_shift.f)
+            # X_12
+            X_stabs[row, TB_edge_index[term]] = Fone
+        end
+        g_shift = simplify(ϕ(edge * g_R2))
+        for term in terms(g_shift.f)
+            # X_14
+            X_stabs[row, LR_edge_index[term]] = Fone
+        end
+
+        g_shift = simplify(ϕ(edge * g_anti_R2))
+        for term in terms(g_shift.f)
+            # Z_12
+            Z_stabs[row, TB_edge_index[term]] = Fone
+        end
+        f_shift = simplify(ϕ(edge * f_anti_R2))
+        for term in terms(f_shift.f)
+            # Z_14
+            Z_stabs[row, LR_edge_index[term]] = Fone
+        end
+        row += 1
+    end
+
+    return CSSCode(X_stabs, Z_stabs)
+end
\ No newline at end of file
diff --git a/src/Quantum/product_codes.jl b/src/Quantum/product_codes.jl
index bda0bd0..40a0faf 100644
--- a/src/Quantum/product_codes.jl
+++ b/src/Quantum/product_codes.jl
@@ -1079,138 +1079,3 @@ function random_homological_product_code(n1::Int, k1::Int, n2::Int, k2::Int)
     d = d1 ⊗ i2 - i1 ⊗ d2
     return CSSCode(d, transpose(d))
 end
-
-"""
-    BivariateBicycleCode(a::MPolyQuoRingElem{FqMPolyRingElem}, b::MPolyQuoRingElem{FqMPolyRingElem})
-
-Return the bivariate bicycle code defined by the residue ring elements `a` and `b`.
-
-# Note
-- This is defined in https://arxiv.org/pdf/2308.07915
-
-# Example
-
-[[360, 12, ≤24]] Bivariate Bicycle Code from Table 3 of [bravyi2024high](@cite).
-
-```jldoctest
-julia> using CodingTheory, Oscar;
-
-julia> S, (x, y) = polynomial_ring(Oscar.Nemo.Native.GF(2), [:x, :y]);
-
-julia> l = 30; m = 6;
-
-julia> R, _ = quo(S, ideal(S, [x^l - 1, y^m - 1]));
-
-julia> a = R(x^9 + y + y^2);
-
-julia> b = R(y^3 + x^25 + x^26);
-
-julia> code = BivariateBicycleCode(a, b);
-
-julia> length(code), dimension(code)
-(360, 12)
-```
-"""
-function BivariateBicycleCode(a::T, b::T) where T <: Union{MPolyQuoRingElem{FqMPolyRingElem}, MPolyQuoRingElem{fpMPolyRingElem}}
-    R = parent(a)
-    R == parent(b) || throw(DomainError("Polynomials must have the same parent."))
-    F = base_ring(base_ring(a))
-    order(F) == 2 || throw(DomainError("This code family is currently only defined over binary fields."))
-    length(symbols(parent(a))) == 2 || throw(DomainError("Polynomials must be over two variables."))
-    g = gens(modulus(R))
-    length(g) == 2 || throw(DomainError("Residue rings must have only two generators."))
-
-    m = -1
-    l = -1
-    for g1 in g
-        exps = collect(exponents(g1))
-        length(exps) == 2 || throw(ArgumentError("Moduli of the incorrect form."))
-        iszero(exps[2]) || throw(ArgumentError("Moduli of the incorrect form."))
-        !iszero(exps[1][1]) && !iszero(exps[1][2]) && throw(ArgumentError("Moduli of the incorrect form."))
-        if iszero(exps[1][1])
-            m = exps[1][2]
-        else
-            l = exps[1][1]
-        end
-    end
-
-    x = matrix(F, [mod1(i + 1, l) == j ? 1 : 0 for i in 1:l, j in 1:l]) ⊗ identity_matrix(F, m)
-    y = identity_matrix(F, l) ⊗ matrix(F, [mod1(i + 1, m) == j ? 1 : 0 for i in 1:m, j in 1:m])
-
-    A = zero_matrix(F, l * m, l * m)
-    for ex in exponents(lift(a))
-        iszero(ex[1]) || iszero(ex[2]) || throw(ArgumentError("Polynomial `a` must not have any `xy` terms"))
-        power, which = findmax(ex)
-        if which == 1
-            A += x^power
-        elseif which == 2
-            A += y^power
-        end
-    end
-
-    B = zero_matrix(F, l * m, l * m)
-    for ex in exponents(lift(b))
-        iszero(ex[1]) || iszero(ex[2]) || throw(ArgumentError("Polynomial `b` must not have any `xy` terms"))
-        power, which = findmax(ex)
-        if which == 1
-            B += x^power
-        elseif which == 2
-            B += y^power
-        end
-    end
-
-
-    return CSSCode(hcat(A, B), hcat(transpose(B), transpose(A)))
-end
-
-"""
-    CoprimeBivariateBicycleCode(a::MPolyQuoRingElem{FqMPolyRingElem}, b::MPolyQuoRingElem{FqMPolyRingElem})
-
-Return the coprime bivariate bicycle code defined by the residue ring elements `a` and `b`.
-
-# Note
-
-- This is defined in https://arxiv.org/pdf/2408.10001.
-"""
-function CoprimeBivariateBicycleCode(a::ResElem, b::ResElem)
-    R = parent(a)
-    S = base_ring(a)
-    R == parent(b) || throw(DomainError("Polynomials must have the same parent."))
-    F = base_ring(S)
-    order(F) == 2 || throw(DomainError("This code family is currently only defined over binary fields."))
-    length(gens(S)) == 1 || throw(DomainError("Polynomials must be over one variable."))
-    f = modulus(R)
-    deg_P = degree(f)
-    f == gen(S)^deg_P - 1 || throw(ArgumentError("Residue ring not of the form π^(l * m) - 1."))
-
-    facs = Nemo.factor(deg_P)
-    length(facs) == 2 || throw(ArgumentError("Residue ring not of the form π^(l * m) - 1."))
-    k = collect(keys(facs.fac))
-    v = collect(values(facs.fac))
-    l = k[1]^v[1]
-    m = k[2]^v[2]
-    # l = 3
-    # m = 7
-    # println("l = $l, m = $m")
-    # actually this check is guaranteed to be true since factor breaks it into its prime
-    # factorization
-    # gcd(l, m) == 1 || throw(ArgumentError("l and m must be coprime"))
-
-    x = matrix(F, [mod1(i + 1, l) == j ? 1 : 0 for i in 1:l, j in 1:l]) ⊗ identity_matrix(F, m)
-    y = identity_matrix(F, l) ⊗ matrix(F, [mod1(i + 1, m) == j ? 1 : 0 for i in 1:m, j in 1:m])
-
-    P = x * y
-    A = zero_matrix(F, deg_P, deg_P)
-    exps = findall(i -> !is_zero(i), collect(coefficients(lift(a)))) .- 1
-    for ex in exps
-        A += P^ex
-    end
-
-    B = zero_matrix(F, deg_P, deg_P)
-    exps = findall(i -> !is_zero(i), collect(coefficients(lift(b)))) .- 1
-    for ex in exps
-        B += P^ex
-    end
-
-    return CSSCode(hcat(A, B), hcat(transpose(B), transpose(A)))
-end
diff --git a/src/Quantum/subsystem_code.jl b/src/Quantum/subsystem_code.jl
index 57b68d3..380a67e 100644
--- a/src/Quantum/subsystem_code.jl
+++ b/src/Quantum/subsystem_code.jl
@@ -170,7 +170,7 @@ function SubsystemCode(G::CTMatrixTypes; char_vec::Union{Vector{zzModRingElem},
             u_bound_bare, _ = _min_wt_row(mat)
 
             # dressed
-            _, mat = _remove_empty(_rref_symp_col_swap(vcat(mat, gauge_ops_mat)), :rows)
+            mat = _remove_empty(_rref_symp_col_swap(vcat(mat, gauge_ops_mat)), :rows)
             u_bound_dressed, _ = _min_wt_row(mat)
             return GraphStateSubsystem(F, n, 0, r, missing, missing, 1, u_bound_bare, 1,
                 u_bound_dressed, stabs, char_vec, signs, missing, false, gauge_ops, gauge_ops_mat,
@@ -178,12 +178,12 @@ function SubsystemCode(G::CTMatrixTypes; char_vec::Union{Vector{zzModRingElem},
         else
             # bare
             # TODO use ! versions throughout after vcats
-            _, mat = _remove_empty(_rref_symp_col_swap(vcat(stabs, bare_logs)), :rows)
+            mat = _remove_empty(_rref_symp_col_swap(vcat(stabs, bare_logs)), :rows)
             anti = hcat(logs_mat[:, n + 1:end], -logs_mat[:, 1:n]) * transpose(_remove_empty(mat, :rows))
             u_bound_bare, _ = minimum(row_wts_symplectic(mat[findall(!iszero(anti[i:i, :]) for i in axes(anti, 1)), :]))
 
             # dressed
-            _, mat = _remove_empty(_rref_symp_col_swap(vcat(mat, gauge_ops_mat)), :rows)
+            mat = _remove_empty(_rref_symp_col_swap(vcat(mat, gauge_ops_mat)), :rows)
             anti = hcat(logs_mat[:, n + 1:end], -logs_mat[:, 1:n]) * transpose(_remove_empty(mat, :rows))
             u_bound_dressed, _ = minimum(row_wts_symplectic(mat[findall(!iszero(anti[i:i, :]) for i in axes(anti, 1)), :]))
             return SubsystemCode(F, n, k, r, missing, missing, 1, u_bound_bare, 1, u_bound_dressed, 
@@ -305,12 +305,12 @@ function SubsystemCode(S::CTMatrixTypes, L::CTMatrixTypes, G::CTMatrixTypes;
             stabs_stand, stand_r, stand_k, P_stand, missing, missing)
     else
         # bare
-        _, mat = _remove_empty(_rref_symp_col_swap(vcat(stabs, bare_logs)), :rows)
+        mat = _remove_empty(_rref_symp_col_swap(vcat(S, bare_logs)), :rows)
         anti = hcat(logs_mat[:, n + 1:end], -logs_mat[:, 1:n]) * transpose(_remove_empty(mat, :rows))
         u_bound_bare, _ = minimum(row_wts_symplectic(mat[findall(!iszero(anti[i:i, :]) for i in axes(anti, 1)), :]))
 
         # dressed
-        _, mat = _remove_empty(_rref_symp_col_swap(vcat(mat, gauge_ops_mat)), :rows)
+        mat = _remove_empty(_rref_symp_col_swap(vcat(mat, gauge_ops_mat)), :rows)
         anti = hcat(logs_mat[:, n + 1:end], -logs_mat[:, 1:n]) * transpose(_remove_empty(mat, :rows))
         u_bound_dressed, _ = minimum(row_wts_symplectic(mat[findall(!iszero(anti[i:i, :]) for i in axes(anti, 1)), :]))
         return SubsystemCode(F, n, k, r, missing, missing, 1, u_bound_bare, 1, u_bound_dressed, S,
@@ -1521,7 +1521,7 @@ function _make_pairs(L::T) where T <: CTMatrixTypes
                     end
                 end
             end
-            iszero(first) && error("Cannot make symplectic basis. Often this is due to the fact that the stabilizers are not maximal and therefore the centralizer still containing part of the isotropic subspace.")
+            iszero(first) && error("Cannot make symplectic basis. Often this is due to the fact that the stabilizers are not maximal and therefore the centralizer still contains part of the isotropic subspace.")
             for c in 2:num_prod
                 if !iszero(prod[first, c])
                     L[c:c, :] += F(prod[first, c]^-1) * L[1:1, :]
@@ -1540,7 +1540,7 @@ function _make_pairs(L::T) where T <: CTMatrixTypes
             # the columns in prod give the commutation relationships between the provided
             # logical operators; they ideally should only consist of {X_1, Z_i} pairs
             # so there should only be one nonzero element in each column
-            prod = hcat(L[:, n + 1:end], -L[:, 1:n]) * transpose(L)
+            prod = hcat(L[:, n + 1:end], L[:, 1:n]) * transpose(L)
             # println("before")
             # display(prod)
             num_prod = ncols(prod)
@@ -1554,13 +1554,13 @@ function _make_pairs(L::T) where T <: CTMatrixTypes
                     end
                 end
             end
-            iszero(first) && error("Cannot make symplectic basis. Often this is due to the fact that the stabilizers are not maximal and therefore the centralizer still containing part of the isotropic subspace.")
+            iszero(first) && error("Cannot make symplectic basis. Often this is due to the fact that the stabilizers are not maximal and therefore the centralizer still contains part of the isotropic subspace.")
             for c in 2:num_prod
                 if !iszero(prod[first, c])
                     L[c:c, :] += L[1:1, :]
                 end
             end
-            prod = hcat(L[:, n + 1:end], -L[:, 1:n]) * transpose(L)
+            prod = hcat(L[:, n + 1:end], L[:, 1:n]) * transpose(L)
             # println("after")
             # display(prod)
             push!(logs, (L[1:1, :], L[first:first, :]))
diff --git a/src/Quantum/types.jl b/src/Quantum/types.jl
index 12dca76..9adbb9a 100644
--- a/src/Quantum/types.jl
+++ b/src/Quantum/types.jl
@@ -21,6 +21,10 @@ abstract type AbstractEASubsystemCode <: AbstractSubsystemCode end
 abstract type AbstractEASubsystemCodeCSS <: AbstractEASubsystemCode end
 abstract type AbstractEAStabilizerCode <: AbstractStabilizerCode end
 abstract type AbstractEAStabilizerCodeCSS <: AbstractEAStabilizerCode end
+abstract type AbstractMonomialCode <: AbstractStabilizerCodeCSS end
+abstract type AbstractBivariateBicycleCode <: AbstractMonomialCode end
+abstract type AbstractGeneralizedToricCode <: AbstractMonomialCode end
+abstract type AbstractGeneralized3DToricCode <: AbstractMonomialCode end
 
 # AbstractQuantumLDPCCode, AbstractQuantumLDPCCSSCode?
 
@@ -336,6 +340,88 @@ mutable struct HypergraphProductCode <: AbstractHypergraphProductCode
       Z_metacheck::Union{CTMatrixTypes, Missing}
 end
 
+#############################
+# Quantum/GeneralizedToricCode.jl
+#############################
+
+
+struct MonomialCode <: AbstractMonomialCode
+      LR::AbstractAlgebra.Generic.LaurentMPolyWrapRing{fpFieldElem, fpMPolyRing}
+      F::CTFieldTypes
+      n::Int
+      f::CTLRPolyElem
+      g::CTLRPolyElem
+end
+
+struct FiniteMonomialCode <: AbstractMonomialCode
+      LR::AbstractAlgebra.Generic.LaurentMPolyWrapRing{fpFieldElem, fpMPolyRing}
+      F::CTFieldTypes
+      n::Int
+      k::Int
+      f::CTLRPolyElem
+      g::CTLRPolyElem
+      a1::Tuple{Int, Int}
+      a2::Tuple{Int, Int}
+end
+
+struct BivariateBicycleCode <: AbstractBivariateBicycleCode
+      R::CTPolyRing
+      F::CTFieldTypes
+      n::Int
+      k::Int
+      f::CTPolyRingElem
+      g::CTPolyRingElem
+      l::Int
+      m::Int
+end
+
+struct CoprimeBivariateBicycleCode <: AbstractBivariateBicycleCode
+      R::CTPolyRing
+      F::CTFieldTypes
+      n::Int
+      k::Int
+      f::CTPolyRingElem
+      g::CTPolyRingElem
+      l::Int
+      m::Int
+end
+
+struct GeneralizedToricCode <: AbstractGeneralizedToricCode
+      LR::AbstractAlgebra.Generic.LaurentMPolyWrapRing{fpFieldElem, fpMPolyRing}
+      F::CTFieldTypes
+      n::Int
+      f::CTLRPolyElem
+      g::CTLRPolyElem
+end
+
+struct FiniteGeneralizedToricCode <: AbstractGeneralizedToricCode
+      LR::AbstractAlgebra.Generic.LaurentMPolyWrapRing{fpFieldElem, fpMPolyRing}
+      F::CTFieldTypes
+      n::Int
+      k::Int
+      f::CTLRPolyElem
+      g::CTLRPolyElem
+      a1::Tuple{Int, Int}
+      a2::Tuple{Int, Int}
+end
+
+struct Generalized3DToricCode <: AbstractGeneralized3DToricCode
+      LR::AbstractAlgebra.Generic.LaurentMPolyWrapRing{fpFieldElem, fpMPolyRing}
+      F::CTFieldTypes
+      f::CTLRPolyElem
+      g::CTLRPolyElem
+end
+
+struct FiniteGeneralized3DToricCode <: AbstractGeneralized3DToricCode
+      LR::AbstractAlgebra.Generic.LaurentMPolyWrapRing{fpFieldElem, fpMPolyRing}
+      F::CTFieldTypes
+      f::CTLRPolyElem
+      g::CTLRPolyElem
+      a1::Tuple{Int, Int}
+      a2::Tuple{Int, Int}
+      l::Int
+end
+
 #############################
           # traits
 #############################
diff --git a/src/automorphisms.jl b/src/automorphisms.jl
new file mode 100644
index 0000000..a66efff
--- /dev/null
+++ b/src/automorphisms.jl
@@ -0,0 +1,36 @@
+# Copyright (c) 2025 Eric Sabo
+# All rights reserved.
+#
+# This source code is licensed under the BSD-style license found in the
+# LICENSE file in the root directory of this source tree.
+
+#############################
+       # Tanner Graphs
+#############################
+
+function Tanner_graph_automorphism_group(H::CTMatrixTypes)
+    nr, nc = size(H)
+    A = vcat(hcat(zero_matrix(base_ring(H), nc, nc), transpose(H)), hcat(H, zero_matrix(base_ring(H), nr, nr)))
+    G = graph_from_adjacency_matrix(Undirected, A)
+    Aut_gens = automorphism_group_generators(G)
+
+    # with the convention above, 1:nc are the var nodes, nc + 1:end are the check nodes
+
+    return Aut_gens
+end
+
+Tanner_graph_automorphism_group(L::AbstractLDPCCode) = Tanner_graph_automorphism_group(L.H)
+
+# I think this needs to be different
+Tanner_graph_automorphism_group(S::AbstractStabilizerCode) = Tanner_graph_automorphism_group(S.stabs)
+
+# these need to pull out the zero's in the symplectic format
+Tanner_graph_automorphism_group_X(S::AbstractStabilizerCSSCode) = Tanner_graph_automorphism_group(S.X_stabs)
+
+Tanner_graph_automorphism_group_Z(S::AbstractStabilizerCSSCode) = Tanner_graph_automorphism_group(S.Z_stabs)
+
+
+
+_foiliated_stabilizer_matrix(H::CodingTheory.CTMatrixTypes, d::Int) = 
+    hcat(identity_matrix(base_ring(H), d) ⊗ H, CodingTheory._rep_pcm_tr(base_ring(H), d) ⊗
+    identity_matrix(base_ring(H), Oscar.nrows(H)))
diff --git a/src/utils.jl b/src/utils.jl
index 2b58ff1..3860a4c 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -315,7 +315,7 @@ function _rref_non_pivot_cols(A::CTMatrixTypes, type::Symbol = :nsp)
     end
 end
 
-function _quotient_space(big::T, small::T, alg::Symbol=:sys_eqs) where T <: CTMatrixTypes
+function _quotient_space(big::T, small::T, alg::Symbol = :sys_eqs) where T <: CTMatrixTypes
     alg ∈ [:VS, :sys_eqs] || throw(ArgumentError("Unknown algorithm type"))
 
     if alg == :VS
@@ -750,7 +750,7 @@ function _rref_col_swap_perm!(A::CTMatrixTypes, row_range::AbstractUnitRange{Int
                     for l in i:nr
                         if !iszero(A[l, k])
                             swap_cols!(A, k, j)
-                            P = P * cperm(sym_group, [k,j])
+                            P = P * cperm(sym_group, [k, j])
                             ind = l
                             break
                         end
@@ -802,7 +802,7 @@ function _rref_col_swap_perm!(A::CTMatrixTypes, row_range::AbstractUnitRange{Int
                     for l in i:nr
                         if !iszero(A[l, k])
                             swap_cols!(A, k, j)
-                            P = P * cperm(sym_group, [k,j])
+                            P = P * cperm(sym_group, [k, j])
                             ind = l
                             break
                         end
@@ -2001,6 +2001,63 @@ function _value_distribution(vals)
     return OrderedDict([(i, count(x -> (x == i), vals)) for i in collect(sort(unique(vals)))])
 end
 
+"""
+    prepare_mult_order_mod(n::Int) -> Tuple{Int,Int,Vector{Tuple{Int,Int}}}
+
+Precompute data for repeated multiplicative-order queries modulo n.
+Returns:
+- prep[1] = n
+- prep[2] = phi = φ(n)
+- prep[3] = pf = prime factorization of φ(n) as (p, e) pairs
+"""
+function prepare_mult_order(n::Int)::Tuple{Int,Int,Vector{Tuple{Int,Int}}}
+    2 <= n || throw(ArgumentError("n must satisfy 2 <= n, got n=$n"))
+    if n > 10_000
+        println("Warning: modulus for multiplicative order is larger than 10000")
+    end
+
+    #could optimize more here if needed. `euler_phi` already has to factor n so this approach is 
+    #doing two calls to `factor` when we could get away with one call.
+    #We would only need to write a variant of euler_phi that takes a `factorization` object as input.
+    #I havent bothered because I expect for the range of n of interest to us the divisor 
+    #search below will be the computational bottleneck rather than the calls to `factor`
+    phi = Int(euler_phi(Nemo.ZZ(n)))
+    pf = Tuple{Int,Int}[(Int(p), Int(e)) for (p, e) in factor(Nemo.ZZ(phi))]
+    return (n, phi, pf)
+end
+
+"""
+    mult_order(a::Int, prep::Tuple{Int,Int,Vector{Tuple{Int,Int}}}) -> Int
+
+Compute ord_n(a) using precomputed tuple from `prepare_mult_order_mod(n)`.
+"""
+function mult_order(a::Int, prep::Tuple{Int,Int,Vector{Tuple{Int,Int}}})::Int
+    n, phi, pf = prep
+    # Note on the validation line below: 
+    # To be pedantic, multiplicative order is only defined when gcd(a,n)=1. 
+    # but its possible we'll run into an application where we're looking at an x (mod n) that's not coprime.
+    # If n is squarefree the sequence x,x^2,... wont include 1 but it is stil periodic. 
+    # We may want to know that period.
+    # Its probably better to write a different function for that
+    gcd(a, n) == 1 || throw(ArgumentError("multiplicative order is only defined when gcd(a,n)=1"))
+
+    a = mod(a, n)
+    a == 1 && return 1
+
+    m = phi
+    for (p, e) in pf
+        for _ in 1:e
+            t = m ÷ p
+            if powermod(a, t, n) == 1
+                m = t
+            else
+                break
+            end
+        end
+    end
+    return m
+end
+
 # #=
 # Example of using the repeated iterator inside of product.
 #
diff --git a/test/Quantum/GeneralizedToricCode_test.jl b/test/Quantum/GeneralizedToricCode_test.jl
new file mode 100644
index 0000000..011fb0a
--- /dev/null
+++ b/test/Quantum/GeneralizedToricCode_test.jl
@@ -0,0 +1,202 @@
+@testitem "Quantum/GeneralizedToricCode.jl" begin
+    using Oscar
+    using CodingTheory
+
+    @testset "Generalized Toric Codes" begin
+        # Table 1 of https://arxiv.org/abs/2503.03827
+        F = Oscar.Nemo.Native.GF(2);
+        R, (x, y) = laurent_polynomial_ring(F, [:x, :y]);
+        f = 1 + x + x * y;
+        g = 1 + y + x * y;
+        a1 = (0, 3);
+        a2 = (2, 1);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 12
+        @test dimension(S) == 4
+        
+        f = 1 + x + y;
+        g = 1 + y + x;
+        a1 = (0, 7);
+        a2 = (1, 2);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 14
+        @test dimension(S) == 6
+
+        f = 1 + x + x * y;
+        g = 1 + y + x * y;
+        a1 = (0, 3);
+        a2 = (3, 0);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 18
+        @test dimension(S) == 4
+                   
+        f = 1 + x + x * y;
+        g = 1 + y + x * y;
+        a1 = (0, 3);
+        a2 = (4, 2);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 24
+        @test dimension(S) == 4
+                   
+        f = 1 + x + x^-1 * y;
+        g = 1 + y + x * y;
+        a1 = (0, 7);
+        a2 = (2, 3);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 28
+        @test dimension(S) == 6
+       
+        f = 1 + x + x^2;
+        g = 1 + y + x^2;
+        a1 = (0, 3);
+        a2 = (5, 1);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 30
+        @test dimension(S) == 4     
+        
+        f = 1 + x + x^-1; 
+        g = 1 + y + y^-1;
+        a1 = (0, 9);
+        a2 = (2, 4);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 36
+        @test dimension(S) == 4 
+                
+        f = 1 + x + x * y;
+        g = 1 + y + x * y^-1;
+        a1 = (0, 7);
+        a2 = (3, 2);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 42
+        @test dimension(S) == 6
+        
+        f = 1 + x + x^2;
+        g = 1 + y + x^2;
+        a1 = (0, 3);
+        a2 = (8, 1);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 48
+        @test dimension(S) == 4
+                   
+        f = 1 + x + x^-1;
+        g = 1 + y + x^3 * y^2;
+        a1 = (0, 3);
+        a2 = (9, 0);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 54
+        @test dimension(S) == 8
+                   
+        f = 1 + x + y^-2;
+        g = 1 + y + x^-2;
+        a1 = (0, 7);
+        a2 = (4, 3);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 56
+        @test dimension(S) == 6
+                   
+        f = 1 + x + y^-2;
+        g = 1 + y + x^2;
+        a1 = (0, 10);
+        a2 = (3, 3);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 60
+        @test dimension(S) == 8
+                   
+        f = 1 + x + x^-1 * y;
+        g = 1 + y + x^-1 * y^-1;
+        a1 = (0, 31);
+        a2 = (1, 13);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 62
+        @test dimension(S) == 10
+                   
+        f = 1 + x + x^-2 * y^-1;
+        g = 1 + y + x^2 * y;
+        a1 = (0, 3);
+        a2 = (11, 2);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 66
+        @test dimension(S) == 4
+                   
+        f = 1 + x + x * y;
+        g = 1 + y + x * y^-1;
+        a1 = (0, 7);
+        a2 = (5, 1);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 70
+        @test dimension(S) == 6
+                   
+        f = 1 + x + x^-1 * y^3;
+        g = 1 + y + x^3 * y^-1;
+        a1 = (0, 12);
+        a2 = (3, 3);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 72
+        @test dimension(S) == 8
+        
+        f = 1 + x + x^-2 * y^-1;
+        g = 1 + y + x^2 * y;
+        a1 = (0, 3);
+        a2 = (13, 1);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 78
+        @test dimension(S) == 4
+                   
+        f = 1 + x + x^-2;
+        g = 1 + y + x^-2 * y^2;
+        a1 = (0, 14);
+        a2 = (3, -6);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 84
+        @test dimension(S) == 6
+                   
+        f = 1 + x + x^-1 * y^-3;
+        g = 1 + y + x^3 * y^-1;
+        a1 = (0, 15);
+        a2 = (3, -6);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 90
+        @test dimension(S) == 8
+                   
+        f = 1 + x + x^-2 * y;
+        g = 1 + y + x * y^-2;
+        a1 = (0, 12);
+        a2 = (4, 2);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 96
+        @test dimension(S) == 4
+
+        f = 1 + x + x^-1 * y^2;
+        g = 1 + y + x^-2 * y^-1;
+        a1 = (0, 7);
+        a2 = (7, 0);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 98
+        @test dimension(S) == 6
+                   
+        f = 1 + x + x^-3 * y;
+        g = 1 + y + x^3 * y^2;
+        a1 = (0, 3);
+        a2 = (17, 2);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 102
+        @test dimension(S) == 4
+                   
+        f = 1 + x + x^-1 * y^-3;
+        g = 1 + y + x^3 * y^-1;
+        a1 = (0, 9);
+        a2 = (6, 0);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 108
+        @test dimension(S) == 8
+                   
+        f = 1 + x + x^-1 * y^3;
+        g = 1 + y + x^3 * y^-1;
+        a1 = (0, 9);
+        a2 = (6, 0);
+        S = CSSCode(FiniteGeneralizedToricCode(f, g, a1, a2))
+        @test length(S) == 108
+        @test dimension(S) == 8
+           
+    end
+end
diff --git a/test/Quantum/product_codes_test.jl b/test/Quantum/product_codes_test.jl
index 1ed9d91..d92a616 100644
--- a/test/Quantum/product_codes_test.jl
+++ b/test/Quantum/product_codes_test.jl
@@ -606,6 +606,96 @@
         @test length(Q) == 196
         @test dimension(Q) == 12
         @test_broken minimum_distance(S) == 8
+
+        # Table 1 of https://arxiv.org/pdf/2407.03973v1
+        # l and m in table are swapped compared to the other papers?
+        # [[108, 16, 6]]
+        l = 6
+        m = 9
+        R, _ = quo(S, ideal(S, [x^l - 1, y^m - 1]))
+        a = R(1 + y + y^3)
+        b = R(y^3 + x^2 + x^4)
+        Q = BivariateBicycleCode(a, b)
+        @test length(Q) == 108
+        @test dimension(Q) == 16
+        @test_broken minimum_distance(S) == 6
+
+        # [[128, 14, 12]]
+        l = 8
+        m = 8
+        R, _ = quo(S, ideal(S, [x^l - 1, y^m - 1]))
+        a = R(x^2 + y + y^3 + y^4)
+        b = R(y^2 + x + x^3 + x^4)
+        Q = BivariateBicycleCode(a, b)
+        @test length(Q) == 128
+        @test dimension(Q) == 14
+        @test_broken minimum_distance(S) == 12
+
+        # [[162, 4, 16]]
+        l = 9
+        m = 9
+        R, _ = quo(S, ideal(S, [x^l - 1, y^m - 1]))
+        a = R(1 + x + y)
+        b = R(x^3 + y + y^2)
+        Q = BivariateBicycleCode(a, b)
+        @test length(Q) == 162
+        @test dimension(Q) == 4
+        @test_broken minimum_distance(S) == 16
+
+        # [[162, 12, 8]]
+        l = 9
+        m = 9
+        R, _ = quo(S, ideal(S, [x^l - 1, y^m - 1]))
+        a = R(1 + x + y^6)
+        b = R(y^3 + x^2 + x^3)
+        Q = BivariateBicycleCode(a, b)
+        @test length(Q) == 162
+        @test dimension(Q) == 12
+        @test_broken minimum_distance(S) == 8
+
+        # [[162, 24, 6]]
+        l = 9
+        m = 9
+        R, _ = quo(S, ideal(S, [x^l - 1, y^m - 1]))
+        a = R(1 + y + y^2)
+        b = R(y^3 + x^3 + x^6)
+        Q = BivariateBicycleCode(a, b)
+        @test length(Q) == 162
+        @test dimension(Q) == 24
+        @test_broken minimum_distance(S) == 6
+
+        # [[270, 8, 18]]
+        l = 9
+        m = 15
+        R, _ = quo(S, ideal(S, [x^l - 1, y^m - 1]))
+        a = R(x^3 + y + y^2)
+        b = R(y^3 + x + x^2)
+        Q = BivariateBicycleCode(a, b)
+        @test length(Q) == 270
+        @test dimension(Q) == 8
+        @test_broken minimum_distance(S) == 18
+
+        # [[98, 6, 12]]
+        l = 7
+        m = 7
+        R, _ = quo(S, ideal(S, [x^l - 1, y^m - 1]))
+        a = R(x + y^3 + y^4)
+        b = R(y + x^3 + x^4)
+        Q = BivariateBicycleCode(a, b)
+        @test length(Q) == 98
+        @test dimension(Q) == 6
+        @test_broken minimum_distance(S) == 12
+
+        # [[162, 8, 12]]
+        l = 9
+        m = 9
+        R, _ = quo(S, ideal(S, [x^l - 1, y^m - 1]))
+        a = R(x^3 + y + y^2)
+        b = R(y + x^3 + x^4)
+        Q = BivariateBicycleCode(a, b)
+        @test length(Q) == 162
+        @test dimension(Q) == 8
+        @test_broken minimum_distance(S) == 12
     end
 
     @testset "CoprimeBivariateBicycleCode" begin