Programa para resolver ecuaciones de congruencia lineal

Una ecuación de congruencia es una congruencia entre expresiones literales, los valores que la satisfacen son sus soluciones. Por ejemplo 13x ≡ 5 (mod 77) es una ecuación de congruencia.

Sus soluciones nos las da x ≡ 30 (mod 77) Es decir, todos los infinitos x que al dividirlos por 77 dan el mismo resto que 30/77 Dicho resto es 30 y todos los x serían múltiplos de 77· + 30: 30, 107, 184, 261…

Si reemplazamos estos x en la ecuación de congruencia lineal inicial veremos que son soluciones:

13*30 ≡ 30 (mod 77) => Resto 5. Así como dan el mismo resto:

13*107 ≡ 30 (mod 77)

13*184 ≡ 30 (mod 77)

Para poder encontrar la solución de una congruencia lineal en primer lugar hay que saber hallar el inverso de a (mod m) Este inverso, a’, es un número tal que a’*a ≡ 1 (mod m)

Por ejemplo, el inverso de 6 (mod 11) es 2: 2*6 ≡ 1 (mod 11)

Teorema: La condición necesaria y suficiente para que dos números a y b sean congruentes respecto al módulo m es que su diferencia sea un múltiplo de m, es decir, m | (a -b)

Por lo tanto en el ejemplo sería 11 | (a’ * 6 – 1) => 6a’ – 1 = 11· , 6a’ = 11· + 1 Como ha de ser un múltiple de 11 + 1 vemos fácilmente que a’ = 2 cumple la condición.

Si para hallar el inverso de 2 (mod 8) hacemos: 8 | (a’ * 2 – 1) Vemos que no hay solución pues 2a’ -1 siempre será un número impar mientras que 8 es par. Por lo tanto no siempre existe un inverso. El siguiente teorema nos permite saber cuando existe el inverso:

La condición necesaria y suficiente para que un entero a posea un inverso módulo m, m > 1, es que mcd(a, m) = 1 Además, si existe tal inverso, ese inverso es único módulo m.

Esto nos resulta familiar con el teorema de Bezout explicado en el artículo anterior. Por ejemplo, para calcular el inverso de 111 (mod 250) primero calculamos mcd(111, 250) = 1 Por lo tanto sabemos que existe un inverso a’ y que habrá un s y un t tal que:

111 *s + 250 * t = 1

Aplicando el algoritmo del artículo obtenemos s = -9 y t = 11 y por lo tanto a’ = -9. Así es:

-9 * 111 ≡ 1 (mod 250)

  • -999 dividido entre 250 da el mismo resto que 1 entre 250, r = 1
  • 250 | (-9 * 111 -1) => 250· = -9 * 111 – 1 => Efectivamente -1000 es múltiplo de 250

Ahora bien, existe un teorema que matiza el anterior:

Existe un entero positivo único a’ menor que m que es un inverso de a (mod m) y cualquier otro inverso de a (mod m) será congruente con a’ (mod m)

Como dijimos en el anterior artículo, el par s y t no es único. La expresión general que nos da todos los coeficientes es:

1 =111(-9 + 250h) + 250(4 – 111h) ∀ h ∈ Z

Para h = 1 tendríamos:

1 = 111*241 + 250*(-107)

Como s = 241 y t = -107 tendremos 111*241 ≡ 1 (mod 250) Por lo tanto es 241 el inverso positivo menor que 250.

Si por un lado tenemos la ecuación de congruencia lineal ax ≡ b (mod m) y esta tiene un a’ ∈ Z tal que a’*a = 1 (mod m) Por otro lado, como cualquier número es congruente consigo mismo módulo m podemos crear la congruencia a’ ≡ a’ (mod m).

Un importante teorema dice:

Sean a, b, c y d enteros y m un entero positivo. Si es ab (mod m) y cd (mod m), entonces también es ac ≡ bd (mod m).

Con este teorema, de las dos congruencias anteriores resulta:

ax * a’ ≡ a’ * b (mod m)

Como a’ * a (mod m) = 1 Queda:

1 * x ≡ a’ * b (mod m)

