Ticket #10765: trac_10765-matrix-misc-routines.patch

File trac_10765-matrix-misc-routines.patch, 47.5 KB (added by rbeezer, 10 years ago)
  • module_list.py

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1297319126 28800
    # Node ID 14e86b60eb9b3b3b816c4127732b4f0f7ac4e643
    # Parent  b995f0056e8f09421b141a642b0ef312807befd1
    10765: Obliterate sage/matrix/misc.py
    
    diff -r b995f0056e8f -r 14e86b60eb9b module_list.py
    a b  
    822822              depends = numpy_depends),
    823823
    824824    Extension('sage.matrix.matrix2',
    825               sources = ['sage/matrix/matrix2.pyx']),
     825              sources = ['sage/matrix/matrix2.pyx'],
     826              libraries=['mpfr']),
    826827   
    827828    Extension('sage.matrix.matrix_complex_double_dense',
    828829              sources = ['sage/matrix/matrix_complex_double_dense.pyx'],
     
    926927              sources = ['sage/matrix/matrix_window_modn_dense.pyx']),
    927928
    928929    Extension('sage.matrix.misc',
    929               sources = ['sage/matrix/misc.pyx'],
    930               libraries=['mpfr','gmp']),
     930              sources = ['sage/matrix/misc.pyx']),
    931931   
    932932    Extension('sage.matrix.strassen',
    933933              sources = ['sage/matrix/strassen.pyx']),
  • sage/matrix/matrix2.pyx

    diff -r b995f0056e8f -r 14e86b60eb9b sage/matrix/matrix2.pyx
    a b  
    2929include "../ext/stdsage.pxi"
    3030include "../ext/python.pxi"
    3131
     32from sage.libs.mpfr cimport *
     33
    3234from sage.misc.randstate cimport randstate, current_randstate
    3335from sage.structure.sequence import Sequence
    3436from sage.combinat.combinat import combinations_iterator
     
    4749import matrix_space
    4850import berlekamp_massey
    4951from sage.modules.free_module_element import is_FreeModuleElement
     52from sage.rings.real_mpfr import  is_RealField
     53from sage.rings.real_mpfr cimport RealNumber
    5054
    5155cdef class Matrix(matrix1.Matrix):
    5256    def _backslash_(self, B):
     
    46364640        self.cache('echelon_form', self)
    46374641        verbose('done with gauss echelon form', tm)
    46384642
     4643    def _echelon_form_multimodular_rational(self, height_guess=None, proof=None):
     4644        r"""
     4645        Returns reduced row-echelon form using a multi-modular algorithm.
     4646        Does not change self.
     4647
     4648        REFERENCE:
     4649
     4650        - Chapter 7 of Stein's "Explicitly Computing Modular Forms".
     4651
     4652        INPUT:
     4653
     4654        -  ``height_guess`` - integer or None
     4655
     4656        -  ``proof`` - boolean (default: None, see
     4657           proof.linear_algebra or sage.structure.proof) Note that the Sage
     4658           global default is proof=True.
     4659
     4660        ALGORITHM:
     4661
     4662        The following is a modular algorithm for computing the echelon
     4663        form.  Define the height of a matrix to be the max of the
     4664        absolute values of the entries.
     4665
     4666        Given Matrix A with n columns (self).
     4667
     4668        0. Rescale input matrix A to have integer entries.  This does
     4669            not change echelon form and makes reduction modulo lots of
     4670            primes significantly easier if there were denominators.
     4671            Henceforth we assume A has integer entries.
     4672
     4673        1. Let c be a guess for the height of the echelon form.  E.g.,
     4674            c=1000, e.g., if matrix is very sparse and application is to
     4675            computing modular symbols.
     4676
     4677        2. Let M = n * c * H(A) + 1,
     4678            where n is the number of columns of A.
     4679
     4680        3. List primes p_1, p_2, ..., such that the product of
     4681            the p_i is at least M.
     4682
     4683        4. Try to compute the rational reconstruction CRT echelon form
     4684            of A mod the product of the p_i.  If rational
     4685            reconstruction fails, compute 1 more echelon forms mod the
     4686            next prime, and attempt again.  Make sure to keep the
     4687            result of CRT on the primes from before, so we don't have
     4688            to do that computation again.  Let E be this matrix.
     4689
     4690        5. Compute the denominator d of E.
     4691            Attempt to prove that result is correct by checking that
     4692
     4693                H(d*E)*ncols(A)*H(A) < (prod of reduction primes)
     4694
     4695            where H denotes the height.   If this fails, do step 4 with
     4696            a few more primes.
     4697
     4698        EXAMPLES:
     4699            sage: A = matrix(QQ, 3, 7, [1..21])
     4700            sage: A._echelon_form_multimodular_rational()
     4701            [ 1  0 -1 -2 -3 -4 -5]
     4702            [ 0  1  2  3  4  5  6]
     4703            [ 0  0  0  0  0  0  0]
     4704
     4705            sage: A = matrix(QQ, 3, 4, [0,0] + [1..9] + [-1/2^20])
     4706            sage: A._echelon_form_multimodular_rational()
     4707            [                1                 0                 0 -10485761/1048576]
     4708            [                0                 1                 0  27262979/4194304]
     4709            [                0                 0                 1                 2]
     4710            sage: A.echelon_form()
     4711            [                1                 0                 0 -10485761/1048576]
     4712            [                0                 1                 0  27262979/4194304]
     4713            [                0                 0                 1                 2]
     4714        """
     4715        from sage.rings.arith import previous_prime, CRT_basis
     4716        from sage.rings.rational_field import is_RationalField
     4717        # next two imports for MAX_MODULUS in the two cases
     4718        import matrix_modn_dense
     4719        import matrix_modn_sparse
     4720
     4721        if not is_RationalField(self.base_ring()):
     4722            raise TypeError('_echelon_form_multimodular_rational() can only be used for matrices over the rationals, not over %s' % self.base_ring())
     4723
     4724        if proof is None:
     4725            from sage.structure.proof.proof import get_flag
     4726            proof = get_flag(proof, "linear_algebra")
     4727
     4728        verbose("Multimodular echelon algorithm on %s x %s matrix"%(self._nrows, self._ncols), caller_name="multimod echelon")
     4729        cdef Matrix E
     4730        if self._nrows == 0 or self._ncols == 0:
     4731            self.cache('in_echelon_form', True)
     4732            self.cache('echelon_form', self)
     4733            self.cache('pivots', [])
     4734            return self
     4735
     4736        B, _ = self._clear_denom()
     4737
     4738        height = self.height()
     4739        if height_guess is None:
     4740            height_guess = 10000000*(height+100)
     4741        tm = verbose("height_guess = %s"%height_guess, level=2, caller_name="multimod echelon")
     4742
     4743        if proof:
     4744            M = self._ncols * height_guess * height  +  1
     4745        else:
     4746            M = height_guess + 1
     4747
     4748        if self.is_sparse():
     4749            p = matrix_modn_sparse.MAX_MODULUS + 1
     4750        else:
     4751            p = matrix_modn_dense.MAX_MODULUS + 1
     4752        X = []
     4753        best_pivots = []
     4754        prod = 1
     4755        problem = 0
     4756        lifts = {}
     4757        while True:
     4758            p = previous_prime(p)
     4759            while prod < M:
     4760                problem = problem + 1
     4761                if problem > 50:
     4762                    verbose("echelon multi-modular possibly not converging?", caller_name="multimod echelon")
     4763                t = verbose("echelon modulo p=%s (%.2f%% done)"%(
     4764                        p, 100*float(len(str(prod))) / len(str(M))), level=2, caller_name="multimod echelon")
     4765
     4766                # We use denoms=False, since we made self integral by calling clear_denom above.
     4767                A = B._mod_int(p)
     4768                t = verbose("time to reduce matrix mod p:",t, level=2, caller_name="multimod echelon")
     4769                A.echelonize()
     4770                t = verbose("time to put reduced matrix in echelon form:",t, level=2, caller_name="multimod echelon")
     4771
     4772                # a worthwhile check / shortcut.
     4773                if self._nrows == self._ncols and len(A.pivots()) == self._nrows:
     4774                    verbose("done: the echelon form mod p is the identity matrix", caller_name="multimod echelon")
     4775                    E = self.parent().identity_matrix()
     4776                    E.cache('pivots', range(self._nrows))
     4777                    E.cache('in_echelon_form', True)
     4778                    self.cache('in_echelon_form', True)
     4779                    self.cache('echelon_form', E)
     4780                    self.cache('pivots', range(self._nrows))
     4781                    return E
     4782
     4783                c = cmp_pivots(best_pivots, A.pivots())
     4784                if c <= 0:
     4785                    best_pivots = A.pivots()
     4786                    X.append(A)
     4787                    prod = prod * p
     4788                else:
     4789                    # do not save A since it is bad.
     4790                    verbose("Excluding this prime (bad pivots).", caller_name="multimod echelon")
     4791                t = verbose("time for pivot compare", t, level=2, caller_name="multimod echelon")
     4792                p = previous_prime(p)
     4793            # Find set of best matrices.
     4794            Y = []
     4795            # recompute product, since may drop bad matrices
     4796            prod = 1
     4797            t = verbose("now comparing pivots and dropping any bad ones", level=2, t=t, caller_name="multimod echelon")
     4798            for i in range(len(X)):
     4799                if cmp_pivots(best_pivots, X[i].pivots()) <= 0:
     4800                    p = X[i].base_ring().order()
     4801                    if not lifts.has_key(p):
     4802                        t0 = verbose("Lifting a good matrix", level=2, caller_name = "multimod echelon")
     4803                        lift = X[i].lift()
     4804                        lifts[p] = (lift, p)
     4805                        verbose("Finished lift", level=2, caller_name= "multimod echelon", t=t0)
     4806                    Y.append(lifts[p])
     4807                    prod = prod * X[i].base_ring().order()
     4808            verbose("finished comparing pivots", level=2, t=t, caller_name="multimod echelon")
     4809            try:
     4810                if len(Y) == 0:
     4811                    raise ValueError("not enough primes")
     4812                t = verbose("start crt linear combination", level=2, caller_name="multimod echelon")
     4813                a = CRT_basis([w[1] for w in Y])
     4814                t = verbose('got crt basis', level=2, t=t, caller_name="multimod echelon")
     4815
     4816                # take the linear combination of the lifts of the elements
     4817                # of Y times coefficients in a
     4818                L = a[0]*(Y[0][0])
     4819                assert Y[0][0].is_sparse() == L.is_sparse()
     4820                for j in range(1,len(Y)):
     4821                    L += a[j]*(Y[j][0])
     4822                verbose("time to take linear combination of matrices over ZZ is",t, level=2, caller_name="multimod echelon")
     4823                t = verbose("now doing rational reconstruction", level=2, caller_name="multimod echelon")
     4824                E = L.rational_reconstruction(prod)
     4825                L = 0  # free memory
     4826                verbose('rational reconstruction completed', t, level=2, caller_name="multimod echelon")
     4827            except ValueError, msg:
     4828                verbose(msg, level=2)
     4829                verbose("Not enough primes to do CRT lift; redoing with several more primes.", level=2, caller_name="multimod echelon")
     4830                M = prod * p*p*p
     4831                continue
     4832
     4833            if not proof:
     4834                verbose("Not checking validity of result (since proof=False).", level=2, caller_name="multimod echelon")
     4835                break
     4836            d   = E.denominator()
     4837            hdE = long((d*E).height())
     4838            if hdE * self.ncols() * height < prod:
     4839                verbose("Validity of result checked.", level=2, caller_name="multimod echelon")
     4840                break
     4841            verbose("Validity failed; trying again with more primes.", level=2, caller_name="multimod echelon")
     4842            M = prod * p*p*p
     4843        #end while
     4844        verbose("total time",tm, level=2, caller_name="multimod echelon")
     4845        self.cache('pivots', best_pivots)
     4846        E.cache('pivots', best_pivots)
     4847        return E
     4848
    46394849    def weak_popov_form(self, ascend=True):
    46404850        """
    46414851        This function computes a weak Popov form of a matrix over a rational
     
    65596769            import misc
    65606770            R = RealField(53, rnd='RNDU')
    65616771            A = self.change_ring(R)
    6562             m1 = misc.hadamard_row_bound_mpfr(A)
     6772            m1 = A._hadamard_row_bound_mpfr()
    65636773            A = A.transpose()
    6564             m2 = misc.hadamard_row_bound_mpfr(A)
     6774            m2 = A._hadamard_row_bound_mpfr()
    65656775            return min(m1, m2)
    6566        
     6776
     6777    def _hadamard_row_bound_mpfr(self):
     6778        r"""
     6779        Given a matrix ``self`` with entries that coerce to RR, compute the row
     6780        Hadamard bound on the determinant.
     6781
     6782        OUTPUT:
     6783
     6784        An integer `n` such that the absolute value of the
     6785        determinant of this matrix is at most `10^n`.
     6786
     6787        EXAMPLES:
     6788
     6789        We create a very large matrix, compute the row Hadamard bound,
     6790        and also compute the row Hadamard bound of the transpose, which
     6791        happens to be sharp. ::
     6792
     6793            sage: a = matrix(ZZ, 2, [2^10000,3^10000,2^50,3^19292])
     6794            sage: (a.change_ring(RR))._hadamard_row_bound_mpfr()
     6795            13976
     6796            sage: len(str(a.det()))
     6797            12215
     6798            sage: (a.transpose().change_ring(RR))._hadamard_row_bound_mpfr()
     6799            12215
     6800
     6801        Note that in the above example using RDF would overflow. ::
     6802
     6803            sage: b = a.change_ring(RDF)
     6804            sage: b._hadamard_row_bound()
     6805            Traceback (most recent call last):
     6806            ...
     6807            OverflowError: cannot convert float infinity to integer
     6808        """
     6809
     6810        if not is_RealField(self.base_ring()):
     6811            raise TypeError("A must have base field an mpfr real field.")
     6812
     6813        cdef RealNumber a, b
     6814        cdef mpfr_t s, d, pr
     6815        cdef Py_ssize_t i, j
     6816
     6817        mpfr_init(s)
     6818        mpfr_init(d)
     6819        mpfr_init(pr)
     6820        mpfr_set_si(d, 0, GMP_RNDU)
     6821
     6822        for i from 0 <= i < self._nrows:
     6823            mpfr_set_si(s, 0, GMP_RNDU)
     6824            for j from 0 <= j < self._ncols:
     6825                a = self.get_unsafe(i, j)
     6826                mpfr_mul(pr, a.value, a.value, GMP_RNDU)
     6827                mpfr_add(s, s, pr, GMP_RNDU)
     6828            mpfr_log10(s, s, GMP_RNDU)
     6829            mpfr_add(d, d, s, GMP_RNDU)
     6830        b = a._new()
     6831        mpfr_set(b.value, d, GMP_RNDU)
     6832        b /= 2
     6833        mpfr_clear(s)
     6834        mpfr_clear(d)
     6835        mpfr_clear(pr)
     6836        return b.ceil()
     6837
    65676838    def find(self,f, indices=False):
    65686839        r"""
    65696840        Find elements in this matrix satisfying the constraints in the
  • sage/matrix/matrix_cyclo_dense.pyx

    diff -r b995f0056e8f -r 14e86b60eb9b sage/matrix/matrix_cyclo_dense.pyx
    a b  
    16501650                verbose("attempting rational reconstruction ...", level=echelon_verbose_level)
    16511651                res = Matrix_cyclo_dense.__new__(Matrix_cyclo_dense, self.parent(),
    16521652                                                 None, None, None)
    1653                 res._matrix = <Matrix_rational_dense>matrix_integer_dense_rational_reconstruction(mat_over_ZZ, prod)
     1653                res._matrix = mat_over_ZZ.rational_reconstruction(prod)
    16541654
    16551655            except ValueError:
    16561656                # If a ValueError is raised here, it means that the
  • sage/matrix/matrix_integer_dense.pyx

    diff -r b995f0056e8f -r 14e86b60eb9b sage/matrix/matrix_integer_dense.pyx
    a b  
    30963096                v.append(str(self.get_unsafe(i,j)))
    30973097        return ' '.join(v)
    30983098
    3099     def rational_reconstruction(self, N):
    3100         """
    3101         Use rational reconstruction to lift self to a matrix over the
     3099    def rational_reconstruction(self, Integer N):
     3100        r"""
     3101        Use rational reconstruction to lift ``self`` to a matrix over the
    31023102        rational numbers (if possible), where we view self as a matrix
    3103         modulo N.
    3104        
     3103        modulo ``N``.  This is done efficiently by assuming there is a
     3104        large common factor dividing the denominators.
     3105
    31053106        INPUT:
    3106        
    3107        
     3107
    31083108        -  ``N`` - an integer
    3109        
    3110        
     3109
    31113110        OUTPUT:
    3112        
    3113        
    3114         -  ``matrix`` - over QQ or raise a ValueError
    3115        
    3116        
    3117         EXAMPLES: We create a random 4x4 matrix over ZZ.
    3118        
    3119         ::
    3120        
     3111
     3112        -  ``matrix`` - over QQ or raise a ``ValueError``
     3113
     3114        EXAMPLES::
     3115
     3116            sage: B = ((matrix(ZZ, 3,4, [1,2,3,-4,7,2,18,3,4,3,4,5])/3)%500).change_ring(ZZ)
     3117            sage: B.rational_reconstruction(500)
     3118            [ 1/3  2/3    1 -4/3]
     3119            [ 7/3  2/3    6    1]
     3120            [ 4/3    1  4/3  5/3]
     3121
     3122        We create a random 4x4 matrix over ZZ.  ::
     3123
    31213124            sage: A = matrix(ZZ, 4, [4, -4, 7, 1, -1, 1, -1, -12, -1, -1, 1, -1, -3, 1, 5, -1])
    3122        
    3123         There isn't a unique rational reconstruction of it::
    3124        
     3125
     3126        There is not a unique rational reconstruction of ``A``.  ::
     3127
    31253128            sage: A.rational_reconstruction(11)
    31263129            Traceback (most recent call last):
    31273130            ...
    31283131            ValueError: Rational reconstruction of 4 (mod 11) does not exist.
    3129        
     3132
    31303133        We throw in a denominator and reduce the matrix modulo 389 - it
    3131         does rationally reconstruct.
    3132        
    3133         ::
    3134        
     3134        does rationally reconstruct.  ::
     3135
    31353136            sage: B = (A/3 % 389).change_ring(ZZ)
    31363137            sage: B.rational_reconstruction(389) == A/3
    31373138            True
    31383139
    31393140        TEST:
    3140        
    3141         Check that ticket #9345 is fixed::
     3141
     3142        Check that ticket #9345 is fixed. ::
    31423143
    31433144            sage: A = random_matrix(ZZ, 3, 3)
    31443145            sage: A.rational_reconstruction(0)
     
    31463147            ...
    31473148            ZeroDivisionError: The modulus cannot be zero
    31483149        """
    3149         import misc
    3150         return misc.matrix_integer_dense_rational_reconstruction(self, N)
     3150        if not N.__nonzero__():
     3151            raise ZeroDivisionError("The modulus cannot be zero")
     3152        cdef Matrix_rational_dense R
     3153        R = Matrix_rational_dense.__new__(Matrix_rational_dense,
     3154                                        self.parent().change_ring(QQ), 0,0,0)
     3155
     3156        cdef mpz_t a, bnd, other_bnd, one, denom
     3157        mpz_init_set_si(denom, 1)
     3158        mpz_init(a)
     3159        mpz_init_set_si(one, 1)
     3160        mpz_init(other_bnd)
     3161
     3162        cdef Integer _bnd
     3163        import math
     3164        _bnd = (N//2).isqrt()
     3165        mpz_init_set(bnd, _bnd.value)
     3166        mpz_sub(other_bnd, N.value, bnd)
     3167
     3168        cdef Py_ssize_t i, j
     3169        cdef int do_it
     3170
     3171        for i from 0 <= i < self._nrows:
     3172            for j from 0 <= j < self._ncols:
     3173                mpz_set(a, self._matrix[i][j])
     3174                if mpz_cmp(denom, one) != 0:
     3175                    mpz_mul(a, a, denom)
     3176                mpz_fdiv_r(a, a, N.value)
     3177                do_it = 0
     3178                if mpz_cmp(a, bnd) <= 0:
     3179                    do_it = 1
     3180                elif mpz_cmp(a, other_bnd) >= 0:
     3181                    mpz_sub(a, a, N.value)
     3182                    do_it = 1
     3183                if do_it:
     3184                    mpz_set(mpq_numref(R._matrix[i][j]), a)
     3185                    if mpz_cmp(denom, one) != 0:
     3186                        mpz_set(mpq_denref(R._matrix[i][j]), denom)
     3187                        mpq_canonicalize(R._matrix[i][j])
     3188                    else:
     3189                        mpz_set_si(mpq_denref(R._matrix[i][j]), 1)
     3190                else:
     3191                    # Otherwise have to do it the hard way
     3192                    mpq_rational_reconstruction(R._matrix[i][j], self._matrix[i][j], N.value)
     3193                    mpz_lcm(denom, denom, mpq_denref(R._matrix[i][j]))
     3194
     3195        mpz_clear(denom)
     3196        mpz_clear(a)
     3197        mpz_clear(one)
     3198        mpz_clear(other_bnd)
     3199        mpz_clear(bnd)
     3200        return R
    31513201
    31523202    def randomize(self, density=1, x=None, y=None, distribution=None, \
    31533203                  nonzero=False):
  • sage/matrix/matrix_integer_sparse.pyx

    diff -r b995f0056e8f -r 14e86b60eb9b sage/matrix/matrix_integer_sparse.pyx
    a b  
    2222include '../modules/binary_search.pxi'
    2323include '../modules/vector_integer_sparse_h.pxi'
    2424include '../modules/vector_integer_sparse_c.pxi'
     25include '../modules/vector_rational_sparse_h.pxi'
     26include '../modules/vector_rational_sparse_c.pxi'
    2527include '../modules/vector_modn_sparse_h.pxi'
    2628include '../modules/vector_modn_sparse_c.pxi'
    2729
     
    3840from sage.rings.integer_ring import ZZ
    3941from sage.rings.finite_rings.integer_mod_ring import IntegerModRing
    4042
     43from matrix_rational_sparse cimport Matrix_rational_sparse
     44
    4145
    4246cdef class Matrix_integer_sparse(matrix_sparse.Matrix_sparse):   
    4347
     
    352356        return res
    353357       
    354358
    355     def rational_reconstruction(self, N):
     359    def rational_reconstruction(self, Integer N):
    356360        """
    357         Use rational reconstruction to lift self to a matrix over the
     361        Use rational reconstruction to lift self to a matrix over the 
    358362        rational numbers (if possible), where we view self as a matrix
    359363        modulo N.
    360364
     
    375379            ...
    376380            ZeroDivisionError: The modulus cannot be zero
    377381        """
    378         import misc
    379         return misc.matrix_integer_sparse_rational_reconstruction(self, N)
     382        from sage.rings.rational_field import QQ
     383
     384        if not N.__nonzero__():
     385            raise ZeroDivisionError("The modulus cannot be zero")
     386        cdef Matrix_rational_sparse R
     387        R = Matrix_rational_sparse.__new__(Matrix_rational_sparse,
     388                                        self.parent().change_ring(QQ), 0,0,0)
     389
     390        cdef mpq_t t
     391        mpq_init(t)
     392
     393        cdef mpz_t a, bnd, other_bnd, one, denom
     394        mpz_init_set_si(denom, 1)
     395        mpz_init(a)
     396        mpz_init_set_si(one, 1)
     397        mpz_init(other_bnd)
     398
     399        cdef Integer _bnd
     400        import math
     401        _bnd = (N//2).isqrt()
     402        mpz_init_set(bnd, _bnd.value)
     403        mpz_sub(other_bnd, N.value, bnd)
     404
     405        cdef Py_ssize_t i, j
     406        cdef int do_it
     407        cdef mpz_vector* A_row
     408        cdef mpq_vector* R_row
     409
     410        for i from 0 <= i < self._nrows:
     411            A_row = &self._matrix[i]
     412            R_row = &R._matrix[i]
     413            reallocate_mpq_vector(R_row, A_row.num_nonzero)
     414            R_row.num_nonzero = A_row.num_nonzero
     415            R_row.degree = A_row.degree
     416            for j from 0 <= j < A_row.num_nonzero:
     417                mpz_set(a, A_row.entries[j])
     418                if mpz_cmp(denom, one) != 0:
     419                    mpz_mul(a, a, denom)
     420                mpz_fdiv_r(a, a, N.value)
     421                do_it = 0
     422                if mpz_cmp(a, bnd) <= 0:
     423                    do_it = 1
     424                elif mpz_cmp(a, other_bnd) >= 0:
     425                    mpz_sub(a, a, N.value)
     426                    do_it = 1
     427                if do_it:
     428                    mpz_set(mpq_numref(t), a)
     429                    if mpz_cmp(denom, one) != 0:
     430                        mpz_set(mpq_denref(t), denom)
     431                        mpq_canonicalize(t)
     432                    else:
     433                        mpz_set_si(mpq_denref(t), 1)
     434                    mpq_set(R_row.entries[j], t)
     435                    R_row.positions[j] = A_row.positions[j]
     436                else:
     437                    # Otherwise have to do it the hard way
     438                    mpq_rational_reconstruction(t, A_row.entries[j], N.value)
     439                    mpq_set(R_row.entries[j], t)
     440                    R_row.positions[j] = A_row.positions[j]
     441                    mpz_lcm(denom, denom, mpq_denref(t))
     442
     443        mpq_clear(t)
     444
     445        mpz_clear(denom)
     446        mpz_clear(a)
     447        mpz_clear(one)
     448        mpz_clear(other_bnd)
     449        mpz_clear(bnd)
     450        return R
    380451
    381452    def _linbox_sparse(self):
    382453        cdef Py_ssize_t i, j
  • sage/matrix/matrix_rational_dense.pyx

    diff -r b995f0056e8f -r 14e86b60eb9b sage/matrix/matrix_rational_dense.pyx
    a b  
    16251625           proof.linear_algebra or sage.structure.proof) Note that the Sage
    16261626           global default is proof=True.
    16271627        """
    1628         import misc
    1629         return misc.matrix_rational_echelon_form_multimodular(self,
    1630                                  height_guess=height_guess, proof=proof)
     1628        return self._echelon_form_multimodular_rational(height_guess=height_guess, proof=proof)
    16311629
    16321630    def _echelon_in_place_classical(self):
    16331631        """
  • sage/matrix/matrix_rational_sparse.pyx

    diff -r b995f0056e8f -r 14e86b60eb9b sage/matrix/matrix_rational_sparse.pyx
    a b  
    560560            height_guess -- integer or None
    561561            proof -- boolean (default: True)
    562562        """
    563         import misc
    564         cdef Matrix E = misc.matrix_rational_echelon_form_multimodular(self,
    565                                  height_guess=height_guess, proof=proof)
     563        cdef Matrix E = self._echelon_form_multimodular_rational(height_guess=height_guess, proof=proof)
    566564        E._parent = self._parent
    567565        return E
    568566
  • sage/matrix/misc.pyx

    diff -r b995f0056e8f -r 14e86b60eb9b sage/matrix/misc.pyx
    a b  
    1 """
     1r"""
    22Misc matrix algorithms
    33
    4 Code goes here mainly when it needs access to the internal structure
    5 of several classes, and we want to avoid circular cimports.
     4All of these routines ahve been moved to other modules.
     5So there is nothing to see here, really.
     6"""
     7# When deprecation period expires, this file can be removed.
     8# Also remove its entry from module_list.py,
     9# as this does not need to be described as an extension there
    610
    7 NOTE: The whole problem of avoiding circular imports -- the reason for
    8 existence of this file -- is now a non-issue, since some bugs in
    9 Cython were fixed.  Probably all this code should be moved into the
    10 relevant classes and this file deleted.
    11 """
    1211
    1312include "../ext/interrupt.pxi"
    1413include "../ext/cdefs.pxi"
    1514include '../ext/stdsage.pxi'
    1615
    17 from sage.libs.mpfr cimport *
    18 
    19 include '../modules/binary_search.pxi'
    20 include '../modules/vector_integer_sparse_h.pxi'
    21 include '../modules/vector_integer_sparse_c.pxi'
    22 include '../modules/vector_rational_sparse_h.pxi'
    23 include '../modules/vector_rational_sparse_c.pxi'
    24 include '../modules/vector_modn_sparse_h.pxi'
    25 include '../modules/vector_modn_sparse_c.pxi'
    26 
    2716cdef extern from "../ext/multi_modular.h":
    2817    ctypedef unsigned long mod_int
    2918    mod_int MOD_INT_OVERFLOW
    3019
    31 
     20# Needed to preserve function signatures
     21from sage.rings.integer  cimport Integer
    3222from matrix0 cimport Matrix
    33 from matrix_modn_dense cimport Matrix_modn_dense
    34 from matrix_modn_sparse cimport Matrix_modn_sparse
    3523from matrix_integer_dense cimport Matrix_integer_dense
    3624from matrix_integer_sparse cimport Matrix_integer_sparse
    37 from matrix_rational_dense cimport Matrix_rational_dense
    38 from matrix_rational_sparse cimport Matrix_rational_sparse
    3925
    40 import matrix_modn_dense
    41 import matrix_modn_sparse
    42 
    43 from sage.rings.integer_ring   import ZZ
    44 from sage.rings.rational_field import QQ
    45 
    46 from sage.rings.integer cimport Integer
    47 from sage.rings.arith import previous_prime, CRT_basis
    48 
    49 from sage.rings.real_mpfr import  is_RealField
    50 from sage.rings.real_mpfr cimport RealNumber
    51 
    52 
    53 from sage.misc.misc import verbose, get_verbose
     26from sage.misc.misc import deprecation
    5427
    5528def matrix_integer_dense_rational_reconstruction(Matrix_integer_dense A, Integer N):
    56     """
    57     Given a matrix over the integers and an integer modulus, do
    58     rational reconstruction on all entries of the matrix, viewed as
    59     numbers mod N.  This is done efficiently by assuming there is a
    60     large common factor dividing the denominators.
     29    r"""
     30    This function is deprecated in favor of
     31    :meth:`sage.matrix.matrix_integer_dense.Matrix_integer_dense.rational_reconstruction` which behaves identically.
    6132
    62     INPUT:
    63         A -- matrix
    64         N -- an integer
     33    EXAMPLE::
    6534
    66     EXAMPLES:
    6735        sage: B = ((matrix(ZZ, 3,4, [1,2,3,-4,7,2,18,3,4,3,4,5])/3)%500).change_ring(ZZ)
    6836        sage: sage.matrix.misc.matrix_integer_dense_rational_reconstruction(B, 500)
    69         [ 1/3  2/3    1 -4/3]
    70         [ 7/3  2/3    6    1]
    71         [ 4/3    1  4/3  5/3]
    72 
    73     TEST:
    74 
    75     Check that ticket #9345 is fixed::
    76 
    77         sage: A = random_matrix(ZZ, 3)
    78         sage: sage.matrix.misc.matrix_integer_dense_rational_reconstruction(A, 0)
    79         Traceback (most recent call last):
    80         ...
    81         ZeroDivisionError: The modulus cannot be zero
     37        doctest:...: DeprecationWarning: (Since Sage Version 4.7) The matrix_integer_dense_rational_reconstruction() function in sage.matrix.misc has been deprecated, use the rational_reconstruction() method from sage.matrix.matrix_integer_dense.Matrix_integer_dense instead
    8238    """
    83     if not N.__nonzero__():
    84         raise ZeroDivisionError("The modulus cannot be zero")
    85     cdef Matrix_rational_dense R
    86     R = Matrix_rational_dense.__new__(Matrix_rational_dense,
    87                                       A.parent().change_ring(QQ), 0,0,0)
    88 
    89     cdef mpz_t a, bnd, other_bnd, one, denom
    90     mpz_init_set_si(denom, 1)
    91     mpz_init(a)
    92     mpz_init_set_si(one, 1)
    93     mpz_init(other_bnd)
    94 
    95     cdef Integer _bnd
    96     import math
    97     _bnd = (N//2).isqrt()
    98     mpz_init_set(bnd, _bnd.value)
    99     mpz_sub(other_bnd, N.value, bnd)
    100 
    101     cdef Py_ssize_t i, j
    102     cdef int do_it
    103    
    104     for i from 0 <= i < A._nrows:
    105         for j from 0 <= j < A._ncols:
    106             mpz_set(a, A._matrix[i][j])
    107             if mpz_cmp(denom, one) != 0:
    108                 mpz_mul(a, a, denom)
    109             mpz_fdiv_r(a, a, N.value)
    110             do_it = 0
    111             if mpz_cmp(a, bnd) <= 0:
    112                 do_it = 1
    113             elif mpz_cmp(a, other_bnd) >= 0:
    114                 mpz_sub(a, a, N.value)
    115                 do_it = 1
    116             if do_it:
    117                 mpz_set(mpq_numref(R._matrix[i][j]), a)
    118                 if mpz_cmp(denom, one) != 0:
    119                     mpz_set(mpq_denref(R._matrix[i][j]), denom)
    120                     mpq_canonicalize(R._matrix[i][j])
    121                 else:
    122                     mpz_set_si(mpq_denref(R._matrix[i][j]), 1)
    123             else:
    124                 # Otherwise have to do it the hard way
    125                 mpq_rational_reconstruction(R._matrix[i][j], A._matrix[i][j], N.value)
    126                 mpz_lcm(denom, denom, mpq_denref(R._matrix[i][j]))
    127            
    128     mpz_clear(denom)
    129     mpz_clear(a)
    130     mpz_clear(one)
    131     mpz_clear(other_bnd)
    132     mpz_clear(bnd)
    133     return R
     39    deprecation('The matrix_integer_dense_rational_reconstruction() function in sage.matrix.misc has been deprecated, use the rational_reconstruction() method from sage.matrix.matrix_integer_dense.Matrix_integer_dense instead', 'Sage Version 4.7')
    13440
    13541def matrix_integer_sparse_rational_reconstruction(Matrix_integer_sparse A, Integer N):
    136     """
    137     Given a sparse matrix over the integers and an integer modulus, do
    138     rational reconstruction on all entries of the matrix, viewed as
    139     numbers mod N.
    140    
    141     EXAMPLES:
     42    r"""
     43    This function is deprecated in favor of
     44    :meth:`sage.matrix.matrix_integer_sparse.Matrix_integer_sparse.rational_reconstruction` which behaves identically.
     45
     46    EXAMPLE::
     47
    14248        sage: A = matrix(ZZ, 3, 4, [(1/3)%500, 2, 3, (-4)%500, 7, 2, 2, 3, 4, 3, 4, (5/7)%500], sparse=True)
    14349        sage: sage.matrix.misc.matrix_integer_sparse_rational_reconstruction(A, 500)
    144         [1/3   2   3  -4]
    145         [  7   2   2   3]
    146         [  4   3   4 5/7]
    147 
    148     TEST:
    149 
    150     Check that ticket #9345 is fixed::
    151 
    152         sage: A = random_matrix(ZZ, 3, sparse=True)
    153         sage: sage.matrix.misc.matrix_integer_sparse_rational_reconstruction(A, 0)
    154         Traceback (most recent call last):
    155         ...
    156         ZeroDivisionError: The modulus cannot be zero
     50        doctest:...: DeprecationWarning: (Since Sage Version 4.7) The matrix_integer_sparse_rational_reconstruction() function in sage.matrix.misc has been deprecated, use the rational_reconstruction() method from sage.matrix.matrix_integer_sparse.Matrix_integer_sparse instead
    15751    """
    158     if not N.__nonzero__():
    159         raise ZeroDivisionError("The modulus cannot be zero")
    160     cdef Matrix_rational_sparse R
    161     R = Matrix_rational_sparse.__new__(Matrix_rational_sparse,
    162                                       A.parent().change_ring(QQ), 0,0,0)
    163    
    164     cdef mpq_t t
    165     mpq_init(t)
    166 
    167     cdef mpz_t a, bnd, other_bnd, one, denom
    168     mpz_init_set_si(denom, 1)
    169     mpz_init(a)
    170     mpz_init_set_si(one, 1)
    171     mpz_init(other_bnd)
    172 
    173     cdef Integer _bnd
    174     import math
    175     _bnd = (N//2).isqrt()
    176     mpz_init_set(bnd, _bnd.value)
    177     mpz_sub(other_bnd, N.value, bnd)
    178 
    179     cdef Py_ssize_t i, j
    180     cdef int do_it
    181     cdef mpz_vector* A_row
    182     cdef mpq_vector* R_row
    183    
    184     for i from 0 <= i < A._nrows:
    185         A_row = &A._matrix[i]
    186         R_row = &R._matrix[i]
    187         reallocate_mpq_vector(R_row, A_row.num_nonzero)
    188         R_row.num_nonzero = A_row.num_nonzero
    189         R_row.degree = A_row.degree
    190         for j from 0 <= j < A_row.num_nonzero:
    191             mpz_set(a, A_row.entries[j])
    192             if mpz_cmp(denom, one) != 0:
    193                 mpz_mul(a, a, denom)
    194             mpz_fdiv_r(a, a, N.value)
    195             do_it = 0
    196             if mpz_cmp(a, bnd) <= 0:
    197                 do_it = 1
    198             elif mpz_cmp(a, other_bnd) >= 0:
    199                 mpz_sub(a, a, N.value)
    200                 do_it = 1
    201             if do_it:
    202                 mpz_set(mpq_numref(t), a)
    203                 if mpz_cmp(denom, one) != 0:
    204                     mpz_set(mpq_denref(t), denom)
    205                     mpq_canonicalize(t)
    206                 else:
    207                     mpz_set_si(mpq_denref(t), 1)
    208                 mpq_set(R_row.entries[j], t)
    209                 R_row.positions[j] = A_row.positions[j]
    210             else:
    211                 # Otherwise have to do it the hard way
    212                 mpq_rational_reconstruction(t, A_row.entries[j], N.value)
    213                 mpq_set(R_row.entries[j], t)
    214                 R_row.positions[j] = A_row.positions[j]
    215                 mpz_lcm(denom, denom, mpq_denref(t))
    216                
    217     mpq_clear(t)
    218 
    219     mpz_clear(denom)
    220     mpz_clear(a)
    221     mpz_clear(one)
    222     mpz_clear(other_bnd)
    223     mpz_clear(bnd)
    224     return R
    225 
     52    deprecation('The matrix_integer_sparse_rational_reconstruction() function in sage.matrix.misc has been deprecated, use the rational_reconstruction() method from sage.matrix.matrix_integer_sparse.Matrix_integer_sparse instead', 'Sage Version 4.7')
    22653
    22754def matrix_rational_echelon_form_multimodular(Matrix self, height_guess=None, proof=None):
    228     """
    229     Returns reduced row-echelon form using a multi-modular
    230     algorithm.  Does not change self.
     55    r"""
     56    This function is deprecated in favor of
     57    :meth:`sage.matrix.matrix2.Matrix._echelon_form_multimodular_rational` which behaves identically.
    23158
    232     REFERENCE: Chapter 7 of Stein's "Explicitly Computing Modular Forms".
     59    EXAMPLE::
    23360
    234     INPUT:
    235         height_guess -- integer or None
    236         proof -- boolean or None (default: None, see proof.linear_algebra or
    237                                sage.structure.proof).
    238                         Note that the global Sage default is proof=True
    239 
    240     ALGORITHM:
    241     The following is a modular algorithm for computing the echelon
    242     form.  Define the height of a matrix to be the max of the
    243     absolute values of the entries.
    244 
    245     Given Matrix A with n columns (self).
    246 
    247      0. Rescale input matrix A to have integer entries.  This does
    248         not change echelon form and makes reduction modulo lots of
    249         primes significantly easier if there were denominators.
    250         Henceforth we assume A has integer entries.
    251 
    252      1. Let c be a guess for the height of the echelon form.  E.g.,
    253         c=1000, e.g., if matrix is very sparse and application is to
    254         computing modular symbols.
    255 
    256      2. Let M = n * c * H(A) + 1,
    257         where n is the number of columns of A.
    258 
    259      3. List primes p_1, p_2, ..., such that the product of
    260         the p_i is at least M.
    261 
    262      4. Try to compute the rational reconstruction CRT echelon form
    263         of A mod the product of the p_i.  If rational
    264         reconstruction fails, compute 1 more echelon forms mod the
    265         next prime, and attempt again.  Make sure to keep the
    266         result of CRT on the primes from before, so we don't have
    267         to do that computation again.  Let E be this matrix.
    268 
    269      5. Compute the denominator d of E.
    270         Attempt to prove that result is correct by checking that
    271 
    272               H(d*E)*ncols(A)*H(A) < (prod of reduction primes)
    273 
    274         where H denotes the height.   If this fails, do step 4 with
    275         a few more primes.
    276 
    277     EXAMPLES:
    27861        sage: A = matrix(QQ, 3, 7, [1..21])
    27962        sage: sage.matrix.misc.matrix_rational_echelon_form_multimodular(A)
    280         [ 1  0 -1 -2 -3 -4 -5]
    281         [ 0  1  2  3  4  5  6]
    282         [ 0  0  0  0  0  0  0]
    283 
    284         sage: A = matrix(QQ, 3, 4, [0,0] + [1..9] + [-1/2^20])
    285         sage: sage.matrix.misc.matrix_rational_echelon_form_multimodular(A)
    286         [                1                 0                 0 -10485761/1048576]
    287         [                0                 1                 0  27262979/4194304]
    288         [                0                 0                 1                 2]
    289         sage: A.echelon_form()
    290         [                1                 0                 0 -10485761/1048576]
    291         [                0                 1                 0  27262979/4194304]
    292         [                0                 0                 1                 2]
     63        doctest:...: DeprecationWarning: (Since Sage Version 4.7) The matrix_rational_echelon_form_multimodular() function in sage.matrix.misc has been deprecated, use the _echelon_form_multimodular_rational() method from sage.matrix.matrix2.Matrix instead
    29364    """
    294 
    295     if proof is None:
    296         from sage.structure.proof.proof import get_flag
    297         proof = get_flag(proof, "linear_algebra")
    298 
    299     verbose("Multimodular echelon algorithm on %s x %s matrix"%(self._nrows, self._ncols), caller_name="multimod echelon")
    300     cdef Matrix E
    301     if self._nrows == 0 or self._ncols == 0:
    302         self.cache('in_echelon_form', True)
    303         self.cache('echelon_form', self)
    304         self.cache('pivots', [])
    305         return self
    306 
    307     B, _ = self._clear_denom()
    308 
    309     height = self.height()
    310     if height_guess is None:
    311         height_guess = 10000000*(height+100)
    312     tm = verbose("height_guess = %s"%height_guess, level=2, caller_name="multimod echelon")
    313 
    314     if proof:
    315         M = self._ncols * height_guess * height  +  1
    316     else:
    317         M = height_guess + 1
    318 
    319     if self.is_sparse():
    320         p = matrix_modn_sparse.MAX_MODULUS + 1
    321     else:
    322         p = matrix_modn_dense.MAX_MODULUS + 1
    323     X = []
    324     best_pivots = []       
    325     prod = 1
    326     problem = 0
    327     lifts = {}
    328     while True:
    329         p = previous_prime(p)
    330         while prod < M:
    331             problem = problem + 1
    332             if problem > 50:
    333                 verbose("echelon multi-modular possibly not converging?", caller_name="multimod echelon")
    334             t = verbose("echelon modulo p=%s (%.2f%% done)"%(
    335                        p, 100*float(len(str(prod))) / len(str(M))), level=2, caller_name="multimod echelon")
    336 
    337             # We use denoms=False, since we made self integral by calling clear_denom above.
    338             A = B._mod_int(p)
    339             t = verbose("time to reduce matrix mod p:",t, level=2, caller_name="multimod echelon")
    340             A.echelonize()
    341             t = verbose("time to put reduced matrix in echelon form:",t, level=2, caller_name="multimod echelon")
    342 
    343             # a worthwhile check / shortcut.
    344             if self._nrows == self._ncols and len(A.pivots()) == self._nrows:
    345                 verbose("done: the echelon form mod p is the identity matrix", caller_name="multimod echelon")
    346                 E = self.parent().identity_matrix()
    347                 E.cache('pivots', range(self._nrows))
    348                 E.cache('in_echelon_form', True)
    349                 self.cache('in_echelon_form', True)
    350                 self.cache('echelon_form', E)
    351                 self.cache('pivots', range(self._nrows))
    352                 return E
    353                
    354             c = cmp_pivots(best_pivots, A.pivots())
    355             if c <= 0:
    356                 best_pivots = A.pivots()
    357                 X.append(A)
    358                 prod = prod * p
    359             else:
    360                 # do not save A since it is bad.
    361                 verbose("Excluding this prime (bad pivots).", caller_name="multimod echelon")
    362             t = verbose("time for pivot compare", t, level=2, caller_name="multimod echelon")
    363             p = previous_prime(p)
    364         # Find set of best matrices.
    365         Y = []
    366         # recompute product, since may drop bad matrices           
    367         prod = 1
    368         t = verbose("now comparing pivots and dropping any bad ones", level=2, t=t, caller_name="multimod echelon")
    369         for i in range(len(X)):               
    370             if cmp_pivots(best_pivots, X[i].pivots()) <= 0:
    371                 p = X[i].base_ring().order()
    372                 if not lifts.has_key(p):
    373                     t0 = verbose("Lifting a good matrix", level=2, caller_name = "multimod echelon")
    374                     lift = X[i].lift()
    375                     lifts[p] = (lift, p)
    376                     verbose("Finished lift", level=2, caller_name= "multimod echelon", t=t0)
    377                 Y.append(lifts[p])
    378                 prod = prod * X[i].base_ring().order()
    379         verbose("finished comparing pivots", level=2, t=t, caller_name="multimod echelon")
    380         try:
    381             if len(Y) == 0:
    382                 raise ValueError("not enough primes")
    383             t = verbose("start crt linear combination", level=2, caller_name="multimod echelon")
    384             a = CRT_basis([w[1] for w in Y])
    385             t = verbose('got crt basis', level=2, t=t, caller_name="multimod echelon")
    386 
    387             # take the linear combination of the lifts of the elements
    388             # of Y times coefficients in a
    389             L = a[0]*(Y[0][0])
    390             assert Y[0][0].is_sparse() == L.is_sparse()
    391             for j in range(1,len(Y)):
    392                 L += a[j]*(Y[j][0])
    393             verbose("time to take linear combination of matrices over ZZ is",t, level=2, caller_name="multimod echelon")
    394             t = verbose("now doing rational reconstruction", level=2, caller_name="multimod echelon")
    395             E = L.rational_reconstruction(prod)
    396             L = 0  # free memory
    397             verbose('rational reconstruction completed', t, level=2, caller_name="multimod echelon")
    398         except ValueError, msg:
    399             verbose(msg, level=2)
    400             verbose("Not enough primes to do CRT lift; redoing with several more primes.", level=2, caller_name="multimod echelon")               
    401             M = prod * p*p*p
    402             continue
    403 
    404         if not proof:
    405             verbose("Not checking validity of result (since proof=False).", level=2, caller_name="multimod echelon")
    406             break
    407         d   = E.denominator()
    408         hdE = long((d*E).height())
    409         if hdE * self.ncols() * height < prod:
    410             verbose("Validity of result checked.", level=2, caller_name="multimod echelon")
    411             break
    412         verbose("Validity failed; trying again with more primes.", level=2, caller_name="multimod echelon")
    413         M = prod * p*p*p
    414     #end while
    415     verbose("total time",tm, level=2, caller_name="multimod echelon")
    416     self.cache('pivots', best_pivots)
    417     E.cache('pivots', best_pivots)
    418     return E
    419        
    420 
    421 ###########################
     65    deprecation('The matrix_rational_echelon_form_multimodular() function in sage.matrix.misc has been deprecated, use the _echelon_form_multimodular_rational() method from sage.matrix.matrix2.Matrix instead', 'Sage Version 4.7')
    42266
    42367def cmp_pivots(x,y):
     68    r"""
     69    This function is deprecated in favor of
     70    :func:`sage.matrix.matrix2.cmp_pivots` which behaves identically.
     71
     72    EXAMPLE::
     73
     74        sage: sage.matrix.misc.cmp_pivots([1,2,3], [4,5,6,7])
     75        doctest:...: DeprecationWarning: (Since Sage Version 4.7) The cmp_pivots() method in sage.matrix.misc has been deprecated, use the identical cmp_pivots() from sage.matrix.matrix2 instead
    42476    """
    425     Compare two sequences of pivot columns.
    426    
    427     If x is shorter than y, return -1, i.e., x < y, "not as good".
    428     If x is longer than y, then x > y, so "better" and return +1.
    429     If the length is the same, then x is better, i.e., x > y
    430     if the entries of x are correspondingly <= those of y with
    431     one being strictly less.
    432        
    433     INPUT:
    434         x, y -- list of integers
     77    deprecation('The cmp_pivots() method in sage.matrix.misc has been deprecated, use the identical cmp_pivots() from sage.matrix.matrix2 instead', 'Sage Version 4.7')
    43578
    436     EXAMPLES:
    437     We illustrate each of the above comparisons.
    438         sage: sage.matrix.misc.cmp_pivots([1,2,3], [4,5,6,7])
    439         -1
    440         sage: sage.matrix.misc.cmp_pivots([1,2,3,5], [4,5,6])
    441         1
    442         sage: sage.matrix.misc.cmp_pivots([1,2,4], [1,2,3])
    443         -1
    444         sage: sage.matrix.misc.cmp_pivots([1,2,3], [1,2,3])
    445         0
    446         sage: sage.matrix.misc.cmp_pivots([1,2,3], [1,2,4])
    447         1
    448     """
    449     if len(x) < len(y):
    450         return -1
    451     if len(x) > len(y):
    452         return 1
    453     if x < y:
    454         return 1
    455     elif x == y:
    456         return 0
    457     else:
    458         return -1
    459    
     79def hadamard_row_bound_mpfr(Matrix A):
     80    r"""
     81    This function is deprecated in favor of
     82    :meth:`sage.matrix.matrix2.Matrix._hadamard_row_bound_mpfr` which behaves identically.
    46083
     84    EXAMPLE::
    46185
    462 #######################################
    463    
    464    
    465 
    466 #######################################
    467 def hadamard_row_bound_mpfr(Matrix A):
    468     """
    469     Given a matrix A with entries that coerce to RR, compute the row
    470     Hadamard bound on the determinant.
    471 
    472     INPUT:
    473         A -- a matrix over RR
    474 
    475     OUTPUT:
    476         integer -- an integer n such that the absolute value of the
    477                    determinant of this matrix is at most $10^n$.
    478 
    479     EXAMPLES:
    480     We create a very large matrix, compute the row Hadamard bound,
    481     and also compute the row Hadamard bound of the transpose, which
    482     happens to be sharp.
    48386        sage: a = matrix(ZZ, 2, [2^10000,3^10000,2^50,3^19292])
    48487        sage: import sage.matrix.misc
    48588        sage: sage.matrix.misc.hadamard_row_bound_mpfr(a.change_ring(RR))
    486         13976
    487         sage: len(str(a.det()))
    488         12215
    489         sage: sage.matrix.misc.hadamard_row_bound_mpfr(a.transpose().change_ring(RR))
    490         12215
     89        doctest:...: DeprecationWarning: (Since Sage Version 4.7) The hadamard_row_bound_mpfr() function in sage.matrix.misc has been deprecated, use the identical _hadamard_row_bound_mpfr() method from sage.matrix.matrix2.Matrix instead
     90    """
     91    deprecation('The hadamard_row_bound_mpfr() function in sage.matrix.misc has been deprecated, use the identical _hadamard_row_bound_mpfr() method from sage.matrix.matrix2.Matrix instead', 'Sage Version 4.7')
    49192
    492     Note that in the above example using RDF would overflow:
    493         sage: b = a.change_ring(RDF)
    494         sage: b._hadamard_row_bound()
    495         Traceback (most recent call last):
    496         ...
    497         OverflowError: cannot convert float infinity to integer
    498     """
    499     if not is_RealField(A.base_ring()):
    500         raise TypeError("A must have base field an mpfr real field.")
    501 
    502     cdef RealNumber a, b
    503     cdef mpfr_t s, d, pr
    504     cdef Py_ssize_t i, j
    505 
    506     mpfr_init(s)
    507     mpfr_init(d)
    508     mpfr_init(pr)
    509     mpfr_set_si(d, 0, GMP_RNDU)
    510    
    511     for i from 0 <= i < A._nrows:
    512         mpfr_set_si(s, 0, GMP_RNDU)
    513         for j from 0 <= j < A._ncols:
    514             a = A.get_unsafe(i, j)
    515             mpfr_mul(pr, a.value, a.value, GMP_RNDU)
    516             mpfr_add(s, s, pr, GMP_RNDU)
    517         mpfr_log10(s, s, GMP_RNDU)
    518         mpfr_add(d, d, s, GMP_RNDU)
    519     b = a._new()
    520     mpfr_set(b.value, d, GMP_RNDU)
    521     b /= 2
    522     mpfr_clear(s)
    523     mpfr_clear(d)
    524     mpfr_clear(pr)
    525     return b.ceil()
    526