Máquinas de Turing en C

Una máquina de Turing es un dispositivo de cálculo lógico que utiliza un input en una o varias cintas que se van moviendo en función de la instrucción que tenga el estado de un autómata, para finalmente obtener un output mediante la reescritura de los datos en la misma cinta. La máquina tiene un cabezal de lectura y este lee el dato que se encuentra en la posición de la cinta. Por buscar una similitud lo mas parecido sería un casette de música, que se rebobina o se adelanta según la canción que nos interese, en la cual podremos borrar o grabar otra canción donde queramos. En el caso de una máquina de Turing de una cinta se define dicha máquina como una 7-tupla de la siguiente manera:

donde:

Símbolo Descripción
Conjunto finito de símbolos que pueden aparecer en la cinta.
Alfabeto de entrada a la máquina un conjunto finito de símbolos (menos el espacio en blanco).
Conjunto finito de estados.
Estado inicial en el que empieza la máquina.
Símbolo denominado blanco y es el único símbolo que se puede repetir un número infinito de veces. Es muy habitual que el símbolo para referirse al blanco sea '' '' o ''.
Conjunto de estados finales de aceptación.
Función de transición donde es un movimiento a la izquierda y es el movimiento a la derecha.

Cada función de transición se escribe así:

Y se interpreta de la siguiente manera:

  1. La máquina se encuentra en el estado
  2. Lee un
  3. Sobreescribe el con una
  4. El cabezal se mueve en la cinta una posición a la derecha (Right)
  5. Finalmente se cambia al estado

Teniendo claros los conceptos mínimos de los que se compone una máquina de Turing se puede diseñar una para saber si una cadena de entrada cumple una expresión lógico-matemática determinada. A continuación algunos ejemplos.

Cadenas con un numero de 'x' seguidas del mismo número de 'y'

La notación de este Lenguaje Independiente del Contexto (LIC) o expresión regular es:

Esta notación representa cualquier cadena compuesta de un número de caracteres seguido del mismo número de caracteres , y que el número de veces que aparecen estos caracteres sea mayor o igual que . Cualquiera de las siguientes cadenas de entrada sería válida:






Y las siguientes serían cadenas no válidas o no aceptadas:






Donde e pueden ser otros símbolos. Para este ejemplo se van a usar los símbolos pertenecientes al conjunto binario, o sea y , y valdrá . Así pues, estos son los datos que conforman la máquina de Turing:





Las funciones de transición serán las siguientes:












Esta lista de funciones de transición también se puede representar mediante la siguiente tabla:

Además de la lista de funciones o la tabla anterior, las funciones de transición también se pueden representar mediante el diagrama de un autómata. El ejemplo de máquina de Turing para el lenguaje se representaría con el siguiente autómata determinista.

A continuación un vídeo que he hecho -como he podido- en el que se puede apreciar el progreso de cada una de las funciones sobre la cinta.

Una vez explicado en detalle cómo funciona una máquina de Turing se puede extrapolar dicho funcionamiento a un programa, por ejemplo en lenguaje de programación C. Existen varias maneras de hacer un programa en C que sea capaz de evaluar si una cadena de entrada contiene un número de seguido del mismo número de , y desde luego hacerlo con una simulación de una máquina de Turing como la que se explica a continuación no es la mejor manera, pero sirve para comprender "la lógica" sobre su funcionamiento, que es realmente el objetivo de este artículo.

El programa tendrá al inicio una serie de variables y constantes. También la declaración de algunos procedimientos que se invocarán durante varias partes del programa, según el flujo de datos que vaya leyendo la máquina. Son los siguientes:

void estado_q0(char*, int);
void estado_q1(char*, int);
void estado_q2(char*, int);
void estado_q3(char*, int);
void stop(bool);
void addBlankL(char *w);
void addBlankR(char *w);
void salidaTrue();
void salidaFalse();
void printLog(char *w, char *ea, int i, char x, char m, char *e);
void stop(bool resultado);
int linea(int s);
void printTablaMT();

Se trata de procedimientos que contienen la lógica de cada uno de los 4 estados definidos en , y además hay un último procedimiento 'void stop(bool)' que servirá para detener la máquina cuando se llegue a un estado de aceptación (en nuestro caso ) o también cuando se llegue a un estado en el que no existen instrucciones para un símbolo leído, como por ejemplo o .

Procedimiento 'void estado_q0(char *w, int i)'

Representa la lógica del estado , el código es el siguiente:

void estado_q0(char *w, int i) {
	estadoActual = "q0";

	if (w[i] == '0') {
		printLog(w, estadoActual, i, 'X', 'R', "q1");
		w[i] = 'X';
		estado_q1(w, i += 1);
	} else if (w[i] == '1' || w[i] == 'X' || w[i] == '_') {
		stop(false);
	} else if (w[i] == 'Y') {
		printLog(w, estadoActual, i, 'Y', 'R', "q3");
		estado_q3(w, i += 1);
	}
}

Ejecuta un condicional que en función del símbolo leído por el cabezal solo continúa si dicho símbolo es o . Si lee cualquier otro símbolo, o sea , o un blanco , la máquina se detiene invocando el procedimiento 'stop(false);' e imprime por pantalla un mensaje indicando que la cadena de entrada no satisface (false) el lenguaje (LIC) definido para esta máquina de Turing. Tanto este procedimiento como los demás del tipo 'estado_qn' reciben dos argumentos, el primero es el puntero a la cadena o string de entrada y el segundo es la posición del cabezal de lectura. Esta posición del cabezal se incrementa si el cabezal se desplaza hacia la derecha (R) y se decrementa si el cabezal se desplaza a la izquierda (L).

Procedimiento 'void estado_q1(char *w, int i)'

Representa la lógica del estado , el código es el siguiente:

void estado_q1(char *w, int i) {
	estadoActual = "q1";

	if (w[i] == '0') {
		printLog(w, estadoActual, i, '0', 'R', "q1");
		estado_q1(w, i += 1);
	} else if (w[i] == '1') {
		printLog(w, estadoActual, i, 'Y', 'L', "q2");
		w[i] = 'Y';
		estado_q2(w, i -= 1);
	} else if (w[i] == 'X' || w[i] == '_') {
		stop(false);
	} else if (w[i] == 'Y') {
		printLog(w, estadoActual, i, 'Y', 'R', "q1");
		estado_q1(w, i += 1);
	}
}

Ejecuta un condicional que hace que el programa continúe solo si el símbolo leído es , o . Así pues, si el símbolo leído es o  la máquina se detendrá.

Procedimiento 'void estado_q2(char *w, int i)'

Representa la lógica del estado , el código es el siguiente:

void estado_q2(char *w, int i) {
	estadoActual = "q2";

	if (w[i] == '0') {
		printLog(w, estadoActual, i, '0', 'L', "q2");
		estado_q2(w, i -= 1);
	} else if (w[i] == '1' || w[i] == '_') {
		stop(false);
	} else if (w[i] == 'X') {
		printLog(w, estadoActual, i, 'X', 'R', "q0");
		estado_q0(w, i += 1);
	} else if (w[i] == 'Y') {
		printLog(w, estadoActual, i, 'Y', 'L', "q2");
		estado_q2(w, i -= 1);
	}
}

Ejecuta un condicional que hace que el programa continúe solo si el símbolo leído es , o . Así pues, si el símbolo leído es o  la máquina se detendrá.

Procedimiento 'void estado_q3(char *w, int i)'

Representa la lógica del estado , el código es el siguiente:

void estado_q3(char *w, int i) {
	estadoActual = "q3";

	if (w[i] == '0' || w[i] == '1' || w[i] == 'X') {
		stop(false);
	} else if (w[i] == 'Y') {
		printLog(w, estadoActual, i, 'Y', 'R', "q3");
		estado_q3(w, i += 1);
	} else if (w[i] == '_') {
		printLog(w, estadoActual, i, '_', 'R', "STOP");
		stop(true);
	}
}

Ejecuta un condicional que hace que el programa continúe solo si el símbolo leído es . Así pues, si el símbolo leído es , , o la máquina se detendrá indicando que no satisface el LIC de entrada, pero si el símbolo leído es un blanco  la máquina se detendrá indicando que la cadena de entrada si satisface el lenguaje mediante esta máquina de Turing.

