Skip to content

Commit c09e1bc

Browse files
author
Release Manager
committed
sagemathgh-40007: Fixed issue in CRT_vectors where moduli are not allowed to be coprime. Fixes sagemath#39158. Fixed the issue by creating functionality for `CRT_basis` if the moduli are not coprime. A more in depth explanation on how `CRT_basis` was expanded can be found in [Here](sagemath#39158). ### 📝 Checklist <!-- Put an `x` in all the boxes that apply. --> - [X] The title is concise and informative. - [X] The description explains in detail what this PR is about. - [X] I have linked a relevant issue or discussion. - [X] I have created tests covering the changes. - [ ] I have updated the documentation and checked the documentation preview. ### ⌛ Dependencies <!-- List all open PRs that this PR logically depends on. For example, --> <!-- - sagemath#12345: short description why this is a dependency --> <!-- - sagemath#34567: ... --> URL: sagemath#40007 Reported by: Noel-Roemmele Reviewer(s): DaveWitteMorris, Noel-Roemmele
2 parents 332d23c + 4c17ea3 commit c09e1bc

File tree

1 file changed

+64
-20
lines changed

1 file changed

+64
-20
lines changed

src/sage/arith/misc.py

Lines changed: 64 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -3638,20 +3638,29 @@ def CRT_list(values, moduli=None):
36383638
return values[0] % moduli[0]
36393639

36403640

3641-
def CRT_basis(moduli):
3641+
def CRT_basis(moduli, *, require_coprime_moduli=True):
36423642
r"""
36433643
Return a CRT basis for the given moduli.
36443644
36453645
INPUT:
36463646
3647-
- ``moduli`` -- list of pairwise coprime moduli `m` which admit an
3647+
- ``moduli`` -- list of moduli `m` which admit an
36483648
extended Euclidean algorithm
36493649
3650+
- ``require_coprime_moduli`` -- boolean (default: ``True``); whether the moduli
3651+
must be pairwise coprime.
3652+
36503653
OUTPUT:
36513654
3652-
- a list of elements `a_i` of the same length as `m` such that
3653-
`a_i` is congruent to 1 modulo `m_i` and to 0 modulo `m_j` for
3654-
`j\not=i`.
3655+
- a list of integers `a_i` of the same length as `m` such that
3656+
if `r` is any list of integers of the same length as `m`, and we
3657+
let `x = \sum r_j a_j`, then `x \equiv r_i \pmod{m_i}` for all `i`
3658+
(if a solution of the system of congruences exists). When the
3659+
moduli are pairwise coprime, this implies that `a_i` is
3660+
congruent to 1 modulo `m_i` and to 0 modulo `m_j` for `j \neq i`.
3661+
3662+
- if ``require_coprime_moduli`` is ``False``, also returns a boolean value
3663+
that is ``True`` if the given moduli are pairwise coprime
36553664
36563665
EXAMPLES::
36573666
@@ -3674,22 +3683,44 @@ def CRT_basis(moduli):
36743683
n = len(moduli)
36753684
if n == 0:
36763685
return []
3677-
M = prod(moduli)
36783686
cs = []
3679-
for m in moduli:
3680-
Mm = M // m
3681-
d, _, v = xgcd(m, Mm)
3682-
if not d.is_one():
3683-
raise ValueError('moduli must be coprime')
3684-
cs.append((v * Mm) % M)
3685-
return cs
3687+
try:
3688+
M = prod(moduli)
3689+
for m in moduli:
3690+
Mm = M // m
3691+
d, _, v = xgcd(m, Mm)
3692+
if not d.is_one():
3693+
raise ValueError('moduli must be coprime')
3694+
cs.append((v * Mm) % M)
3695+
if require_coprime_moduli:
3696+
return cs
3697+
# also return a boolean flag to report that the moduli are coprime
3698+
return [cs, True]
3699+
except ValueError:
3700+
if require_coprime_moduli:
3701+
raise
3702+
e = [1]
3703+
M_i = moduli[0]
3704+
for i in range(1, n):
3705+
m_i = moduli[i]
3706+
d_i = gcd(M_i, m_i)
3707+
e_i = CRT(0, 1, M_i // d_i, m_i // d_i)
3708+
e.append(e_i)
3709+
M_i = M_i.lcm(m_i)
3710+
partial_prod_table = [1]
3711+
for i in range(1, n):
3712+
partial_prod_table.append((1 - e[-i]) * partial_prod_table[-1])
3713+
for i in range(n):
3714+
cs.append(e[i] * partial_prod_table[-i-1])
3715+
# also return a boolean flag to report that the moduli are not coprime
3716+
return [cs, False]
36863717

36873718

36883719
def CRT_vectors(X, moduli):
36893720
r"""
36903721
Vector form of the Chinese Remainder Theorem: given a list of integer
3691-
vectors `v_i` and a list of coprime moduli `m_i`, find a vector `w` such
3692-
that `w = v_i \pmod m_i` for all `i`.
3722+
vectors `v_i` and a list of moduli `m_i`, find a vector `w` such
3723+
that `w = v_i \pmod{m_i}` for all `i`.
36933724
36943725
This is more efficient than applying :func:`CRT` to each entry.
36953726
@@ -3708,18 +3739,31 @@ def CRT_vectors(X, moduli):
37083739
37093740
sage: CRT_vectors([vector(ZZ, [2,3,1]), Sequence([1,7,8], ZZ)], [8,9]) # needs sage.modules
37103741
[10, 43, 17]
3742+
3743+
``CRT_vectors`` also works for some non-coprime moduli::
3744+
3745+
sage: CRT_vectors([[6],[0]],[10, 4])
3746+
[16]
3747+
sage: CRT_vectors([[6],[0]],[10, 10])
3748+
Traceback (most recent call last):
3749+
...
3750+
ValueError: solution does not exist
37113751
"""
37123752
# First find the CRT basis:
37133753
if not X or len(X[0]) == 0:
37143754
return []
37153755
n = len(X)
37163756
if n != len(moduli):
37173757
raise ValueError("number of moduli must equal length of X")
3718-
a = CRT_basis(moduli)
3719-
modulus = prod(moduli)
3720-
return [sum(a[i] * X[i][j] for i in range(n)) % modulus
3721-
for j in range(len(X[0]))]
3722-
3758+
res = CRT_basis(moduli, require_coprime_moduli=False)
3759+
a = res[0]
3760+
modulus = LCM_list(moduli)
3761+
candidate = [sum(a[i] * X[i][j] for i in range(n)) % modulus
3762+
for j in range(len(X[0]))]
3763+
if not res[1] and any((X[i][j] - candidate[j]) % moduli[i] != 0 for i in range(n)
3764+
for j in range(len(X[i]))):
3765+
raise ValueError("solution does not exist")
3766+
return candidate
37233767

37243768
def binomial(x, m, **kwds):
37253769
r"""

0 commit comments

Comments
 (0)