Si existe el inverso modular, es decir, si m y a son primos relativos, o lo que es lo mismo, si mcd(m, a) = 1, todo x congruente con a’ * b (mod m) será solución a la ecuación de congruencia.

Vamos a calcular la solución de 3x ≡ 5 (mod 13):

El programa del artículo anterior, para a y m, es decir, 3 y 13 nos devuelve:

m.c.d. = 1

Bezout s = -4

Bezout t = 1

Como a y m son primos entre si, sin más cálculos podemos ver que el inverso de 3 (mod 13) es 9, pero procedamos como hemos visto. Si el inverso es -4 necesitamos encontrar el positivo, busquemos para ello otro par s y t con h = 1:

1 = 3*(-4 + 13*1) + 13*(1 – 3*1) = 3*9 + 13*(-2)

Por esta vía también obtenemos a’ = 9. Todo x congruente con 9*5 (mod 13) será solución. Como 45 (mod 13) = 6, todos los múltiples de 13· + 6 serán las soluciones a 3x ≡ 5 (mod 13). Por lo tanto, x ≡ 6 (mod 13) es una congruencia equivalente más simple.

Vamos a calcular la solución de 5x ≡ 2 (mod 6):

El programa nos devuelve:

m.c.d. =  1
Bezout s =  -1
Bezout t =  1

Como a y m son primos entre si, tendremos que el inverso a’ es -1. Busquemos el positivo:

1 = 5 * (-1 + 6*1) + 6*(1 – 5*1) = 5 * 5 + 6 * -4

El inverso positivo menor que m es 5.

5*2 (mod 6), 10 (mod 6), x ≡ 10 (mod 6) Los infinitos x congruentes con 10 módulo 6 serán solución de 5x ≡ 2 (mod 6). Es decir, todos los múltiples 6· + 4 Por lo tanto, es equivalente x ≡ 4 (mod 6)

Cuando a y m no son primos entre si, es decir mcd(a, m) > 1 y b es un múltiplo del mcd, es decir mcd | b, la congruencia tiene solución. Si queremos entender el proceso para hallar las soluciones y no limitarnos a aplicar fórmulas será necesario introducir los conjuntos de restos y las clases de equivalencia, pero para ello será necesario un nuevo artículo.

El siguiente programa en Python calcula las ecuaciones de congruencia que se han tratado en este artículo.

#! /usr/bin/env python
# -*- coding: utf-8 -*-
# Resuelve ecuaciones de congruencia tipo ax ≡ b (mod m)

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])
m = long(argv[3])
#Para las soluciones no necesitaremos t pero lo recogemos para evitar el error.
mcd, s, t = extendedEuclideanAlgorithm(a, m)
if mcd == 1:
    inverse = (s * b) % m
    print "x ≡ {0} (mod {1})" . format(inverse, m)
elif b%mcd == 0:
    print "Ver siguiente artículo"
else:
    print "No tiene solución"

Para resolver la ecuación 3x ≡ 5 (mod 13) haríamos:

vic@LESBIAN:~/mates$ ./inversoModulo.py 3 5 13
x ≡ 6 (mod 13)

Añadido el 11/03/2018:

Casi un año después de su publicación, observo en Google Analytics que este artículo tiene un porcentaje de rebote muy elevado. Tal vez este sea debido a que la mayoría de visitantes estén buscando simplemente un programa que resuelva congruencias lineales. Si es tu caso, para ejecutar el programa que aparece al final de esta entrada, puedes instalar fácilmente el intérprete de Python 2 en tu ordenador, es muy ligero y existe para todas las plataformas. Como alternativa, dispones de intérpretes en línea (1, 2) donde puedes copiar y pegar el código; después deberás modificarlo ligeramente porque los parámetros no serán argumentos pasados a través de la línea de comandos.


Cálculo del mcd y la identidad de Bezout

Editado el 11/03/2018:

Casi un año después de su publicación, observo en Google Analytics que este artículo tiene un porcentaje de rebote muy elevado. Tal vez este sea debido a que la mayoría de visitantes estén buscando simplemente un programa que calcule el mcd de dos números y sus coeficientes de Bezout. Si es tu caso, puedes instalar fácilmente el intérprete de Python 2 en tu ordenador, es muy ligero y existe para todas las plataformas. Como alternativa, dispones de intérpretes en línea (1, 2) donde puedes copiar y pegar el código; después deberás modificarlo ligeramente porque los parámetros no serán argumentos pasados a través de la línea de comandos.


