1
1
"""
2
2
Misc matrix algorithms
3
-
4
- Code goes here mainly when it needs access to the internal structure
5
- of several classes, and we want to avoid circular cimports.
6
-
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
3
"""
12
4
13
5
from cysignals.signals cimport sig_check
14
6
15
- cimport sage.rings.abc
16
-
17
7
from sage.arith.misc import CRT_basis, previous_prime
18
8
from sage.arith.rational_reconstruction cimport mpq_rational_reconstruction
19
9
from sage.data_structures.binary_search cimport *
20
10
from sage.ext.mod_int cimport *
21
- from sage.libs.flint.fmpq cimport fmpq_set_mpq, fmpq_canonicalise
22
- from sage.libs.flint.fmpq_mat cimport fmpq_mat_entry_num, fmpq_mat_entry_den, fmpq_mat_entry
23
- from sage.libs.flint.fmpz cimport fmpz_set_mpz, fmpz_one
24
11
from sage.libs.gmp.mpq cimport *
25
12
from sage.libs.gmp.mpz cimport *
26
- from sage.libs.mpfr cimport *
13
+ from sage.misc.lazy_import import LazyImport
27
14
from sage.misc.verbose import verbose
28
15
from sage.modules.vector_integer_sparse cimport *
29
16
from sage.modules.vector_modn_sparse cimport *
30
17
from sage.modules.vector_rational_sparse cimport *
31
18
from sage.rings.integer cimport Integer
32
19
from sage.rings.rational_field import QQ
33
- from sage.rings.real_mpfr cimport RealNumber
34
20
35
21
from .matrix0 cimport Matrix
36
- from .matrix_integer_dense cimport Matrix_integer_dense
37
22
from .matrix_integer_sparse cimport Matrix_integer_sparse
38
- from .matrix_rational_dense cimport Matrix_rational_dense
39
23
from .matrix_rational_sparse cimport Matrix_rational_sparse
40
24
41
-
42
- def matrix_integer_dense_rational_reconstruction (Matrix_integer_dense A , Integer N ):
43
- """
44
- Given a matrix over the integers and an integer modulus, do
45
- rational reconstruction on all entries of the matrix, viewed as
46
- numbers mod N. This is done efficiently by assuming there is a
47
- large common factor dividing the denominators.
48
-
49
- INPUT:
50
-
51
- A -- matrix
52
- N -- an integer
53
-
54
- EXAMPLES::
55
-
56
- sage: B = ((matrix(ZZ, 3,4, [1,2,3,-4,7,2,18,3,4,3,4,5])/3)%500).change_ring(ZZ)
57
- sage: sage.matrix.misc.matrix_integer_dense_rational_reconstruction(B, 500)
58
- [ 1/3 2/3 1 -4/3]
59
- [ 7/3 2/3 6 1]
60
- [ 4/3 1 4/3 5/3]
61
-
62
- TESTS:
63
-
64
- Check that :trac:`9345` is fixed::
65
-
66
- sage: A = random_matrix(ZZ, 3)
67
- sage: sage.matrix.misc.matrix_integer_dense_rational_reconstruction(A, 0)
68
- Traceback (most recent call last):
69
- ...
70
- ZeroDivisionError: The modulus cannot be zero
71
- """
72
- if not N:
73
- raise ZeroDivisionError (" The modulus cannot be zero" )
74
- cdef Matrix_rational_dense R
75
- R = Matrix_rational_dense.__new__ (Matrix_rational_dense,
76
- A.parent().change_ring(QQ), 0 ,0 ,0 )
77
-
78
- cdef mpz_t a, bnd, other_bnd, denom, tmp
79
- cdef mpq_t qtmp
80
- cdef Integer _bnd
81
- cdef Py_ssize_t i, j
82
- cdef int do_it
83
-
84
- mpz_init_set_si(denom, 1 )
85
- mpz_init(a)
86
- mpz_init(tmp)
87
- mpz_init(other_bnd)
88
- mpq_init(qtmp)
89
-
90
- _bnd = (N// 2 ).isqrt()
91
- mpz_init_set(bnd, _bnd.value)
92
- mpz_sub(other_bnd, N.value, bnd)
93
-
94
- for i in range (A._nrows):
95
- for j in range (A._ncols):
96
- sig_check()
97
- A.get_unsafe_mpz(i, j, a)
98
- if mpz_cmp_ui(denom, 1 ) != 0 :
99
- mpz_mul(a, a, denom)
100
- mpz_fdiv_r(a, a, N.value)
101
- do_it = 0
102
- if mpz_cmp(a, bnd) <= 0 :
103
- do_it = 1
104
- elif mpz_cmp(a, other_bnd) >= 0 :
105
- mpz_sub(a, a, N.value)
106
- do_it = 1
107
- if do_it:
108
- fmpz_set_mpz(fmpq_mat_entry_num(R._matrix, i, j), a)
109
- if mpz_cmp_ui(denom, 1 ) != 0 :
110
- fmpz_set_mpz(fmpq_mat_entry_den(R._matrix, i, j), denom)
111
- fmpq_canonicalise(fmpq_mat_entry(R._matrix, i, j))
112
- else :
113
- fmpz_one(fmpq_mat_entry_den(R._matrix, i, j))
114
- else :
115
- # Otherwise have to do it the hard way
116
- A.get_unsafe_mpz(i, j, tmp)
117
- mpq_rational_reconstruction(qtmp, tmp, N.value)
118
- mpz_lcm(denom, denom, mpq_denref(qtmp))
119
- fmpq_set_mpq(fmpq_mat_entry(R._matrix, i, j), qtmp)
120
-
121
- mpz_clear(denom)
122
- mpz_clear(a)
123
- mpz_clear(tmp)
124
- mpz_clear(other_bnd)
125
- mpz_clear(bnd)
126
- mpq_clear(qtmp)
127
-
128
- return R
25
+ matrix_integer_dense_rational_reconstruction = \
26
+ LazyImport(' sage.matrix.misc_flint' , ' matrix_integer_dense_rational_reconstruction' ,
27
+ deprecation = 35758 )
28
+ hadamard_row_bound_mpfr = \
29
+ LazyImport(' sage.matrix.misc_mpfr' , ' hadamard_row_bound_mpfr' ,
30
+ deprecation = 35758 )
129
31
130
32
131
33
def matrix_integer_sparse_rational_reconstruction (Matrix_integer_sparse A , Integer N ):
132
- """
34
+ r """
133
35
Given a sparse matrix over the integers and an integer modulus, do
134
36
rational reconstruction on all entries of the matrix, viewed as
135
- numbers mod N .
37
+ numbers mod `N` .
136
38
137
39
EXAMPLES::
138
40
139
41
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)
140
- sage: sage.matrix.misc.matrix_integer_sparse_rational_reconstruction(A, 500)
42
+ sage: from sage. matrix. misc import matrix_integer_sparse_rational_reconstruction
43
+ sage: matrix_integer_sparse_rational_reconstruction( A, 500)
141
44
[1/3 2 3 -4 ]
142
45
[ 7 2 2 3 ]
143
46
[ 4 3 4 5/7 ]
@@ -424,32 +327,33 @@ def matrix_rational_echelon_form_multimodular(Matrix self, height_guess=None, pr
424
327
425
328
426
329
def cmp_pivots (x , y ):
427
- """
330
+ r """
428
331
Compare two sequences of pivot columns.
429
332
430
- If x is shorter than y , return -1 , i.e., x < y, "not as good".
431
- If x is longer than y , then x > y, so "better" and return +1 .
432
- If the length is the same, then x is better, i.e., x > y
433
- if the entries of x are correspondingly <= those of y with
333
+ If `x` is shorter than `y` , return `-1` , i. e. , ` x < y` , "not as good".
334
+ If `x` is longer than `y` , then ` x > y` , so "better" and return ` + 1` .
335
+ If the length is the same, then `x` is better, i. e. , ` x > y`
336
+ if the entries of `x` are correspondingly ` \l eq` those of `y` with
434
337
one being strictly less.
435
338
436
339
INPUT:
437
340
438
- - x, y -- lists or tuples of integers
341
+ - ``x``, ``y`` -- lists or tuples of integers
439
342
440
343
EXAMPLES:
441
344
442
345
We illustrate each of the above comparisons. ::
443
346
444
- sage: sage.matrix.misc.cmp_pivots([1,2,3], [4,5,6,7])
347
+ sage: from sage. matrix. misc import cmp_pivots
348
+ sage: cmp_pivots( [1,2,3 ], [4,5,6,7 ])
445
349
-1
446
- sage: sage.matrix.misc. cmp_pivots([1,2,3,5], [4,5,6])
350
+ sage: cmp_pivots( [1,2,3,5 ], [4,5,6 ])
447
351
1
448
- sage: sage.matrix.misc. cmp_pivots([1,2,4], [1,2,3])
352
+ sage: cmp_pivots( [1,2,4 ], [1,2,3 ])
449
353
-1
450
- sage: sage.matrix.misc. cmp_pivots([1,2,3], [1,2,3])
354
+ sage: cmp_pivots( [1,2,3 ], [1,2,3 ])
451
355
0
452
- sage: sage.matrix.misc. cmp_pivots([1,2,3], [1,2,4])
356
+ sage: cmp_pivots( [1,2,3 ], [1,2,4 ])
453
357
1
454
358
"""
455
359
x = tuple (x)
@@ -464,70 +368,3 @@ def cmp_pivots(x, y):
464
368
return 0
465
369
else :
466
370
return - 1
467
-
468
-
469
- def hadamard_row_bound_mpfr (Matrix A ):
470
- """
471
- Given a matrix A with entries that coerce to RR, compute the row
472
- Hadamard bound on the determinant.
473
-
474
- INPUT:
475
-
476
- A -- a matrix over RR
477
-
478
- OUTPUT:
479
-
480
- integer -- an integer n such that the absolute value of the
481
- determinant of this matrix is at most $10^n$.
482
-
483
- EXAMPLES:
484
-
485
- We create a very large matrix, compute the row Hadamard bound,
486
- and also compute the row Hadamard bound of the transpose, which
487
- happens to be sharp. ::
488
-
489
- sage: a = matrix(ZZ, 2, [2^10000,3^10000,2^50,3^19292])
490
- sage: import sage.matrix.misc
491
- sage: sage.matrix.misc.hadamard_row_bound_mpfr(a.change_ring(RR))
492
- 13976
493
- sage: len(str(a.det()))
494
- 12215
495
- sage: sage.matrix.misc.hadamard_row_bound_mpfr(a.transpose().change_ring(RR))
496
- 12215
497
-
498
- Note that in the above example using RDF would overflow::
499
-
500
- sage: b = a.change_ring(RDF)
501
- sage: b._hadamard_row_bound()
502
- Traceback (most recent call last):
503
- ...
504
- OverflowError: cannot convert float infinity to integer
505
- """
506
- if not isinstance (A.base_ring(), sage.rings.abc.RealField):
507
- raise TypeError (" A must have base field an mpfr real field." )
508
-
509
- cdef RealNumber a, b
510
- cdef mpfr_t s, d, pr
511
- cdef Py_ssize_t i, j
512
-
513
- mpfr_init(s)
514
- mpfr_init(d)
515
- mpfr_init(pr)
516
- mpfr_set_si(d, 0 , MPFR_RNDU)
517
-
518
- for i in range (A._nrows):
519
- mpfr_set_si(s, 0 , MPFR_RNDU)
520
- for j in range (A._ncols):
521
- sig_check()
522
- a = A.get_unsafe(i, j)
523
- mpfr_mul(pr, a.value, a.value, MPFR_RNDU)
524
- mpfr_add(s, s, pr, MPFR_RNDU)
525
- mpfr_log10(s, s, MPFR_RNDU)
526
- mpfr_add(d, d, s, MPFR_RNDU)
527
- b = a._new()
528
- mpfr_set(b.value, d, MPFR_RNDU)
529
- b /= 2
530
- mpfr_clear(s)
531
- mpfr_clear(d)
532
- mpfr_clear(pr)
533
- return b.ceil()
0 commit comments