Programa completo

A continuación el programa completo con comentarios en el código.

/*
 * turinMachine v0.05
 * Copyleft - 2017  Javier Dominguez Gomez
 * Written by Javier Dominguez Gomez <jdg@member.fsf.org>
 * GnuPG Key: D6648E2B
 * Madrid, Spain
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 * Compilation: gcc -std=gnu99 -O0 -g3 -Wall -c -o turing.o "turing.c"
 *              gcc -o turing turing.o
 *
 * Usage:       ./turing
 *
 * Case {x^{n}y^{n} | n>=0}
 * +--------------------------------------------------------------------+
 * |        0           1           X           Y           _           |
 * +--------------------------------------------------------------------+
 * | q0     (X,R,q1)    -           -           (Y,R,q3)    (_,STOP,q4) |
 * | q1     (0,R,q1)    (Y,L,q2)    -           (Y,R,q1)    -           |
 * | q2     (0,L,q2)    (1,L,q2)    (X,R,q0)    (Y,L,q2)    -           |
 * | q3     -           -           -           (Y,R,q3)    (_,STOP,q4) |
 * | q4     -           -           -           -           (STOP)      |
 * +--------------------------------------------------------------------+
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <stdbool.h>

/* Longitud maxima + 1. */
#define MAX_LONG 256

/* Declaracion de variables globales. */
char *estadoActual;
char *funcion = "0^{n}1^{n} | n >= 0";
int i, longitudW;

void estado_q0(char*, int);
void estado_q1(char*, int);
void estado_q2(char*, int);
void estado_q3(char*, int);
void stop(bool);
void addBlankL(char *w);
void addBlankR(char *w);
void salidaTrue();
void salidaFalse();
void printLog(char *w, char *ea, int i, char x, char m, char *e);
void stop(bool resultado);
int linea(int s);
void printTablaMT();

int main(int argc, char *argv[]) {

	/* Se asigna memoria para la cadena w. */
	char *w = malloc(MAX_LONG);

	printf("Introduce una cadena: ");

	/* Guarda la cadena w, con limite de longitud. */
	fgets(w, MAX_LONG, stdin);

	/* Se obtiene la logitud de la cadena de entrada,
	 * sin tener en cuanta el salto de linea. */
	longitudW = strlen(w) - 1;

	/* Elimina salto de linea del final, si existiera. */
	if ((longitudW + 1 > 0) && (w[longitudW] == '\n')) {
		w[longitudW] = '\0';
	}

	/* Primera condicion es que la cadena debe ser par */
	if (strlen(w) % 2 != 0) {
		salidaFalse();
	} else if(strlen(w) == 0){
		/* Acepta la cadena vacia. */
		salidaTrue();
	} else {
		printTablaMT();

		/* Se incorpora un blanco por la izquierda
		 * y otro por la derecha de la cadena w. */
		addBlankL(w);
		addBlankR(w);

		/* Le digo a la maquina de Turing que empiece
		 * por el caracter w[1] en el estado q0 */
		estado_q0(w, 1);
	}

	/* Se libera la memoria reservada para la daena w. */
	free(w);

	return 0;
}

void addBlankL(char *w) {
	/* Recorro el array de derecha a izquierda y desplazo
	 * los caracteres una posicion a la derecha */
	for (i = strlen(w); i > 0; i--) {
		w[i] = w[i - 1];
	}

	/* Concateno la cadena w con un blanco por la izquierda. */
	w[i] = '_';
}

void addBlankR(char *w) {
	/* Concateno la cadena w con un blanco por la derecha. */
	strcat(w, "_");
}

void salidaTrue() {
	printf("\n\nEnhorabuena! se cumple que %s\n\n", funcion);
}

void salidaFalse() {
	printf("\n\nNo se cumple que %s\n\n", funcion);
}

void printLog(char *w, char *ea, int i, char x, char m, char *e) {
	printf(
			"\nw es: \"%s\". Estoy en %s, leo \"%c\", lo cambio por \"%c\", me muevo a la %c y me voy al estado %s.",
			w, ea, w[i], x, m, e);
}

void estado_q0(char *w, int i) {
	estadoActual = "q0";

	// Si leo 0
	if (w[i] == '0') {

		printLog(w, estadoActual, i, 'X', 'R', "q1");

		// lo cambio por X
		w[i] = 'X';

		// me muevo a la derecha (R) y me voy al estado q1.
		estado_q1(w, i += 1);

		// Si leo 1, X o un blanco (_)
	} else if (w[i] == '1' || w[i] == 'X' || w[i] == '_') {

		// se para la maquina de Turing con un resultado false.
		stop(false);

		// Si leo Y
	} else if (w[i] == 'Y') {

		printLog(w, estadoActual, i, 'Y', 'R', "q3");

		// dejo Y sin cambiar, me muevo a la derecha (R) y me voy al estado q3.
		estado_q3(w, i += 1);
	}
}

void estado_q1(char *w, int i) {
	estadoActual = "q1";

	// Si leo 0
	if (w[i] == '0') {

		printLog(w, estadoActual, i, '0', 'R', "q1");

		// dejo 0 sin cambiar, me muevo a la derecha (R) y sigo en el estado q1.
		estado_q1(w, i += 1);

		// Si leo 1
	} else if (w[i] == '1') {

		printLog(w, estadoActual, i, 'Y', 'L', "q2");

		// lo cambio por Y
		w[i] = 'Y';

		// me muevo a la izquierda (L) y me voy al estado q2.
		estado_q2(w, i -= 1);

		// Si leo X o un blanco (_)
	} else if (w[i] == 'X' || w[i] == '_') {

		// se para la maquina de Turing con un resultado false.
		stop(false);

		// Si leo Y
	} else if (w[i] == 'Y') {

		printLog(w, estadoActual, i, 'Y', 'R', "q1");

		// dejo Y sin cambiar, me muevo a la derecha (R) y me voy al estado q1.
		estado_q1(w, i += 1);
	}
}

void estado_q2(char *w, int i) {
	estadoActual = "q2";

	// Si leo 0
	if (w[i] == '0') {

		printLog(w, estadoActual, i, '0', 'L', "q2");

		// dejo 0 sin cambiar, me muevo a la izquierda (L) y sigo en el estado q2.
		estado_q2(w, i -= 1);

		// Si leo 1 o un blanco (_)
	} else if (w[i] == '1' || w[i] == '_') {

		// se para la maquina de Turing con un resultado false.
		stop(false);

		// Si leo X
	} else if (w[i] == 'X') {

		printLog(w, estadoActual, i, 'X', 'R', "q0");

		// dejo X sin cambiar, me muevo a la derecha (R) y me voy al estado q0.
		estado_q0(w, i += 1);

		// Si leo Y
	} else if (w[i] == 'Y') {

		printLog(w, estadoActual, i, 'Y', 'L', "q2");

		// dejo Y sin cambiar, me muevo a la izquierda (L) y me voy al estado q2.
		estado_q2(w, i -= 1);
	}
}

void estado_q3(char *w, int i) {
	estadoActual = "q3";

	// Si leo 0, 1 o X
	if (w[i] == '0' || w[i] == '1' || w[i] == 'X') {

		// se para la maquina de Turing con un resultado false.
		stop(false);

		// Si leo Y
	} else if (w[i] == 'Y') {

		printLog(w, estadoActual, i, 'Y', 'R', "q3");

		// dejo Y sin cambiar, me muevo a la derecha (R) y sigo en el estado q3.
		estado_q3(w, i += 1);

		// Si leo un blanco (_)
	} else if (w[i] == '_') {

		printLog(w, estadoActual, i, '_', 'R', "STOP");

		// se para la maquina de Turing con un resultado true.
		stop(true);
	}
}

void stop(bool resultado) {
	if (resultado) {
		salidaTrue();
	} else {
		salidaFalse();
	}
}

int linea(int s) {
	fprintf(stdout,"+");
	for (int i = 0; i <= s; i++) {
		fprintf(stdout,"-");
	}
	fprintf(stdout,"+\n");
	return 0;
}

