En el artículo anterior expliqué cómo se resuelven las ecuaciones diofánticas y su relación con las ecuaciones de congruencia. En la presente entrada vamos a ver un programa en Python que las resuelve.

#! /usr/bin/env python
# -*- coding: utf-8 -*-
# Resuelve ecuaciones diofánticas tipo ax + by = c

import sys
from sys import argv

def extendedEuclideanAlgorithm(old_r, r):
    negative = False
    s, old_t = 0, 0
    old_s, t = 1, 1

    if (r < 0):
        r = abs(r)
        negative = True
        
    while r > 0:
        q = old_r / r
        #MCD:
        r, old_r = old_r - q * r, r
        #Coeficiente s:
        s, old_s = old_s - q * s, s
        #Coeficiente t:
        t, old_t = old_t - q * t, t
        
    if negative:
        old_t = old_t * -1
        
    return old_r, old_s, old_t

a = long(argv[1])
b = long(argv[2])
c = long(argv[3])

mcd, s, t = extendedEuclideanAlgorithm(a, b)
if c % mcd == 0:
    a1, b1, c1 = -a / mcd, b / mcd, c / mcd
    x1, y1 = s * c1, t * c1
    print "x = {0}{1:+d}k" . format(x1, b1)
    print "y = {0}{1:+d}k" . format(y1, a1)
else:
    print "No tiene solución"

Para calcuar 23x -4y = 11 hacemos:

vic@LESBIAN:~/mates$ ./diofanticas.py 23 -4 11
x = -11-4k
y = -66-23k

Aplicando el concepto de clase de equivalencia, tal y como se explica en el anterior artículo, podemos computar la forma paramétrica más simplificada como muestra el siguiente programa:

#! /usr/bin/env python
# -*- coding: utf-8 -*-
# Resuelve ecuaciones diofánticas tipo ax + by = c
# El primero par x e y es la forma parmétrica más simplificada
 
import sys
from sys import argv

def extendedEuclideanAlgorithm(old_r, r):
    negative = False
    s, old_t = 0, 0
    old_s, t = 1, 1

if (r < 0):
    r = abs(r)
    negative = True

while r > 0:
    q = old_r / r
    #MCD:
    r, old_r = old_r - q * r, r
    #Coeficiente s:
    s, old_s = old_s - q * s, s
    #Coeficiente t:
    t, old_t = old_t - q * t, t

if negative:
    old_t = old_t * -1

return old_r, old_s, old_t

a = long(argv[1])
b = long(argv[2])
c = long(argv[3])

mcd, s, t = extendedEuclideanAlgorithm(a, b)
if c % mcd == 0:
    a1, b1, c1 = a / mcd, b / mcd, c / mcd
    x1, y1 = s * c1, t * c1
    # Uso abs() pues Python no hace la división Euclídea con cociente negativo
    equivClass = x1 % abs(b1)
    print "x = {0}{1:+d}k" . format(equivClass, b1)
    print "y = {0}{1:+d}k" . format((c1 - (a1 * equivClass)) / b1, -a1)
    print "x = {0}{1:+d}k" . format(x1, b1)
    print "y = {0}{1:+d}k" . format(y1, -a1)
else:
    print "No tiene solución"

Este es el resultado de diferentes ejecuciones:

vic@LESBIAN:~/mates$ ./diofanticas.py 4 7 29
x = 2+7k
y = 3-4k
x = 58+7k
y = -29-4k

vic@LESBIAN:~/mates$ ./diofanticas.py 23 -4 11
x = 1-4k
y = 3-23k
x = -11-4k
y = -66-23k

El uso de la función abs(), que nos devuelve el valor absoluto, es debido a que Python no hace la división Euclídea cuando el cociente es negativo. Los dos pasos:

• q = old_r / r
• old_r = old_r – q * r

Se podrían unificar con la función divmod(), pero me parece menos claro para quien no conoce el lenguaje, siendo los dos pasos más parecidos al pseudocódigo.