Programa para resolver ecuaciones diofánticas

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.

Deja un comentario

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *

Este sitio usa Akismet para reducir el spam. Aprende cómo se procesan los datos de tus comentarios.