diff --git a/CHANGELOG.md b/CHANGELOG.md
index fc74781..d1efef8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,21 @@ All notable changes to QuestBase.jl will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [Unreleased]
+
+### Fixed
+
+- `simplify_exp_products` left `exp(a)*exp(a)` alone. A product stores a repeated factor as a power, so the term arrives as the single factor `exp(a)^2` and the `isexp` test skipped it.
+- `collapse_pythagorean` took the first squared trig factor it found in a summand, so a summand carrying two of them, such as `cos(3t)^2*sin(t)^2*a`, failed to pair up with its partner. Every way of splitting a summand is now a candidate.
+- `declare_variable` bound its result inside `QuestBase` with `@eval`. Nothing read that binding, and evaluating into a closed module prevents a downstream package from precompiling a workload that declares variables.
+
+### Performance
+
+- `fraction_free_linear_solve` splits the system into the independent subsystems its sparsity pattern allows and eliminates each on its own. The coefficient matrix of an n-harmonic ansatz is n uncoupled 2x2 blocks; solved whole, every solution came out over the determinant of the full matrix. Over six representative systems `rearrange_standard` drops from 134.2 MiB to 10.1 MiB and from 0.114 s to 0.028 s.
+- Sums and products are rebuilt with SymbolicUtils' n-ary `add_worker`/`mul_worker` instead of being folded with `sum`/`prod`/`+=`, which is quadratic in the number of terms.
+- `trig_reduce` no longer calls `expand_all`. Its extra `Postwalk(expand_exp_power)` was the most expensive step of the averaging, and the following `simplify_exp_products` already normalises `exp(a)^n` at every node it reaches.
+- `is_rearranged` decides whether a derivative appears on the left-hand side by walking the expression tree rather than rendering both sides to strings.
+
 ## [0.5.0]
 
 ### Breaking
diff --git a/src/HarmonicEquation.jl b/src/HarmonicEquation.jl
index 2900425..65be1fd 100644
--- a/src/HarmonicEquation.jl
+++ b/src/HarmonicEquation.jl
@@ -236,8 +236,8 @@ function is_rearranged(eom::HarmonicEquation)
 
     HB_bool = isequal(rhs, dvar)
     hopf_bool = in("Hopf", getfield.(eom.variables, :type))
-    MF_bool =
-        !any([occursin(str1, str2) for str1 in string.(dvar) for str2 in string.(lhs)])
+    # no derivative anywhere on the left-hand side, asked of the tree rather than `string(lhs)`
+    MF_bool = !any(occurs_in(dv, l) for dv in dvar, l in lhs)
 
     # Hopf-containing equations or MF equation are arranged by construstion
     return HB_bool || hopf_bool || MF_bool
diff --git a/src/QuestBase.jl b/src/QuestBase.jl
index 57a6300..bb832e7 100644
--- a/src/QuestBase.jl
+++ b/src/QuestBase.jl
@@ -19,7 +19,10 @@ using SymbolicUtils:
     add_with_div,
     is_literal_number,
     unwrap_const,
-    unwrap
+    unwrap,
+    vartype,
+    add_worker,
+    mul_worker
 
 using Symbolics:
     Symbolics,
diff --git a/src/Symbolics/Symbolics_utils.jl b/src/Symbolics/Symbolics_utils.jl
index 0fb77bd..83e2b96 100644
--- a/src/Symbolics/Symbolics_utils.jl
+++ b/src/Symbolics/Symbolics_utils.jl
@@ -16,7 +16,10 @@ function expand_fraction(x::BasicSymbolic)
         args = arguments(x)
         num = SymbolicUtils.expand(args[1])
         if isadd(num)
-            sum(expand_fraction(arg / args[2]) for arg in arguments(num))
+            add_worker(
+                vartype(typeof(x)),
+                map(arg -> expand_fraction(arg / args[2]), arguments(num)),
+            )
         else
             x
         end
@@ -26,12 +29,19 @@ function expand_fraction(x::BasicSymbolic)
 end
 expand_fraction(x::Num) = Num(expand_fraction(unwrap(x)))
 
