RDF[] for libsingular, delegate ideal comparison and .divides() to singular

src/sage/libs/singular/decl.pxd

@@ -990,11 +990,34 @@ cdef extern from "singular/coeffs/coeffs.h":
 
     number *ndCopyMap(number *, const n_Procs_s* src,const n_Procs_s* dst)
 
+    ctypedef struct LongComplexInfo:
+        short float_len
+        short float_len2
+        const char* par_name
+
 cdef extern from "singular/coeffs/rmodulo2m.h":
 
     #init 2^m from a long
     number *nr2mMapZp(number *,const n_Procs_s* src,const n_Procs_s* dst)
 
+cdef extern from "singular/coeffs/shortfl.h":
+    """
+    static inline SI_FLOAT sage_nrFloat(number n) { // copy from shortfl.cc to allow inlining
+        SI_FLOAT f = 0;
+        memcpy(&f, &n, sizeof(f) < sizeof(n) ? sizeof(f) : sizeof(n));
+        return f;
+    }
+
+    static inline number sage_nrInit(SI_FLOAT f) {
+        number n = 0;
+        memcpy(&n, &f, sizeof(f) < sizeof(n) ? sizeof(f) : sizeof(n));
+        return n;
+    }
+    """
+    ctypedef double SI_FLOAT  # actually might be double or float
+    SI_FLOAT nrFloat "sage_nrFloat" (number *n)
+    number *sage_nrInit(SI_FLOAT)
+
 cdef extern from "singular/kernel/maps/gen_maps.h":
 
     # mapping from p in r1 by i2 to r2
@@ -1013,6 +1036,20 @@ cdef extern from "singular/polys/ext_fields/algext.h":
 
     nMapFunc naSetMap(const n_Procs_s* src, const n_Procs_s* dst)
 
+cdef extern from "singular/coeffs/mpr_complex.h":
+    cdef cppclass gmp_float:
+        gmp_float(double v)
+        double to_double "operator double"() const
+
+    cdef cppclass gmp_complex:
+        gmp_complex(double re, double im)
+        gmp_float real() const
+        gmp_float imag() const
+
+    char *floatToStr(const gmp_float & r, const unsigned int oprec)
+    char *complexToStr(gmp_complex & c, const unsigned int oprec, const n_Procs_s* src)
+
+
 cdef extern from "singular/coeffs/rmodulon.h":
 
     # init ZmodN from GMP

src/sage/libs/singular/ring.pyx

@@ -26,11 +26,10 @@ from sage.libs.singular.decl cimport rChangeCurrRing, rComplete, rDelete, idInit
 from sage.libs.singular.decl cimport omAlloc0, omStrDup, omAlloc
 from sage.libs.singular.decl cimport ringorder_dp, ringorder_Dp, ringorder_lp, ringorder_ip, ringorder_ds, ringorder_Ds, ringorder_ls, ringorder_M, ringorder_c, ringorder_C, ringorder_wp, ringorder_Wp, ringorder_ws, ringorder_Ws, ringorder_a, rRingOrder_t
 from sage.libs.singular.decl cimport prCopyR
-from sage.libs.singular.decl cimport n_unknown, n_algExt, n_transExt, n_Z, n_Zn,  n_Znm, n_Z2m
-from sage.libs.singular.decl cimport n_coeffType
+from sage.libs.singular.decl cimport n_unknown, n_R, n_algExt, n_transExt, n_long_C, n_Z, n_Zn, n_Znm, n_Z2m
+from sage.libs.singular.decl cimport n_coeffType, LongComplexInfo
 from sage.libs.singular.decl cimport rDefault, GFInfo, ZnmInfo, nInitChar, AlgExtInfo, TransExtInfo
 
-
 from sage.rings.integer cimport Integer
 from sage.rings.integer_ring cimport IntegerRing_class
 import sage.rings.abc
@@ -327,6 +326,7 @@ cdef ring *singular_ring_new(base_ring, n, names, term_order) except NULL:
     cdef AlgExtInfo extParam
     cdef TransExtInfo trextParam
     cdef n_coeffType _type = n_unknown
+    cdef LongComplexInfo info
 
     #cdef cfInitCharProc myfunctionptr;
 
@@ -510,6 +510,17 @@ cdef ring *singular_ring_new(base_ring, n, names, term_order) except NULL:
         _cf = nInitChar( n_Z, NULL) # integer coefficient ring
         _ring = rDefault (_cf, nvars, _names, nblcks, _order, _block0, _block1, _wvhdl)
 
