diff --git a/src/module_list.py b/src/module_list.py index d7b5f459cf0..17df5141a79 100644 --- a/src/module_list.py +++ b/src/module_list.py @@ -795,6 +795,10 @@ def uname_specific(name, value, alternative): Extension('sage.matrix.matrix2', sources = ['sage/matrix/matrix2.pyx']), + Extension("sage.matrix.matrix_real_ball_dense", + ["sage/matrix/matrix_real_ball_dense.pyx"], + libraries=[arb_dylib_name]), + Extension("sage.matrix.matrix_complex_ball_dense", ["sage/matrix/matrix_complex_ball_dense.pyx"], libraries=[arb_dylib_name]), diff --git a/src/sage/libs/arb/arb_mat.pxd b/src/sage/libs/arb/arb_mat.pxd new file mode 100644 index 00000000000..a975baf6287 --- /dev/null +++ b/src/sage/libs/arb/arb_mat.pxd @@ -0,0 +1,10 @@ +# distutils: depends = arb_mat.h + +from sage.libs.arb.types cimport arb_t, arb_mat_t + +# arb_mat.h +cdef extern from "arb_wrap.h": + arb_t arb_mat_entry(arb_mat_t mat, unsigned long i, unsigned long j) + void arb_mat_init(arb_mat_t mat, long r, long c) + void arb_mat_clear(arb_mat_t mat) + bint arb_mat_solve(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, long prec) diff --git a/src/sage/libs/arb/types.pxd b/src/sage/libs/arb/types.pxd index 14ae24512c9..9f707b7eaeb 100644 --- a/src/sage/libs/arb/types.pxd +++ b/src/sage/libs/arb/types.pxd @@ -1,4 +1,4 @@ -# distutils: depends = mag.h arf.h arb.h acb.h acb_mat.h acb_poly.h acb_calc.h +# distutils: depends = mag.h arf.h arb.h acb.h arb_mat.h acb_mat.h acb_poly.h acb_calc.h # mag.h cdef extern from "arb_wrap.h": @@ -39,6 +39,12 @@ cdef extern from "arb_wrap.h": ctypedef acb_struct * acb_ptr ctypedef const acb_struct * acb_srcptr +# arb_mat.h +cdef extern from "arb_wrap.h": + ctypedef struct arb_mat_struct: + pass + ctypedef arb_mat_struct[1] arb_mat_t + # acb_mat.h cdef extern from "arb_wrap.h": ctypedef struct acb_mat_struct: diff --git a/src/sage/matrix/matrix_real_ball_dense.pxd b/src/sage/matrix/matrix_real_ball_dense.pxd new file mode 100644 index 00000000000..1e9b7dd2b2f --- /dev/null +++ b/src/sage/matrix/matrix_real_ball_dense.pxd @@ -0,0 +1,5 @@ +from sage.libs.arb.types cimport arb_mat_t +from .matrix_dense cimport Matrix_dense + +cdef class Matrix_real_ball_dense(Matrix_dense): + cdef arb_mat_t value diff --git a/src/sage/matrix/matrix_real_ball_dense.pyx b/src/sage/matrix/matrix_real_ball_dense.pyx new file mode 100644 index 00000000000..8e05c1eb4ff --- /dev/null +++ b/src/sage/matrix/matrix_real_ball_dense.pyx @@ -0,0 +1,199 @@ +r""" +Arbitrary precision :class:`RealBallField` matrices. + +This module contains :class:`Matrix_real_ball_dense`, a minimal +subclass of :class:`Matrix_dense` that can solving linear interval +systems. +""" +from cysignals.signals cimport sig_on, sig_str, sig_off +from sage.libs.arb.arb cimport * +from sage.libs.arb.arb_mat cimport * +from .args cimport SparseEntry, MatrixArgs_init +from sage.structure.parent cimport Parent +from sage.rings.real_arb cimport (RealBall) + +cdef class Matrix_real_ball_dense(Matrix_dense): + r""" + A matrix whose entries live in a `RealBallField`. At the moment, + the only nontrivial method it implements + :meth:`_solve_right_nonsingular_square`, to avoid the issues that + arise when the naive superclass implementation is used to perform + interval arithmetic. This is part of :trac:`29729`. + + The underlying matrix type is ``arb_mat`` of the Arb library. + + TESTS: + + Ensure that we can create a matrix over a field with a random + number of bits of precision, to catch corner cases:: + + sage: set_random_seed() + sage: bits = ZZ.random_element(2,1000) + sage: R = RealBallField(bits) + sage: nrows = ZZ.random_element(100) + sage: ncols = ZZ.random_element(100) + sage: matrix.random(R, nrows, ncols).base_ring() + Real ball field with... + + """ + def __cinit__(self): + r""" + Allocate memory for and initialize my underlying Arb matrix C + structure. + """ + sig_str("Arb exception") + arb_mat_init(self.value, self._nrows, self._ncols) + sig_off() + + + def __dealloc__(self): + """ + Free the memory allocated for my underlying Arb matrix C structure. + + TESTS: + + Ensure that nothing outrageous happens if we call ``del`` on + an instance of this class:: + + sage: set_random_seed() + sage: bits = ZZ.random_element(2,100) + sage: R = RealBallField(bits) + sage: nrows = ZZ.random_element(10) + sage: ncols = ZZ.random_element(10) + sage: A = matrix.random(R, nrows, ncols).base_ring() + sage: del(A) + sage: A + Traceback (most recent call last): + ... + NameError: name 'A' is not defined + + """ + arb_mat_clear(self.value) + + + def __init__(self, parent, entries=None, copy=None, bint coerce=True): + r""" + Initialize this dense matrix over a real ball field. + + INPUT: + + - ``parent`` -- a matrix space over a real ball field + + - ``entries`` -- the entries of this matrix, see :func:`matrix` + + - ``copy`` -- ignored (for backwards compatibility) + + - ``coerce`` -- if ``False``, assume without checking that the + entries lie in the base ring (this can crash SageMath if + you're wrong) + + EXAMPLES: + + We create a two-by-two matrix with the default precision:: + + sage: from sage.matrix.matrix_real_ball_dense import * + sage: MS = MatrixSpace(RBF, 2, 2) + sage: Matrix_real_ball_dense(MS, range(4)) + [ 0 1.000000000000000] + [2.000000000000000 3.000000000000000] + + With ``coerce=False``, nonsense (or worse...) can result:: + + sage: from sage.matrix.matrix_real_ball_dense import * + sage: MS = MatrixSpace(RBF, 1, 1) + sage: Matrix_real_ball_dense(MS, [pi], coerce=False) + [[+/- inf]] + + """ + ma = MatrixArgs_init(parent, entries) + cdef RealBall z + for t in ma.iter(coerce, True): + se = t + z = se.entry + arb_set(arb_mat_entry(self.value, se.i, se.j), z.value) + + + cdef set_unsafe(self, Py_ssize_t i, Py_ssize_t j, object x): + r""" + Set entry ``i,j1` of this matrix to ``x``, unsafely. + + The new value ``x`` must be of type ``RealBall`` and + ``i,j`` must be within bounds. + + INPUT: + + - ``i`` -- row index + + - ``j`` -- column index + + - ``x`` -- a ``RealBall``; the value to set self[i,j] to. + + TESTS: + + The unsafe setter is indirectly used by the safe one: + + sage: A = matrix(RBF, 1, 1, 0); A + [0] + sage: A[0,0] = e; A + [[2.718281828459045 +/- 5.35e-16]] + + """ + arb_set(arb_mat_entry(self.value, i, j), ( x).value) + + + cdef get_unsafe(self, Py_ssize_t i, Py_ssize_t j): + """ + Get entry ``i,j`` of this matrix, unsafely. + + The indices ``i,j`` must be within bounds. + + INPUT: + + - ``i`` -- row index + + - ``j`` -- column index + + .. NOTE:: + + The unsafe getter is indirectly tested by the safe getter that + is used everywhere (in particular, to display a matrix). + + """ + cdef RealBall z = RealBall.__new__(RealBall) + z._parent = self._base_ring + arb_set(z.value, arb_mat_entry(self.value, i, j)) + return z + + + def _solve_right_nonsingular_square(self, + Matrix_real_ball_dense rhs, + check_rank=None): + r""" + Override the superclass method with a square/nonsingular + algorithm (from Arb) tailored to interval systems. + + EXAMPLES: + + Numerical noise should not trick SageMath into thinking that + this system is nonsingular (:trac:`29729`):: + + sage: A = matrix(RBF, [[2/3, 1], [2/5, 3/5]]) + sage: b = vector(RBF, [1, 1]) + sage: A.solve_right(b) + Traceback (most recent call last): + ... + ValueError: unable to invert this matrix + + """ + result_space = self.matrix_space(self._nrows, rhs._ncols) + cdef Matrix_real_ball_dense res + res = Matrix_real_ball_dense(result_space, coerce=False) + soln_prec = min(self._base_ring._prec, rhs._base_ring._prec) + + sig_on() + success = arb_mat_solve(res.value, self.value, rhs.value, soln_prec) + sig_off() + if success: + return res + else: + raise ValueError("unable to invert this matrix") diff --git a/src/sage/matrix/matrix_space.py b/src/sage/matrix/matrix_space.py index acb2d53f210..0d437bdbf5c 100644 --- a/src/sage/matrix/matrix_space.py +++ b/src/sage/matrix/matrix_space.py @@ -185,6 +185,8 @@ def get_matrix_class(R, nrows, ncols, sparse, implementation): sage: type(matrix(GF(16007), 2, range(4))) + sage: type(matrix(RBF, 2, range(4))) + sage: type(matrix(CBF, 2, range(4))) sage: type(matrix(GF(2), 2, range(4))) @@ -245,6 +247,11 @@ def get_matrix_class(R, nrows, ncols, sparse, implementation): from . import matrix_symbolic_dense return matrix_symbolic_dense.Matrix_symbolic_dense + from sage.rings.real_arb import RealBallField + if isinstance(R, RealBallField): + from . import matrix_real_ball_dense + return matrix_real_ball_dense.Matrix_real_ball_dense + from sage.rings.complex_arb import ComplexBallField if isinstance(R, ComplexBallField): from . import matrix_complex_ball_dense