Adding new functions rank_interpolation, rank_interpolation_mc and evaluation_points.

src/Fields.jl

@@ -47,3 +47,32 @@ function canonical_unit(x::FieldElement)
   iszero(x) && return one(x)
   return x
 end
+
+@doc raw"""
+    evaluation_points(K::Field, n::Int) -> Vector{FieldElem}
+
+If $K$ contains at least $n$ elements, this function returns a vector $v$ with $n$ distinct elements from $K$. Otherwise it returns an empty vector. 
+"""
+function evaluation_points(K::Field, n::Int)
+   @req n > 0 "n must be positive."
+   if is_finite(K)
+      v = elem_type(K)[]
+      if order(K) < n
+         return v
+      else
+         append!(v, Iterators.take(K, n))
+      end
+   else
+      @assert is_zero(characteristic(K)) "This should not happen. Please open a bug report."
+      v = Vector{elem_type(K)}(undef, n)
+      v[1] = K(div(-n, 2))
+      for i = 2:n
+         v[i] = v[i-1] + 1
+      end
+   end
+   return v
+end
+
+function ext_of_degree(K::Field, n::Int)
+   error("Cannot extend field")
+end

src/Matrix.jl

@@ -2873,6 +2873,219 @@ function rank(M::MatrixElem{T}) where {T <: FieldElement}
    return lu!(P, A)
 end
 