-"Apply a function f on every member of a sum or a product"
+"The vartype every expression here is built on, for callers with no `BasicSymbolic` to read it off."
+const NUM_VARTYPE = vartype(fieldtype(Num, 1))
+
+"""
+Apply a function f on every member of a sum or a product.
+
+Rebuilt with `add_worker`/`mul_worker`; folding with `sum` is `O(n²)`.
+"""
 function _apply_termwise(f, x::BasicSymbolic)
     if isadd(x)
-        sum(f(arg) for arg in arguments(x))
+        add_worker(vartype(typeof(x)), map(f, arguments(x)))
     elseif ismul(x)
-        prod(f(arg) for arg in arguments(x))
+        mul_worker(vartype(typeof(x)), map(f, arguments(x)))
     elseif isdiv(x)
         args = arguments(x)
         _apply_termwise(f, args[1]) / _apply_termwise(f, args[2])
@@ -107,9 +117,9 @@ get_independent(x, t::Num) = x
 
 function get_independent(x::BasicSymbolic, t::Num)
     if isadd(x)
-        sum(get_independent(arg, t) for arg in arguments(x))
+        add_worker(vartype(typeof(x)), map(arg -> get_independent(arg, t), arguments(x)))
     elseif ismul(x)
-        prod(get_independent(arg, t) for arg in arguments(x))
+        mul_worker(vartype(typeof(x)), map(arg -> get_independent(arg, t), arguments(x)))
     elseif isdiv(x)
         args = arguments(x)
         !is_function(args[2], t) ? get_independent(args[1], t) / args[2] : 0
@@ -131,7 +141,14 @@ function get_all_terms(x::Equation)
 end
 function _get_all_terms(x::BasicSymbolic)
     if isadd(x)
-        vcat([_get_all_terms(arg) for arg in arguments(x)]...)
+        # not `vcat(parts...)`: splatting a runtime-length collection dispatches dynamically
+        terms = Any[]
+        for arg in arguments(x)
+            part = _get_all_terms(arg)
+            # the fallback method returns a bare non-symbolic, not a list of one
+            part isa AbstractVector ? append!(terms, part) : push!(terms, part)
+        end
+        terms
     elseif ismul(x)
         arguments(x)
     elseif isdiv(x)
