Default to exact arithmetic (#28)

* Default to exact arithmetic

* Bump minor

* Fix

Author: blegat

.github/workflows/ci.yml

@@ -17,7 +17,7 @@ jobs:
           - version: '1'
             os: ubuntu-latest
             arch: x64
-          - version: '1.0'
+          - version: '1.6'
             os: ubuntu-latest
             arch: x64
           - version: '1'

Project.toml

@@ -1,7 +1,7 @@
 name = "SemialgebraicSets"
 uuid = "8e049039-38e8-557d-ae3a-bc521ccf6204"
 repo = "https://github.com/JuliaAlgebra/SemialgebraicSets.jl.git"
-version = "0.2.5"
+version = "0.3.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -12,7 +12,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 [compat]
 MultivariatePolynomials = "0.4.2"
 MutableArithmetics = "0.3.1, 1"
-julia = "1"
+julia = "1.6"
 
 [extras]
 DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"

src/ideal.jl

@@ -6,6 +6,13 @@ struct EmptyPolynomialIdeal <: AbstractPolynomialIdeal end
 ideal(p::FullSpace, algo=defaultgröbnerbasisalgorithm(p)) = EmptyPolynomialIdeal()
 Base.rem(p::AbstractPolynomialLike, I::EmptyPolynomialIdeal) = p
 
+promote_for_division(::Type{T}) where {T} = T
+promote_for_division(::Type{T}) where {T<:Integer} = Rational{big(T)}
+promote_for_division(::Type{Rational{T}}) where {T} = Rational{big(T)}
+function promote_for_division(::Type{Complex{T}}) where {T}
+    return Complex{promote_for_division(T)}
+end
+
 mutable struct PolynomialIdeal{T, PT<:APL{T}, A<:AbstractGröbnerBasisAlgorithm} <: AbstractPolynomialIdeal
     p::Vector{PT}
     gröbnerbasis::Bool
@@ -15,7 +22,7 @@ function PolynomialIdeal{T, PT}(p::Vector{PT}, algo::A) where {T, PT<:APL{T}, A<
     PolynomialIdeal{T, PT, A}(p, false, algo)
 end
 function PolynomialIdeal(p::Vector{PT}, algo) where {T, PT<:APL{T}}
-    S = Base.promote_op(/, T, T)
+    S = promote_for_division(T)
     PolynomialIdeal{S, polynomialtype(PT, S)}(AbstractVector{polynomialtype(PT, S)}(p), algo)
 end
 function PolynomialIdeal{T, PT}() where {T, PT<:APL{T}}

src/schur.jl

@@ -18,20 +18,18 @@ function clusterordschur(M::AbstractMatrix{<:Real}, ɛ)
         # See bug JuliaLang/julia#...
         return Matrix{float(eltype(M))}(undef, 0, 0), Vector{Int}[]
     else
-        _clusterordschur(M, ɛ)
+        return _clusterordschur(M, ɛ)
     end
 end
+function _clusterordschur(M::AbstractMatrix{BigFloat}, ɛ)
+    Z, clusters = _clusterordschur(Float64.(M), ɛ)
+    return convert(Matrix{BigFloat}, Z), clusters
+end
 function _clusterordschur(M::AbstractMatrix{<:Real}, ɛ)
     # M = Z * T * Z' and "values" gives the eigenvalues
-    if VERSION >= v"0.7-"
-        sf = schur(M)
-        Z = sf.Z
-        v = sf.values
-    else
-        sf = schurfact(M)
-        Z = sf[:Z]
-        v = sf[:values]
-    end
+    sf = LinearAlgebra.schur(M)
+    Z = sf.Z
+    v = sf.values
     # documentation says that the error on the eigenvalues is ɛ * norm(T) / conditionnumber
     nT = norm(sf.T)
 
@@ -101,5 +99,5 @@ function _clusterordschur(M::AbstractMatrix{<:Real}, ɛ)
             break
         end
     end
-    Z, clusters
+    return Z, clusters
 end

src/variety.jl

@@ -11,7 +11,7 @@ defaultalgebraicsetlibrary(p::Vector{<:APL}, solveroralgo...) = defaultalgebraic
 mutable struct AlgebraicSet{T, PT<:APL{T}, A, S<:AbstractAlgebraicSolver, U} <: AbstractAlgebraicSet
     I::PolynomialIdeal{T, PT, A}
     projective::Bool
-    elements::Vector{Vector{T}}
+    elements::Vector{Vector{U}}
     elementscomputed::Bool
     iszerodimensional::Bool
     solver::S

test/macro.jl

@@ -4,24 +4,24 @@ struct DummySolver <: SemialgebraicSets.AbstractAlgebraicSolver end
     Mod.@polyvar x y
     @test isa(FullSpace(), FullSpace)
     V = @set x * y == 1
-    @test V isa AlgebraicSet
+    @test V isa AlgebraicSet{Rational{BigInt}}
     @test_throws ArgumentError addinequality!(V, x*y)
     @testset "Basic" begin
         S = @set x - y == 0 && x^2*y >= 1
         addequality!(S, x^2 - y)
         addinequality!(S, x + y - 1)
-        # Algebraic set forces `Float64`
-        @test S isa BasicSemialgebraicSet{Float64}
+        # Algebraic set forces `Rational{BigInt}`
+        @test S isa BasicSemialgebraicSet{Rational{BigInt}}
         @test S == basicsemialgebraicset(S.V, S.p)
-        @test sprint(show, S) == "{ (x, y) | x - y = 0, x^2 - y = 0, x^2*y - 1.0 ≥ 0, x + y - 1.0 ≥ 0 }"
-        @test sprint(show, MIME"text/plain"(), S) == "Basic semialgebraic Set defined by 2 equalities\n x - y = 0\n x^2 - y = 0\n2 inequalities\n x^2*y - 1.0 ≥ 0\n x + y - 1.0 ≥ 0\n"
-        @test S.V isa AlgebraicSet{Float64}
+        @test sprint(show, S) == "{ (x, y) | x - y = 0, x^2 - y = 0, x^2*y - 1//1 ≥ 0, x + y - 1//1 ≥ 0 }"
+        @test sprint(show, MIME"text/plain"(), S) == "Basic semialgebraic Set defined by 2 equalities\n x - y = 0\n x^2 - y = 0\n2 inequalities\n x^2*y - 1//1 ≥ 0\n x + y - 1//1 ≥ 0\n"
+        @test S.V isa AlgebraicSet{Rational{BigInt}}
         @test sprint(show, S.V) == "{ (x, y) | x - y = 0, x^2 - y = 0 }"
         @test sprint(show, MIME"text/plain"(), S.V) == "Algebraic Set defined by 2 equalities\n x - y = 0\n x^2 - y = 0\n"
-        @test S === MultivariatePolynomials.changecoefficienttype(S, Float64)
-        @test S.V === MultivariatePolynomials.changecoefficienttype(S.V, Float64)
+        @test S === MultivariatePolynomials.changecoefficienttype(S, Rational{BigInt})
+        @test S.V === MultivariatePolynomials.changecoefficienttype(S.V, Rational{BigInt})
         @test S.V.I === convert(typeof(S.V.I), S.V.I)
-        @test BasicSemialgebraicSet{Int, polynomialtype(x, Int)}() isa BasicSemialgebraicSet{Float64, polynomialtype(x, Float64)}
+        @test BasicSemialgebraicSet{Int, polynomialtype(x, Int)}() isa BasicSemialgebraicSet{Rational{BigInt}, polynomialtype(x, Rational{BigInt})}
         @test Int32(2)*x^2*y isa MultivariatePolynomials.AbstractTerm{Int32}
         Sf = MultivariatePolynomials.changecoefficienttype(S, Float32)
         @test Sf isa BasicSemialgebraicSet{Float32}
@@ -125,8 +125,8 @@ Algebraic Set defined by no equality
                   @set(x ≤ y && x == 1 && 2 + y == x),
                   @set(2 + y == x && x ≤ y && x == 1),
                   @set(2 + y == x && x ≤ y && x == 1)]
-            @test S isa BasicSemialgebraicSet{Float64}
-            @test S.V isa AlgebraicSet{Float64}
+            @test S isa BasicSemialgebraicSet{Rational{BigInt}}
+            @test S.V isa AlgebraicSet{Rational{BigInt}}
         end
 
         solver = DummySolver()
@@ -136,8 +136,8 @@ Algebraic Set defined by no equality
                   @set(x ≤ y && x == 1 && 2 + y == x, solver),
                   @set(2 + y == x && x ≤ y && x == 1, solver),
                   @set(2 + y == x && x ≤ y && x == 1, solver)]
-            @test S isa BasicSemialgebraicSet{Float64}
-            @test S.V isa AlgebraicSet{Float64}
+            @test S isa BasicSemialgebraicSet{Rational{BigInt}}
+            @test S.V isa AlgebraicSet{Rational{BigInt}}
             @test S.V.solver === solver
         end
     end

test/solve.jl

@@ -55,11 +55,7 @@ solver = ReorderedSchurMultiplicationMatricesSolver(sqrt(eps(Float64)), Mersenne
     Mod.@polyvar x y z
     V = @set x == y
     @test !iszerodimensional(V)
-    @static if VERSION >= v"0.7-"
-        @test_throws ErrorException iterate(V)
-    else
-        @test_throws ErrorException start(V)
-    end
+    @test_throws ErrorException iterate(V)
     @test_throws ErrorException length(V)
     V = @set 4x^2 == -5x && 3x^3 == 0 solver
     @test V.solver.solver === solver