This issue tracker has been migrated to GitHub, and is currently read-only.
For more information, see the GitHub FAQs in the Python's Developer Guide.

Author Dennis Sweeney
Recipients Dennis Sweeney
Date 2020-02-16.21:40:38
SpamBayes Score -1.0
Marked as misclassified Yes
Message-id <1581889239.13.0.672135436609.issue39657@roundup.psfhosted.org>
In-reply-to
Content
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)
History
Date User Action Args
2020-02-16 21:40:39Dennis Sweeneysetrecipients: + Dennis Sweeney
2020-02-16 21:40:39Dennis Sweeneysetmessageid: <1581889239.13.0.672135436609.issue39657@roundup.psfhosted.org>
2020-02-16 21:40:39Dennis Sweeneylinkissue39657 messages
2020-02-16 21:40:38Dennis Sweeneycreate