@@ -166,6 +183,22 @@ is_function(f, var) = unwrap(var) in get_variables(f)
 """
 $(TYPEDSIGNATURES)
 
+Return true if `needle` occurs as a subexpression of `haystack`.
+
+Structural, not `occursin(string(needle), string(haystack))`: pretty-printing costs a full
+render and normalises, reorders and elides.
+"""
+occurs_in(needle, haystack) = _occurs_in(unwrap(needle), unwrap(haystack))
+
+function _occurs_in(target, expr)
+    isequal(expr, target) && return true
+    expr isa BasicSymbolic && SymbolicUtils.iscall(expr) || return false
+    return any(arg -> _occurs_in(target, arg), arguments(expr))
+end
+
+"""
+$(TYPEDSIGNATURES)
+
 Return the subexpressions of `x` in which a variable of `vars` (or a derivative of one)
 appears inside an operation other than an addition, a multiplication or a power with a
 non-negative integer exponent. An empty result means `x` is a polynomial in `vars`.
diff --git a/src/Symbolics/exponentials.jl b/src/Symbolics/exponentials.jl
index b680de7..cca1c27 100644
--- a/src/Symbolics/exponentials.jl
+++ b/src/Symbolics/exponentials.jl
@@ -1,12 +1,28 @@
 "Returns true if expr is an exponential"
 isexp(expr) = isterm(expr) && operation(expr) === exp
 
+"""
+The argument `a` of an exponential factor, reading `exp(a)` and `exp(a)^n` alike (the latter
+as `a*n`); `nothing` otherwise.
+
+A `Mul` stores repeated factors as a power, so `exp(a)*exp(a)` is held as `exp(a)^2` and a
+plain [`isexp`](@ref) test misses it.
+"""
+function _exp_argument(expr)
+    isexp(expr) && return arguments(expr)[1]
+    if expr isa BasicSymbolic && ispow(expr)
+        base, exponent = arguments(expr)
+        isexp(base) && return arguments(base)[1] * exponent
+    end
+    return nothing
+end
+
 "Expand powers of exponential such that exp(x)^n => exp(x*n)"
 function expand_exp_power(expr::BasicSymbolic)
     if isadd(expr)
-        sum(expand_exp_power(arg) for arg in arguments(expr))
+        add_worker(vartype(typeof(expr)), map(expand_exp_power, arguments(expr)))
     elseif ismul(expr)
-        prod(expand_exp_power(arg) for arg in arguments(expr))
+        mul_worker(vartype(typeof(expr)), map(expand_exp_power, arguments(expr)))
     elseif ispow(expr) && isexp(arguments(expr)[1])
         exp(arguments(arguments(expr)[1])[1] * arguments(expr)[2])
     else
@@ -26,7 +42,7 @@ function simplify_exp_products(expr::BasicSymbolic)
     elseif ismul(expr)
         _simplify_exp_products_mul(expr)
     elseif ispow(expr) && isexp(arguments(expr)[1])
-        # exp(x)^n => exp(x*n) — fixes bug where exp powers were left unexpanded
+        # exp(x)^n => exp(x*n), which is otherwise left unexpanded here
         exp(arguments(arguments(expr)[1])[1] * arguments(expr)[2])
     else
         expr
@@ -35,10 +51,12 @@ end
 
 function _simplify_exp_products_mul(expr)
     args = arguments(expr)
-    ind = findall(isexp, args)
+    exp_args = map(_exp_argument, args)
+    ind = findall(!isnothing, exp_args)
     rest_ind = setdiff(1:length(args), ind)
-    rest = isempty(rest_ind) ? 1 : prod(args[rest_ind])
-    total = isempty(ind) ? 0 : sum(arguments(args[i])[1] for i in ind)
+    vtype = vartype(typeof(expr))
+    rest = isempty(rest_ind) ? 1 : mul_worker(vtype, args[rest_ind])
+    total = isempty(ind) ? 0 : add_worker(vtype, exp_args[ind])
     if is_literal_number(total)
         iszero(unwrap_const(total)) && return rest
     end
diff --git a/src/Symbolics/fourier.jl b/src/Symbolics/fourier.jl
index 60579e9..b398db7 100644
--- a/src/Symbolics/fourier.jl
+++ b/src/Symbolics/fourier.jl
@@ -25,7 +25,10 @@ function trig_reduce(x)
     x = add_div(x) # a/b + c/d = (ad + bc)/bd
     x = expand(x) # open all brackets
     x = trig_to_exp(x)
-    x = expand_all(x) # expand products of exponentials
+    # Not `expand_all`: its extra `Postwalk(expand_exp_power)` is the most expensive step of
+    # the averaging, and `simplify_exp_products` below already normalises `exp(a)^n` at every
+    # node it reaches.
+    x = expand(x) # expand products of exponentials
     x = simplify_exp_products(x) # simplify products of exps
     x = exp_to_trig(x)
     x = Num(simplify_complex(expand(x)))
diff --git a/src/Symbolics/linear_solve.jl b/src/Symbolics/linear_solve.jl
index f4c3cf5..3021abf 100644
--- a/src/Symbolics/linear_solve.jl
+++ b/src/Symbolics/linear_solve.jl
@@ -83,6 +83,10 @@ The coefficients are converted together into SymbolicUtils' sparse polynomial
 representation. Bareiss elimination then performs only exact polynomial divisions
 and produces one determinant denominator per solution, instead of the nested
 fractions produced by symbolic LU.
+
+The system is first split into independent subsystems (see [`_system_blocks`](@ref)) and each
+is eliminated on its own: for a harmonic ansatz that is `n` separate `2 × 2` solves instead
+of one `2n × 2n`.
 """
 function fraction_free_linear_solve(equations, variables)
     coefficients, offsets, islinear = Symbolics.linear_expansion(equations, variables)
@@ -92,7 +96,78 @@ function fraction_free_linear_solve(equations, variables)
     rows == columns == length(offsets) ||
         throw(DimensionMismatch("fraction-free linear solve requires a square system"))
 
