[New-bugs-announce] [issue39657] Bezout and Chinese Remainder Theorem in the Standard Library?

[Dennis Sweeney]

New submission from Dennis Sweeney <sweeney.dennis650 at gmail.com>:

Should something like the following go in the standard library, most likely in the math module? I know I had to use such a thing before pow(a, -1, b) worked, but Bezout is more general. And many of the easy stackoverflow implementations of CRT congruence-combining neglect the case where the divisors are not coprime, so that's an easy thing to miss. 

def bezout(a, b):
    """
    Given integers a and b, return a tuple (x, y, g),
    where x*a + y*b == g == gcd(a, b).
    """
    # Apply the Extended Euclidean Algorithm:
    # use the normal Euclidean Algorithm on the RHS
    # of the equations
    #     u1*a + v1*b == r1
    #     u2*a + v2*b == r2
    # But carry the LHS along for the ride.
    u1, v1, r1 = 1, 0, a
    u2, v2, r2 = 0, 1, b

    while r2:
        q = r1 // r2
        u1, u2 = u2, u1-q*u2
        v1, v2 = v2, v1-q*v2
        r1, r2 = r2, r1-q*r2
        assert u1*a + v1*b == r1
        assert u2*a + v2*b == r2

    if r1 < 0:
        u1, v1, r1 = -u1, -v1, -r1

    # a_coefficient, b_coefficient, gcd
    return (u1, v1, r1)

def crt(cong1, cong2):
    """
    Apply the Chinese Remainder Theorem:
        If there are any integers x such that 
        x == a1 (mod n1) and x == a2 (mod n2),
        then there are integers a and n such that the
        above congruences both hold iff x == a (mod n)
    Given two compatible congruences (a1, n1), (a2, n2),
    return a single congruence (a, n) that is equivalent
    to both of the given congruences at the same time.
    
    Not all congruences are compatible. For example, there
    are no solutions to x == 1 (mod 2) and x == 2 (mod 4).
    For congruences (a1, n1), (a2, n2) to be compatible, it
    is sufficient, but not necessary, that gcd(n1, n2) == 1.
    """
    a1, n1 = cong1
    a2, n2 = cong2
    c1, c2, g = bezout(n1, n2)
    assert n1*c1 + n2*c2 == g
    
    if (a1 - a2) % g != 0:
        raise ValueError(f"Incompatible congruences {cong1} and {cong2}.")
    
    lcm = n1 // g * n2
    rem = (a1*c2*n2 + a2*c1*n1)//g
    return rem % lcm, lcm

assert crt((1,4),(2,3)) == (5, 12)
assert crt((1,6),(7,4)) == (7, 12)

----------
components: Library (Lib)
messages: 362106
nosy: Dennis Sweeney
priority: normal
severity: normal
status: open
title: Bezout and Chinese Remainder Theorem in the Standard Library?
type: enhancement
versions: Python 3.9

_______________________________________
Python tracker <report at bugs.python.org>
<https://bugs.python.org/issue39657>
_______________________________________