void printTablaMT(){
	printf("\nMaquina de Turing para %s\n\n", funcion);
	linea(62);
	printf("|%-5.5s      %-5.5s      %-5.5s      %-5.5s      %-5.5s   %-10.10s |\n", " ", "0", "1", "X", "Y", "    _");
	linea(62);
	printf("|%-5.5s   %-8.8s   %-8.8s   %-8.8s   %-8.8s   %-10.10s |\n", " q0", "(X,R,q1)", "   -", "   -", "(Y,R,q3)", "    _");
	printf("|%-5.5s   %-8.8s   %-8.8s   %-8.8s   %-8.8s   %-10.10s |\n", " q1", "(0,R,q1)", "(Y,L,q2)", "   -", "(Y,R,q1)", "    _");
	printf("|%-5.5s   %-8.8s   %-8.8s   %-8.8s   %-8.8s   %-10.10s |\n", " q2", "(0,L,q2)", "   -", "(X,R,q0)", "(Y,L,q2)", "    _");
	printf("|%-5.5s   %-8.8s   %-8.8s   %-8.8s   %-8.8s   %-10.10s |\n", " q3", "   -", "   -", "   -", "(Y,R,q3)", "(#,R,STOP)");
	linea(62);
}

Por último aclarar que las máquinas de Turing, que no dejan de ser un autómata, actualmente en Ciencias de la Computación se utilizan para crear compiladores de algunos lenguajes de programación. También para desarrollo de inteligencia artificial y redes neuronales. Así pues, no tiene mucho sentido hacer un programa como el de este artículo para otro fin que no sea el de explicar el funcionamiento forense de una máquina de Turing.

Espero que os haya gustado. Como siempre, cualquier pregunta, sugerencia o corrección será bien venida.

Operaciones con matrices en C

En este artículo voy a explicar cómo se realizan algunas de las operaciones mas básicas con matrices y su implementación en un programa C que te calculará automáticamente el determinante, la matriz traspuesta, la matriz adjunta y la matriz invertida (si existiera). Para el ejemplo voy a usar una matriz cuadrada de orden 3. Pero antes refresquemos la memoria recordando qué es cada cosa y poniendo un fragmento de código para ver representada la explicación en C.

Matriz inicial de orden 3

Tenemos una matriz de entrada, por ejemplo:

La primera parte del programa comienza con un bucle que nos pide introducir los datos de nuestra matriz cuadrada de orden 3:

Matriz matriz = malloc(sizeof(struct MatrizStruct));

for (i = 0; i < 3; i++) {
	for (j = 0; j < 3; j++) {
		printf("Numero %d,%d: ", i + 1, j + 1);
		scanf("%d", &matriz->matrix_a[i][j]);
	}
}

./Matrices
Numero 1,1: 3
Numero 1,2: 5
Numero 1,3: 1
Numero 2,1: 2
Numero 2,2: -1
Numero 2,3: 0
Numero 3,1: -1
Numero 3,2: 3
Numero 3,3: 1

Los datos se van a ir asignando a las direcciones de memoria correspondientes de cada coordenada en un array bidimensional que previamente he declarado en un struct, pensado para almacenar toda la información sobre la matriz.

struct MatrizStruct {
	int matrix_a[3][3];    //Este es el Array inicial
	int matrix_aT[3][3];
	int matrix_aAdj[3][3];
	int matrix_aInv[3][3];
};

Antes de todo esto, en un header de C he declarado Matriz como un tipo de dato (typedef). De este modo podré utilizar y modificar toda la estructura de la que se compone cada matriz.

struct MatrizStruct;
typedef struct MatrizStruct* Matriz;

Una vez que ya tenemos la matriz de entrada ya podemos hacer operaciones con ella.

Determinante

El determinante de una matriz (cuadrada) lo designaremos por  o bien por , aunque la primera forma es la mas utilizada. Se trata del número que se le asocia tras realizar el siguiente cálculo:

Cuyo desarrollo es idéntico al que nos proporciona la conocida Regla de Sarrus:

Para obtener el determinante de la matriz de entrada he utilizado la siguiente función a la que es necesario pasarle por argumento la estructura Matriz. Esta calcula el determinante y devuelve un dato de tipo int.

int determinante(Matriz m) {
   return (m->matrix_a[0][0] * m->matrix_a[1][1] * m->matrix_a[2][2])
   + (m->matrix_a[2][0] * m->matrix_a[0][1] * m->matrix_a[1][2])
   + (m->matrix_a[0][2] * m->matrix_a[1][0] * m->matrix_a[2][1])
   - (m->matrix_a[2][0] * m->matrix_a[1][1] * m->matrix_a[0][2])
   - (m->matrix_a[0][0] * m->matrix_a[1][2] * m->matrix_a[2][1])
   - (m->matrix_a[0][1] * m->matrix_a[1][0] * m->matrix_a[2][2]);
}

En este ejemplo, el determinante de la matriz de entrada es el siguiente:

Matriz traspuesta

La matriz traspuesta de se designa y es una matriz en la que los elementos que formaban una fila en pasan a ser los elementos de una columna en . Véase en la siguiente notación:

Por ejemplo:

Para obtener la matriz traspuesta de la matriz de entrada he utilizado la siguiente función a la que también es necesario pasarle por argumento la estructura Matriz. Esta función no solo modifica el array del struct si no que también devuelve un dato de tipo Matriz. Esto nos puede venir muy bien para otro tipo de operaciones que no se tratan en este programa.

Matriz traspuesta(Matriz m) {
	m->matrix_aT[0][0] = m->matrix_a[0][0];
	m->matrix_aT[0][1] = m->matrix_a[1][0];
	m->matrix_aT[0][2] = m->matrix_a[2][0];
	m->matrix_aT[1][0] = m->matrix_a[0][1];
	m->matrix_aT[1][1] = m->matrix_a[1][1];
	m->matrix_aT[1][2] = m->matrix_a[2][1];
	m->matrix_aT[2][0] = m->matrix_a[0][2];
	m->matrix_aT[2][1] = m->matrix_a[1][2];
	m->matrix_aT[2][2] = m->matrix_a[2][2];

	return m;
}

Matriz adjunta

Se denomina matriz adjunta de o a la matriz obtenida al sustituir los elementos por sus adjuntos , es decir a la matriz formada por los adjuntos de los elementos. Por ejemplo:

El cálculo sería el siguiente:

Así pues, la matriz adjunta de es:

Al igual que con la función para obtener la matriz traspuesta, la matriz adjunta también se obtiene con una función a a que pasaremos como argumento la estructura Matriz, a continuación se calculará cada elemento del array bidimensional perteneciente al struct y se asignarán los valores para esta nueva matriz. Esta función también tiene un return de tipo Matriz.

Matriz adjunta(Matriz m) {
	m->matrix_aAdj[0][0] = m->matrix_a[1][1] * m->matrix_a[2][2] - m->matrix_a[2][1] * m->matrix_a[1][2];
	m->matrix_aAdj[1][0] = -(m->matrix_a[1][0] * m->matrix_a[2][2] - m->matrix_a[2][0] * m->matrix_a[1][2]);
	m->matrix_aAdj[2][0] = m->matrix_a[1][0] * m->matrix_a[2][1] - m->matrix_a[2][0] * m->matrix_a[1][1];
	m->matrix_aAdj[0][1] = -(m->matrix_a[0][1] * m->matrix_a[2][2] - m->matrix_a[2][1] * m->matrix_a[0][2]);
	m->matrix_aAdj[1][1] = m->matrix_a[0][0] * m->matrix_a[2][2] - m->matrix_a[2][0] * m->matrix_a[0][2];
	m->matrix_aAdj[2][1] = -(m->matrix_a[0][0] * m->matrix_a[2][1] - m->matrix_a[2][0] * m->matrix_a[0][1]);
	m->matrix_aAdj[0][2] = m->matrix_a[0][1] * m->matrix_a[1][2] - m->matrix_a[1][1] * m->matrix_a[0][2];
	m->matrix_aAdj[1][2] = -(m->matrix_a[0][0] * m->matrix_a[1][2] - m->matrix_a[1][0] * m->matrix_a[0][2]);
	m->matrix_aAdj[2][2] = m->matrix_a[0][0] * m->matrix_a[1][1] - m->matrix_a[1][0] * m->matrix_a[0][1];

	return m;
}