+@doc raw"""
+    rank_interpolation(M::MatrixElem{T}) where {T <: PolyRingElem} -> Int
+    rank_interpolation(M::MatrixElem{T}) where {T <: MPolyRingElem} -> Int
+    rank_interpolation(M::MatrixElem{T}) where {T <: AbstractAlgebra.Generic.RationalFunctionFieldElem} -> Int
+
+Returns the rank of $A$ using an interpolation-like method. 
+
+# Examples
+
+```jldoctest
+julia> Qx, x = polynomial_ring(QQ, :x);
+
+julia> AbstractAlgebra.rank_interpolation(matrix(Qx, 2, 2, [x 2x; x^2 1]))
+2
+
+julia> Qy, y = rational_function_field(QQ, :y);
+
+julia> AbstractAlgebra.rank_interpolation(matrix(Qy, 2, 2, [1//y -y^2; 3 2-y]))
+2
+```
+"""
+function rank_interpolation(M::MatElem{T}) where {T <: PolyRing}
+   n = nrows(M)
+   m = ncols(M)
+   if is_zero(n) || is_zero(m)
+      return 0
+   end
+   Kx = base_ring(M)
+   K = base_ring(Kx)
+   #The maximum degree of det(M') is calculated where M' is an arbitrary quadratic submatrix of M.
+   min_ = min(n, m)
+   if (min_ == n)
+      maxdetdeg = sum(maximum(degree(M[i, j]) for j in 1:m) for i in 1:n)
+   else 
+      maxdetdeg = sum(maximum(degree(M[i, j]) for i in 1:n) for j in 1:m)
+   end
+   r = 0
+   eval_set = evaluation_points(K, maxdetdeg+1)
+   if !is_empty(eval_set)
+      #rank(M) is calculated by computing the rank of det_deg+1 matrices, evaluated in an element of eval_set, respectively.
+      for elem in eval_set
+         M_eval = map_entries(p -> evaluate(p, elem), M)
+         r = max(rank_interpolation(M_eval), r)
+         if r == min_
+            break
+         end
+      end
+      return r
+   else
+      #function get_eval_set returns an empty set if and only if K is a finite field and order(K) < det_deg+1 holds. 
+      #In this case a field extension L of K such that order(L) >= det_deg+1 is constructed.
+      #d is the smallest natural number such that order(K)^d >= det_deg+1.
+      d = clog(order(K), ZZ(maxdetdeg+1))
+      L, l = ext_of_degree(K, d)
+      Lx, _ = polynomial_ring(L, var(Kx)) 
+      #The given matrix M is embedded into the space of matrices over Lx
+      A = matrix(Lx, n, m, [map_coefficients(l, M[i, j]; parent = Lx) for i in 1:n, j in 1:m])
+      return rank_interpolation(A)
+   end
+end
+
+function rank_interpolation(M::MatElem{T}) where {T <: MPolyRing}
+   n = nrows(M)
+   m = ncols(M)
+   if is_zero(n) || is_zero(m)
+      return 0
+   end
+   Kx = base_ring(M)
+   K = base_ring(Kx)
+   num_vars = number_of_variables(Kx)
+   #The maximum degree of det(M') is calculated where M' is an arbitrary quadratic submatrix of M.
+   min_ = min(n, m)
+   if min_ == n
+      maxdetdeg = [sum(maximum(degree(M[i, j], k) for j in 1:m) for i in 1:n) for k in 1:num_vars]
+   else
+      maxdetdeg = [sum(maximum(degree(M[i, j], k) for i in 1:n) for j in 1:m) for k in 1:num_vars] 
+   end
+   r = 0
+   eval_set = Vector{Vector{elem_type(K)}}(undef, num_vars)
+   for i = 1:num_vars
+      eval_set[i] = evaluation_points(K, maxdetdeg[i]+1)
+      if is_empty(eval_set[i])
+         #function get_eval_set returns an empty set if and only if K is a finite field and order(K) < det_deg+1 holds. 
+         #In this case a field extension L of K such that order(L) >= det_deg+1 is constructed.
+         #d is the smallest natural number such that order(K)^d >= det_deg+1.
+         d = clog(order(K), ZZ(maximum(maxdetdeg)+1))
+         L, l = ext_of_degree(K, d)
+         Lx, _ = polynomial_ring(L, symbols(Kx)) 
+         #The given matrix M is embedded into the space of matrices over Lx
+         A = matrix(Lx, n, m, [map_coefficients(l, M[i, j]; parent = Lx) for i in 1:n, j in 1:m])
+         return rank_interpolation(A)
+      end
+   end
+   #rank(M) is calculated by computing the rank of (det_deg+1)^number_of_variables(Kx) matrices, evaluated in elements of eval_set.
+   for tup in ProductIterator(eval_set)
+      M_eval = map_entries(p -> evaluate(p, tup), M)
+      r = max(rank_interpolation(M_eval), r)
+      if r == min_
+         break
+      end
+   end
+   return r
+end
+
+function rank_interpolation(M::MatElem{<: RingElem})
+   return rank(M)
+end
+
+@doc raw"""
+    rank_interpolation_mc(M::MatrixElem{T}, ϵ::Float64) where {T <: PolyRingElem} -> Int
+    rank_interpolation_mc(M::MatrixElem{T}, ϵ::Float64) where {T <: MPolyRingElem} -> Int
+    rank_interpolation_mc(M::MatrixElem{T}, ϵ::Float64) where {T <: AbstractAlgebra.Generic.RationalFunctionFieldElem} -> Int
+
+Returns the rank of $A$ with error probability < $ϵ$ using an interpolation-like method.
+
+# Examples
+
+```jldoctest
+julia> Qx, x = polynomial_ring(QQ, :x);
+
+julia> AbstractAlgebra.rank_interpolation_mc(matrix(Qx, 2, 2, [x 2x; x^2 1]), 0.01)
+2
+
+julia> Qy, y = rational_function_field(QQ, :y);
+
+julia> AbstractAlgebra.rank_interpolation_mc(matrix(Qy, 2, 2, [1//y -y^2; 3 2-y]), 0.00001)
+2
+```
+"""
+function rank_interpolation_mc(M::MatElem{T}, err::Float64) where {T <: PolyRingElem}
+   n = nrows(M)
+   m = ncols(M)
+   if is_zero(n) || is_zero(m)
+      return 0
+   end
+   Kx = base_ring(M)
+   K = base_ring(Kx)
+   min_ = min(n, m)
+   #The maximum degree of det(M') is calculated where M' is an arbitrary quadratic submatrix of M.
+   maxdetdeg = min_*maximum(degree(M[i, j]) for i in 1:n, j in 1:n)
+   S = evaluation_points(K, 10*maxdetdeg)
+   #k is the minimum amount of evaluations of M needed to compute the correct rank of M with error probability < err
+   k = ceil(Base.log(10, 1/err))
+   r = 0
+   if !is_empty(S)
+      #rank(M) is calculated by computing the rank of k matrices, evaluated in elements of eval_set, and taking the maximum
+      for _ = 1:k
+         a = rand(S)
+         M_eval = map_entries(p -> evaluate(p, a), M)
+         r = max(rank(M_eval), r)
+         if r == min_
+            return r
+         end
+      end
+   else
+      #function get_eval_set returns an empty set if and only if K is a finite field and order(K) < det_deg+1 holds. 
+      #In this case a field extension L of K such that order(L) >= det_deg+1 is constructed.
+      #d is the smallest natural number such that order(K)^d >= det_deg+1.
+      d = clog(order(K), ZZ(maxdetdeg*10))
+      L, l = ext_of_degree(K, d)
+      Lx, _ = polynomial_ring(L, var(Kx)) 
+      #The given matrix M is embedded into the space of matrices over Lx
+      A = matrix(Lx, n, m, [map_coefficients(l, M[i, j]; parent = Lx) for i in 1:n, j in 1:m])
+      return rank_interpolation_mc(A, err)
+   end
+   return r
+end
+
+function rank_interpolation_mc(M::MatElem{T}, err::Float64) where {T <: MPolyRingElem}
+   n = nrows(M)
+   m = ncols(M)
+   if is_zero(n) || is_zero(m)
+      return 0
+   end
+   Kx = base_ring(M)
+   K = base_ring(Kx)
+   num_vars = number_of_variables(Kx)
+   min_ = min(n, m)
+   #The maximum degree of det(M') is calculated where M' is an arbitrary quadratic submatrix of M.
+   maxdetdeg = min_*maximum(degree(M[i, j], k) for i in 1:n, j in 1:m, k in 1:num_vars)
+   S = evaluation_points(K, 10*maxdetdeg)
+   #k is the minimum amount of evaluations of M needed to compute the correct rank of M with error probability < err
+   k = ceil(Base.log(10, 1/err))
+   r = 0
+   if !is_empty(S)
+      #rank(M) is calculated by computing the rank of k matrices, evaluated in elements of eval_set, and taking the maximum
+      for _ = 1:k
+         a = Vector{elem_type(K)}(undef, num_vars)
+         for i = 1:num_vars
+            a[i] = rand(S)
+         end
+         M_eval = map_entries(p -> evaluate(p, a), M)
+         r = max(rank(M_eval), r)
+         if r == min_
+            return r
+         end
+      end
+      return r
+   else
+      #function get_eval_set returns an empty set if and only if K is a finite field and order(K) < det_deg*s holds. 
+      #In this case a field extension L of K is constructed such that order(L) >= det_deg*s.
+      #d is the smallest natural number such that order(K)^d >= det_deg*s.
+      d = clog(order(K), ZZ(maxdetdeg*10))
+      L, l = ext_of_degree(K, d)
+      Lx, _ = polynomial_ring(L, symbols(Kx)) 
+      #The given matrix M is embedded into the space of matrices over Lx
+      A = matrix(Lx, n, m, [map_coefficients(l, M[i, j]; parent = Lx) for i in 1:n, j in 1:m])
+      return rank_interpolation_mc(A, err)
+   end
+end
+
+
+
 ###############################################################################
 #
 #   Linear solving