He hecho un breve programa en Python que calcula el máximo común divisor (mcd) de dos números pasados como parámetros y sus coeficientes de Bezout. El teorema de Bezout nos dice que dados dos números a y b, cuyo mcd(a, b) = d existen dos números enteros s y t tales que:

a*s + b*t = d

Estos dos números s y t son los llamados coeficientes de Bezout. Este par s y t no tiene por qué ser único, pueden haber otros s y t que cumplan la igualdad. La identidad de Bezout no es una mera curiosidad matemática sino que tiene muchas aplicaciones: desde calcular el inverso de a módulo m, lo que permite resolver ecuaciones de congruencia lineal, hasta resolver ecuaciones diofánticas.

El máximo común divisor se puede calcular a mano como nos enseñaron en la escuela mediante el algoritmo de Euclides. Los coeficientes de Bezout los podemos calcular a mano mediante el algoritmo extendido de Euclides, pero a la hora de programar es más fácil implementar el siguiente algoritmo explicado en la Wikipedia.

Todos estos valores se calculan mediante la división euclídea y esta impone la condición 0 >= r < |d| Es decir, el resto debe ser mayor o igual a 0 y menor que el valor absoluto del divisor. Como expliqué en un artículo anterior, Python no hace la división euclídea cuando el divisor es negativo, por lo tanto en el código hago una pequeña manipulación al signo del segundo número para obtener los coeficientes de Bezout con los signos correctos.

#! /usr/bin/env python
# -*- coding: utf-8 -*-

import sys
from sys import argv

old_r = long(argv[1])
r = long(argv[2])
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
    
print "m.c.d. = ", old_r
print "Bezout s = ", old_s
print "Bezout t = ", old_t

Ejemplo de su ejecución desde la terminal:

vic@LESBIAN:~/mates$ ./bezoutYmcd.py 23 -4
m.c.d. =  1
Bezout s =  -1
Bezout t =  -6

Si convertimos este programa en una función tendremos una herramienta útil para calcular el inverso de a módulo m, resolver ecuaciones de congruencia lineal y ecuaciones diofánticas.

Simplificación de expresiones booleanas mediante álgebra de Boole

George Boole

George Boole, el patrón.

Cuando estamos repasando código podemos encontrarnos expresiones booleanas en condiciones lógicas más complejas de lo necesario que dificultan entender lo que el código hace. Esto puede ser debido a que la persona no tiene los conocimientos necesarios para hacerlo mejor (en esta profesión uno se encuentra hasta biólogos sin mayor interés en la profesión que el salario) o al constreñimiento que sufrió el programador cuando desarrolló el código que estamos repasando. Por lo tanto, puede ser nuestro propio código el que incurra en esos errores que ahora, con más calma, descubrimos y nos preguntamos en qué estaríamos pensando cuando hicimos esos whiles y esos ifs.

Para reducir el total de términos de una expresión booleana existen técnicas como los mapas de Karnaugh o el método de Quine-McCluskey pero en programación generalmente es suficiente con las leyes del álgebra de Boole. Teniendo en cuenta que la condición lógica AND es el producto booleano, OR es la suma y NOT el complemento (que se puede representar con los símbolos    o ¬), las leyes son: Sigue leyendo

Un algoritmo de generación de números pseudoaleatorios

números aleatorios

En la nota a pie del artículo La privacidad no existe citaba un artículo de la Wikipedia acerca de cómo generar números aleatorios mediante software. En esta entrada voy a poner un ejemplo sencillo desarrollado en Python. Obviamente Python tiene sus propias funciones para ello, pero nada nos impide desarrollar nuestra propia función.

En primer lugar veamos la teoría que hay detrás. Estos algoritmos lo que realmente generan son números pseudoaleatorios porque usan procedimientos sistemáticos. El método más empleado es el de congruencia lineal. Se eligen cuatro números enteros: un módulo m, un multiplicador a, el incremento c y la semilla x0 Estos números deben ser números enteros positivos mayores de 0. El artículo de la Wikipedia también impone que  a, c y x0 deben ser menores que m pero no encuentro ninguna razón para ello. De hecho, al final el artículo entra en contradicción al afirmar que el estándar POSIX para C define a = 1103515245 y m = 32768 Con lo que a > m

