sidef-scripts/Math/linear_diophantine_equation_invmod_search.sf

#!/usr/bin/ruby

# Find the smallest solution in positive integers to the equation:
#
#   a*x - b*y = n
#
# where a,b,n are given and n is an intermediate value in Euclid's GCD(a,b) algorithm.

# See also:
#   https://en.wikipedia.org/wiki/Diophantine_equation
#   https://mathworld.wolfram.com/DiophantineEquation.html

func modular_inverse_search (a, n, z) {

    var (u, w) = (1, 0)
    var (q, r) = (0, 0)

    var c = n

    while (c != 0) {
        (q, r) = divmod(a, c)
        (a, c) = (c, r)
        (u, w) = (w, u - q*w)
        break if (a == z)
    }

    u += n if (u < 0)

    return u
}

func solve(a, b, c) {

    var x = modular_inverse_search(a, b, c)
    var y = ((a*x - c) / b)

    return (x, y)
}

var tests = [
    [79, 23, 1],
    [97, 43, 1],
    [43, 97, 1],
    [55, 28, 1],
    [42, 22, 2],
    [79, 23, 10],
]

for a,b,n in tests {
    var (x, y) = solve(a, b, n)

    assert(x.is_int)
    assert(b.is_int)

    assert_eq(a*x - b*y, n)
    printf("#{a}*x - #{b}*y = %2s  -->  (x, y) = (%2s, %2s)\n", n, x, y)
}

say "\n>> Extra tests:"

var a = 43*97
var b = 41*57

for n in ([1, 7, 8, 23, 31, 147, 178, 325, 503, 1834]) {
    var (x, y) = solve(a, b, n)

    assert(x.is_int)
    assert(y.is_int)

    assert_eq(a*x - b*y, n)
    printf("#{a}*x - #{b}*y = %4s  -->  (x, y) = (%5s, %5s)\n", n, x, y)
}

__END__
79*x - 23*y =  1  -->  (x, y) = ( 7, 24)
97*x - 43*y =  1  -->  (x, y) = ( 4,  9)
43*x - 97*y =  1  -->  (x, y) = (88, 39)
55*x - 28*y =  1  -->  (x, y) = (27, 53)
42*x - 22*y =  2  -->  (x, y) = (21, 40)
79*x - 23*y = 10  -->  (x, y) = ( 1,  3)

>> Extra tests:
4171*x - 2337*y =    1  -->  (x, y) = ( 2035,  3632)
4171*x - 2337*y =    7  -->  (x, y) = (  223,   398)
4171*x - 2337*y =    8  -->  (x, y) = ( 2258,  4030)
4171*x - 2337*y =   23  -->  (x, y) = (   65,   116)
4171*x - 2337*y =   31  -->  (x, y) = ( 2323,  4146)
4171*x - 2337*y =  147  -->  (x, y) = (    9,    16)
4171*x - 2337*y =  178  -->  (x, y) = ( 2332,  4162)
4171*x - 2337*y =  325  -->  (x, y) = (    4,     7)
4171*x - 2337*y =  503  -->  (x, y) = ( 2336,  4169)
4171*x - 2337*y = 1834  -->  (x, y) = (    1,     1)