fix

Project.toml

@@ -1,22 +1,26 @@
 name = "QuestBase"
 uuid = "7e80f742-43d6-403d-a9ea-981410111d43"
 authors = ["Orjan Ameye <orjan.ameye@hotmail.com>", "Jan Kosata <kosataj@phys.ethz.ch>", "Javier del Pino <jdelpino@phys.ethz.ch>"]
-version = "0.4.0"
+version = "0.4.1"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
 DocStringExtensions = "0.9.4"
+DynamicPolynomials = "0.6"
 SymbolicUtils = "4"
 Symbolics = "7"
 julia = "1.12"
 Random = "1.10"
 LinearAlgebra = "1.10"
+MultivariatePolynomials = "0.5"
 Test = "1.10"
 OrderedCollections = "1.8, 2"
 Aqua = "0.8.11"

src/DifferentialEquation.jl

@@ -182,9 +182,16 @@ the left-hand sides. Uses symbolic linear solving to determine the right-hand si
 the equations in place.
 """
 function rearrange!(eom::DifferentialEquation, new_lhs::Vector{Num})
-    soln = Symbolics.symbolic_linear_solve(
-        get_equations(eom), new_lhs; simplify=false, check=true
-    )
+    soln = try
+        fraction_free_linear_solve(get_equations(eom), new_lhs)
+    catch error
+        error isa BareissFailure || rethrow()
+        Num.(
+            Symbolics.symbolic_linear_solve(
+                get_equations(eom), new_lhs; simplify=false, check=true
+            ),
+        )
+    end
     # Use original variable keys (not extracted from new_lhs, as v7's get_variables
     # treats derivatives as variables)
     eom.equations = OrderedDict(zip(collect(keys(eom.equations)), new_lhs .~ soln))

src/HarmonicEquation.jl

@@ -171,15 +171,16 @@ end
 
 "Rearrange an equation system such that the field equations is equal to the vector specified in new_lhs"
 function rearrange!(eom::HarmonicEquation, new_rhs::Vector{Num})
-    soln = Symbolics.symbolic_linear_solve(
-        eom.equations, new_rhs; simplify=false, check=true
-    )
-    # `symbolic_linear_solve` returns nested fractions whose denominators hold the
-    # coefficient determinant; for a trigonometric ansatz it contains identities like
-    # cos(ωt)² + sin(ωt)² = 1 that SymbolicUtils 4's `simplify` no longer collapses.
-    # Flatten the nest (polynomial-level, no AC matching) and reduce the denominator,
-    # so no consumer works with superficially time-dependent denominators.
-    soln = reduce_denominator.(Symbolics.simplify_fractions.(Num.(soln)))
+    soln = try
+        fraction_free_linear_solve(eom.equations, new_rhs)
+    catch error
+        error isa BareissFailure || rethrow()
+        fallback = Symbolics.symbolic_linear_solve(
+            eom.equations, new_rhs; simplify=false, check=true
+        )
+        Symbolics.simplify_fractions.(Num.(fallback))
+    end
+    soln = reduce_denominator.(soln)
     eom.equations = soln .~ new_rhs
     return nothing
 end

src/QuestBase.jl

@@ -1,6 +1,8 @@
 module QuestBase
 
 using DocStringExtensions
+import DynamicPolynomials as DP
+import MultivariatePolynomials as MP
 using OrderedCollections: OrderedCollections, OrderedDict, OrderedSet
 using LinearAlgebra: LinearAlgebra
 
@@ -40,6 +42,7 @@ include("Symbolics/Symbolics_utils.jl")
 include("Symbolics/exponentials.jl")
 include("Symbolics/fourier.jl")
 include("Symbolics/drop_powers.jl")
+include("Symbolics/linear_solve.jl")
 include("DifferentialEquation.jl")
 include("Variables.jl")
 include("HarmonicVariable.jl")

src/Symbolics/linear_solve.jl

@@ -0,0 +1,183 @@
+struct BareissFailure <: Exception
+    message::String
+end
+
+const _BAREISS_VARIABLE = only(DP.@polyvar __questbase_bareiss_variable)
+
+function Base.showerror(io::IO, error::BareissFailure)
+    return print(io, error.message)
+end
+
+"""
+    fraction_free_linear_solve(equations, variables)
+
+Solve a square symbolic linear system using pivoted Bareiss–Jordan elimination.
+
+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.
+"""
+function fraction_free_linear_solve(equations, variables)
+    coefficients, offsets, islinear = Symbolics.linear_expansion(equations, variables)
+    islinear || throw(ArgumentError("the equation system is not linear in the variables"))
+
+    rows, columns = size(coefficients)
+    rows == columns == length(offsets) ||
+        throw(DimensionMismatch("fraction-free linear solve requires a square system"))
+
+    return _bareiss_jordan_solve(coefficients, -offsets)
+end
+
+function _bareiss_jordan_solve(coefficients, right_hand_side)
+    n = size(coefficients, 1)
+    symbolic_type = eltype(coefficients)
+    poly_variable_type = typeof(_BAREISS_VARIABLE)
+    poly_to_symbolic = Dict{poly_variable_type,symbolic_type}()
+    symbolic_to_poly = Dict{symbolic_type,poly_variable_type}()
+    augmented = Matrix{Any}(undef, n, n + 1)
+
+    for column in 1:n, row in 1:n
+        augmented[row, column] = _to_polynomial(
+            coefficients[row, column], symbolic_to_poly, poly_to_symbolic
+        )
+    end
+    for row in 1:n
+        augmented[row, n + 1] = _to_polynomial(
+            right_hand_side[row], symbolic_to_poly, poly_to_symbolic
+        )
+    end
+
+    previous_pivot = 1
+    for pivot_index in 1:n
+        pivot_row = _find_bareiss_pivot(augmented, pivot_index, n)
+        isnothing(pivot_row) && throw(LinearAlgebra.SingularException(pivot_index))
+        if pivot_row != pivot_index
+            for column in axes(augmented, 2)
+                swapped = augmented[pivot_index, column]
+                augmented[pivot_index, column] = augmented[pivot_row, column]
+                augmented[pivot_row, column] = swapped
+            end
+        end
+
+        old = copy(augmented)
+        pivot = old[pivot_index, pivot_index]
+        for row in 1:n
+            row == pivot_index && continue
+            for column in axes(augmented, 2)
+                column == pivot_index && continue
+                numerator =
+                    pivot * old[row, column] -
+                    old[row, pivot_index] * old[pivot_index, column]
+                augmented[row, column] = if pivot_index == 1
+                    numerator
+                else
+                    _exact_polynomial_division(numerator, previous_pivot)
+                end
+            end
+            augmented[row, pivot_index] = 0
+        end
+        previous_pivot = pivot
+    end
+
+    solution = Vector{Num}(undef, n)
+    for row in 1:n
+        numerator_poly = augmented[row, n + 1]
+        denominator_poly = augmented[row, row]
+        leading_coefficient = if denominator_poly isa Number
+            denominator_poly
+        else
+            MP.leading_coefficient(denominator_poly)
+        end
+        if leading_coefficient isa Real && leading_coefficient < 0
+            numerator_poly = -numerator_poly
+            denominator_poly = -denominator_poly
+        end
+        numerator = _from_polynomial(numerator_poly, poly_to_symbolic)
+        denominator = _from_polynomial(denominator_poly, poly_to_symbolic)
+        solution[row] = numerator / denominator
+    end
+    return solution
+end
+
+function _to_polynomial(expression::Num, symbolic_to_poly, poly_to_symbolic)
+    return _to_polynomial(unwrap(expression), symbolic_to_poly, poly_to_symbolic)
+end
+_to_polynomial(expression::Number, _, _) = expression
+function _to_polynomial(expression::BasicSymbolic, symbolic_to_poly, poly_to_symbolic)
+    if SymbolicUtils.isconst(expression)
+        return unwrap_const(expression)
+    elseif isadd(expression)
+        result = 0
+        for argument in arguments(expression)
+            result += _to_polynomial(argument, symbolic_to_poly, poly_to_symbolic)
+        end
+        return result
+    elseif ismul(expression)
+        result = 1
+        for argument in arguments(expression)
+            result *= _to_polynomial(argument, symbolic_to_poly, poly_to_symbolic)
+        end
+        return result
+    elseif ispow(expression)
+        base, exponent = arguments(expression)
+        if SymbolicUtils.isconst(exponent)
+            value = unwrap_const(exponent)
+            if value isa Integer && value >= 0
+                return _to_polynomial(base, symbolic_to_poly, poly_to_symbolic)^value
+            end
+        end
+    end
+
+    variable = get!(symbolic_to_poly, expression) do
+        name = Symbol(:questbase_atom_, hash(expression))
+        return MP.similar_variable(_BAREISS_VARIABLE, name)
+    end
+    get!(poly_to_symbolic, variable, expression)
+    return variable
+end
+
+_from_polynomial(value::Number, _) = Num(value)
+function _from_polynomial(variable::typeof(_BAREISS_VARIABLE), poly_to_symbolic)
+    return Num(poly_to_symbolic[variable])
+end
+function _from_polynomial(polynomial, poly_to_symbolic)
+    variables = MP.variables(polynomial)
+    result = Num(0)
+    for term in MP.terms(polynomial)
+        expression = Num(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
+        end
+        result += expression
+    end
+    return result
+end
+
+function _find_bareiss_pivot(augmented, pivot_index, n)
+    best_row = nothing
+    best_size = typemax(Int)
+    for row in pivot_index:n
+        candidate = augmented[row, pivot_index]
+        iszero(candidate) && continue
+        candidate_size = candidate isa Number ? 1 : length(MP.terms(candidate))
+        if candidate_size < best_size
+            best_row = row
+            best_size = candidate_size
+        end
+    end
+    return best_row
+end
+
+function _exact_polynomial_division(numerator, denominator)
+    quotient = try
+        MP.div_multiple(numerator, denominator)
+    catch error
+        throw(BareissFailure("Bareiss exact division failed: $(sprint(showerror, error))"))
+    end
+    iszero(quotient * denominator - numerator) ||
+        throw(BareissFailure("Bareiss division produced a nonzero remainder"))
+    return quotient
+end

test/DifferentialEquations.jl

@@ -87,6 +87,34 @@ end
     @eqtest new_eq != diff_eq
 end
 
+@testset "State-dependent derivative mass matrix" begin
+    @variables α β μ γx γy kx ky Fx Fy
+    ẍ = d(x, t, 2)
+    ÿ = d(y, t, 2)
+    ẋ = d(x, t)
+    ẏ = d(y, t)
+
+    m11 = 1 + α * x^2
+    m12 = μ * x
+    m21 = μ * y
+    m22 = 1 + β * y^2
+    force_x = Fx * cos(ω * t) - γx * ẋ - kx * x
+    force_y = Fy * sin(ω * t) - γy * ẏ - ky * y
+    equations = [m11 * ẍ + m12 * ÿ ~ force_x, m21 * ẍ + m22 * ÿ ~ force_y]
+    system = DifferentialEquation(equations, [x, y])
+
+    rearrange_standard!(system)
+    rearranged = collect(values(system.equations))
+    determinant = m11 * m22 - m12 * m21
+
+    @eqtest rearranged[1].lhs == ẍ
+    @eqtest rearranged[2].lhs == ÿ
+    @eqtest rearranged[1].rhs ==
+        Symbolics.expand(m22 * force_x - m12 * force_y) / Symbolics.expand(determinant)
+    @eqtest rearranged[2].rhs ==
+        Symbolics.expand(m11 * force_y - m21 * force_x) / Symbolics.expand(determinant)
+end
+
 @testset "Error Cases" begin
     @testset "Equation{Vector}" begin
         # define equation of motion

test/symbolics.jl

@@ -173,12 +173,52 @@ end
     @eqtest reduce_denominator(a / (ω^2 + 1)) == a / (ω^2 + 1)
 end
 
+@testset "fraction-free linear solve" begin
+    using LinearAlgebra
+    using QuestBase: fraction_free_linear_solve
+
+    @variables a b c x y z
+    equations = [2x + y - z ~ a, x + 3y + z ~ b, 2x - y + 4z ~ c]
+    solution = fraction_free_linear_solve(equations, [x, y, z])
+    @eqtest solution[1] == (13a - 3b + 4c) / 31
+    @eqtest solution[2] == (-2a + 10b - 3c) / 31
+    @eqtest solution[3] == (-7a + 4b + 5c) / 31
+
+    # A zero diagonal requires a row pivot.
+    pivoted = fraction_free_linear_solve([y ~ a, x + y ~ b], [x, y])
+    @eqtest pivoted[1] == b - a
+    @eqtest pivoted[2] == a
+
+    @test_throws LinearAlgebra.SingularException fraction_free_linear_solve(
+        [x + y ~ a, 2x + 2y ~ b], [x, y]
+    )
+
+    # Compare arbitrary nonsingular integer systems against exact rational arithmetic.
+    rng = Random.MersenneTwister(0x5eed)
+    for dimension in 2:5
+        matrix = rand(rng, -4:4, dimension, dimension)
+        while iszero(det(matrix))
+            matrix = rand(rng, -4:4, dimension, dimension)
+        end
+        rhs = rand(rng, -4:4, dimension)
+        symbolic_variables = only(@variables q[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
-    # `symbolic_linear_solve` returns nested fractions whose denominators hold the
-    # coefficient determinant; for a trigonometric ansatz it contains cos² + sin² = 1.
-    # `rearrange!` must flatten and collapse it, otherwise every consumer works with
-    # superficially time-dependent denominators (this stalled the Krylov-Bogoliubov
-    # order-2 averaging for hours).
+    # Fraction-free elimination produces one determinant denominator instead of nested
+    # LU fractions. For a trigonometric ansatz that determinant contains
+    # cos² + sin² = 1, which `rearrange!` must collapse before downstream averaging.
     using QuestBase: HarmonicEquation, HarmonicVariable, DifferentialEquation
     using QuestBase: rearrange!, d
     @variables ω t T u1(T) v1(T) x(t)