Matriz inversa

Solo tienen matriz inversa aquellas matrices que su determinante es distinto de cero, en cuyo caso su matriz inversa es igual a la traspuesta de la matriz adjunta dividida por el determinante de . Esta sería la notación:

Una notación mas reducida seria:

Para calcular la matriz inversa también hay una función que recibe la estructura Matriz como argumento. Se asignan los valores correspondientes a la matriz traspuesta de la adjunta y finalmente tiene un return de la estructura de tipo Matriz ya modificada. Solo se invocará a esta función en el caso de que el determinante sea distinto de .

Matriz inversa(Matriz m) {
	m->matrix_aInv[0][0] = m->matrix_aAdj[0][0];
	m->matrix_aInv[1][0] = m->matrix_aAdj[1][0];
	m->matrix_aInv[2][0] = m->matrix_aAdj[2][0];
	m->matrix_aInv[0][1] = m->matrix_aAdj[0][1];
	m->matrix_aInv[1][1] = m->matrix_aAdj[1][1];
	m->matrix_aInv[2][1] = m->matrix_aAdj[2][1];
	m->matrix_aInv[0][2] = m->matrix_aAdj[0][2];
	m->matrix_aInv[1][2] = m->matrix_aAdj[1][2];
	m->matrix_aInv[2][2] = m->matrix_aAdj[2][2];

	return m;
}

Pero ojo, Esta función solo asigna los valores de la matriz traspuesta de la adjunta. No se generan fracciones por que el programa las ejecutaría y en ocasiones daría como resultado números decimales muy largos y no quedaría presentable, así pues, he decidido que se almacenen como numerador de la fracción que mas adelante se imprimirá por pantalla añadiéndole su denominador, que será el determinante de la matriz inicial .

printf("\nMatriz inversa:\n");
	/* Solo se calculara la inversa si el determinante es != 0 */
	if(determinante(m) != 0){
		inversa(m);

		/* Imprime la matriz inversa. */
		for (i = 0; i < 3; i++) {
			for (j = 0; j < 3; j++) {
				printf("%3d/%d", m->matrix_aInv[i][j], determinante(m));
			}
			printf("\n");
		}
	} else {
		printf("No tiene invertida por que su determinante es 0.\n");
	}

Programa completo

Si queréis probar este modesto programa en vuestra máquina aquí os dejo el contenido de los dos archivos que necesitáis para compilarlo. Se admiten correcciones y mejoras que seguro que tiene un montón.

Archivo matriz.h
/*
 * matriz.h v0.01
 * Copyleft - 2017  Javier Dominguez Gomez
 * Written by Javier Dominguez Gomez <jdg@member.fsf.org>
 * GnuPG Key: 6ECD1616
 * Madrid, Spain
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#ifndef MATRIZ_H_
#define MATRIZ_H_

struct MatrizStruct;
typedef struct MatrizStruct* Matriz;

void printDatos(Matriz m);
int determinante(Matriz m);
Matriz traspuesta(Matriz m);
Matriz adjunta(Matriz m);
Matriz inversa(Matriz m);

#endif /* MATRIZ_H_ */

 

Archivo matriz.c
/*
 * matriz.c v0.01
 * Copyleft - 2017  Javier Dominguez Gomez
 * Written by Javier Dominguez Gomez <jdg@member.fsf.org>
 * GnuPG Key: 6ECD1616
 * Madrid, Spain
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 * Compilation: gcc -O0 -g3 -Wall -c -MF"matriz.d" -MT"matriz.o" -o "matriz.o" "matriz.c"
 *              gcc  -o "matriz" ./matriz.o
 *
 * Usage:       ./matriz
 */

#include "matriz.h"

#include <stdio.h>
#include <stdlib.h>

int i, j;

struct MatrizStruct {
	int matrix_a[3][3];
	int matrix_aT[3][3];
	int matrix_aAdj[3][3];
	int matrix_aInv[3][3];
};

int determinante(Matriz m) {
	return (m->matrix_a[0][0] * m->matrix_a[1][1] * m->matrix_a[2][2])
			+ (m->matrix_a[2][0] * m->matrix_a[0][1] * m->matrix_a[1][2])
			+ (m->matrix_a[0][2] * m->matrix_a[1][0] * m->matrix_a[2][1])
			- (m->matrix_a[2][0] * m->matrix_a[1][1] * m->matrix_a[0][2])
			- (m->matrix_a[0][0] * m->matrix_a[1][2] * m->matrix_a[2][1])
			- (m->matrix_a[0][1] * m->matrix_a[1][0] * m->matrix_a[2][2]);
}

Matriz traspuesta(Matriz m) {
	m->matrix_aT[0][0] = m->matrix_a[0][0];
	m->matrix_aT[0][1] = m->matrix_a[1][0];
	m->matrix_aT[0][2] = m->matrix_a[2][0];
	m->matrix_aT[1][0] = m->matrix_a[0][1];
	m->matrix_aT[1][1] = m->matrix_a[1][1];
	m->matrix_aT[1][2] = m->matrix_a[2][1];
	m->matrix_aT[2][0] = m->matrix_a[0][2];
	m->matrix_aT[2][1] = m->matrix_a[1][2];
	m->matrix_aT[2][2] = m->matrix_a[2][2];

	return m;
}

Matriz adjunta(Matriz m) {
	m->matrix_aAdj[0][0] = m->matrix_a[1][1] * m->matrix_a[2][2] - m->matrix_a[2][1] * m->matrix_a[1][2];
	m->matrix_aAdj[1][0] = -(m->matrix_a[1][0] * m->matrix_a[2][2] - m->matrix_a[2][0] * m->matrix_a[1][2]);
	m->matrix_aAdj[2][0] = m->matrix_a[1][0] * m->matrix_a[2][1] - m->matrix_a[2][0] * m->matrix_a[1][1];
	m->matrix_aAdj[0][1] = -(m->matrix_a[0][1] * m->matrix_a[2][2] - m->matrix_a[2][1] * m->matrix_a[0][2]);
	m->matrix_aAdj[1][1] = m->matrix_a[0][0] * m->matrix_a[2][2] - m->matrix_a[2][0] * m->matrix_a[0][2];
	m->matrix_aAdj[2][1] = -(m->matrix_a[0][0] * m->matrix_a[2][1] - m->matrix_a[2][0] * m->matrix_a[0][1]);
	m->matrix_aAdj[0][2] = m->matrix_a[0][1] * m->matrix_a[1][2] - m->matrix_a[1][1] * m->matrix_a[0][2];
	m->matrix_aAdj[1][2] = -(m->matrix_a[0][0] * m->matrix_a[1][2] - m->matrix_a[1][0] * m->matrix_a[0][2]);
	m->matrix_aAdj[2][2] = m->matrix_a[0][0] * m->matrix_a[1][1] - m->matrix_a[1][0] * m->matrix_a[0][1];

	return m;
}

Matriz inversa(Matriz m) {
	/* Una matriz solo es invertible si su determinante es distinto de 0. */

	/* La matriz inversa es igual a (1/|A|)*(Adj(A))^t */

	/* Construyo solo la matriz traspuesta de la adjunta. Cada
	 * elemento de la matriz sera el numerador de una fraccion que
	 * tendra por denominador el determinite de la matriz de entrada. */
	m->matrix_aInv[0][0] = m->matrix_aAdj[0][0];
	m->matrix_aInv[1][0] = m->matrix_aAdj[1][0];
	m->matrix_aInv[2][0] = m->matrix_aAdj[2][0];
	m->matrix_aInv[0][1] = m->matrix_aAdj[0][1];
	m->matrix_aInv[1][1] = m->matrix_aAdj[1][1];
	m->matrix_aInv[2][1] = m->matrix_aAdj[2][1];
	m->matrix_aInv[0][2] = m->matrix_aAdj[0][2];
	m->matrix_aInv[1][2] = m->matrix_aAdj[1][2];
	m->matrix_aInv[2][2] = m->matrix_aAdj[2][2];

	return m;
}

