Merge pull request #385 from SciML/small-unimportant-fix

Let each numerator have its own denominator

src/RationalFunctionFields/IdealMQS.jl

@@ -20,17 +20,13 @@ https://mediatum.ub.tum.de/doc/685465/685465.pdf
 mutable struct IdealMQS{T} <: AbstractBlackboxIdeal
     funcs_den_nums::Vector{Vector{T}}
     den_lcm_orig::T
+    # Indices of pivot elements for each component of funcs_den_nums 
+    pivots_indices::Vector{Int}
     parent_ring_param::Generic.MPolyRing{T}
     # Numerators and denominators over QQ
-    # NB: this lists may have different length
     nums_qq::Vector{T}
     dens_qq::Vector{T}
     sat_qq::T
-    # dens_indices[i] is a pair of form (from, to).
-    # Denominator as index i corresponds to numerators at indices [from..to]
-    dens_indices::Vector{Tuple{Int, Int}}
-    # Indices of pivot elements for each component of funcs_den_nums 
-    pivots_indices::Vector{Int}
     den_lcm::T
     sat_var_index::Int
     # Numerators and denominators over GF. 
@@ -72,7 +68,6 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal
         # degree and length. Such element will serve as a normalizing term for
         # the component
         funcs_den_nums = map(plist -> filter(!iszero, plist), funcs_den_nums)
-        # funcs_den_nums = filter(plist -> length(plist) > 1, funcs_den_nums)
         @assert !isempty(funcs_den_nums) "All elements of the ideal are zero"
         pivots =
             map(plist -> findmin(p -> (total_degree(p), length(p)), plist), funcs_den_nums)
@@ -113,14 +108,9 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal
         den_lcm_sat = parent_ring_change(den_lcm, R_sat)
         sat_qq = den_lcm_sat * t_sat - 1
         # We construct the array of numerators nums_qq and the array of
-        # denominators dens_qq. Since there are usually more unique numerators
-        # than denominators, we condense the array dens_qq and produce an
-        # additional array dens_indices. The element dens_indices[i] tells us
-        # the indices of the numerators that correspond to a given denominator
-        # dens_qq[i]
+        # denominators dens_qq of the same length
         nums_qq = empty(funcs_den_nums[1])
         dens_qq = empty(funcs_den_nums[1])
-        dens_indices = Vector{Tuple{Int, Int}}()
         for i in 1:length(funcs_den_nums)
             plist = funcs_den_nums[i]
             den = plist[pivots_indices[i]]
@@ -130,8 +120,6 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal
                 matching = :byindex,
                 shift = Int(sat_var_index == 1),
             )
-            push!(dens_qq, den)
-            push!(dens_indices, (length(nums_qq) + 1, length(nums_qq) + length(plist) - 1))
             for j in 1:length(plist)
                 j == pivots_indices[i] && continue
                 num = plist[j]
@@ -141,23 +129,25 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal
                     matching = :byindex,
                     shift = Int(sat_var_index == 1),
                 )
-                push!(nums_qq, num)
+                G = gcd(num, den)
+                _num, _den = num / G, den / G
+                push!(nums_qq, _num)
+                push!(dens_qq, _den)
             end
         end
         parent_ring_param, _ = polynomial_ring(ring, varnames, internal_ordering = ordering)
         @debug "Constructed MQS ideal in $R_sat with $(length(nums_qq) + 1) elements"
-        @assert length(pivots_indices) == length(dens_indices) == length(dens_qq)
         @assert length(pivots_indices) == length(funcs_den_nums)