Una semilla válida la podemos obtener del reloj interno de la CPU: el RTC. En el ejemplo usaré el método time() de la librería time. Éste nos da el tiempo que ha pasado desde el tiempo Unix: el 1 de Enero de 1970. Como se puede ver, devuelve un número en coma flotante:

>>> time.time()
1487416084.126488

Como necesitamos un entero lo convertimos, el número de segundos será suficiente para el ejemplo. Sigue leyendo

La operación módulo en teoría de números y cómo la implementan diferentes lenguajes

En informática, la operación módulo, dados dos números, D el dividendo y d el divisor, nos da el resto r de la división entera de D por d, siendo d≠0 y |r|<|d|. Aunque no parece gran cosa, aun menos teniendo en cuenta que hace décadas que las CPUs incorporan coprocesadores matemáticos capaces de hacer a gran velocidad operaciones aparentemente más relevantes como raíces o logaritmos, la realidad es que es una operación fundamental en informática: Desde la generación de tablas de hash y números aleatorios hasta la criptografía, la generación de checksums y el álgebra computacional.

Veamos como Python resuelve esta operación. Abrimos la consola de Python y mediante el operador % le pedimos el módulo con D=27 y d=13:

>>> 27%13
 1

Así es: 27 = 13 * 2 + 1

Lo mismo nos devuelve PHP:

echo 27%13;
1

Nada nuevo bajo el sol, matemáticas básicas. Probemos esto ahora:

>>> -27%13
12

Así es: -27 = 13 * (-3) + 12

Cambiando los signos de D y d:

>>> 27%-13
-12

¡¿Comorll?!?

El resultado esperado es 1:

27 = -13 * (-2) + 1

¿Por qué da -12? ¿Y PHP qué opina de todo esto?

echo 27%-13;
1

Buen chico. ¿Será que un lenguaje tan usado y prestigioso como Python no sabe dividir? ¿Por qué entonces lo escogen en tantas facultades de todo el mundo para la enseñanza? La NASA también lo usa, ¿han escogido un lenguaje que no sabe ni calcular el resto de una división entera?

Siguiendo con PHP:

echo -27%13;
-1

Habíamos visto antes que el resultado era 12 ¿PHP anda también perdido? La definición de la Wikipedia tampoco ayuda a esclarecer que está pasando. En mi opinión, esa definición es inexacta pues dice que la operación módulo calcula el resto de la división euclídea y en la división euclídea se impone la condición r ≥ 0 Cuando no se impone esa condición, y no lo hice en la definición de la operación módulo que di al iniciar esta entrada, las divisiones tendrán dos cocientes y dos residuos posibles: por defecto y por exceso. Mientras que con r ≥ 0 sólo habrá un resultado posible. Si hay un cociente q habrá otro cociente q’ cuyo valor será q + 1 y cada cociente tendrá su respectivo residuo r y r’ Con D=-27 y d=13 tendremos q=-3 y r=12 por un lado y por otro q=-2 y r=-1

-27 = 13 * (-3) + 12

-27 = 13 * (-2) + (-1)

Por lo tanto ambos lenguajes están calculando correctamente. En teoría de números siempre se escoge r ≥ 0 y algunos lenguajes, como Pascal, Scheme o Stata así lo hacen también, para estos lenguajes sí aplica la definición de la Wikipedia, pero la mayoría y los más usados van a escoger la división por exceso o por defecto. Algunos, como PHP, C o C++ escogerán siempre para el resto el signo del dividendo D mientras que otros como Python o Mathematica el del divisor d. Aparentemente para conseguir esta igualdad de signo, usarán la división por exceso o por defecto según sea necesario, aunque lo que realmente hacen es una operación más rápida. Los que usan el mismo signo que el dividendo van realmente rápido, truncan la división:

r = D – d * truncar(D / d)