-    return _bareiss_jordan_solve(coefficients, -offsets)
+    right_hand_side = -offsets
+    blocks = _system_blocks(coefficients)
+    isnothing(blocks) && return _bareiss_jordan_solve(coefficients, right_hand_side)
+
+    solution = Vector{Num}(undef, rows)
+    for (block_rows, block_columns) in blocks
+        solution[block_columns] = _bareiss_jordan_solve(
+            coefficients[block_rows, block_columns], right_hand_side[block_rows]
+        )
+    end
+    return solution
+end
+
+"Is `entry` zero on its face, so the row and column it joins are unrelated?"
+function _is_structural_zero(entry)
+    value = unwrap_const(unwrap(entry))
+    return value isa Number && iszero(value)
+end
+
+"""
+Split a system into the independent subsystems its sparsity pattern allows.
+
+The coefficient matrix of an `n`-harmonic ansatz is `n` uncoupled `2×2` blocks. Solving it
+whole is worse than wasteful: a Bareiss entry carries the determinant of everything
+eliminated so far, so the solution comes out over the determinant of the full matrix rather
+than the one block determinant that survives.
+
+Blocks are the connected components of the bipartite graph joining row `r` to column `c`
+where the coefficient is not [`_is_structural_zero`](@ref). Returns `(rows, columns)` index
+pairs, or `nothing` when there is one component or an unbalanced one (structurally
+singular, left to the full solve to reject).
+"""
+function _system_blocks(coefficients)
+    n = size(coefficients, 1)
+    # union-find over rows `1:n` and columns `n .+ (1:n)`
+    parent = collect(1:(2n))
+    function representative(index)
+        root = index
+        while parent[root] != root
+            root = parent[root]
+        end
+        while parent[index] != root  # path compression
+            parent[index], index = root, parent[index]
+        end
+        return root
+    end
+    for column in 1:n, row in 1:n
+        _is_structural_zero(coefficients[row, column]) && continue
+        row_root, column_root = representative(row), representative(n + column)
+        row_root == column_root || (parent[row_root] = column_root)
+    end
+
+    block_of_root = Dict{Int,Int}()
+    block_rows = Vector{Int}[]
+    for row in 1:n
+        block = get!(block_of_root, representative(row)) do
+            push!(block_rows, Int[])
+            return length(block_rows)
+        end
+        push!(block_rows[block], row)
+    end
+    length(block_rows) > 1 || return nothing
+
+    block_columns = [Int[] for _ in block_rows]
+    for column in 1:n
+        block = get(block_of_root, representative(n + column), 0)
+        iszero(block) && return nothing  # a column no row reaches: structurally singular
+        push!(block_columns[block], column)
+    end
+    all(length.(block_rows) .== length.(block_columns)) || return nothing
+
+    return collect(zip(block_rows, block_columns))
 end
 
 function _bareiss_jordan_solve(coefficients, right_hand_side)
@@ -217,19 +292,31 @@ _from_polynomial(value::Number, _) = Num(value)
 function _from_polynomial(variable::typeof(_BAREISS_VARIABLE), poly_to_symbolic)
     return Num(poly_to_symbolic[variable])
 end
+"""
+Rebuild a symbolic expression from a Bareiss polynomial.
+
+A solved entry carries thousands of monomials, so the sum and the per-monomial product go to
+`add_worker`/`mul_worker` whole rather than through `+=`/`*=`. See [`_apply_termwise`](@ref).
+"""
 function _from_polynomial(polynomial, poly_to_symbolic)
     variables = MP.variables(polynomial)
-    result = Num(0)
+    addends = Any[]
+    factors = Any[]
     for term in MP.terms(polynomial)
-        expression = Num(MP.coefficient(term))
+        empty!(factors)
+        push!(factors, MP.coefficient(term))
         exponents = MP.exponents(MP.monomial(term))
         for (variable, exponent) in zip(variables, exponents)
             iszero(exponent) && continue
-            expression *= Num(poly_to_symbolic[variable])^exponent
+            push!(factors, unwrap(poly_to_symbolic[variable])^exponent)
         end
-        result += expression
+        push!(
+            addends,
+            length(factors) == 1 ? only(factors) : mul_worker(NUM_VARTYPE, copy(factors)),
+        )
     end