+    elif isinstance(base_ring, sage.rings.abc.RealDoubleField):
+        _cf = nInitChar(n_R, NULL)
+        _ring = rDefault(_cf, nvars, _names, nblcks, _order, _block0, _block1, _wvhdl)
+
+    elif isinstance(base_ring, sage.rings.abc.ComplexDoubleField):
+        info.float_len = 15
+        info.float_len2 = 0
+        info.par_name = "I"
+        _cf = nInitChar(n_long_C, <void *>&info)
+        _ring = rDefault(_cf, nvars, _names, nblcks, _order, _block0, _block1, _wvhdl)
+
     elif isinstance(base_ring, sage.rings.abc.IntegerModRing):
 
         ch = base_ring.characteristic()

src/sage/libs/singular/singular.pyx

@@ -40,6 +40,10 @@ from sage.rings.finite_rings.finite_field_base import FiniteField
 from sage.rings.polynomial.polynomial_ring import PolynomialRing_field
 from sage.rings.fraction_field import FractionField_generic
 
+from sage.rings.real_double cimport RealDoubleElement
+from sage.rings.complex_double cimport ComplexDoubleElement
+import sage.rings.abc
+
 from sage.rings.finite_rings.finite_field_prime_modn import FiniteField_prime_modn
 from sage.rings.finite_rings.finite_field_givaro import FiniteField_givaro
 from sage.rings.finite_rings.finite_field_ntl_gf2e import FiniteField_ntl_gf2e
@@ -1599,6 +1603,9 @@ cdef object si2sa(number *n, ring *_ring, object base):
 
     An element of ``base``
     """
+    cdef gmp_float* f
+    cdef gmp_complex* c
+    cdef char* s
     if isinstance(base, FiniteField_prime_modn) and _ring.cf.type == n_Zp:
         return base(_ring.cf.cfInt(n, _ring.cf))
 
@@ -1629,6 +1636,17 @@ cdef object si2sa(number *n, ring *_ring, object base):
     elif isinstance(base, IntegerModRing_generic):
         return si2sa_ZZmod(n, _ring, base)
 
+    elif isinstance(base, sage.rings.abc.RealDoubleField) and _ring.cf.type == n_R:
+        return base(nrFloat(n))
+
+    elif isinstance(base, sage.rings.abc.RealDoubleField) and _ring.cf.type == n_long_R:
+        f = <gmp_float*>n
+        return base(f[0].to_double())
+
+    elif isinstance(base, sage.rings.abc.ComplexDoubleField) and _ring.cf.type == n_long_C:
+        c = <gmp_complex*>n
+        return base(c[0].real().to_double(), c[0].imag().to_double())
+
     else:
         raise ValueError("cannot convert from SINGULAR number")
 
@@ -1648,6 +1666,9 @@ cdef number *sa2si(Element elem, ring * _ring) noexcept:
     a (pointer to) a singular number
     """
     cdef int i = 0
+    cdef ComplexDoubleElement z
+    cdef RealDoubleElement re
+    cdef RealDoubleElement im
 
     if isinstance(elem._parent, FiniteField_prime_modn) and _ring.cf.type == n_Zp:
         return n_Init(int(elem),_ring.cf)
@@ -1671,6 +1692,15 @@ cdef number *sa2si(Element elem, ring * _ring) noexcept:
         return sa2si_NF(elem, _ring)
     elif isinstance(elem._parent, IntegerModRing_generic):
         return sa2si_ZZmod(elem, _ring)
+    elif isinstance(elem._parent, sage.rings.abc.RealDoubleField) and _ring.cf.type == n_R:
+        return sage_nrInit((<RealDoubleElement>elem)._value)
+    elif isinstance(elem._parent, sage.rings.abc.RealDoubleField) and _ring.cf.type == n_long_R:
+        return <number*><void*>new gmp_float((<RealDoubleElement>elem)._value)
+    elif isinstance(elem._parent, sage.rings.abc.ComplexDoubleField) and _ring.cf.type == n_long_C:
+        z = <ComplexDoubleElement>elem
+        re = <RealDoubleElement>z.real()
+        im = <RealDoubleElement>z.imag()
+        return <number*><void*>new gmp_complex(re._value, im._value)
     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)

src/sage/rings/ideal.py

@@ -367,24 +367,29 @@ def _richcmp_(self, other, op):
             sage: I == J
             True
 
-        TESTS:
-
-        Check that the example from :issue:`37409` raises an error
-        rather than silently returning an incorrect result::
+        TESTS::
 
             sage: R.<x> = ZZ[]
             sage: I = R.ideal(1-2*x,2)
-            sage: I.is_trivial()  # not implemented -- see #37409
+            sage: I.is_trivial()  # needs sage.libs.singular
+            True
+            sage: R.ideal(x, 2) <= R.ideal(1)  # needs sage.libs.singular
             True