src/generic/RationalFunctionField.jl

@@ -683,3 +683,76 @@ end
 function is_squarefree(f::Union{PolyRingElem{T}, MPolyRingElem{T}}) where {T <: RationalFunctionFieldElem}
   is_squarefree(map_coefficients(data, f; cached = false))
 end
+
+################################################################################
+#
+#  Generation of finite subset and rank interpolation
+#
+################################################################################
+
+function evaluation_points(K::RationalFunctionField, n::Int)
+   @req n > 0 "n must be positive."
+   F = base_ring(K)
+   v = evaluation_points(F, n)
+   if is_empty(v)
+      v = elem_type(K)[]
+      base_v = evaluation_points(F, Int(order(F)))
+      d = clog(order(F), ZZ(n))
+      genKpowers = powers(gen(K), d - 1)
+      for coeffs in ProductIterator(base_v, d; inplace=true)
+         push!(v, sum(coeffs[j] * genKpowers[j] for j in 1:d))
+         if length(v) == n
+            break
+         end
+      end
+   end
+   return v
+end
+
+function rank_interpolation(M::MatrixElem{<: RationalFunctionFieldElem})
+   n = nrows(M)
+   m = ncols(M)
+   if is_zero(n) || is_zero(m)
+      return 0
+   end
+   Kx = base_ring(Generic.underlying_fraction_field(base_ring(M)))
+   A = deepcopy(M)
+   #All rows/columns of M are multiplied with polynomials such that all entries of M have denominator 1.
+   #Then M can be considered as a matrix over a polynomial ring. 
+   for i = 1:n
+      for j = 1:m
+         if !is_one(denominator(A[i, j]))
+            if n < m
+               multiply_column!(A, denominator(A[i, j]), j)
+            else
+               multiply_row!(A, denominator(A[i, j]), i)
+            end
+         end
+      end
+   end
+   return rank_interpolation(matrix(Kx, n, m, [numerator(A[i, j]) for i in 1:n, j in 1:m]))
+end
+
+function rank_interpolation_mc(M::MatrixElem{<: RationalFunctionFieldElem}, err::Float64)
+   n = nrows(M)
+   m = ncols(M)
+   if is_zero(n) || is_zero(m)
+      return 0
+   end
+   Kx = base_ring(Generic.underlying_fraction_field(base_ring(M)))
+   A = deepcopy(M)
+   #All rows/columns of M are multiplied with polynomials such that all entries of M have denominator 1.
+   #Then M can be considered as a matrix over some polynomial ring.
+   for i = 1:n
+      for j = 1:m
+         if is_one(denominator(A[i, j]))
+            if n < m
+               multiply_column!(A, denominator(A[i, j]), j)
+            else
+               multiply_row!(A, denominator(A[i, j]), i)
+            end
+         end
+      end
+   end
+   return rank_interpolation_mc(matrix(Kx, n, m, [numerator(A[i, j]) for i in 1:n, j in 1:m]), err)
+end