-    return result
+    isempty(addends) && return Num(0)
+    return Num(length(addends) == 1 ? only(addends) : add_worker(NUM_VARTYPE, addends))
 end
 
 function _find_bareiss_pivot(augmented, pivot_index, n)
diff --git a/src/Symbolics/pythagorean.jl b/src/Symbolics/pythagorean.jl
index a9095a6..6088996 100644
--- a/src/Symbolics/pythagorean.jl
+++ b/src/Symbolics/pythagorean.jl
@@ -34,37 +34,43 @@ function _is_squared(x)
     return value isa Number && value == 2
 end
 
+const _TrigSplit = Tuple{Any,Any,Num}
+
 """
-Split a summand into `(cos-or-sin, argument, remaining factors)` when it carries a squared
-trig factor, `nothing` otherwise.
+Every way to read a summand as `(cos-or-sin, argument, remaining factors)` by pulling out one
+squared trig factor; empty when it carries none.
+
+A summand can hold more than one squared trig factor, as `cos(3θ)^2*sin(θ)^2*a` does, and
+only one of its splits pairs up with a partner elsewhere in the sum. Committing to the first
+factor found leaves `cos(3θ)^2*sin(θ)^2*a + sin(3θ)^2*sin(θ)^2*a` unpaired, so every split
+has to stay a candidate.
 """
-function _squared_trig_term(summand)
+function _squared_trig_terms(summand)
     if summand isa BasicSymbolic && ispow(summand)
         base, exponent = arguments(summand)
-        _is_squared(exponent) || return nothing
+        _is_squared(exponent) || return _TrigSplit[]
         trig = _trig_operation(base)
-        isnothing(trig) && return nothing
-        return trig[1], trig[2], Num(1)
+        isnothing(trig) && return _TrigSplit[]
+        return _TrigSplit[(trig[1], trig[2], Num(1))]
     elseif summand isa BasicSymbolic && ismul(summand)
-        trig = nothing
-        rest = Num(1)
-        for factor in arguments(summand)
-            if isnothing(trig) && factor isa BasicSymbolic && ispow(factor)
-                base, exponent = arguments(factor)
-                if _is_squared(exponent)
-                    candidate = _trig_operation(base)
-                    if !isnothing(candidate)
-                        trig = candidate
-                        continue
-                    end
-                end
+        factors = collect(arguments(summand))
+        splits = _TrigSplit[]
+        for (i, factor) in enumerate(factors)
+            (factor isa BasicSymbolic && ispow(factor)) || continue
+            base, exponent = arguments(factor)
+            _is_squared(exponent) || continue
+            trig = _trig_operation(base)
+            isnothing(trig) && continue
+
+            rest = Num(1)
+            for (j, other) in enumerate(factors)
+                j == i || (rest *= wrap(other))
             end
-            rest *= wrap(factor)
+            push!(splits, (trig[1], trig[2], rest))
         end
-        isnothing(trig) && return nothing
-        return trig[1], trig[2], rest
+        return splits
     end
-    return nothing
+    return _TrigSplit[]
 end
 
 "Pair up squared sines and cosines among the summands of a single `Add` node."
@@ -72,36 +78,40 @@ function _collapse_pythagorean_node(x)
     (x isa BasicSymbolic && isadd(x)) || return x
     summands = collect(arguments(x))
     length(summands) > 1 || return x
-    parts = map(_squared_trig_term, summands)
-    any(!isnothing, parts) || return x
+    parts = map(_squared_trig_terms, summands)
+    any(!isempty, parts) || return x
 
     taken = falses(length(summands))
-    result = Num(0)
+    # collected, not accumulated with `+=`: see `_apply_termwise`
+    kept = Any[]
     collapsed_any = false
     for i in eachindex(summands)
         taken[i] && continue
-        part = parts[i]
-        if !isnothing(part)
-            operation_i, argument_i, coefficient_i = part
+        matched = false
+        for (operation_i, argument_i, coefficient_i) in parts[i]
             partner_operation = operation_i === cos ? sin : cos
+            matches_partner(split) =
+                split[1] === partner_operation &&
+                isequal(split[2], argument_i) &&
+                isequal(split[3], coefficient_i)
             partner = findfirst(eachindex(summands)) do j
                 j == i && return false
                 taken[j] && return false
