def egcd(a, b):
if not (a and b):
return (1 if a else 0, b if b else 0, a + b)
d, r = divmod(a, b)
x1, y1, g = egcd(b, r)
return (y1, x1 - y1 * d, g)
def solve(a, b, c):
"""
return (x, y) s.t. ax + by = c
x >= 0, y >= 0
"""
x1, y1, g = egcd(a, b)
if c % g != 0:
raise ValueError("Impossible")
x = x1 * c / g
y = y1 * c / g
swapped = False
if x < 0:
swapped = True
x, y = y, x
a, b = b, a
if y < 0:
u = a / g
v = b / g
k = x / v
x -= k * v
y += k * u
if y < 0:
raise ValueError("Impossible for positive x, y")
return (y, x) if swapped else (x, y)
print solve(18, 15, 198)
CgpkZWYgZWdjZChhLCBiKToKICAgIGlmIG5vdCAoYSBhbmQgYik6CiAgICAgICAgcmV0dXJuICgxIGlmIGEgZWxzZSAwLCBiIGlmIGIgZWxzZSAwLCBhICsgYikKICAgIGQsIHIgPSBkaXZtb2QoYSwgYikKICAgIHgxLCB5MSwgZyA9IGVnY2QoYiwgcikKICAgIHJldHVybiAoeTEsIHgxIC0geTEgKiBkLCBnKQoKZGVmIHNvbHZlKGEsIGIsIGMpOgogICAgIiIiCiAgICByZXR1cm4gKHgsIHkpIHMudC4gYXggKyBieSA9IGMKICAgIHggPj0gMCwgeSA+PSAwCiAgICAiIiIKICAgIHgxLCB5MSwgZyA9IGVnY2QoYSwgYikKICAgIGlmIGMgJSBnICE9IDA6CiAgICAgICAgcmFpc2UgVmFsdWVFcnJvcigiSW1wb3NzaWJsZSIpCgogICAgeCA9IHgxICogYyAvIGcKICAgIHkgPSB5MSAqIGMgLyBnCgogICAgc3dhcHBlZCA9IEZhbHNlCgogICAgaWYgeCA8IDA6CiAgICAgICAgc3dhcHBlZCA9IFRydWUKICAgICAgICB4LCB5ID0geSwgeAogICAgICAgIGEsIGIgPSBiLCBhIAoKICAgIGlmIHkgPCAwOgogICAgICAgIHUgPSBhIC8gZwogICAgICAgIHYgPSBiIC8gZwoKICAgICAgICBrID0geCAvIHYKCiAgICAgICAgeCAtPSBrICogdgogICAgICAgIHkgKz0gayAqIHUKICAgICAgICBpZiB5IDwgMDoKICAgICAgICAJcmFpc2UgVmFsdWVFcnJvcigiSW1wb3NzaWJsZSBmb3IgcG9zaXRpdmUgeCwgeSIpCiAgICByZXR1cm4gKHksIHgpIGlmIHN3YXBwZWQgZWxzZSAoeCwgeSkKCgpwcmludCBzb2x2ZSgxOCwgMTUsICAxOTgpCg==