-            sage: I.is_trivial()
-            Traceback (most recent call last):
-            ...
-            NotImplementedError: ideal comparison in Univariate Polynomial Ring in x over Integer Ring is not implemented
         """
         if self.is_zero():
             return rich_to_bool(op, other.is_zero() - 1)
         if other.is_zero():
             return rich_to_bool(op, 1)  # self.is_zero() is already False
+
+        from sage.rings.integer_ring import ZZ
+        if self.ring().base_ring() is ZZ:
+            #assert self.ring().implementation() != "singular"
+            from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
+            Rs = PolynomialRing(ZZ, self.ring().variable_names(), implementation="singular")
+            Is = Rs.ideal(self.gens())
+            Js = Rs.ideal(other.gens())
+            return Is.__richcmp__(Js, op)
+
         S = set(self.gens())
         T = set(other.gens())
         if S == T:

src/sage/rings/polynomial/multi_polynomial_libsingular.pyx

@@ -170,8 +170,6 @@ AUTHORS:
 #   * pNext and pIter don't need currRing
 #   * p_Normalize apparently needs currRing
 
-from warnings import warn
-
 from libc.limits cimport INT_MAX
 
 from cpython.object cimport Py_NE
@@ -4687,6 +4685,38 @@ cdef class MPolynomial_libsingular(MPolynomial_libsingular_base):
 
         A :exc:`ValueError` exception is raised if ``g (== self)`` does not belong to ``I``.
 
+        .. WARNING::
+
+            This method is unreliable when the base ring is inexact::
+
+                sage: R.<x> = PolynomialRing(RDF, implementation="singular")
+                sage: g = 5*x^3 + x - 7; m = x^4 - 12*x + 13; R(1).lift((g, m))
+                Traceback (most recent call last):
+                ...
+                ValueError: polynomial is not in the ideal
+                sage: g.inverse_mod(m)
+                Traceback (most recent call last):
+                ...
+                ArithmeticError: element is non-invertible
+                sage: inverse_mod(g, m)
+                Traceback (most recent call last):
+                ...
+                ArithmeticError: element is non-invertible
+
+            While the generic implementation is fine::
+
+                sage: R.<x> = PolynomialRing(RDF)
+                sage: g = 5*x^3 + x - 7; m = x^4 - 12*x + 13; inverse_mod(g, m)  # abs tol 1e-14
+                -0.03196361250430942*x^3 - 0.03832697590106863*x^2 - 0.046305090023464945*x + 0.3464796877258926
+                sage: g.inverse_mod(m)  # abs tol 1e-14
+                -0.03196361250430942*x^3 - 0.03832697590106863*x^2 - 0.046305090023464945*x + 0.3464796877258926
+
+            But this is because of Singular:
+
+            .. code-block:: bash
+
+                Singular -q -c 'ring r = real,(x),lp; poly f = 5*x3 + x - 7; poly g = x4 - 12*x + 13; ideal M = f,g; ideal G = std(M); ideal SM = f; matrix T = lift(G,SM); T; quit;'
+
         INPUT:
 
         - ``I`` -- an ideal in ``self.parent()`` or tuple of generators of that ideal
@@ -4803,9 +4833,10 @@ cdef class MPolynomial_libsingular(MPolynomial_libsingular_base):
             finally:
                 s = check_error()
             if s:
+                msg = "polynomial is not in the ideal"
                 if s != ('2nd module does not lie in the first',):
-                    warn(f'unexpected error from singular: {s}')
-                raise ValueError("polynomial is not in the ideal")
+                    msg += f'; unexpected error from singular: {s}'
+                raise ValueError(msg)
 
             l = [new_MP(parent, pTakeOutComp(&res.m[i], 1))
                  for i in range(IDELEMS(res)) for _ in range(IDELEMS(_I))]

src/sage/rings/polynomial/polynomial_element.pyx

@@ -1669,13 +1669,13 @@ cdef class Polynomial(CommutativePolynomial):
         """
         from sage.rings.polynomial.polynomial_singular_interface import can_convert_to_singular
         P = a.parent()
-        if can_convert_to_singular(P):
+        if P.is_exact() and can_convert_to_singular(P):
             from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
             try:
                 R = PolynomialRing(P.base_ring(), P.variable_names(), implementation="singular")
             except NotImplementedError:
-                # PolynomialRing over RDF/CDF etc. are still not implemented in libsingular
-                # (in particular singular_ring_new) even though they are implemented in _do_singular_init_
+                # singular(PolynomialRing(Frac(ZZ["u", "v", "w"]), "x")) works, but
+                # PolynomialRing(Frac(ZZ["u", "v", "w"]), "x", implementation="singular") fails
                 pass
             else:
                 return P(R(a).inverse_mod(R.ideal(m)))
