diff --git a/src/sage/libs/singular/decl.pxd b/src/sage/libs/singular/decl.pxd index 3b4a2581f85..076a884756f 100644 --- a/src/sage/libs/singular/decl.pxd +++ b/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 diff --git a/src/sage/libs/singular/ring.pyx b/src/sage/libs/singular/ring.pyx index eae7c1a7549..501f6b53647 100644 --- a/src/sage/libs/singular/ring.pyx +++ b/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, &info) + _ring = rDefault(_cf, nvars, _names, nblcks, _order, _block0, _block1, _wvhdl) + elif isinstance(base_ring, sage.rings.abc.IntegerModRing): ch = base_ring.characteristic() diff --git a/src/sage/libs/singular/singular.pyx b/src/sage/libs/singular/singular.pyx index 3665030f014..4354c0359a7 100644 --- a/src/sage/libs/singular/singular.pyx +++ b/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 = n + return base(f[0].to_double()) + + elif isinstance(base, sage.rings.abc.ComplexDoubleField) and _ring.cf.type == n_long_C: + c = 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((elem)._value) + elif isinstance(elem._parent, sage.rings.abc.RealDoubleField) and _ring.cf.type == n_long_R: + return new gmp_float((elem)._value) + elif isinstance(elem._parent, sage.rings.abc.ComplexDoubleField) and _ring.cf.type == n_long_C: + z = elem + re = z.real() + im = z.imag() + return 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) diff --git a/src/sage/rings/ideal.py b/src/sage/rings/ideal.py index eacdda6dd2f..752ad05c35f 100644 --- a/src/sage/rings/ideal.py +++ b/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. = 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: diff --git a/src/sage/rings/polynomial/multi_polynomial_libsingular.pyx b/src/sage/rings/polynomial/multi_polynomial_libsingular.pyx index 2ace29d079f..86dbcfd02aa 100644 --- a/src/sage/rings/polynomial/multi_polynomial_libsingular.pyx +++ b/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. = 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. = 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))] diff --git a/src/sage/rings/polynomial/polynomial_element.pyx b/src/sage/rings/polynomial/polynomial_element.pyx index 917dece9552..d040cc19858 100644 --- a/src/sage/rings/polynomial/polynomial_element.pyx +++ b/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. = 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 diff --git a/src/sage/rings/polynomial/polynomial_ring_constructor.py b/src/sage/rings/polynomial/polynomial_ring_constructor.py index aab9d9393b6..910a46c5e22 100644 --- a/src/sage/rings/polynomial/polynomial_ring_constructor.py +++ b/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. = PolynomialRing(RDF, implementation='singular') + sage: x^2 + 1 + x^2 + (1.000e + 00) + sage: C. = 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