-                isnothing(parts[j]) && return false
-                return parts[j][1] === partner_operation &&
-                       isequal(parts[j][2], argument_i) &&
-                       isequal(parts[j][3], coefficient_i)
+                return any(matches_partner, parts[j])
             end
             if !isnothing(partner)
                 taken[i] = taken[partner] = true
-                result += coefficient_i
+                push!(kept, unwrap(coefficient_i))
                 collapsed_any = true
-                continue
+                matched = true
+                break
             end
         end
+        matched && continue
         taken[i] = true
-        result += wrap(summands[i])
+        push!(kept, summands[i])
     end
     collapsed_any || return x
-    return unwrap(result)
+    isempty(kept) && return unwrap(Num(0))
+    return length(kept) == 1 ? unwrap(Num(only(kept))) : add_worker(NUM_VARTYPE, kept)
 end
diff --git a/src/Variables.jl b/src/Variables.jl
index 5ddf1b3..c4f9484 100644
--- a/src/Variables.jl
+++ b/src/Variables.jl
@@ -4,22 +4,37 @@ function d(f::Num, x::Num, deg=1)::Num
 end
 d(funcs::Vector{Num}, x::Num, deg=1) = Num[d(f, x, deg) for f in funcs]
 
-"Declare a variable in the the currently active Module namespace"
-function declare_variable(name::String)
-    var_sym = Symbol(name)
-    @eval($(var_sym) = first(Symbolics.@variables $var_sym))
-    return eval(var_sym)
-end
+"""
+The symtype `Symbolics.@variables f(t)` gives its variable. Spelled out in full because
+`Symbolics.variable` needs a concrete type and `FnType{Tuple,Real}` is still a `UnionAll`.
+"""
+const _FnTypeReal = SymbolicUtils.FnType{Tuple,Real,Nothing}
+
+"""
+$(TYPEDSIGNATURES)
+
+Declare a symbolic variable named `name`.
+
+`Symbolics.@variables` needs the name as a literal, so names only known at runtime (`u1`,
+`v2`, the bracket-free copies [`_remove_brackets`](@ref) makes) have to go through
+`Symbolics.variable` instead.
+
+This used to build the variable with `@eval` and bind the result inside this module. Nothing
+ever read that binding, since it lands in `QuestBase` rather than in the caller's namespace,
+and evaluating into a closed module makes the package impossible to precompile with a
+workload: `@compile_workload` runs the very code that creates these variables.
+"""
+declare_variable(name::String) = Symbolics.variable(Symbol(name))
 
 declare_variable(x::Num) = declare_variable(string(x))
 
-"Declare a variable that is a function of another variable in the Module namespace"
+"""
+$(TYPEDSIGNATURES)
+
+Declare a variable named `name` that is a function of `independent_variable`.
+"""
 function declare_variable(name::String, independent_variable::Num)
-    # independent_variable = declare_variable(independent_variable) convert string into Num
-    var_sym = Symbol(name)
-    new_var = Symbolics.@variables $var_sym(independent_variable)
-    @eval($(var_sym) = first($new_var)) # store the variable under "name" in this namespace
-    return eval(var_sym)
+    return Symbolics.variable(Symbol(name); T=_FnTypeReal)(independent_variable)
 end
 
 "Return the name of a variable (excluding independent variables)"
diff --git a/test/symbolics.jl b/test/symbolics.jl
index b12cb4a..152757d 100644
--- a/test/symbolics.jl
+++ b/test/symbolics.jl
@@ -23,6 +23,15 @@ end
     @eqtest simplify_exp_products(exp(a) * exp(b)) == exp(a + b)
     @eqtest simplify_exp_products(exp(3a) * exp(4b)) == exp(3a + 4b)
     @eqtest simplify_exp_products(im * exp(3a) * exp(4b)) == im * exp(3a + 4b)