void printDatos(Matriz m) {

	/* Imprime la matriz de entrada. */
	printf("\nMatriz:\n");
	for (i = 0; i < 3; i++) {
		for (j = 0; j < 3; j++) {
			printf("%3d", m->matrix_a[i][j]);
		}
		printf("\n");
	}

	/* Imprime su determinante. */
	printf("\nEl determinante es: %d\n", determinante(m));

	/* Imprime la matriz traspuesta. */
	printf("\nMatriz traspuesta:\n");
	for (i = 0; i < 3; i++) {
		for (j = 0; j < 3; j++) {
			printf("%3d", m->matrix_aT[i][j]);
		}
		printf("\n");
	}

	/* Imprime la matriz adjunta. */
	printf("\nMatriz adjunta:\n");
	for (i = 0; i < 3; i++) {
		for (j = 0; j < 3; j++) {
			printf("%3d", m->matrix_aAdj[i][j]);
		}
		printf("\n");
	}

	printf("\nMatriz inversa:\n");
	/* Solo se calculara la inversa si el determinante es distinto de cero. */
	if(determinante(m) != 0){
		inversa(m);

		/* Imprime la matriz inversa. */
		for (i = 0; i < 3; i++) {
			for (j = 0; j < 3; j++) {
				printf("%3d/%d", m->matrix_aInv[i][j], determinante(m));
			}
			printf("\n");
		}
	} else {
		printf("No tiene invertida por que su determinante es 0.\n");
	}
}

int main(void) {
	Matriz matriz = malloc(sizeof(struct MatrizStruct));

	/* Pide datos de la matriz 3X3 */
	for (i = 0; i < 3; i++) {
		for (j = 0; j < 3; j++) {
			printf("Numero %d,%d: ", i + 1, j + 1);
			scanf("%d", &matriz->matrix_a[i][j]);
		}
	}

	/* Se calcula todo */
	determinante(matriz);
	traspuesta(matriz);
	adjunta(matriz);

	/* Se muestran todos los datos por pantalla. */
	printDatos(matriz);

	return 0;
}

Para compilarlo solo tenéis que guardar estos dos archivos en una carpeta y ejecutar lo siguiente:

gcc -O0 -g3 -Wall -c -MF"matriz.d" -MT"matriz.o" -o "matriz.o" "matriz.c"
gcc  -o "matriz" ./matriz.o

Una simulación de su ejecución sería la siguiente:

 

Un millón de dígitos de Pi

En la siguiente imagen que he generado se pueden apreciar el primer millón de dígitos de Pi. Se trata de una foto MUY grande (22.8MB) a la que podéis acceder y descargar:

http://ingenieriainversa.org/1millon_de_digitos_de_Pi.png

Para obtener este primer millón de dígitos de Pi he utilizado el programa C que emplea el Algoritmo de Chudnovsky del que ya hablé en este otro post.

Constante de Kaprekar

El número 6174 es un número realmente misterioso. A simple vista, puede no parecer tan obvio, pero estamos a punto de ver que cualquier persona que sepa restar también puede descubrir el misterio que hace al número 6174 tan especial. En en año 1949 el matemático Dattatreya Ramachandra Kaprekar (Devlali, India) ideó un proceso que hoy en día conocemos como la Operación de Kaprekar. Esta operación está sujeta a las siguientes reglas:

Esto es:

  • El número inicial ha de ser una cadena de dígitos que formen un numero natural.
  • La longitud de ha de ser de 4 dígitos.
  • Cada uno de los dígitos ha de ser un número natural comprendido entre el y el , ambos inclusive.
  •  No puede haber tres dígitos seguidos iguales.

Por ejemplo, aquí algunos ejemplos de números aceptados y no aceptados conforme a las reglas de arriba:

Aceptados No aceptados
1234 1118
7391 4333
8033 5555
2764 0000

La operación es muy sencilla, consiste en elegir un número que cumpla la reglas ya descritas anteriormente. Por ejemplo el número . A continuación, reorganizar los dígitos para obtener el número más grande y el más pequeño que estos dígitos pueden hacer.

Número mas grande posible
Número mas pequeño posible

Después, restar al número más grande el número más pequeño para obtener un nuevo número.

Si se repite una y otra vez la operación con cada nuevo número, se termina obteniendo la constante de Kaprekar, o sea el número .




El motivo por el que no se puede escoger un número con tres dígitos consecutivos iguales es por que o bien en la primera resta (cuatro dígitos seguidos iguales) o en la segunda resta (tres dígitos seguidos iguales) termina dando resto . Por ejemplo:


Al final, da igual el número que se escoja (siempre que cumpla las reglas anteriores), al final siempre devuelve el número .

A continuación dejo el código fuente de un programa en C que hace la Operación de Keprekar sobre un número introducido.

/*
 * constanteKaprekar v0.01
 * Copyleft - 2016  Javier Dominguez Gomez
 * Written by Javier Dominguez Gomez <jdg@member.fsf.org>
 * GnuPG Key: 6ECD1616
 * Madrid, Spain
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 * Compilation: gcc -std=gnu99 -O0 -g3 -Wall -c -o constanteKaprekar.o "constanteKaprekar.c"
 *              gcc -o constanteKaprekar constanteKaprekar.o
 *
 * Usage:       ./constanteKaprekar
 */

#include <stdio.h>
#include <stdbool.h>

/* Procedimiento para ordenar dos numeros. */
void swap(int *x, int *y) {
	if (*x < *y) {
		int aux;
		aux = *x;
		*x = *y;
		*y = aux;
	}
}

/* Funcion que devuelve la longitud de un int. */
int LongitudInt(int n) {
	int i = 1;
	while (n > 9) {
		i++;
		n /= 10;
	}
	return i;
}

/* Funcion que devuelve un int con sus digitos en orden descendente. */
int OrdenacionDescendenteInt(int n) {
	int i, n_ordenado = 0;
	for (i = 9; i >= 0; i--) {
		int aux = n;
		while (aux > 0) {
			int digit = aux % 10;
			if (digit == i) {
				n_ordenado *= 10;
				n_ordenado += digit;
			}
			aux /= 10;
		}
	}
	return n_ordenado;
}

/* Funcion que devuelve un int con sus digitos en orden ascendente. */
int OrdenacionAscendenteInt(int n) {
	int i, n_ordenado = 0;
	for (i = 0; i <= 9; i++) {
		int aux = n;
		while (aux > 0) {
			int digit = aux % 10;
			if (digit == i) {
				n_ordenado *= 10;
				n_ordenado += digit;
			}
			aux /= 10;
		}
	}
	return n_ordenado;
}

/* Funcion que comprueba si el numero introducido tiene 4 digitos. */
static bool CompruebaDigitos(int n) {
	if (LongitudInt(n) < 4) {
		printf("Longitud del numero %d invalida. Tiene que tener 4 digitos.\n", n);
		return false;
	} else {
		return true;
	}
}

int main() {
	int n, x, y, kaprekar;

	do {
		printf("Numero de 4 digitos: ");
		scanf("%d", &n);
	} while (CompruebaDigitos(n) == false);

	x = OrdenacionDescendenteInt(n);
	y = OrdenacionAscendenteInt(n);
	kaprekar = x - y;

	printf("%d - %d = %d\n", x, y, kaprekar);

	do {
		if (kaprekar == 0) {
			printf("Numero %d invalido por resto 0.\n", n);
			break;
		}
		x = OrdenacionDescendenteInt(kaprekar);
		y = OrdenacionAscendenteInt(kaprekar);
		kaprekar = x - y;
		printf("%d - %d = %d\n", x, y, kaprekar);
	} while (kaprekar != 6174);

	return 0;
}

Una vez compilado este es su funcionamiento:

./constanteKaprekar
Numero de 4 digitos: 9581
9851 - 1589 = 8262
8622 - 2268 = 6354
6543 - 3456 = 3087
8730 - 378 = 8352
8532 - 2358 = 6174

Cálculo de números primos en C

Un número primo es un número natural que solo es divisible (con resto cero) por si mismo y por . Por ejemplo, el número no es primo por que se puede dividir entre , y , sin embargo el número si es primo, ya que solo se puede dividir entre y si mismo. Veamos el detalle:

4 % 1 = 0 (true)
4 % 2 = 0 (true)
4 % 3 = 1 (false)
4 % 4 = 0 (true)

