fix: handle bug in remove_denominators

Author: AayushSabharwal

ext/CatalystHomotopyContinuationExtension.jl

@@ -6,7 +6,7 @@ import DynamicPolynomials
 import ModelingToolkit as MT
 import HomotopyContinuation as HC
 import Setfield: @set
-import Symbolics: unwrap, wrap, Rewriters, symtype, issym
+import Symbolics: unwrap, wrap, Rewriters, symtype, issym, maketerm, BasicSymbolic, metadata
 using Symbolics: iscall
 
 # Creates and exports hc_steady_states function.

ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl

@@ -59,7 +59,7 @@ function steady_state_polynomial(rs::ReactionSystem, ps, u0)
     eqs_pars_funcs = vcat(equations(ns), conservedequations(rs))
     eqs = map(eq -> substitute(eq.rhs - eq.lhs, p_dict), eqs_pars_funcs)
     eqs_intexp = make_int_exps.(eqs)
-    ss_poly = Catalyst.to_multivariate_poly(remove_denominators.(eqs_intexp))
+    ss_poly = Catalyst.to_multivariate_poly(remove_denominators.(common_denominator.(eqs_intexp)))
     return poly_type_convert(ss_poly)
 end
 
@@ -78,6 +78,45 @@ function ___make_int_exps(expr)
     end
 end
 
+function common_denominator(expr)
+    iscall(expr) || return expr
+    if operation(expr) == /
+        num, den = arguments(expr)
+        num = common_denominator(num)
+        den = common_denominator(den)
+        return num / den
+    end
+    if operation(expr) == +
+        num = 0
+        den = 1
+        for arg in arguments(expr)
+            arg = common_denominator(arg)
+            if iscall(arg) && operation(arg) == /
+                n, d = arguments(arg)
+            else
+                n = arg
+                d = 1
+            end
+            num = num * d + den * n
+            den *= d
+        end
+        return num / den
+    end
+    if operation(expr) == ^
+        base, pow = arguments(expr)
+        base = common_denominator(base)
+        if iscall(base) && operation(base) == /
+            num, den = arguments(base)
+        else
+            num, den = base, 1
+        end
+        num ^= pow
+        den ^= pow
+        return num / den
+    end
+    return maketerm(BasicSymbolic, operation(expr), map(common_denominator, arguments(expr)), metadata(expr))
+end
+
 # If the input is a fraction, removes the denominator.
 function remove_denominators(expr)
     s_expr = simplify_fractions(expr)