|
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=99999) |
| 28 | +hadamard_row_bound_mpfr = \ |
| 29 | + LazyImport('sage.matrix.misc_mpfr', 'hadamard_row_bound_mpfr', |
| 30 | + deprecation=99999) |
129 | 31 |
|
130 | 32 |
|
131 | 33 | def matrix_integer_sparse_rational_reconstruction(Matrix_integer_sparse A, Integer N):
|
@@ -464,70 +366,3 @@ def cmp_pivots(x, y):
|
464 | 366 | return 0
|
465 | 367 | else:
|
466 | 368 | 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