+
+    # a product stores a repeated factor as a power, so `exp(a)*exp(a)` never reaches the
+    # merge as two separate `exp` factors. It arrives as the single factor `exp(a)^2`.
+    @eqtest simplify_exp_products(exp(a) * exp(a)) == exp(2a)
+    @eqtest simplify_exp_products(b * exp(a) * exp(a)) == b * exp(2a)
+    @eqtest simplify_exp_products(exp(a)^2 * exp(b)) == exp(2a + b)
+    # and the exponents must cancel to nothing, as they do for two plain factors
+    @eqtest simplify_exp_products(exp(a)^2 * exp(-2a) * b) == b
+    @eqtest simplify_exp_products(exp(a)^2 * exp(-a)^2) == 1
 end
 
 @testset "euler" begin
@@ -249,6 +258,76 @@ end
     end
 end
 
+@testset "fraction-free linear solve splits independent blocks" begin
+    using LinearAlgebra
+    using QuestBase: fraction_free_linear_solve, _system_blocks
+
+    @variables a b c x y z w p q
+
+    # Two 2x2 blocks, interleaved so the decomposition cannot rely on contiguous indices.
+    # Each solution carries only its own block's determinant: solved whole, every one of
+    # them would come out over the product of both.
+    equations = [p * x + q * z ~ a, x - z ~ b, p * y + w ~ c, y - w ~ 0]
+    solution = fraction_free_linear_solve(equations, [x, y, z, w])
+    @eqtest solution[1] == (a + b * q) / (p + q)
+    @eqtest solution[3] == (a - b * p) / (p + q)
+    @eqtest solution[2] == c / (1 + p)
+    @eqtest solution[4] == c / (1 + p)
+
+    # Transcendental coefficients do not change the decomposition, and the trig-free block
+    # is not dragged through the rotation block's determinant.
+    trig = fraction_free_linear_solve(
+        [cos(p) * x - sin(p) * z ~ a, sin(p) * x + cos(p) * z ~ b, y ~ c], [x, y, z]
+    )
+    @eqtest trig[2] == c
+    # Put the rotation on a rational point of the unit circle so the residual is exact.
+    numeric = Dict(cos(p) => 3 // 5, sin(p) => 4 // 5, a => 1, b => 2)
+    @eqtest substitute(cos(p) * trig[1] - sin(p) * trig[3], numeric) == 1
+    @eqtest substitute(sin(p) * trig[1] + cos(p) * trig[3], numeric) == 2
+
+    # Whole matrix nonzero: one component, so there is nothing to split and the solve runs
+    # as before.
+    @test isnothing(_system_blocks(Num[1 1; 1 1]))
+
+    # A structurally singular pattern has an unbalanced component, and is left to the full
+    # solve to reject rather than being silently split.
+    @test isnothing(_system_blocks(Num[1 1; 0 0]))
+    @test_throws LinearAlgebra.SingularException fraction_free_linear_solve(
+        [x + y ~ a, 0 * x + 0 * y ~ b], [x, y]
+    )
+
+    # Random block-diagonal integer systems against exact rational arithmetic.
+    rng = Random.MersenneTwister(0xb10c)
+    for _ in 1:8
+        sizes = rand(rng, 1:3, 3)
+        dimension = sum(sizes)
+        matrix = zeros(Int, dimension, dimension)
+        offset = 0
+        for size in sizes
+            block = rand(rng, -4:4, size, size)
+            while iszero(det(block))
+                block = rand(rng, -4:4, size, size)
+            end
+            matrix[(offset + 1):(offset + size), (offset + 1):(offset + size)] = block
+            offset += size
+        end
+        permutation = Random.randperm(rng, dimension)
+        matrix = matrix[permutation, :]
+        rhs = rand(rng, -4:4, dimension)
+        symbolic_variables = only(@variables r[1:dimension])
+        random_equations = [
+            sum(
+                matrix[row, column] * symbolic_variables[column] for column in 1:dimension
+            ) ~ rhs[row] for row in 1:dimension
+        ]
+        actual = fraction_free_linear_solve(random_equations, symbolic_variables)
+        expected = Rational{Int}.(matrix) \ Rational{Int}.(rhs)
+        for index in eachindex(actual)
+            @eqtest actual[index] == expected[index]
+        end
+    end
+end
+
 @testset "rearrange! reduces determinant denominators" begin
     # Fraction-free elimination produces one determinant denominator instead of nested
     # LU fractions. For a trigonometric ansatz that determinant contains