sumiya11 / Groebner.jl

Groebner bases in (almost) pure Julia
https://sumiya11.github.io/Groebner.jl/
GNU General Public License v2.0
66 stars 13 forks source link

Allow polynomial coefficient fields to be finite_field where field degree > 1; and also algebraic field extensions of QQ #162

Open xyz333math opened 1 week ago

xyz333math commented 1 week ago

Currently the only finite fields allowed seem to be the Nemo fields of degree one, for example GF(p) where p is a prime integer; I think that GF is the same as finite_field in Nemo.

I did a quick experiment to find out that function io_convert_polynomials_to_ir(polynomials, options::KeywordArguments) in ./input_output/AbstractAlgebra.jl is the function that checks coeffient types and either allows or disallows the field type. E.g. (suppose p=13), it will allow finite_field(13) but disallow finite_field(13, 3, "a") which are of degree 1 and 3 respectively.

Initially I was going to volunteer to try to code this change myself. But after having looked at the code briefly, I came to the conclusion that it is not at all an easy fix, especially for someone coming in cold, who is not familiar with the code. I cloned the project and played around with the code a bit. It's not obvious what is going on in the functions io_convert_polynomials_to_ir and io_convert_ir_to_polynomials which I believe would be one of the necessary changes to incorporating finite_field of degree > 1. I assume that the naming "_irto" and "_toir" means converting to and from internal representation? It appears that maybe the internal representation of an element of a GF(p) field is just a UINT. And if that is the case then I suppose to represent an element of a finite field of of degree 3 say would require an array of length 3 of UINTS. I.e. elements of finite_field(p, 1) are basically 1-dimensional vectors (i.e. a scalar) and elements of finite_field(p, m) are m-dimensional vectors. Furthermore if there are matrices whose entries are coefficients of monomials, which are scalars for degree 1 fields, then that would be quite a bit trickier for coeffients that are m-dimensional vectors. On the otherhand, for you guys that implemented this package, it goes without saying that it would be much easier for you to add this feature than someone coming in cold. Let me know your thoughts.

sumiya11 commented 6 days ago

Thanks for raising this issue !

Implementing support for finite fields of degree > 1 is possible in principle, and I agree that this would be a nice feature.

So that we are on the same page, I mean support for something like:

using Nemo, Groebner
K, a = finite_field(7, 2, "a");
R, (x,y) = K["x","y"];
groebner([x + a*y]);

Implementation-wise, I think Julia generic types simplify our life:

https://github.com/sumiya11/Groebner.jl/blob/545ec1d4e92d5ce6e1c922a5c7a52472cf8b0202/src/f4/matrix.jl#L30

So, it seems feasible, but indeed would require some work.

I assume that the naming "ir_to" and "to_ir" means converting to and from internal representation?

You are close to the truth :^). Here "ir" stands for intermediate representation. The full pipeline is:

input polynomials -> ir -> internal -> F4 -> internal -> ir -> output polynomials.

(The functions names are io_convert_polynomials_to_ir -> ir_convert_ir_to_internal -> f4! -> ir_convert_internal_to_ir-> io_convert_ir_to_polynomials).

It appears that maybe the internal representation of an element of a GF(p) field is just a UINT. And if that is the case then I suppose to represent an element of a finite field of of degree 3 say would require an array of length 3 of UINTS.

Yes, I think this is one possible approach. Another possibly simpler option is to wrap elements of F_q^n in a struct (say inside of io_convert_polynomials_to_ir), forward field operations to existing implementation (which is say somewhere in Nemo or Oscar), and use this struct in F4:

struct MyFqCoeff
    actual_coeff
end

import Base: +, *
+(x::MyFqCoeff, y::MyFqCoeff) = MyFqCoeff(x.actual_coeff + y.actual_coeff)
*(x::MyFqCoeff, y::MyFqCoeff) = MyFqCoeff(x.actual_coeff * y.actual_coeff)
# etc

K, a = finite_field(7, 2, "a");
c1 = MyFqCoeff(K(5))
c1*c1

In this case, we need to teach io_convert_polynomials_to_ir and friends to recognize MyFqCoeff.

Furthermore if there are matrices whose entries are coefficients of monomials, which are scalars for degree 1 fields, then that would be quite a bit trickier for coeffients that are m-dimensional vectors.

By the way, Groebner.jl already works with vectors of scalars, though for a different purpose. Here, Coeff defines all supported coefficient types in F4. You can notice that it supports CompositeCoeffZp which is a tuple of several integers.

https://github.com/sumiya11/Groebner.jl/blob/545ec1d4e92d5ce6e1c922a5c7a52472cf8b0202/src/types.jl#L25

sumiya11 commented 6 days ago

@xyz333math In the next several weeks I won't be able to work on this. If you decide to give it a try, I could try to answer your questions.

Out of curiosity, what is the application where you use GB over F_q^n ?

xyz333math commented 6 days ago

