RDF[] for libsingular, delegate ideal comparison to singular

Author: user202729

src/sage/libs/singular/ring.pyx

@@ -26,10 +26,16 @@ 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_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
 from sage.libs.singular.decl cimport rDefault, GFInfo, ZnmInfo, nInitChar, AlgExtInfo, TransExtInfo
 
+cdef extern from "coeffs/coeffs.h":
+    ctypedef struct LongComplexInfo:
+        short float_len
+        short float_len2
+        const char* par_name
+
 
 from sage.rings.integer cimport Integer
 from sage.rings.integer_ring cimport IntegerRing_class
@@ -327,6 +333,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 +517,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,16 +40,37 @@ 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
 from sage.libs.gmp.all cimport *
 
+cdef extern from "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)
+
 from sage.cpython.string import FS_ENCODING
 from sage.cpython.string cimport str_to_bytes, char_to_str, bytes_to_str
 
 from sage.rings.polynomial.multi_polynomial_libsingular cimport MPolynomialRing_libsingular
 
+ctypedef union _nf_real:
+    double f
+    number* n
+
 ctypedef struct fraction "fractionObject":
     poly *numerator
     poly *denominator
@@ -1599,6 +1620,10 @@ cdef object si2sa(number *n, ring *_ring, object base):
 
     An element of ``base``
     """
+    cdef gmp_float* f
+    cdef gmp_complex* c
+    cdef _nf_real u
+    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 +1654,18 @@ 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:
+        u.n = n
+        return base(u.f)
+
+    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 +1685,10 @@ cdef number *sa2si(Element elem, ring * _ring) noexcept:
     a (pointer to) a singular number
     """
     cdef int i = 0
+    cdef _nf_real u
+    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 +1712,16 @@ 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:
+        u.f = (<RealDoubleElement>elem)._value
+        return u.n
+    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,22 @@ 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))
+                sage: g.inverse_mod(m)
+                sage: inverse_mod(g, m)
+                sage: R.<x> = PolynomialRing(RDF)
+                sage: g = 5*x^3 + x - 7; m = x^4 - 12*x + 13; inverse_mod(g, m)
+                sage: g.inverse_mod(m)
+
+            But this is because of Singular::
+
+                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 +4817,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,16 +1669,10 @@ 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_
-                pass
-            else:
-                return P(R(a).inverse_mod(R.ideal(m)))
+            R = PolynomialRing(P.base_ring(), P.variable_names(), implementation="singular")
+            return P(R(a).inverse_mod(R.ideal(m)))
         from sage.rings.ideal import Ideal_generic
         if isinstance(m, Ideal_generic):
             v = m.gens_reduced()
@@ -11588,6 +11582,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 +11620,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 +11640,13 @@ 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 can_convert_to_singular(P):
+                from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
+                R = PolynomialRing(P.base_ring(), P.variable_names(), implementation="singular")
+                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,52 @@ 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: R.<x> = PolynomialRing(RDF, implementation='singular')                    # needs sage.libs.singular
+        sage: x^2 + 1
+        x^2 + (1.000e + 00)
+
+        sage: S.<x> = PolynomialRing(CDF, implementation='singular')                    # needs sage.libs.singular
+        sage: x^2 + 1
+        x^2 + 1
+
+    The ring printing via libSingular matches the output of the Singular
+    interface::
+
+        sage: from sage.libs.singular.function import singular_function               # needs sage.libs.singular
+        sage: R.<x> = PolynomialRing(RDF, implementation="singular")                  # needs sage.libs.singular
+        sage: C.<x> = PolynomialRing(CDF, implementation="singular")                  # needs sage.libs.singular
+        sage: print(singular_function("print")(R))                                    # needs sage.libs.singular
+        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))                                    # needs sage.libs.singular
+        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))                                                      # needs sage.libs.singular
+        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))                                                      # needs sage.libs.singular
+        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::