123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293 |
- #!/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)
|