Merge pull request #48 from rprebet/output-interval

Add the possibility to output isolating intervals

Author: ederc

src/algorithms/solvers.jl

@@ -1,6 +1,6 @@
 import msolve_jll: libmsolve
 
-export real_solutions, rational_solutions, rational_parametrization
+export inter_solutions, real_solutions, rational_solutions, rational_parametrization
 
 @doc Markdown.doc"""
     _get_rational_parametrization(nr::Int32, lens::Vector{Int32}, cfs::Ptr{BigInt})
@@ -134,6 +134,7 @@ function _core_msolve(
         C,x = polynomial_ring(QQ,"x")
         I.rat_param = RationalParametrization(Symbol[], ZZRingElem[], C(-1), C(-1), PolyRingElem[])
         I.real_sols = QQFieldElem[]
+        I.inter_sols = Vector{QQFieldElem}[]
         return I.rat_param, I.real_sols
     end
     [nterms += jl_len[i] for i=1:jl_ld]
@@ -152,6 +153,7 @@ function _core_msolve(
                                           rat_param[3], rat_param[4])
     if jl_nb_sols == 0
         I.real_sols = QQFieldElem[]
+        I.inter_sols = Vector{QQFieldElem}[]
         return rat_param, Vector{QQFieldElem}[]
     end
 
@@ -164,25 +166,32 @@ function _core_msolve(
     # end
 
     #= solutions are returned as intervals, i.e. a minimum and a maximum entry for
-     = the numerator and denominator; thus we sum up and divide by  =#
+     = the numerator and denominator; we also sum up and divide by 2 =#
+
+    inter_solutions = Vector{Vector{Vector{QQFieldElem}}}(undef, jl_nb_sols)
     solutions = Vector{Vector{QQFieldElem}}(undef, jl_nb_sols)
 
     len = 2*jl_nb_sols*nr_vars
     i = 1
     k = 1
     while i <= len
         j = 1
+        inter_tmp = Vector{Vector{QQFieldElem}}(undef, nr_vars)
         tmp = Vector{QQFieldElem}(undef, nr_vars)
         while j <= nr_vars
-            tmp[perm[j]]  = QQFieldElem(unsafe_load(jl_sols_num, i)) >> Int64(unsafe_load(jl_sols_den, i))
-            tmp[perm[j]] += QQFieldElem(unsafe_load(jl_sols_num, i+1)) >> Int64(unsafe_load(jl_sols_den, i+1))
-            tmp[perm[j]] = tmp[perm[j]] >> 1
+            inter_tmp_coeff = Vector{QQFieldElem}(undef, 2)
+            inter_tmp_coeff[1] = QQFieldElem(unsafe_load(jl_sols_num, i)) >> Int64(unsafe_load(jl_sols_den, i))
+            inter_tmp_coeff[2] = QQFieldElem(unsafe_load(jl_sols_num, i+1)) >> Int64(unsafe_load(jl_sols_den, i+1))
+            inter_tmp[perm[j]] = inter_tmp_coeff
+            tmp[perm[j]] = sum(inter_tmp_coeff) >> 1
             i += 2
             j += 1
         end
+        inter_solutions[k] = inter_tmp
         solutions[k] = tmp
         k += 1
     end
+    I.inter_sols = inter_solutions
     I.real_sols = solutions
 
     ccall((:free_msolve_julia_result_data, libmsolve), Nothing ,
@@ -401,3 +410,62 @@ function real_solutions(
 
     return I.real_sols
 end
+
+@doc Markdown.doc"""
+    inter_solutions(I::Ideal{T} where T <: MPolyRingElem, <keyword arguments>)
+
+Given an ideal `I` with a finite solution set over the complex numbers, return
+the intervals with a given precision (default 32 bits), containing the roots of the ideal.
+
+**Note**: At the moment only QQ is supported as ground field. If the dimension of the ideal
+is greater than zero an empty array is returned.
+
+# Arguments
+- `I::Ideal{T} where T <: MPolyRingElem`: input generators.
+- `initial_hts::Int=17`: initial hash table size `log_2`.
+- `nr_thrds::Int=1`: number of threads for parallel linear algebra.
+- `max_nr_pairs::Int=0`: maximal number of pairs per matrix, only bounded by minimal degree if `0`.
+- `la_option::Int=2`: linear algebra option: exact sparse-dense (`1`), exact sparse (`2`, default), probabilistic sparse-dense (`42`), probabilistic sparse(`44`).
+- `info_level::Int=0`: info level printout: off (`0`, default), summary (`1`), detailed (`2`).
+- `precision::Int=32`: bit precision for the computed solutions.
+
+# Examples
+```jldoctest
+julia> using AlgebraicSolving
+
+julia> R,(x1,x2,x3) = polynomial_ring(QQ, ["x1","x2","x3"])
+(Multivariate polynomial ring in 3 variables over QQ, QQMPolyRingElem[x1, x2, x3])
+
+julia> I = Ideal([x1+2*x2+2*x3-1, x1^2+2*x2^2+2*x3^2-x1, 2*x1*x2+2*x2*x3-x2])
+QQMPolyRingElem[x1 + 2*x2 + 2*x3 - 1, x1^2 - x1 + 2*x2^2 + 2*x3^2, 2*x1*x2 + 2*x2*x3 - x2]
+
+julia> inter_solutions(I)
+4-element Vector{Vector{Vector{QQFieldElem}}}:
+ [[1354207349//2147483648, 2708414699//4294967296], [1354207349//4294967296, 677103675//2147483648], [-355532291286331190123863132844989723573//2722258935367507707706996859454145691648, -1422129165145324760495452531379958894291//10889035741470030830827987437816582766592]]
+ [[1, 1], [0, 0], [0, 0]]
+ [[972985841//4294967296, 486492921//2147483648], [60811615//536870912, 486492921//4294967296], [93053303113782607831679264095876317083//340282366920938463463374607431768211456, 372213212455130431326717056383505268333//1361129467683753853853498429727072845824]]
+ [[1431655765//4294967296, 715827883//2147483648], [-1//4294967296, 1//4294967296], [1814839290245005138471331239636097127765//5444517870735015415413993718908291383296, 907419645122502569235665619818048563883//2722258935367507707706996859454145691648]]
+```
+"""
+function inter_solutions(
+    I::Ideal{T} where T <: MPolyRingElem;     # input generators
+    initial_hts::Int=17,                  # hash table size, default 2^17
+    nr_thrds::Int=1,                      # number of threads
+    max_nr_pairs::Int=0,                  # number of pairs maximally chosen
+                                          # in symbolic preprocessing
+    la_option::Int=2,                     # linear algebra option
+    info_level::Int=0,                    # info level for print outs
+    precision::Int=32                     # precision of the solution set
+    )
+
+isdefined(I, :inter_sols) ||
+_core_msolve(I,
+             initial_hts = initial_hts,
+             nr_thrds = nr_thrds,
+             max_nr_pairs = max_nr_pairs,
+             la_option = la_option,
+             info_level = info_level,
+             precision = precision)
+
+return I.inter_sols
+end

src/types.jl

@@ -29,6 +29,7 @@ mutable struct Ideal{T <: MPolyRingElem}
     gens::Vector{T}
     dim::Int
     gb::Dict{Int, Vector{T}}
+    inter_sols::Vector{Vector{Vector{QQFieldElem}}}
     real_sols::Vector{Vector{QQFieldElem}}
     rat_sols::Vector{Vector{QQFieldElem}}
     rat_param::RationalParametrization

test/algorithms/solvers.jl

@@ -12,9 +12,18 @@
         [1611414365//8589934592, 673053615//8589934592, 632173751//8589934592, 708759148891639684402860468800417934359477//2787593149816327892691964784081045188247552],
         [2863311531//8589934592, 0, 0, 14518714321960041107770649917088777022123//43556142965880123323311949751266331066368]
                             ]
+    inter_sols = Vector{Vector{QQFieldElem}}[
+        [[2431274387//4294967296, 607818597//1073741824], [320390731//2147483648, 640781463//4294967296], [1097534103//4294967296, 137191763//536870912], [-16357136847304379909836796974131395209669//87112285931760246646623899502532662132736, -4089284211826094977459199243532848802417//21778071482940061661655974875633165533184]],
+        [[1889817751//4294967296, 236227219//536870912], [1319238065//4294967296, 659619033//2147483648], [113559211//1073741824, 454236845//4294967296], [-46316921746739551124221118038944540401681//348449143727040986586495598010130648530944, -2894807609171221945263819877434033775105//21778071482940061661655974875633165533184]],
+        [[1, 1], [0, 0], [0, 0], [0, 0]],
+        [[3205239737//4294967296, 1602619869//2147483648], [501382663//2147483648, 1002765327//4294967296], [-792885089//4294967296, -24777659//134217728], [3478667982170137579718084770105695880909725//44601490397061246283071436545296723011960832, 27829343857361100637744678160845567047277801//356811923176489970264571492362373784095686656]],
+        [[402853591//2147483648, 805707183//4294967296], [336526807//4294967296, 42065851//536870912], [316086875//4294967296, 79021719//1073741824], [177189787222909921100715117200104483589869//696898287454081973172991196020261297061888, 354379574445819842201430234400208967179739//1393796574908163946345982392040522594123776]],
+        [[1431655765//4294967296, 715827883//2147483648], [-1//4294967296, 1//4294967296], [-1//4294967296, 1//4294967296], [7259357160980020553885324958544388511061//21778071482940061661655974875633165533184, 3629678580490010276942662479272194255531//10889035741470030830827987437816582766592]]
+                                            ]
     rat_sols = Vector{QQFieldElem}[[49, 0, 0, 0], [49//3, 0, 0, 1//3]]
 
     @test sols == real_solutions(I)
+    @test inter_sols == inter_solutions(I)
     @test rat_sols == rational_solutions(I)
     @test I.real_sols == real_solutions(I)
 
@@ -45,10 +54,12 @@
 
     I = Ideal([x1^2-x2, x1*x3-x4, x2*x4-12, x4^3-x3^2])
     real_solutions(I)
+    inter_solutions(I)
     @test I.rat_param.vars == Symbol[]
 
     I = Ideal([x1^2-x2, x1*x3, x2-12])
 	@test_throws ErrorException real_solutions(I)
+    @test_throws ErrorException inter_solutions(I)
 	@test_throws ErrorException rational_solutions(I)
 
     # check variable permutation