Merge pull request #538 from JuliaHomotopyContinuation/feature/suppor…

…t-sqrt Support sqrt in modelkit
JuliaHomotopyContinuation · Jun 6, 2023 · 41150ee · 41150ee
2 parents 20c7c4d + 3fc1304
commit 41150ee
Show file tree

Hide file tree

Showing 4 changed files with 26 additions and 15 deletions.
diff --git a/src/DoubleDouble.jl b/src/DoubleDouble.jl
@@ -467,9 +467,10 @@ Base.abs(a::DoubleF64) = a.hi < 0.0 ? -a : a
 Base.eps(::DoubleF64) = double_eps
 Base.eps(::Type{DoubleF64}) = double_eps # 2^-104
 
+Base.floatmin(::Type{DoubleF64}) = 2.0041683600089728e-292 # = 2^(-1022 + 53)
 # Base.realmin(::Type{DoubleF64}) = 2.0041683600089728e-292 # = 2^(-1022 + 53)
-# Base.realmin(::Type{DoubleF64}) = 2.0041683600089728e-292 # = 2^(-1022 + 53)
-# Base.realmax(::Type{DoubleF64}) = DoubleF64(1.79769313486231570815e+308, 9.97920154767359795037e+291);
+Base.floatmax(::Type{DoubleF64}) =
+    DoubleF64(1.79769313486231570815e+308, 9.97920154767359795037e+291)
 # Base.realmax(::Type{DoubleF64}) = DoubleF64(1.79769313486231570815e+308, 9.97920154767359795037e+291);
 
 Base.isnan(a::DoubleF64) = isnan(a.hi) || isnan(a.lo)

diff --git a/src/model_kit/intermediate_representation.jl b/src/model_kit/intermediate_representation.jl
@@ -233,9 +233,9 @@ function pow!(ir::IntermediateRepresentation, @nospecialize(a), @nospecialize(b)
         return add_op!(ir, OP_POW_INT, a, k)
     elseif k == 1 // 2
         return add_op!(ir, OP_SQRT, a)
+    else
+        error("Cannot handle exponent: $b")
     end
-
-    return add_op!(ir, OP_POW, a, b)
 end
 
 function split_off_minus_one(ex)
@@ -379,16 +379,14 @@ function split_into_num_denom!(ir, ex, cse, pse)
         if class(e) == :Pow
             x, k_ex = args(e)
             xref = expr_to_ir_statements!(ir, x, cse, pse)
-            if class(k_ex) == :Integer
-                k = convert(Int32, k_ex)
-                if k < 0
-                    push!(denoms, pow!(ir, xref, -k))
-                else
-                    push!(nums, pow!(ir, xref, k))
-                end
+            k = to_number(k_ex)
+            if k isa Basic
+                throw(ExprError("Cannot handle non-constant exponents"))
+            end
+            if k < 0
+                push!(denoms, pow!(ir, xref, -k))
             else
-                kref = expr_to_ir_statements!(ir, k_ex, cse, pse)
-                push!(nums, add_op!(ir, OP_POW, xref, kref))
+                push!(nums, pow!(ir, xref, k))
             end
         else
             push!(nums, expr_to_ir_statements!(ir, e, cse, pse))

diff --git a/src/model_kit/taylor.jl b/src/model_kit/taylor.jl
@@ -399,8 +399,13 @@ end
     taylor_impl(K, M - 1) do list, D
         v₀ = add_op!(list, OP_SQRT, D[:x, 0])
         d = mul!(list, 2, v₀)
-        ids = [v₀]
-        for k = 1:K
+        ids = []
+        push!(ids, v₀)
+        if K >= 1
+            v₁ = div!(list, D[:x, 1], d)
+            push!(ids, v₁)
+        end
+        for k = 2:K
             s = nothing
             for j = 1:(k-1)
                 s = muladd!(list, ids[j+1], ids[k-j+1], s)

diff --git a/test/test_systems.jl b/test/test_systems.jl
@@ -994,6 +994,12 @@ function small_rational()
     )
 end
 
+function sqrt_parameters()
+    @var x y a b
+
+    System([sqrt(a + b) * x^2 - y, (x * y + a - sqrt(b))^2 - 3], [x, y], [a, b])
+end
+
 TEST_SYSTEM_COLLECTION = [
     ("moments3", moments3()),
     ("fano_quintic", fano_quintic()),
@@ -1008,4 +1014,5 @@ TEST_SYSTEM_COLLECTION = [
     ("steiner", steiner()),
     ("four_bar", four_bar()),
     ("small_rational", small_rational()),
+    ("sqrt_parameters", sqrt_parameters()),
 ]