fix: close the quote on the version string

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

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

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,

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`.

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

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)))

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)