La función truncar simplemente extermina todos los decimales, o dicho más fino: se redondea al número entero inmediato más cercano a 0. Esto es lo que hace PHP o C:

-1 = -27 – [13*truncar(-27 / 13)] = -27 – [13*truncar(-2.07…)] = -27 – [13*(-2)] = -27 – (-26)

Por otro lado, los que usan el signo del divisor d, usan la función floor (suelo). Esta función proporciona el número entero más pequeño inmediato o igual al número que se le pasa. De manera que por ejemplo Π quedaría en 3 y -Π en -4. Está función se expresa mediante ⌊⌋: floor(x) = ⌊x⌋

r = D – d * ⌊D / d⌋

Esto es lo que Python hace:

12 = -27 – 13[⌊-27 / 13⌋] = -27 – 13[⌊-2.07…⌋] = -27 -13[-3] = -27 + 39

Con la primera extrañeza que nos encontramos, cuando pusimos un divisor negativo, ahora sabemos qué hace Python:

-12 = 27 – (-13)*[⌊27 / -13⌋] = 27 – (-13)⌊-2.07…⌋ = 27 – (-13)*(-3) = 27 – 39

Lenguajes respetuosos con la teoría, como Pascal, hacen esto:

Si d > 0 => q = ⌊D / d⌋

Si d < 0 => q = ⌈D / d⌉

Es decir, si el divisor es negativo usan la función ceiling (techo) Dicha función da el número entero más grande inmediato o igual al número que se le pasa. De manera que Π quedaría en 4 y -Π en -3. Está función se expresa mediante ⌈⌉: ceiling(x) =⌈x⌉

Hay que tener en cuenta que no es esto lo que hacen los lenguajes que cogen el signo del divisor, es decir, no es lo mismo ⌊x⌋ que truncar(), no en el caso de números negativos. Volviendo al ejemplo D = 27 y d = -13 que en Python dio r = -12 mientras que PHP r = 1:

-12 ≠ 27 – (-13*⌈27 / -13⌉) = 27 – (-13*(-2)) = 27 – 26 = 1

Es lo que nos devolvió PHP, porque para números negativos truncar(x) equivale a ⌈x⌉ Para los positivos, en cambio, truncar(x) equivale a ⌊x⌋, por eso para D y d positivos todo lenguaje da el mismo r. Las diferencias surgen con los negativos. Veamos el caso en que ambos son negativos, D = -27 y d = -13. PHP usa truncar(), que en este caso sí equivale a ⌈⌉ :

-27 – (-13* truncar(-27 / -13) = -27 – (-13*truncar(2.07)) = -27 – (-13*2) = -27 + 26 = -1

>>> -27%-13
-1
echo -27%-13;
-1

Mientras que ⌊⌋, que ningún lenguaje la usa, daría:

-27 – (-13*⌊-27 / -13⌋) = -27 – (-13*(⌊2.07⌋) = -27 – (-13)*3 = -27 + 39 = 12

Como q = (D – r) / d Acaba resultando una ecuación sin solución entera:

(-27 – 12) / -13

Por lo tanto r ≠ 12

En la definición de división entera que di al principio introduje una «trampa» para que los resultados de PHP, Python o C fueran correctos: añadí la condición |r|<|d| En teoría de números, la división entera se define sólo para el conjunto de los naturales y la única condición es d ≠ 0 En cambio, para el conjunto de los números enteros (que incluye los naturales más los negativos), se usa la división euclídea, la cual añade la condición r ≥ 0 con el objetivo de obtener unos q y r únicos.

En mi opinión, la razón principal por la que algunos lenguajes usan la división entera también para números negativos es obtener un conjunto de r más amplio. Sin la condición r ≥ 0 se obtienen también r negativos. Con ello, por ejemplo, se consigue reducir el número de colisiones en una tabla de hash. Secundariamente, porque es más rápido el cálculo: truncar es más rápido que redondear. Sacrifican purismo matemático por fines prácticos. Además, el programador puede seguir obteniendo el resto de la división euclídea por otras vías. Puede programar su propia función con ceil() y floor(), o con floor() y truncate() si viene de microsegundos. También es más que probable que el lenguaje tenga una función para obtener dicho resto. Por ejemplo C tiene la función fmod en la librería math.h