Use libsingular for polynomial rings over function fields

[Simon King]

Singular has polynomial rings over function fields. However, Sage makes no use of it:

sage: P.<a> = QQ[]
sage: F = P.fraction_field()
sage: type(F['x','y'])
<class 'sage.rings.polynomial.multi_polynomial_ring.MPolynomialRing_polydict_domain_with_category'>
sage: F = FunctionField(QQ,'a')
sage: type(F['x','y'])
<class 'sage.rings.polynomial.multi_polynomial_ring.MPolynomialRing_polydict_domain_with_category'>

Update:

This works now:

sage: F = PolynomialRing(QQ,'a').fraction_field()
sage: R.<x,y> = F[]
sage: F.inject_variables()
Defining a
sage: I = R.ideal([a*x+5*y^2, (1+a)/(1-a)*x^3-3*y*x])
sage: I.groebner_basis()
[x^3 + (3*a - 3)/(a + 1)*x*y, y^2 + (a)/5*x]

}}}

[Simon King]

It is not a big deal to create the MPolynomialRing_libsingular with parameters. However, I need to find out how to create a Singular number with parameters. Martin, can you give me a pointer to a relevant file of the Singular sources?

[Martin Albrecht]

Looks like you want kernel/longtrans.h:

/*
* ABSTRACT:   numbers in transcendental field extensions,
              i.e., in rational function fields
*/

[Simon King]

Replying to malb:

Looks like you want kernel/longtrans.h:

/*
* ABSTRACT:   numbers in transcendental field extensions,
              i.e., in rational function fields
*/

Thank you, Martin!

So far, I only found number (*nPar)(int i). If I understand correctly, it returns the i-th parameter as a number. But this would probably be cumbersome to use for transforming an element of a Sage function field resp. of a Sage quotient field of polynomial ring.

[Simon King]

Seems like one actually needs to construct a number by starting with nPar(i) (which is a function pointer to ntPar, IIUC) and conversion of integers, and then add, multiply and divide. Or is there a way to directly convert a polynomial in variables a,b,c into a number in parameters a,b,c? After all, there is this globally accessible ring nacRing hosting numerator and denominator of a function field element.

[Martin Albrecht]

I'd assume there should be a way to convert a polynomial to a number as these numbers are polynomials internally as far as I know. Maybe ask on libsingular-devel?

[Simon King]

Miguel asked me to push the basic steps that I did for this ticket.

All what works is to set up polynomial rings over "fraction fields of polynomial rings" as libsingular polynomial rings with parameters. What is missing is everything else. In particular: the element constructor must be modified so that elements of quotient fields of polynomial rings are converted into libsingular numbers.

---

New commits:

​d15c0cc

Use libsingular for polynomial ring with parameters.

[Miguel Marco]

Ping!

I am trying (once again) to work on this. I also got to the point of creating the ring, and am stuck in the point of translating the coefficients.