src/generic/imports.jl

@@ -130,6 +130,7 @@ import ..AbstractAlgebra: divrem
 import ..AbstractAlgebra: domain
 import ..AbstractAlgebra: elem_type
 import ..AbstractAlgebra: evaluate
+import ..AbstractAlgebra: evaluation_points
 import ..AbstractAlgebra: exp
 import ..AbstractAlgebra: exponent_vector
 import ..AbstractAlgebra: exponent_vectors
@@ -204,6 +205,8 @@ import ..AbstractAlgebra: pretty_sort
 import ..AbstractAlgebra: primpart
 import ..AbstractAlgebra: promote_rule
 import ..AbstractAlgebra: pseudorem
+import ..AbstractAlgebra: rank_interpolation
+import ..AbstractAlgebra: rank_interpolation_mc
 import ..AbstractAlgebra: reduce!
 import ..AbstractAlgebra: reduced_form
 import ..AbstractAlgebra: remove

test/Matrix-test.jl

@@ -214,3 +214,46 @@ end
   @test M * QQ(1) == QQ(1) * M == QQ.(M)
   @test N * ZZ(1) == ZZ(1) * N == QQ.(N)
 end
+
+@testset "Rank Interpolation" begin
+  Fx, x = polynomial_ring(GF(29), :x)
+  Qy, y = polynomial_ring(QQ, 3, :y)
+  Qz, z = rational_function_field(QQ, :z)
+
+  A = rand(matrix_space(Fx, 2, 2), 0:5)
+  B = rand(matrix_space(Qy, 2, 2), 1:3, 1:4, 1:2)
+  C = rand(matrix_space(Qz, 2, 2), 0:5, 1:5)
+
+  @test rank(A) == AbstractAlgebra.rank_interpolation(A)
+  @test rank(A) >= AbstractAlgebra.rank_interpolation_mc(A, 0.01)
+  @test rank(B) == AbstractAlgebra.rank_interpolation(B)
+  @test rank(B) >= AbstractAlgebra.rank_interpolation_mc(B, 0.01)
+  @test rank(C) == AbstractAlgebra.rank_interpolation(C)
+  @test rank(C) >= AbstractAlgebra.rank_interpolation_mc(C, 0.01)
+
+  multiply_row!(A, 0, 1)
+  multiply_row!(B, 0, 1)
+  multiply_row!(C, 0, 1)
+
+  @test rank(A) == AbstractAlgebra.rank_interpolation(A, 0.01)
+  @test rank(A) >= AbstractAlgebra.rank_interpolation_mc(A)
+  @test rank(B) == AbstractAlgebra.rank_interpolation(B, 0.01)
+  @test rank(B) >= AbstractAlgebra.rank_interpolation_mc(B)
+  @test rank(C) == AbstractAlgebra.rank_interpolation(C, 0.01)
+  @test rank(C) >= AbstractAlgebra.rank_interpolation_mc(C)
+
+  A = rand(matrix_space(Fx, 0, 2), 0:5)
+  B = rand(matrix_space(Qy, 0, 2), 1:3, 1:4, 1:2)
+  C = rand(matrix_space(Qz, 2, 0), 0:5, 1:5)
+
+  @test rank(A) == AbstractAlgebra.rank_interpolation(A)
+  @test rank(A) >= AbstractAlgebra.rank_interpolation_mc(A, 0.01)
+  @test rank(B) == AbstractAlgebra.rank_interpolation(B)
+  @test rank(B) >= AbstractAlgebra.rank_interpolation_mc(B, 0.01)
+  @test rank(C) == AbstractAlgebra.rank_interpolation(C)
+  @test rank(C) >= AbstractAlgebra.rank_interpolation_mc(C, 0.01)
+
+  # A = rand(matrix_space(Fx, 10, 10), 0:5)
+  # @test_throws AbstractAlgebra.rank_interpolation(A)
+  # @test_throws AbstractAlgebra.rank_interpolation_mc(A, 0.01)
+end