support need: coefficients in algebraic extensions of finite fields (elements in L):

K, a = finite_field(p, m, "a")  # where p is a prime int and m is any positive int
Kx, x = polynomial_ring(K, "x")
L, b = finite_field(f(x), "b")  # where f(x) is an irreducible poly in Kx
# exmpl element in L:  3 + 2*a^2 - a*b^3

so L is a Nemo field and I believe the type is FqField. Although I think L here would be isomorphic to some finite_field(q, n, "c") I think it may have a different structure in julia code. If we need to convert elements in L to an intermediate form like an array of Int, then it might be necessary to programmatically derive what q, n are, depending on what f(x) is, for finite_field(q, n, "c"); it'd be nice if there is a way to avoid that. Can we simply let f4 matrix entries be elements of L directly (i.e. type FqField)? I think the key to this whole thing is to figure out what the intermediate and internal representations of FqField elements need to be.

Since L is a FqField, then elements of L should automatically implement all the necessary arithmetic operations, so we would add FqField to the Union of types for Coeff.

My application: I'm currently a math PhD student at UTA (Texas) doing research that requires usage of multivariate polynomials over finite fields.

I appreciate your offer to try to answer my questions. I may give the coding a try, I haven't quite decided yet. If I do, then there is no harm if I create a git branch right?
How about this:

  1. I'll create a branch
  2. add a unit test
  3. make code changes as far as I can get
  4. check it in and have you look at it at your convenience.

This has the potential for us to knock this out quickly.

sumiya11 commented 6 days ago

Thanks for the clarification.

It seemed quite interesting so I went ahead and here is my attempt: https://github.com/sumiya11/Groebner.jl/pull/163.

The following works now:

using Nemo, Groebner
K, a = Nemo.finite_field(3, 2, "a")
R, (X, Y) = K["X", "Y"]
groebner([a*X - Y, X*Y - 1]) == [Y^2 + 2*a, X + (2*a + 1)*Y]

Your example with multiple extensions should naturally work as well, but I did not test yet.

The approach is simple: when the coefficients are GF(p) or QQ, nothing is changed. For all other fields, wrap the coefficients in a custom CoeffGeneric struct and passes it to F4 directly. Whence, intermediate and internal representation are now allowed to contain CoeffGeneric coefficients.

Can we simply let f4 matrix entries be elements of L directly (i.e. type FqField)?

Right, I also decided to do it almost like this.

My application: I'm currently a math PhD student at UTA (Texas) doing research that requires usage of multivariate polynomials over finite fields.

That's great !

For reference:

https://github.com/sumiya11/Groebner.jl/pull/163/files#diff-bd89d3371814e68223de34429fdae327c61bb942a94b0e4e4e32567519e0ca55R45-R69

https://github.com/sumiya11/Groebner.jl/pull/163/files#diff-fced79f727a11041cb1e1f1ae828388223e615e68fc63e377f452f6a58415229R157-R191

https://github.com/sumiya11/Groebner.jl/pull/163/files#diff-8df34a15ef488c0828d958d846300f5f2aab3d5b68f526537bd8f5aa0e588307R28-R42

xyz333math commented 6 days ago

That is outstanding! Thanks for doing that. I tried it out and it works!

Below is the function I wrote to test it. Since I'm using Oscar (not a Groebner.jl dep) I have this function in a separate project. It certainly appears like it correctly computes the groebner basis. I noticed the default monomial ordering is Lex. The only issue I noticed: it prints out Warning messages: Warning: Groebner.jl does not have a native implementation for the given field

I'm very excited about this and I anticipate using it very heavily; so over time I could probably write a ton of unit/integrated/functional tests for you. Also I could compare results with Oscar (but Oscar doesn't support degree>1 finite fields yet which is why I tried Groebner.jl in the first place) Now in my case, I came from a heavy java dev background, and so I was annoyed by the clumsy default unit test toolset for julia, and thus I never really properly incorporated the julia way of doing unit testing into my projects. So you'd have to guide me on how you want me to add tests to Groebner.jl.

using Oscar
using Groebner

function testJLGroebnerFinFldCoef1()
    K, a = Oscar.finite_field(13, 3, "a")
    Ks, s = Oscar.polynomial_ring(K, "s")
    # we know that s^2 + 11 is irreducible in K[s].
    # we extend K to a splitting field for s^2 + 11
    L, b = finite_field(s^2 + 11, "b")
    PolRng, (x, y, z) = Oscar.polynomial_ring(L, ["x", "y", "z"])
    l1 = one(L)
    c1 = 2*a + b
    c2 = a + l1
    c3 = b+L(7)
    c4 = a*b + L(5)

    p = c1*x+y+z
    q = l1*x + c3*y^2 + c2*z^3
    r = l1*x^3 + c4*y^3 + l1*z^3

    ideal1 = [p,q,r]
    println("input ideal:")
    for g in ideal1
        println(g)
    end
    basis1 = groebner(ideal1, ordering=Groebner.DegLex())
    println("groebner basis (deglex monom ordering):")
    for g in basis1
        println(g)
    end
    basis2 = groebner(ideal1, ordering=Groebner.Lex())
    println("groebner basis (lex monom ordering):")
    for g in basis2
        println(g)
    end