I have started with conversion from sage coefficients to singular ones (the sa2si function in singular.pyx; but for some reason i don't understand, the generated .c code does not import the definitions from transext.h.

[Miguel Marco]

Also, i am not sure if i am dealing with the pointers correctly. It could use some review freom someone with more c/cython experience.

[Nils Bruin]

It's a little easier if the branch you're talking about is actually attached.

The cython docs for "cimport extern" is here:

​https://cython.readthedocs.io/en/latest/src/userguide/external_C_code.html#referencing-c-header-files

It actually does say there that the extern declaration should let cython generate an #include "singular/polys/ext_fields/transext.h". However, cython with take it on blind fate that the cython declarations given in the block itself somehow correspond to something that makes sense in C, because cython will not generate any declarations to define the types you specify: the "extern" promises that this is somehow resolved (by the given include or something else). So you should probably check that the generated C file indeed has an #include included. Otherwise, if you made a subtle typo compared with what's in C, the error would only arise in the C compiler, so you should probably carefully read the error that the C compiler generates and check that your spelling of everything matches up. Hope this helps?

[Miguel Marco]

Thanks for the answer!

I did check that: the generated .cpp file has this line

#include "singular/polys/ext_fields/transext.h"

and the file SAGE_LOCAL/local/include/singular/polys/ext_fields/transext.h has the following declarations:

struct fractionObject
{
  poly numerator;
  poly denominator;
  int complexity;
};

typedef struct fractionObject * fraction;

So I still have no idea about what could be going on.

One possible guess is that the generated file is a .cpp, that is, no pure C, but C++. Maybe some problem with scopes?

[Miguel Marco]

I think I finally got it. This branch seems to work.

Still needs some polishing (documentation, cleanup, and check that we are handling memory correctly... which means lots of testing).

So please be welcome to take a look and stress-test it.

[Miguel Marco]

New commits:

​3de10e6	Hacky solution for sa2si (coefficients that fit in int size only)
​2d89757	Use mpz to create coefficients
​64f4f57	Working implementation. Still needs cleanup, documentation, and checking for memory leaks
​3d3419f	Fixed case with just one parameter
​136a89a	Fix problem with several parameters introduced by previous fix

[Miguel Marco]

I found a subtle bug:

sage: F = PolynomialRing(QQ,'t').fraction_field()
sage: FA = FreeAlgebra(F, 6, 'x1,x2,x3,x4,x5,x6')
sage: N = FA.g_algebra({FA.gen(j)*FA.gen(i):-FA.gen(i)*FA.gen(j) for i in range(5) for j in range(i+1,6)})
sage: I = N.ideal([g^2 for g in N.gens()],side='twosided')
sage: N.inject_variables()
Defining x1, x2, x3, x4, x5, x6
sage: I.reduce(x1*x2*x3 + x2^2*x4)
x1*x2*x3 + x2^2*x4

Note that it doesn't cancel out the square of x2. It is very subtle, since:

sage: I.reduce(x1*x2*x3)
x1*x2*x3
sage: I.reduce(x2^2*x4)
0

that is, the summands are correctly reduced. Moreover

sage: I.reduce(x3*x2*x6 + x2^2*x4)
-x2*x3*x6
sage: I.reduce(x1*x2*x4 + x2^2*x4)
x1*x2*x4 + x2^2*x4
sage: I.reduce(x1*x2*x5 + x2^2*x4)
x1*x2*x5
sage: I.reduce(x1*x2*x6 + x2^2*x4)
x1*x2*x6
sage: I.reduce(x3*x2*x6 + x2^2*x4)
-x2*x3*x6
sage: I.reduce(x3*x2*x5 + x2^2*x4)
-x2*x3*x5

The bug only appears with a very precise set of possibilities for the summands.

It feels that it would be very hard to investigate this bug, but it doesn't appear in the previously affected fields:

sage: FA = FreeAlgebra(QQ, 6, 'x1,x2,x3,x4,x5,x6')
sage: N = FA.g_algebra({FA.gen(j)*FA.gen(i):-FA.gen(i)*FA.gen(j) for i in range(5) for j in range(i+1,6)})
sage: I = N.ideal([g^2 for g in N.gens()],side='twosided')
sage: N.inject_variables()
Defining x1, x2, x3, x4, x5, x6
sage: I.reduce(x1*x2*x3 + x2^2*x4)
x1*x2*x3

So this code doesn't really introduce a bug in previous functionalities.

---

New commits:

​6271a79	Hack in plural.pyx to construct elements from numbers (through the base field)
​2989430	Make coercion to base ring in plural less agresive

[Miguel Marco]

Replying to mmarco:

A hint might be that it doesn't happen whoth more that one parameter. So it might be related to how we treat that case (it is treated differently since by default univariate polynomial rings in sage are not wrapped around singular).

[Miguel Marco]

Replying to mmarco:

Replying to mmarco:

A hint might be that it doesn't happen whoth more that one parameter. So it might be related to how we treat that case (it is treated differently since by default univariate polynomial rings in sage are not wrapped around singular).

Sorry, this is not correct. The bug appears with several parameters as well.

[Dima Pasechnik]

We saw this "not reducing x²" thing in #27508 - a ticket that shows that we don't fully understand how linsingular works, and the only documentation on it is basically source code, or whatever you can extract from the upstream author...

[Miguel Marco]

Thanks that was actually very helpful!

The solution to #27508 (to force noTailReduction = False) was only applied to GroebnerStrategy. We just have to apply it to NCGroebnerStrategy as well

---

New commits:

​50ac936

Force tail reduction in NCGroebnerStrategy

[Miguel Marco]

I am doing the final touches before it is ready, but running the whole test suite I got the following error:

sage -t --warn-long 67.8 --random-seed=0 src/sage/schemes/generic/morphism.py
**********************************************************************
File "src/sage/schemes/generic/morphism.py", line 1605, in sage.schemes.generic.morphism.SchemeMorphism_polynomial.specialization
Failed example:
    f.specialization({alpha:5,beta:10})
Exception raised:
    Traceback (most recent call last):
      File "/home/mmarco/sage/local/lib/python3.9/site-packages/sage/doctest/forker.py", line 694, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/home/mmarco/sage/local/lib/python3.9/site-packages/sage/doctest/forker.py", line 1088, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.schemes.generic.morphism.SchemeMorphism_polynomial.specialization[32]>", line 1, in <module>
        f.specialization({alpha:Integer(5),beta:Integer(10)})
      File "/home/mmarco/sage/local/lib/python3.9/site-packages/sage/dynamics/arithmetic_dynamics/generic_ds.py", line 334, in specialization
        F = self.as_scheme_morphism().specialization(D, phi, homset)
      File "/home/mmarco/sage/local/lib/python3.9/site-packages/sage/schemes/generic/morphism.py", line 1618, in specialization
        phi = FractionSpecializationMorphism(self[0].parent(), D)
      File "/home/mmarco/sage/local/lib/python3.9/site-packages/sage/rings/polynomial/flatten.py", line 663, in __init__
        self._specialization = SpecializationMorphism(domain.base(), D)
      File "/home/mmarco/sage/local/lib/python3.9/site-packages/sage/rings/polynomial/flatten.py", line 540, in __init__
        raise NameError("argument " + str(var) + " is not a generator anywhere in the polynomial tower")
    NameError: argument (alpha) is not a generator anywhere in the polynomial tower
**********************************************************************
File "src/sage/schemes/generic/morphism.py", line 1609, in sage.schemes.generic.morphism.SchemeMorphism_polynomial.specialization
Failed example:
    f.specialization({alpha:5}).specialization({beta:10}) == f.specialization({alpha:5,beta:10})
Exception raised:
    Traceback (most recent call last):
      File "/home/mmarco/sage/local/lib/python3.9/site-packages/sage/doctest/forker.py", line 694, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/home/mmarco/sage/local/lib/python3.9/site-packages/sage/doctest/forker.py", line 1088, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.schemes.generic.morphism.SchemeMorphism_polynomial.specialization[33]>", line 1, in <module>
        f.specialization({alpha:Integer(5)}).specialization({beta:Integer(10)}) == f.specialization({alpha:Integer(5),beta:Integer(10)})
      File "/home/mmarco/sage/local/lib/python3.9/site-packages/sage/dynamics/arithmetic_dynamics/generic_ds.py", line 334, in specialization
        F = self.as_scheme_morphism().specialization(D, phi, homset)
      File "/home/mmarco/sage/local/lib/python3.9/site-packages/sage/schemes/generic/morphism.py", line 1618, in specialization
        phi = FractionSpecializationMorphism(self[0].parent(), D)
      File "/home/mmarco/sage/local/lib/python3.9/site-packages/sage/rings/polynomial/flatten.py", line 663, in __init__
        self._specialization = SpecializationMorphism(domain.base(), D)
      File "/home/mmarco/sage/local/lib/python3.9/site-packages/sage/rings/polynomial/flatten.py", line 540, in __init__
        raise NameError("argument " + str(var) + " is not a generator anywhere in the polynomial tower")
    NameError: argument (alpha) is not a generator anywhere in the polynomial tower
**********************************************************************
1 item had failures:
   2 of  35 in sage.schemes.generic.morphism.SchemeMorphism_polynomial.specialization
    [476 tests, 2 failures, 1.00 s]

Which actually comes from the behaviour of FlatteningMorphism? in src/sage/rings/polynomial/flatten.py

For some reason, it did work with the previous (pexpect) singular interface, but breaks with this branch. I will try to sort it out, but I am not sure it does make any mathematical sense to have this kind of morphisms.

I mean: If I understand correctly, the flattening morphism is suposed to be constructed from, say the ring QQ(a,b)[x,y] and the dictionary {a:3}, and will return a "morphism" from QQ(a,b)[x,y] to QQ(b)[x,y], given by substituting a=3.

But that is not defined for all elements of QQ(a,b)[x,y]:

sage: from sage.rings.polynomial.flatten import SpecializationMorphism
sage: F = PolynomialRing(QQ,'a,b').fraction_field()
sage: F.inject_variables()
Defining a, b
R.<x,y> = F[]
sage: phi = SpecializationMorphism(R, {a:3})
sage: phi
phi = SpecializationMorphism(R, {a:3})

phi

Specialization morphism:
  From: Multivariate Polynomial Ring in x, y over Fraction Field of Multivariate Polynomial Ring in a, b over Rational Field
  To:   Multivariate Polynomial Ring in x, y over Fraction Field of Univariate Polynomial Ring in b over Rational Field
sage: el = (1/(3-a))*x
sage: el.parent()
Multivariate Polynomial Ring in x, y over Fraction Field of Multivariate Polynomial Ring in a, b over Rational Field
sage: phi(el)    #   booom!!
---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
...
ZeroDivisionError: fraction field element division by zero

So, do you think this is ok to have (even if it is not really a morphism), or should we just eliminate it?

---

New commits:

​7ba8fc5	Doctest the tail reduction fix
​328dcc7	Make sure we are in the correct libsingular ring when creating the object
​57d9d7d	Cleanup variables
​e57f349	Make sure we free memory during coefficient creation

[Miguel Marco]

Regardles of what we do with the previous "morphisms", this last patch solves the underlying issues (related to the fact that the same variable was represented differently depending on weather it was considered in the base field or in the polynomial ring).

All doctests pass in my machine, so I think it is ready for review.

---

New commits:

​d2e4ff7

Make sure single variables are represented consistently, and fix doctests broken because of the change in representation

[Ben Hutz]

Weighing on the goals of the Flattening-Specialization framework.

The idea is that given a family of objects under some parameterization. It is incredibly useful to have a way to specialize to a specific member of that family (by choosing a value of one of the parameters). This is obtained in the .specialization() functions available for various objects. It is not at all surprising the you can choose parameter values for which the specialized member is not defined. The error returned to the user should be understandable in this case.

To implement this well, a number of structures needed to be created. The Flatenning structures take a stacked polynomial ring tower such as K[a,b][c,d] and create the polynomial ring K[a,b,c,d]. A value can then be substituted in for one of the variables as a specialization and the tower rebuilt with Unflattening.

So Flattening doesn't go from QQ[a,b][x,y] to QQ(b)[x,y] as mentioned above. Rather it just goes from Q[a,b][x,y] to QQ[a,b,x,y]. (How function fields are dealt with in specialization is a little more complicated, but essentially flattening the poly ring and specializing the numerator/denominator, then rebuilding correctly).

Looks like this issues mentioned were already resolved.

---

New commits:

​58f528c

Clean unnecessary parenthesis in single variables.

[Miguel Marco]

Yes, I found that the reason that the previous fail was that the flattening process depended on the string representation of variables, and libsingular represents parameters in polynomial rinms with parenthesis. Changing that representation fixed that problem.

There was another problem caused by the fact that univariate polynomials are usually a sifferent class than multivariate ones. To make use of libsingular over rings with one parameter, a "multivariate" polynomial ring with just one variable has to be created. Then some subtlety in the way morphisms work ended giving the fraction field element constructor a list as an input, instead of an element. Adding that option to the element constructor fixed that second problem.

I see the usefulness of the specialization morphisms. My concerns is more a conceptual one, about them being no "real" morphisms in the mathematical sense when fraction fields are involved, since a morphism should be defined by all values of the domain. But it might be that the usefulness overwheights the mathematical inacuracy

[Nils Bruin]

Replying to mmarco:

I see the usefulness of the specialization morphisms. My concerns is more a conceptual one, about them being no "real" morphisms in the mathematical sense when fraction fields are involved, since a morphism should be defined by all values of the domain. But it might be that the usefulness overwheights the mathematical inacuracy

The domain of definition for such specializations would be some relevant local ring. I think it will be easier to handle these using a partial homomorphism as is done now. Perhaps a doctest and a comment on that these maps may be partial.

[Miguel Marco]

Ok then. Maybe that comment could be handled in a different ticket (maybe we could even try to implement localization rings)?.

This one already has plenty of code to review.

[Miguel Marco]

New commits:

​b45fc7e	Clean intermediate coefficients too
​4a02d25	Improve documentation

[Miguel Marco]

New commits:

​13e30ad	Implemented transcendental extensions over prime fields
​cb3070f	Add documentation for prime fields
​5555f49	Fix doctests for prime fields
​ecbcaf9	Use cfPower instead of repeated multiplications

[Miguel Marco]

I found that there was a branch problem with the branch name starting with a slash. I changed it and now git trac seems to work ok.

Sorry to those that had problems trying to check it out.

[Travis Scrimshaw]

Some minor things:

• Not that it matters too much since they are cdef functions, but in si2sa_QQ and similar methods, the INPUT: and OUTPUT: blocks should not be indented, and it should be

    - ``_ ring`` -- a (pointer to) a singular ring, in whose coefficient field
      lives ``n``
  

• It would be nice to get rid of some of the blanklines between blocks in the documentation. For example, in sa2si_transext_FF.
• These parentheses are extraneous: if (nMapFuncPtr is NULL):.
• Is it possible to combine these if statements?

  +    elif (isinstance(elem._parent, FractionField_generic)
  +          and isinstance(elem._parent.base(), (MPolynomialRing_libsingular, PolynomialRing_field))
  +        if isinstance(elem._parent.base().base_ring(), RationalField):
               return sa2si_transext_QQ(elem, _ring)
  +        elif isinstance(elem._parent.base().base_ring(), FiniteField_prime_modn):
  +            return sa2si_transext_FF(elem, _ring)
  -    elif isinstance(elem._parent, FractionField_generic) and isinstance(elem._parent.base(), (MPolynomialRing_libsingular, PolynomialRing_field)) and isinstance(elem._parent.base().base_ring(), RationalField):
  -        return sa2si_transext_QQ(elem, _ring)
  -
  -    elif isinstance(elem._parent, FractionField_generic) and isinstance(elem._parent.base(), (MPolynomialRing_libsingular, PolynomialRing_field)) and isinstance(elem._parent.base().base_ring(), FiniteField_prime_modn):
  -        return sa2si_transext_FF(elem, _ring)
  

  They are checking all but the last condition, and I don't think a failure of the last one. Then remove the else: (but keeping the raise ValueError of course). This makes it more clear about how the code should work.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​9424a24

Some minor reviewer's changes

[Miguel Marco]

Thanks for the review.

[Travis Scrimshaw]

Setting the author so the patchbot will run.

[Miguel Marco]

Linter complains about an unused variable in a method for matrix groups, which is rewriten to use this interface in #19391.

Anyways, in the current implementation, that assignation is necesary to create the corresponding singular ring, even if the python variable is not used later. The linter is just not smart enough in this case.

[Travis Scrimshaw]

There are some failures reported by the patchbot:

sage -t --long --random-seed=65607106580737169923545380617831998347 src/sage/dynamics/arithmetic_dynamics/projective_ds.py  # 1 doctest failed
sage -t --long --random-seed=65607106580737169923545380617831998347 src/sage/combinat/root_system/non_symmetric_macdonald_polynomials.py  # 1 doctest failed
sage -t --long --random-seed=65607106580737169923545380617831998347 src/sage/interfaces/expect.py  # 2 doctests failed

[Miguel Marco]

I can't reproduce the error in pexpect or projective_ds.py (I suspect somehow the file i am testing is different from the one in the patchbot).

About the error in macdonald polynomials, I found that it is due to one function creating polynomials in x0, x1... and the other on x1, x2 ... Will look into it.

[Miguel Marco]

By rebasing to the newest develop version, I can reproduce the error in `src/sage/dynamics/arithmetic_dynamics/projective_ds.py. It has to do with code added in #31994.

I am trying to debug it, but it is quite complicated. Adding the people involved in that ticket to cc for help.

[Alexander Galarraga]

I've traced the error to the affine_patch function in src/sage/schemes/projective/projective_subscheme.py. The error is caused by code equivalent to the following

sage: T.<c,d> = QQ[]
sage: F = FractionField(T)
sage: R.<x,y,z> = F[]
sage: Q = F['x', 'z']
sage: f = R(d*z^2 + c*y*z^2)
sage: f((Q.gens()[0], 1, Q.gens()[1]))
TypeError

My usual workaround here is to create a dictionary and attempt f.subs(), however,

f.subs({x: Q.gens()[0], y: 1, z:Q.gens()[1]})

crashes Sage.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​1a94d13	Merge branch 'develop' into plural
​2704119	Fix creating polynomial from strings in the case of field with several generators
​84a04d1	Fix variable indices in doctest
​d412c77	Fix crash in .subs()

[Miguel Marco]

Thanks for the hint. I could pinpoint the origin of both problems. This should solve it.

I still can't reproduce the error in expect.py. Can somebody confirm it happens?

[Miguel Marco]

One of the patchbot passes all tests, and complains only abut the unused variable (that is actually needed to prevent gargage collection). The other patchbots have what seem to be unrelated problems.

Is it ok to merge this in this state? or should I do some dummy operation with the variable to make pyflakes happy?

[Travis Scrimshaw]

Then it is a spurious error. I wouldn't worry about the pyflakes error. Thank you.

[Volker Braun]

On 32-bit:

**********************************************************************
File "src/sage/rings/polynomial/polydict.pyx", line 1628, in sage.rings.polynomial.polydict.ETuple.eadd
Failed example:
    y^(2^32)
Expected:
    Traceback (most recent call last):
    ...
    OverflowError: exponent overflow (...)
Got:
    <BLANKLINE>
    Traceback (most recent call last):
      File "/var/lib/buildbot/slave/sage_git/build/local/var/lib/sage/venv-python3.9.7/lib/python3.9/site-packages/sage/doctest/forker.py", line 694, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/var/lib/buildbot/slave/sage_git/build/local/var/lib/sage/venv-python3.9.7/lib/python3.9/site-packages/sage/doctest/forker.py", line 1088, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.rings.polynomial.polydict.ETuple.eadd[6]>", line 1, in <module>
        y**(Integer(2)**Integer(32))
      File "sage/rings/polynomial/multi_polynomial_libsingular.pyx", line 2467, in sage.rings.polynomial.multi_polynomial_libsingular.MPolynomial_libsingular.__pow__ (build/cythonized/sage/rings/polynomial/multi_polynomial_libsingular.cpp:23166)
        singular_polynomial_pow(&_p, self._poly, exp, _ring)
    OverflowError: Python int too large to convert to C unsigned long
**********************************************************************
1 item had failures:
   1 of   8 in sage.rings.polynomial.polydict.ETuple.eadd
    [276 tests, 1 failure, 0.26 s]
----------------------------------------------------------------------
sage -t --long --random-seed=307417638528231788267817062701073757912 src/sage/rings/polynomial/polydict.pyx  # 1 doctest failed
----------------------------------------------------------------------

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​6ce2bed

Adapt overflow doctest to 32 bits

[Travis Scrimshaw]

LGTM.