Como se puede ver, el número no solo es divisible con resto cero entre y él mismo, también lo es entre , por lo tanto no es un número primo. Veamos ahora el detalle con el número :

11 % 1  = 0 (true)
11 % 2  = 1 (false)
11 % 3  = 2 (false)
11 % 4  = 3 (false)
11 % 5  = 1 (false)
11 % 6  = 5 (false)
11 % 7  = 4 (false)
11 % 8  = 3 (false)
11 % 9  = 2 (false)
11 % 10 = 1 (false)
11 % 11 = 0 (true)

En el caso del número si podemos afirmar que es un número primo ya que solo es divisible por y por él mismo.

Ahora que ya sabemos cómo podemos identificar un número primo solo queda hacer la lógica en un programa C en el que indicaremos como valor para la variable limite la cantidad de números a analizar, para este ejemplo he puesto .

#include <stdio.h>

int main(int argc, char **argv) {
	int i, limite = 1000, primo, j;

	for (i = 2; i <= limite; i++) {
		for (j = 2; j < i; j++) {
			primo = 1;
			if (i % j == 0) {
				primo = 0;
				break;
			}
		}
		if (primo != 0) printf("%d ", i);
	}
	return 0;
}

Y tras compilarlo, al ejecutarlo nos devuelve lo siguiente:

./numerosPrimos
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659 661 673 677 683 691 701 709 719 727 733 739 743 751 757 761 769 773 787 797 809 811 821 823 827 829 839 853 857 859 863 877 881 883 887 907 911 919 929 937 941 947 953 967 971 977 983 991 997

 

Para obtener más cantidad de números primos solo hay que modificar el valor de la variable limite.

Cálculo aproximado de Pi

A lo largo de la historia, varias personas han intentado hallar una fórmula que sirva para obtener una aproximación lo mas exacta posible del número Pi. Algunas fórmulas son muy precisas, pero todas tienen un margen de error cada cierto número elevado de decimales. Por suerte hoy disponemos de computadoras que permiten realizar cálculos a enormes velocidades.

Pi

A finales del año 2010 los matemáticos Alexander Yee y Shigeru Kondo utilizaron el Algoritmo de Chudnovsky como método para calcular los decimales de Pi y una computadora casera con 96 GB de memoria RAM para conseguir superar el record que ya poseían, y que actualmente es de 10 Billones de decimales. Para ello tardaron 371 días con la computadora trabajando a su máximo rendimiento. Toda una proeza, sin duda.

La expresión matemática del Algoritmo de Chudnovsky es la siguiente:

Aquí he hecho un Hack con el algoritmo en C que devuelve como resultado un valor escandalosamente aproximado de Pi con 16 dígitos:

#include <stdio.h>
#include <math.h>

//Funcion que devuelve el factorial de un numero
double fact(double n) {
 if(n == 0) return 1;
 else return n * fact(n-1);
}

int main(void) {

 double k, pi;

 //Algoritmo de Chudnovsky
 for(k=0.0; k<=1.0; k++) {
  pi += fact((6.0 * k)) * (13591409.0 + (545140134.0 * k)) / ((fact(3.0 * k) * pow(fact(k),3.0)) * pow(-640320.0,(3.0 * k)));
 }

 pi = 426880 * sqrt(10005)/pi;

 printf("%.16f\n", pi);

 return 0;
}

El resultado de su ejecución sería el siguiente:

./chudnovsky
3.1415926535897931

Pero esto no es suficiente. Si lo que queremos es obtener, por ejemplo, un millón de decimales o más tendremos que usar un método un poco diferente, y para ello utilizaremos la librería GMP de la Free Software Foundation. Esta librería nos proporcionará una serie de funciones y procedimientos aritméticos para nuestros cálculos, así que para poder probarlo es necesario instalar la librería en nuestro sistema antes de nada. Este es un ejemplo* de código C que tras ser compilado y ejecutado nos devolverá tantos dígitos de Pi como le hayamos indicado como argumento:

/*
* Compute pi to a certain number of decimal digits, and print it.
*
*   gcc -O2 -Wall -o chudnovsky chudnovsky.c -lgmp
*
* WARNING: This is a demonstration program only, is not optimized for
* speed, and should not be used for serious work!
*
* The Chudnovsky Algorithm:
*                               _____
*                     426880 * /10005
*  pi = ---------------------------------------------
*         _inf_
*         \     (6*k)! * (13591409 + 545140134 * k)
*          \    -----------------------------------
*          /     (3*k)! * (k!)^3 * (-640320)^(3*k)
*         /____
*          k=0
*
* http://en.wikipedia.org/wiki/Pi#Rapidly_convergent_series
*
* First million digits: http://www.piday.org/million.php
*
* Copyright (c) 2012 Brian "Beej Jorgensen" Hall <beej@beej.us>
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject to
* the following conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
* CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
* TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
* SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gmp.h>

// how many to display if the user doesn't specify:
#define DEFAULT_DIGITS 60

// how many decimal digits the algorithm generates per iteration:
#define DIGITS_PER_ITERATION 14.1816474627254776555

/**
 * Compute pi to the specified number of decimal digits using the
 * Chudnovsky Algorithm.
 *
 * http://en.wikipedia.org/wiki/Pi#Rapidly_convergent_series
 *
 * NOTE: this function returns a malloc()'d string!
 *
 * @param digits number of decimal digits to compute
 *
 * @return a malloc'd string result (with no decimal marker)
 */
char *chudnovsky(unsigned long digits)
{
	mpf_t result, con, A, B, F, sum;
	mpz_t a, b, c, d, e;
	char *output;
	mp_exp_t exp;
	double bits_per_digit;

	unsigned long int k, threek;
	unsigned long iterations = (digits/DIGITS_PER_ITERATION)+1;
	unsigned long precision_bits;

	// roughly compute how many bits of precision we need for
	// this many digit:
	bits_per_digit = 3.32192809488736234789; // log2(10)
	precision_bits = (digits * bits_per_digit) + 1;

	mpf_set_default_prec(precision_bits);

	// allocate GMP variables
	mpf_inits(result, con, A, B, F, sum, NULL);
	mpz_inits(a, b, c, d, e, NULL);

	mpf_set_ui(sum, 0); // sum already zero at this point, so just FYI

	// first the constant sqrt part
	mpf_sqrt_ui(con, 10005);
	mpf_mul_ui(con, con, 426880);

	// now the fun bit
	for (k = 0; k < iterations; k++) {
		threek = 3*k;

		mpz_fac_ui(a, 6*k);  // (6k)!

		mpz_set_ui(b, 545140134); // 13591409 + 545140134k
		mpz_mul_ui(b, b, k);
		mpz_add_ui(b, b, 13591409);

		mpz_fac_ui(c, threek);  // (3k)!

		mpz_fac_ui(d, k);  // (k!)^3
		mpz_pow_ui(d, d, 3);

		mpz_ui_pow_ui(e, 640320, threek); // -640320^(3k)
		if ((threek&1) == 1) { mpz_neg(e, e); }

		// numerator (in A)
		mpz_mul(a, a, b);
		mpf_set_z(A, a);

		// denominator (in B)
		mpz_mul(c, c, d);
		mpz_mul(c, c, e);
		mpf_set_z(B, c);

		// result
		mpf_div(F, A, B);

		// add on to sum
		mpf_add(sum, sum, F);
	}

	// final calculations (solve for pi)
	mpf_ui_div(sum, 1, sum); // invert result
	mpf_mul(sum, sum, con); // multiply by constant sqrt part

	// get result base-10 in a string:
	output = mpf_get_str(NULL, &exp, 10, digits, sum); // calls malloc()

	// free GMP variables
	mpf_clears(result, con, A, B, F, sum, NULL);
	mpz_clears(a, b, c, d, e, NULL);

	return output;
}

/**
 * Print a usage message and exit
 */
void usage_exit(void)
{
	fprintf(stderr, "usage: chudnovsky [digits]\n");
	exit(1);
}

/**
 * MAIN
 *
 * See usage_exit() for usage.
 */