@@ -11588,6 +11588,12 @@ cdef class Polynomial(CommutativePolynomial):
             sage: q.divides(p * q)
             True
             sage: R.<x> = Zmod(6)[]
+            sage: f = x - 2
+            sage: g = 3
+            sage: R.ideal(f*g) <= R.ideal(f)
+            True
+            sage: f.divides(f*g)
+            True
             sage: p = 4*x + 3
             sage: q = 5*x**2 + x + 2
             sage: q.divides(p)
@@ -11620,9 +11626,7 @@ cdef class Polynomial(CommutativePolynomial):
             sage: p = 4*x + 3
             sage: q = 2*x**2 + x + 2
             sage: p.divides(q)
-            Traceback (most recent call last):
-            ...
-            NotImplementedError: divisibility test only implemented for polynomials over an integral domain unless obvious non divisibility of leading terms
+            False
         """
         if p.is_zero():
             return True          # everything divides 0
@@ -11642,7 +11646,17 @@ cdef class Polynomial(CommutativePolynomial):
             return False
 
         if not self.base_ring().is_integral_domain():
-            raise NotImplementedError("divisibility test only implemented for polynomials over an integral domain unless obvious non divisibility of leading terms")
+            from sage.rings.polynomial.polynomial_singular_interface import can_convert_to_singular
+            P = self.parent()
+            if P.is_exact() and can_convert_to_singular(P):
+                from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
+                try:
+                    R = PolynomialRing(P.base_ring(), P.variable_names(), implementation="singular")
+                except NotImplementedError:
+                    pass
+                else:
+                    return bool(R(self).divides(R(p)))
+            raise NotImplementedError("divisibility test only implemented for polynomials over an integral domain unless Singular can be used")
 
         try:
             return (p % self).is_zero()      # if quo_rem is defined

src/sage/rings/polynomial/polynomial_ring_constructor.py

@@ -488,6 +488,51 @@ def PolynomialRing(base_ring, *args, **kwds):
         ...
         NotImplementedError: polynomials over Real Field with 53 bits of precision are not supported in Singular
 
+    Polynomial rings over double-precision real and complex fields are supported
+    by the ``'singular'`` implementation::
+
+        sage: # needs sage.libs.singular
+        sage: R.<x> = PolynomialRing(RDF, implementation='singular')
+        sage: x^2 + 1
+        x^2 + (1.000e + 00)
+        sage: C.<x> = PolynomialRing(CDF, implementation='singular')
+        sage: x^2 + 1
+        x^2 + 1
+
+    The ring printing via libsingular matches the output of the Singular
+    interface::
+
+        sage: # needs sage.libs.singular
+        sage: from sage.libs.singular.function import singular_function
+        sage: print(singular_function("print")(R))
+        polynomial ring, over a field, global ordering
+        // coefficients: Float() considered as a field
+        // number of vars : 1
+        //        block   1 : ordering dp
+        //                  : names    x
+        //        block   2 : ordering C
+        sage: print(singular_function("print")(C))
+        polynomial ring, over a field, global ordering
+        // coefficients: real[I](complex:15 digits, additional 0 digits)/(I^2+1) considered as a field
+        // number of vars : 1
+        //        block   1 : ordering dp
+        //                  : names    x
+        //        block   2 : ordering C
+        sage: print(singular(R))
+        polynomial ring, over a field, global ordering
+        // coefficients: Float() considered as a field
+        // number of vars : 1
+        //        block   1 : ordering dp
+        //                  : names    x
+        //        block   2 : ordering C
+        sage: print(singular(C))
+        polynomial ring, over a field, global ordering
+        // coefficients: real[I](complex:15 digits, additional 0 digits)/(I^2+1) considered as a field
+        // number of vars : 1
+        //        block   1 : ordering dp
+        //                  : names    x
+        //        block   2 : ordering C
+
     The following corner case used to result in a warning message from
     ``libSingular``, and the generators of the resulting polynomial
     ring were not zero::
@@ -854,6 +899,9 @@ def _multi_variate(base_ring, names, sparse=None, order='degrevlex', implementat
     # yield the same implementation. We need this for caching.
     implementation_names = set([implementation])
 
+    if implementation is None and isinstance(base_ring, (sage.rings.abc.RealDoubleField, sage.rings.abc.ComplexDoubleField)):
+        implementation = "generic"  # singular has some issues with RDF/CDF, do not make singular the default
+
     if implementation is None or implementation == "singular":
         try:
             from sage.rings.polynomial.multi_polynomial_libsingular import MPolynomialRing_libsingular