end
sumiya11 commented 6 days ago

I noticed the default monomial ordering is Lex.

Right. The default in Oscar.polynomial_ring is :lex, so Groebner inherits it. To override, this should work: Oscar.polynomial_ring(..., internal_ordering=:degrevlex) (or groebner(..., ordering=...) as you did.

The only issue I noticed: it prints out Warning messages: Warning: Groebner.jl does not have a native implementation for the given field

For now the warning signals that generic coefficients are used. Groebner is a low-level package that targets performance, and, the generic fallback may be not as fast as one might expect in some cases.

Now in my case, I came from a heavy java dev background, and so I was annoyed by the clumsy default unit test toolset for julia, and thus I never really properly incorporated the julia way of doing unit testing into my projects. So you'd have to guide me on how you want me to add tests to Groebner.jl.

It would be amazing if you can contribute tests for towers of fields :). Most tests in Groebner are simple black box tests: test that groebner(sys) == res for some predefined sys, res. Once I merge that PR (in several days I think), you are welcome to open a PR and add test cases, say, here https://github.com/sumiya11/Groebner.jl/blob/5741d237aef69f2c44635e5d0b5fa843ccddaeda/test/input_output/Nemo.jl#L28

I am not too sure about adding Oscar dependency in tests (it is a bit large and does not work on Windows). We could use Nemo.

sumiya11 commented 6 days ago

But after having looked at the code briefly, I came to the conclusion that it is not at all an easy fix, especially for someone coming in cold, who is not familiar with the code. I cloned the project and played around with the code a bit. It's not obvious what is going on in the functions io_convert_polynomials_to_ir and io_convert_ir_to_polynomials which I believe would be one of the necessary changes to incorporating finite_field of degree > 1.

BTW: if you would have comments on how to improve the clarity of some code in Groebner please let me know.

sumiya11 commented 5 days ago

Most tests in Groebner are simple black box tests: test that groebner(sys) == res for some predefined sys, res. Once I merge that PR (in several days I think), you are welcome to open a PR and add test cases, say, here

But of course if you have other ideas on how to better test this, feel free to propose it.

xyz333math commented 4 days ago

In my opinion you should not add Oscar as a dependency for the reasons you mentioned above and: better to keep your package dependencies as small and lean as possible.
For finite_fields, I believe that Oscar uses Nemo implementation.

I will make an effort to add some tests for tower of fields (create a new testset "Nemo.jl, finite-fields"). While developing a test, is there a way to run a single test or a single testset without having to run the entire runtests.jl? Does TestSetExtensions support anything like that?

sumiya11 commented 4 days ago

While developing a test, is there a way to run a single test or a single testset without having to run the entire runtests.jl?

In Groebner.jl, no easy way I think. My workflow has been: write code - test in REPL - put in test file - run runtests.jl and go for a coffee.

In general, there is https://github.com/JuliaTesting/ReTest.jl and https://discourse.julialang.org/t/prerelease-of-new-testing-framework-and-test-run-ui-in-vs-code/86355/52, which can run tests individually, but I have not tried it.

sumiya11 commented 4 days ago

Does TestSetExtensions support anything like that?

In fact it was only used for @includetests. I replaced it with simple includes while I am at it.

sumiya11 commented 4 days ago

Merged #163.

xyz333math commented 1 day ago

Hello, I've got another similar item for you, to implement support for field extensions in general. E.g. the code below as you can see uses the rationals extended by the sqrt(2).

R, s = polynomial_ring(QQ, "s")
K, q = number_field(s^2 - 2, "q")
PolRng, (x, y, z) = polynomial_ring(K, ["x", "y", "z"])

idl = [
    x^2+y+z,
    x^2 + q*y^2 + (q+1)*z^3,
    x^3 + 3*y^3 + z^3
]

basis1 = groebner(idl, ordering=Groebner.DegLex())
println("groebner basis (deglex monom ordering):")
for g in basis1
    println(g)
end

running this script results in an error in function io_extract_ring(polynomials) in file Groebner.jl/src/input_output/AbstractAlgebra.jl throws an error here:

# Detect that K is either AbstractAlgebra.QQ or Nemo.QQ
base = AbstractAlgebra.base_ring(K)
if !(AbstractAlgebra.base_ring(base) == Union{})
        ground = :generic
end
sumiya11 commented 1 day ago

Could you try the master branch ?

xyz333math commented 17 hours ago

Ok it works! Sorry, I didn't realize you had made further changes.

sumiya11 commented 16 hours ago

Great ! Thanks for checking