int main(int argc, char **argv)
{
	char *pi, *endptr;
	long digits;

	switch (argc) {
		case 1:
			digits = DEFAULT_DIGITS;
			break;

		case 2:
			digits = strtol(argv[1], &endptr, 10);
			if (*endptr != '\0') { usage_exit(); }
			break;

		default:
			usage_exit();
	}

	if (digits < 1) { usage_exit(); }

	pi = chudnovsky(digits);

	// since there's no decimal point in the result, we'll print the
	// first digit, then the rest of it, with the expectation that the
	// decimal will appear after "3", as per usual:
	printf("%.1s.%s\n", pi, pi+1);

	// chudnovsky() malloc()s the result string, so let's be proper:
	free(pi);

	return 0;
}

Por ejemplo, si queremos saber los 1000 primeros dígitos de Pi tendríamos que ejecutarlo de la siguiente manera:

./chudnovsky 1000

 Y nos devuelve Pi con los siguientes 1000 dígitos:

3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642019

Como curiosidad, si obtenemos decimales de Pi podremos observar que en la posición aparece la cadena... . Alucinante, ¿verdad? Aquí se puede ver un fragmento de los últimos dígitos:

Cadena 44899 en la posición 44899 de los decimales de Pi.
Cadena 44899 en la posición 44899 de los decimales de Pi.

IMPORTANTE: No te recomiendo que pruebes a calcular un número de decimales de Pi muy elevado (por ejemplo un millón), ya que es posible que tu CPU comience a calentarse considerablemente. Este era un truco muy viejo que se usaba para "romper" el ordenador de alguien hace muchos años, y consistía en meter en el arranque un algoritmo que calculase decimales de Pi como si no hubiera mañana en un sumatorio de infinitas iteraciones.

*Agradecimientos a Brian "Beej Jorgensen" Hall.

Curiosidades de la simetría numérica

¿Quien dijo que las matemáticas eran aburridas? Son lo suficientemente complejas como para estar en "casi todos" los rincones del Universo, pero también tienen algunas fórmulas sencillas que dan como resultado unas series numéricas que forman patrones simétricos, capicúas, progresiones geométricas y un montón de cosas más. A continuación dejo unos ejemplos matemáticos y su desarrollo en C para que lo pruebe el que quiera:

Ejemplo 1
#include <stdio.h>

int concat();

int main(void) {
	int num = 1;

	for (int i = 1; i <= 9; i++) {
		printf("%9d * 8 + %d = %d\n", num, i, num*8+i);
		num = concat(num, i+1);
	}
}

int concat(int x, int y) {
	int temp = y;
	while (y != 0) {
		x *= 10;
		y /= 10;
	}
	return x + temp;
}

El resultado de la ejecución de este primer ejemplo será la siguiente:

        1 x 8 + 1 = 9
       12 x 8 + 2 = 98
      123 x 8 + 3 = 987
     1234 x 8 + 4 = 9876
    12345 x 8 + 5 = 98765
   123456 x 8 + 6 = 987654
  1234567 x 8 + 7 = 9876543
 12345678 x 8 + 8 = 98765432
123456789 x 8 + 9 = 987654321
Ejemplo 2
#include <stdio.h>

int concat();

int main(void) {
	int num = 1;

	for (int i = 1; i <= 9; i++) {
		printf("%9d * 9 + %2d = %d\n", num, i+1, num*9+i+1);
		num = concat(num, i + 1);
	}
}

int concat(int x, int y) {
	int temp = y;
	while (y != 0) {
		x *= 10;
		y /= 10;
	}
	return x + temp;
}

El resultado de la ejecución de este segundo ejemplo será la siguiente:

        1 x 9 +  2 = 11
       12 x 9 +  3 = 111
      123 x 9 +  4 = 1111
     1234 x 9 +  5 = 11111
    12345 x 9 +  6 = 111111
   123456 x 9 +  7 = 1111111
  1234567 x 9 +  8 = 11111111
 12345678 x 9 +  9 = 111111111
123456789 x 9 + 10 = 1111111111
Ejemplo 3
#include <stdio.h>
#include <math.h>

void f(int x) {
	printf("f(%2d) = ((10^%2d)-1)/9 = %.0f\n", x, x, (pow(10, x) - 1) / 9);
}

int main(int argc, char **argv) {
	int x;
	printf("f(x) = ((10^x)-1)/9\n\n");
	for (x = 1; x <= 20; x++) f(x);
	return 0;
}

El resultado de la ejecución de este tercer ejemplo será la siguiente:

f(x) = ((10^x)-1)/9

f( 1) = ((10^ 1) - 1) / 9 = 1
f( 2) = ((10^ 2) - 1) / 9 = 11
f( 3) = ((10^ 3) - 1) / 9 = 111
f( 4) = ((10^ 4) - 1) / 9 = 1111
f( 5) = ((10^ 5) - 1) / 9 = 11111
f( 6) = ((10^ 6) - 1) / 9 = 111111
f( 7) = ((10^ 7) - 1) / 9 = 1111111
f( 8) = ((10^ 8) - 1) / 9 = 11111111
f( 9) = ((10^ 9) - 1) / 9 = 111111111
f(10) = ((10^10) - 1) / 9 = 1111111111
f(11) = ((10^11) - 1) / 9 = 11111111111
f(12) = ((10^12) - 1) / 9 = 111111111111
f(13) = ((10^13) - 1) / 9 = 1111111111111
f(14) = ((10^14) - 1) / 9 = 11111111111111
f(15) = ((10^15) - 1) / 9 = 111111111111111
f(16) = ((10^16) - 1) / 9 = 1111111111111111
f(17) = ((10^17) - 1) / 9 = 11111111111111111
f(18) = ((10^18) - 1) / 9 = 111111111111111111
f(19) = ((10^19) - 1) / 9 = 1111111111111111111
f(20) = ((10^20) - 1) / 9 = 11111111111111111111
Ejemplo 4
#include <stdio.h>

int concat();

int main(void) {
	int num = 9;

	for (int i = 9; i >= 2; i--) {
		printf("%8d * 9 + %d = %d\n", num, i-2, num*9+i-2);
		num = concat(num, i-1);
	}
}

int concat(int x, int y) {
	int temp = y;
	while (y != 0) {
		x *= 10;
		y /= 10;
	}
	return x + temp;
}

El resultado de la ejecución de este tercer ejemplo será la siguiente:

       9 x 9 + 7 = 88
      98 x 9 + 6 = 888
     987 x 9 + 5 = 8888
    9876 x 9 + 4 = 88888
   98765 x 9 + 3 = 888888
  987654 x 9 + 2 = 8888888
 9876543 x 9 + 1 = 88888888
98765432 x 9 + 0 = 888888888
Ejemplo 5
#include <stdio.h>

int main(void) {
	double num = 1;

	for (int i = 1; i <= 9; i++) {
		printf("%9.0f * %-9.0f = %0.f\n", num, num, num*num);
		num = (num*10)+1;
	}
}

El resultado de la ejecución de este cuarto ejemplo será la siguiente:

        1 x 1         = 1
       11 x 11        = 121
      111 x 111       = 12321
     1111 x 1111      = 1234321
    11111 x 11111     = 123454321
   111111 x 111111    = 12345654321
  1111111 x 1111111   = 1234567654321
 11111111 x 11111111  = 123456787654321
111111111 x 111111111 = 12345678987654321

Sorprendente, ¿verdad?

Obtención de la función canónica a partir de la tabla de la verdad

Ya vimos en un post anterior sobre el álgebra de Boole los principios mas básicos de su aplicación en sistemas digitales. En este nuevo post toca explicar los que son las funciones canónicas y cómo podemos obtenerlas a partir de la tabla de la verdad de una función lógica. Esto nos va a venir muy bien para simplificar circuitos que inicialmente pueden ser enormes y muy aparatosos. Así pues, se define término canónico de una función lógica a todo producto o suma en el que aparecen todas las variables en su forma directa o complementada . Estas son las dos formas canónicas minterm y maxterm:

1ª forma canónica minterm suma de productos canónicos.
2ª forma canónica maxterm producto de sumas canónicas.

OBTENCIÓN A PARTIR DE LA TABLA DE LA VERDAD

Supongamos que tenemos la siguiente tabla de la verdad para una función lógica de tres variables , y :