+        @assert length(nums_qq) == length(dens_qq)
 
         new{elem_type(R_sat)}(
             funcs_den_nums,
             den_lcm_orig,
+            pivots_indices,
             parent_ring_param,
             nums_qq,
             dens_qq,
             sat_qq,
-            dens_indices,
-            pivots_indices,
             den_lcm,
             sat_var_index,
             Dict(),
@@ -212,11 +202,8 @@ function fractionfree_generators_raw(mqs::IdealMQS)
     polys = Vector{elem_type(big_ring)}(undef, length(nums_qq) + 1)
     @inbounds for i in 1:length(dens_qq)
         den_y, den_x = dens_y[i], dens_x[i]
-        span = mqs.dens_indices[i]
-        for j in span[1]:span[2]
-            num_y, num_x = nums_y[j], nums_x[j]
-            polys[j] = num_y * den_x - den_y * num_x
-        end
+        num_y, num_x = nums_y[i], nums_x[i]
+        polys[i] = num_y * den_x - den_y * num_x
     end
     polys[end] = sat_y
     main_var_indices = 1:(length(varnames) + 1)
@@ -256,9 +243,9 @@ function ParamPunPam.specialize_mod_p(
     @assert haskey(mqs.cached_nums_gf, K_1)
     nums_gf, dens_gf, sat_gf =
         mqs.cached_nums_gf[K_1], mqs.cached_dens_gf[K_1], mqs.cached_sat_gf[K_1]
-    dens_indices = mqs.dens_indices
     K_2 = base_ring(nums_gf[1])
     @assert K_1 == K_2
+    @assert length(nums_gf) == length(dens_gf)
     @assert length(point) == nvars(ParamPunPam.parent_params(mqs))
     # +1 actual variable because of the saturation!
     @assert length(point) + 1 == nvars(parent(nums_gf[1]))
@@ -269,11 +256,8 @@ function ParamPunPam.specialize_mod_p(
     @inbounds for i in 1:length(dens_gf_spec)
         den, den_spec = dens_gf[i], dens_gf_spec[i]
         iszero(den_spec) && __throw_unlucky_evaluation("Point: $point")
-        span = dens_indices[i]
-        for j in span[1]:span[2]
-            num, num_spec = nums_gf[j], nums_gf_spec[j]
-            polys[j] = num * den_spec - den * num_spec
-        end
+        num, num_spec = nums_gf[i], nums_gf_spec[i]
+        polys[i] = num * den_spec - den * num_spec
     end
     polys[end] = sat_gf
     if !saturated
@@ -285,7 +269,6 @@ end
 function specialize(mqs::IdealMQS, point::Vector{Nemo.QQFieldElem}; saturated = true)
     @debug "Evaluating MQS ideal over QQ at $point"
     nums_qq, dens_qq, sat_qq = mqs.nums_qq, mqs.dens_qq, mqs.sat_qq
-    dens_indices = mqs.dens_indices
     K = base_ring(mqs)
     @assert length(point) == nvars(ParamPunPam.parent_params(mqs))
     # +1 actual variable because of the saturation!
@@ -297,11 +280,8 @@ function specialize(mqs::IdealMQS, point::Vector{Nemo.QQFieldElem}; saturated =
     @inbounds for i in 1:length(dens_qq_spec)
         den, den_spec = dens_qq[i], dens_qq_spec[i]
         iszero(den_spec) && __throw_unlucky_evaluation("Point: $point")
-        span = dens_indices[i]
-        for j in span[1]:span[2]
-            num, num_spec = nums_qq[j], nums_qq_spec[j]
-            polys[j] = num * den_spec - den * num_spec
-        end
+        num, num_spec = nums_qq[i], nums_qq_spec[i]
+        polys[i] = num * den_spec - den * num_spec
     end
     polys[end] = sat_qq
     if !saturated

test/RationalFunctionFields/raw_generators.jl

@@ -47,14 +47,14 @@ import Groebner
     gb = Groebner.groebner(eqs, ordering = Groebner.DegRevLex())
     # GB is linear
     @test length(gb) == length(gens(parent(eqs[1])))
-    expected = 10e6
+    expected = 9202476
     str = join(map(string, eqs), ",")
     @info "" length(str)
     @test abs(length(str) - expected) / expected * 100 < 5
 
     # Part 2: over Q
     eqs = StructuralIdentifiability.fractionfree_generators_raw(rff.mqs)[1]
-    expected = 26e6
+    expected = 21486079
     str = join(map(string, eqs), ",")
     @info "" length(str)
     @test abs(length(str) - expected) / expected * 100 < 5