Término minterm Término maxterm
0 0 0 0 0 0
1 1 0 0 1 1
2 2 0 1 0 1
3 3 0 1 1 0
4 4 1 0 0 0
5 5 1 0 1 1
6 6 1 1 0 1
7 7 1 1 1 1

Minterms: Se toman las salidas que son "" y se expresa como suma de términos producto en los que las variables que son "" se expresan como literales y las que son "" como invertidas


Maxterms: Se toman las salidas que son "" y se expresa como producto de términos suma en los que las variables que son "" se expresan como literales y las que son "" como invertidas.


PASAR DE LA 1ª FORMA CANÓNICA A LA 2ª FORMA CANÓNICA

1. Se saca la función minterm invertida con los términos que son .

2. Se hace la inversa de la función aplicando la Ley de Morgan a los términos canónicos.


3. Se obtiene directamente cambiando los términos minúscula por mayúscula.

PASAR DE LA 2ª FORMA CANÓNICA A LA 1ª FORMA CANÓNICA

1. Se representa la función invertida, tomando los términos maxterm que no aparecen.

2. Se hace la inversa de la función aplicando la Ley de Morgan a los términos canónicos.


3. Se obtiene directamente cambiando los términos mayúscula por minúscula.

EJEMPLOS

En este ejemplo veremos qué sucede cuando uno de los términos no es canónico. ¿Que significa que "no es canónico"? Supongamos que tenemos que hallar la 2ª forma canónica de la siguiente función lógica:

Aquí parecen dos términos, y ambos han de ser obligatoriamente canónicos, pero el primero no lo es por que aparece una sola. Para que sea un término canónico han da parecer todas las variables en cada término, que en este caso son dos, y . Para hacerlo canónico podemos aplicar algunas leyes del álgebra de Boole. Por ejemplo:

En este momento ya tenemos lo tres términos de forma canónica, o lo que es lo mismo, la función tiene dos variables y y en cada término aparecen las dos variables. Para terminar, con la tabla de la verdad de esta función lógica de dos variables podremos obtener las dos funciones canónicas:

Tabla de la verdad:

Término minterm Término maxterm
0 0 0 0 1
1 1 0 1 1
2 2 1 0 0
3 3 1 1 0

Minterm:

Maxterm:

Álgebra de Boole

George Boole (1815 - 1864) desarrolló una herramienta matemática para solucionar problemas en función de dos únicos estados, el estado SI y el estado NO. Posteriormente se utilizaría para el estudio de los primeros computadores, George Booleen los que se pudo comprobar que para realizar grandes y complejos cómputos era mas factible usar la electrónica digital en vez de la analógica. Hoy lo seguimos utilizando prácticamente para todo. Puede parecer superfluo o poco relevante, pero en realidad se trata de la la base fundamental para entender las exigencias computacionales del procesamiento digital de la información, vamos ¡casi nada!

La aplicación del álgebra de Boole en computadores es del tipo binario, lo que quiere decir que el estado de un elemento del circuito lógico viene representado por una variable que puede valer "1" o "0". Un ejemplo que se utiliza siempre en sistemas digitales es un interruptor y como base una carga, de modo que si el interruptor está abierto hay tensión y si está cerrado no. Es decir, SI/NO, 1/0, verdadero/falso.

Por poner un ejemplo práctico, podremos decir "" cuando una lámpara está apagada y "" cuando está encendida.

Hoy en día tenemos sistemas de computadoras muy complicados en los que por muy complejos que estos sean al final internamente tienen un montón de dispositivos que o están abiertos o están cerrados. Lo complicado es que son muchos, y como son muchos pues nos veremos obligados a realizar funciones matemáticas con todos esos dispositivos para hacer funcionar un sistema digital.

En matemáticas, una función no es otra cosa mas que una representación de una serie de variables de entrada relacionadas ente si con las operaciones aritméticas correspondientes, de forma que el resultado de esa combinación nos va a dar un resultado de salida.

En el álgebra de Boole es exactamente igual, solo que el resultado de salida va a ser siempre un "" o un "". Por ejemplo, una función de tres variables de entrada se podría representar de la siguiente manera:

Para poder facilitar el cálculo de los valores de salida utilizaremos una "tabla de la verdad". Se trata de una tabla que recoge todas las combinaciones de las variables de entrada y los valores que toman las salidas, por ejemplo, para la función de tres variables del ejemplo la tabla de verdad sería:

0 0 0 0
0 0 1 0
0 1 0 0
0 1 1 1
1 0 0 0
1 0 1 0
1 1 0 1
1 1 1 1

En este caso desconocemos cuál es la función , pero si sabemos que salidas tiene. De momento eso no es importante, solo es un ejemplo. Ahora bien, la salida de para cada caso dependerá de la función.

A continuación se detallan algunos ejemplos de operaciones que se podrán realizar en el álgebra de Boole con puertas lógicas.

OPERACIONES EN EL ÁLGEBRA DE BOOLE

  1. Unión o adición: 
    Es la suma de las variables, en lógica equivalente a la función .OR BooleEn un circuito con interruptores podríamos representar esta operación de la siguiente forma (interruptores en paralelo):
    Unión o adiciónEn este caso, para que la salida sea y se encienda la bombilla bastaría con que uno de los dos (o los dos) interruptores o valgan . Por lo tanto, su tabla de verdad sería:
    0 0 0
    0 1 1
    1 0 1
    1 1 1
  2. Intersección o producto:
    Es el producto de las variables, en lógica equivalente a la función .AND BooleEn un circuito con interruptores podríamos representar esta operación de la siguiente forma (interruptores en serie):
    Intersección o productoEn este caso, para que la salida sea y se encienda la bombilla bastaría con que los dos interruptores y valgan . Por lo tanto, su tabla de verdad sería:
    0 0 0
    0 1 0
    1 0 0
    1 1 1
  3. Complementación o inversión:
    Esta operación es simplemente la inversión del valor de la variable a la que se le aplique, por ejemplo o .

Tabla de verdad completa con todas las operaciones:

0 0 0 0 1 1
0 1 1 1 1 0
1 0 1 0 0 1
1 1 1 1 0 0

Leyes fundamentales

A continuación unas leyes muy básicas y fáciles de recordar que nos van a servir principalmente para simplificar y de este modo para hacer circuitos más reducidos, por lo tanto más eficaces:









Conmutativa:
Asociativa:
Distributiva:
Absorción:
Morgan:
Teorema de Shannon:

Composición de circuitos

Podemos crear un circuito a partir de una función y también al revés, por ejemplo, el circuito resultante de la siguiente función sería:


Circuito ejemplo Boole

En los próximos artículos veremos la obtención de las diferentes funciones canónicas a partir de una tabla de la verdad como las que ya hemos visto aquí y también la simplificación mediante mapas de Karnaugh.

Conversor de Binario a Decimal en Java

En marzo escribí una entrada en la que explicaba la notación matemática para convertir un número binario a decimal. Aunque Java no es de mis lenguajes preferidos, quise hacer la implementación de una clase en este lenguaje para compararlo con el mismo método en otros de mas bajo nivel como C o Ensamblador. Sin mas dilación, aquí va el código:

import java.util.Scanner;
 
public class BinarioDecimal {
 
	//Atributos
	private static String binario;
	private static int longitud;
	private static int decimal;
 
	//Constructor vacío
	public BinarioDecimal() {}
 
	//Programa principal
	public static void main(String[] args) {
		binario = getStringInput("Binario: ");
		longitud = binario.length();
		int n = longitud;
		for (int i = 0; i < longitud; i++) {
			n--;
			int b = Character.getNumericValue(binario.charAt(i));
			decimal += b*(int)Math.pow(2,n);
		}
		System.out.printf("Decimal: %d", decimal);
	}
 
	//Metodo para leer la entrada por teclado
	public static String getStringInput(String prompt) {
		Scanner in = new Scanner(System.in);
		System.out.print(prompt);
		String input = in.nextLine();
		in.close();
		return input;
	}
}

He añadido alguna cosa de más, como la entrada de datos mediante la clase Scanner de java.util para hacerlo interactivo, aunque lo importante es que aquí se puede apreciar la fórmula matemática de conversión: