Java Para Cálculos Numéricos

A veces se utiliza Java para hacer cálculos numéricos intensivos a pesar de su fama de ser lento. Muchas veces la rapidez de cálculo es un aspecto importante, pero casi nunca el único a la hora de elegir el lenguaje para implementar estos cálculos. Java tiene muchas cosas atractivas para este tipo de aplicaciones como el hecho de generar programas portables, de tener muy buen soporte para programación paralela y tener «garbage collector» entre otros. Pero, ¿qué tan lento o rápido es realmente para este tipo de cálculos?

Para tener una idea decidí comparar Java con C que es el lenguaje rápido por excelencia aunque en este campo muchos consideran a Fortran como en rey de la computación numérica.

Para la comparación tomé el algoritmo de Chudnovsky, que calcula el número PI con la cantidad de dígitos que se quiera, mientras se disponga de memoria suficiente para almacenarlos.

Para el cálculo en C usé la biblioteca GMP que permite hacer cálculos con precisión arbitraria y se autoproclama como la implementación más rápida del cálculo de PI.

En Java usé la biblioteca Apfloat que también ofrece la posibilidad de hacer cálculos con precisión arbitraria y que promete ser más rápida que la implementación de BigDecimal de Java para números grandes.

En esta era de computadoras con múltiples núcleos, opté por realizar el cálculo utilizando todo el paralelismo disponible y para ello utilicé dos implementaciones paralelizadas del algoritmo en cuestión. Apfloat ya trae una implementación paralela del algoritmo de Chudnovsky como ejemplo de uso y la implementación de C está basada en OpenMP.

  • Algoritmo de Chudnovsky en C
  • Algoritmo de Chudnovsky en Java
  • Opciones de compilación de C: GCC: 4.7.2  -fopenmp -Wall -O2 -lgmp -lm
  • Versión de GMP 5.0.2
  • Versión de Apfloat 1.7.1
  • Versión Java: OpenJDK 1.7.0u15
  • Hardware: AMD Phenom X4 9550 4 GB RAM
  • Sistema operativo: Ubuntu 12.10 64 bits
  • Memoria disponible para Java 2 GB

Parámetros usados para Apfloat

  • builderFactory=org.apfloat.internal.LongBuilderFactory
  • defaultRadix=10
  • cacheL1Size=131072
  • cacheL2Size=524288
  • cacheBurst=64
  • memoryTreshold=402653184
  • sharedMemoryTreshold=268435456
  • blockSize=1048576
  • numberOfProcessors=4
  • filePath=
  • fileInitialValue=0
  • fileSuffix=.ap
  • cleanupAtExit=true

Resultados

Dígitos Apfloat (segundos)
GMP (segundos) Apfloat/GMP
100.000 0,435 0,041 10,61
1.000.000 3,987 0,595 6,70
10.000.000 67,821 8,044 8,43
100.000.000 745,732 111,824 6,67

Los resultados muestran que en esta prueba la versión Java llega a ser 6 veces más lenta que GMP.

Como Apfloat utiliza almacenamiento en disco cuando los números son más grandes que un cierto tamaño (configurable), elegir la configuración adecuada lleva un poco de pruebas. La configuración por omisión no está del todo optimizada para computadoras actuales con varios gigas de RAM.

Código de la Prueba

El código C es el citado más arriba.

Java

package calculatepi;

import java.io.File;
import java.io.IOException;
import java.io.PrintWriter;
import java.io.Writer;
import java.util.ArrayList;
import java.util.List;
import java.util.logging.Level;
import java.util.logging.Logger;
import org.apfloat.Apfloat;
import org.apfloat.ApfloatRuntimeException;
import org.apfloat.samples.Operation;
import org.apfloat.samples.Pi;
import org.apfloat.samples.PiParallel;

public class App {

    public static void main(String[] args) {
        final int max = Integer.parseInt(args[0]);
        Writer nullWritter = new Writer() {
            @Override
            public void write(char[] chars, int i, int i1) throws IOException {
            }

            @Override
            public void flush() throws IOException {
            }

            @Override
            public void close() throws IOException {
            }
        };

        try (
                PrintWriter out = new PrintWriter(nullWritter);
                PrintWriter err = new PrintWriter(new File("/home/user/apfloatPierr.txt"))) {
            Pi.setOut(out);
            Pi.setErr(err);
            final Operation warmUpOp = new PiParallel.ParallelChudnovskyPiCalculator(1000000, 10);

            final List<Operation> ops = new ArrayList<>();

            for (int i = 100000; i <= max; i = i * 10) {
                ops.add(new PiParallel.ParallelChudnovskyPiCalculator(i, 10));
            }

            for (int i = 0; i < 5; i++) {
                PiParallel.run(1000000, 10, warmUpOp);
            }
            int digits = 100000;
            for (Operation op : ops) {
                err.println("====================================");
                PiParallel.run(digits, 10, op);
                digits = digits * 10;
            }

        } catch (IOException | ApfloatRuntimeException ex) {
            Logger.getLogger(App.class.getName()).log(Level.SEVERE, null, ex);
        }

    }
}

Bash

#!/bin/sh

HEAP=2g
OUTFILE=apfloatPierr.txt
rm $OUTFILE

echo "Java..."

java -Xmx$HEAP -Xms$HEAP -jar /home/user/CalculatePi.jar 100000000
sleep 1

echo "C..."
echo "100000" >> $OUTFILE
./pi 100000 0 4 1> /dev/null 2>> $OUTFILE
echo "1000000" >> $OUTFILE
./pi 1000000 0 4 1> /dev/null 2>> $OUTFILE
echo "10000000" >> $OUTFILE
./pi 10000000 0 4 1> /dev/null 2>> $OUTFILE
echo "100000000" >> $OUTFILE
./pi 100000000 0 4 1> /dev/null 2>> $OUTFILE

Entorno de Desarrollo de C/C++ para Windows

Instalación de Cygwin

Cygwin es una colección de herramientas de GNU para Windows. Entre las herramientas de desarrollo nos interesa el compilador GCC, el debugger GBD y demás herramientas necesarias para compilar como MAKE.

Ir a http://www.cygwin.com/ y bajar el instalador de Cygwin (http://www.cygwin.com/setup.exe).
cygwin-1

Ejecutar setup.exe y seguir la instalación.


Si se está en una computadora con conexión a Internet, se debe elegir la opción “Install from Internet”.

En la carpeta llamada Root Directory se van a instalar las herramientas que seleccionemos.

En la carpeta seleccionada quedarán los paquetes que descarguemos para una futura reinstalación o para copiarlos a otra computadora e instalarlos en ella.


Se debe elegir desde dónde descargar los paquetes. Si se elige un servidor muy lento o que no funciona, se puede cambiar por otro.

Muchos paquetes ya están seleccionados, pero además se deben marcar para instalar los siguientes de la categoría devel: automake1.9, cppcheck, gcc4, gdb, make, libmpfr4




Configurar Path

Finalizada la instalación se debe agregar la carpeta donde están las herramientas al Path de Windows. Yendo a Panel de ControlSystem and SecuritySystem y eligiendo Advanced system settings.

En la solapa Advanced, seleccionar el botón Environment Variables….

Se debe crear o editar una variable de usuario seleccionado en botón New… o Edit… debajo del recuadro User variables for….
Si la variable Path ya existe, se la edita y si no existe, se la crea nueva.

En caso de crearla, se le pone el nombre Path y como valor el nombre de la carpeta elegida para instalar Cygwin. Si la variable Path ya existía, simplemente al modificar su valor agregar al final un punto y coma y luego el nombre completo de la carpeta: C:\path\anterior;C:\Otro\Path\anterior;C:\cygwin\bin

Instalación de un Editor

Para escribir nuestros programas sirve cualquier editor que esté hecho para programar. Típicamente estos editores resaltan la sintaxis del lenguaje y permiten mostrar números de línea, entre otras cosas. Gedit es uno de ellos. El bloc de notas de Windows no es recomendable.

Está disponible para Windows.

Su instalación es muy sencilla.







Uso de las Herramientas

Cygwin crea un entorno UNIX completo dentro de Windows. Este incluye, además de las herramientas de programación que instalamos, un shell. Para usar el compilador y las herramientas de Cygwin debemos hacerlo mediante ese shell al que accedemos ejecutando Cygwin Terminal. El directorio “activo” inicialmente al abrir Cygwin Terminal está ubicado dentro de la carpeta donde lo instalamos. Podemos situar nuestros programas ahí o en cualquier otra parte y movernos hacia esa otra carpeta usando los comandos del shell. En la imagen se observa un programa guardado en la carpeta “c”.

Para compilar y ejecutar el programa que escribimos, debemos abrir Cygwin Terminal. No es lo mismo abrir el shell de Windows “cmd” ya que Windows no soporta algunas de las características que tiene UNIX como los links simbólicos y es probable que tengamos errores extraños al ejecutar las herramientas de Cygwin.

En este ejemplo para compilar y ejecutar el programa que escribimos y guardamos debemos ejecutar estos comandos:

  • cd c
  • gcc -o helloworld -Wall -std=c99 helloworld.c
  • ./helloworld

Instalación sin acceso a Internet

Para instalar en una computadora sin conexión a Internet se debe primero ejecutar el instalador en otra donde sí se disponga de conexión y utilizar la opción
Download Without Installing. Luego proseguir con la instalación. El instalador va a descargar los paquetes seleccionados en una carpeta local. Esa carpeta puede ser copiada a otra computadora que no tenga conexión a Internet.

En la computadora sin conexión a Internet se debe ejecutar el instalador y seleccionar la opción Install From Local Directory.

Luego se debe indicar en qué carpeta se copiaron los paquetes descargados anteriormente y proseguir con la instalación.

Aritmética Binaria y Decimal

Por Qué Nunca Usar float y double para Representar Montos de Dinero

Cuando aprendemos a programar una de las primeras lecciones es que se usan diferentes tipos de datos para representar números enteros y números no enteros. Una variable entera no puede tener el valor 0,5. En Pascal se usan los tipos Integer y Real y en C están por ejemplo los tipos int, long, float y double.

Supongamos que tenemos que hacer un sistema de permita administrar sumas de dinero; más específicamente, deudas y pagos. Puede ser para un banco, una tarjeta de crédito, un almacén o cualquier aplicación comercial. Para dicho sistema tenemos que hacer una función en C muy sencilla que reciba la deuda actual y un arreglo de pagos que realizó el deudor y que devuelva la deuda actualizada una vez descontados los pagos. Este podría ser el prototipo:

double actualizarDeuda(double deuda, double pagos[], int cantPagos);

Como la deuda puede ser con centavos, no podemos usar int para representarla entonces usamos double.

Lo que hace la función es restar de la deuda cada elemento del arreglo pagos y devolver el monto actualizado.

double actualizarDeuda(double deuda, double pagos[], int cantPagos){
    double deudaActualizada = deuda;
    int i;
    for(i = 0; i<cantPagos;i++){
        deudaActualizada -= pagos[i];
    }
    return deudaActualizada;
}

Vamos a constatar que la función hace lo esperado con una prueba.

Supongamos que un deudor tenía una deuda de $308,21 y realizó 4 pagos por $8,01; $0,20; $299,99 y $0,01. Si bien no es sencillo pagar en efectivo montos como $8,01, es un caso que puede darse perfectamente en transacciones electrónicas.

double deuda = 308.21d;
double pagos[4] = {8.01d, 0.20d, 299.99d, 0.01d};

En teoría la deuda está saldada así que la función debería devolvernos que la deuda es $0,00.

double miDeuda = actualizarDeuda(deuda, pagos, 4);
if(miDeuda == 0.0d){
    printf("Libre deuda.\n");
} else {
    printf("Al Veraz.\n");
}

Sin embargo la prueba imprime Al Veraz. es decir que nuestro sistema sigue registrando una deuda a pesar de que está saldada.

La nueva deuda no es $0,00, pero casi. El valor es $-0,00000000000000909516. Es decir, nuestra función nos calculó que el deudor pagó 9 milbillonésimas de peso más de lo que debía.

Esto nos hace pensar que hay un error de redondeo así que decidimos usar la función round.

miDeuda = round(actualizarDeuda(deuda, pagos, 4));
if(miDeuda == 0.0d){
    printf("Redondeando, logro el libre deuda.\n");
} else {
    printf("Redondeando, igual voy al Veraz.\n");
}

Ahora sí obtenemos el resultado Redondeando, logro el libre deuda lo que es en principio correcto.

A pesar de que en muchos sistemas el problema de la precisión está “solucionado” de esa manera, no se está atacando la causa del problema sino una de sus consecuencias. Como se observa a continuación, redondear el resultado puede llevarnos a dar por saldada una deuda que todavía no está totalmente paga.

pagos[2]-=0.4d;
miDeuda = round(actualizarDeuda(deuda, pagos, 4));
if(miDeuda == 0.0d){
    printf("Con redondeo, debiendo $0,40 igual logro el libre deuda.\n");
} else {
    printf("A pesar de redondear, como debo voy al Veraz.\n");
}

Si descontamos $0,40 de uno de los pagos, está claro que no deberíamos extender el libre deuda. Sin embargo, obtenemos el cartel Con redondeo, debiendo $0,40 igual logro el libre deuda lo que es incorrecto y le hace perder dinero quien nos encargó hacer el sistema. Es una situación complicada para nosotros ya que desde el punto de vista del acreedor, es una atrocidad.

Este problema lo tenemos en cualquier lenguaje de programación y también lo tenemos en muchas aplicaciones muy usadas para hacer cálculos con montos de dinero.

arit-excelLa única solución que ofrecen al problema en el caso de esa planilla de cálculo tan famosa es “redondear” (http://support.microsoft.com/kb/214118/es).

La causa del problema está en la representación binaria de los números de coma flotante. Como no todos los números son representables se tiene que buscar el número más cercano al especificado y eso implica un error en la representación. Por ejemplo el número 0,11 no es representable como double (porque en binario es periódico) y es aproximado con 0,11000000000000000056.

Los números de coma flotante están definidos por el estándar IEEE-754 que data de 1985. A partir de 2008, el estándar incorporó además de los números de binarios de coma flotante que ya existían (que usan base 2 como los float y double) nuevos números decimales de coma flotante (que usan base 10).

Esa es la solución al problema de la precisión que estamos teniendo. Vamos a reescribir la función de esta forma:

_Decimal64 actualizarDeuda(_Decimal64 deuda, _Decimal64 pagos[],int cantPagos){
    _Decimal64 deudaActualizada = deuda;
    int i;
    for(i = 0; i<cantPagos;i++){
        deudaActualizada -= pagos[i];
    }
    return deudaActualizada;
}

Usamos el tipo _Decimal64 en lugar de double. El resto es igual.

Y ahora sin aplicar ningún redondeo, las cuentas simplemente dan el resultado esperado.

_Decimal64 deuda = 308.21dd; 
/* el sufijo dd es el que corresponde al tipo _Decimal64 */
_Decimal64 pagos[4] = {8.01dd, 0.20dd, 299.99dd, 0.01dd};
_Decimal64 miDeuda = actualizarDeuda(deuda, pagos, 4);
if(miDeuda == 0.0dd){
    printf("Con _Decimal64, obtengo en libre deuda normalmente.\n");
} else {
    printf("Con _Decimal64, al Veraz.\n");
}

El soporte de decimales de coma flotante está disponible a partir de GCC 4.2, pero no está disponible para todas las plataformas. En particular en Windows con Cygwin obtenemos un error que dice “error: decimal floating point not supported for this target”.

Otros compiladores sí lo soportan en Windows como por ejemplo el de Intel. http://software.intel.com/en-us/articles/using-decimal-floating-point-with-intel-c-compiler/.

Más información

  1. http://en.wikipedia.org/wiki/IEEE_754-2008
  2. http://speleotrove.com/decimal/

Programar en C con Garbage Collection

De los 20 lenguajes mejor ubicados en el índice TIOBE, 14 tienen alguna forma de “Garbage Collection“. Esos 14 lenguajes representan más de un 47% de peso en el índice. Los únicos lenguajes “importantes” que no tienen recolección automática de basura son C, C++ y Objective-C (Delphi, Pascal y Ada son minoritarios).

El caso de Objective-C es especial ya que desde la versión 2.0 tiene un garbage collector, pero es opcional y no está en la versión que se usa para desarrollar aplicaciones para el iPhone.

En C y C++ se puede usar el garbage collector de Boehm-Demers-Weiser. Básicamente en lugar de llamar a malloc, llamas a GC_MALLOC y los llamados a GC_FREE son opcionales.

En Ubuntu, debemos instalar libgc

 sudo apt-get install libgc-dev libgc1c2

y al compilar debemos usar la opción -lgc para indicarle al linker que use libgc.

gcc -std=c99 -Wall  -o programa *.c -lgc

La función scanf y la lectura de la entrada de teclado

El objetivo es explicar por qué nunca hay que usar la función scanf para leer la entrada del teclado.

Cuando se lee la entrada generada por el teclado se deben tener en cuenta algunos aspectos:

  1. lo leído es generado por un ser humano que como tal se equivoca;
  2. el programa no accede al teclado mismo sino que es el sistema operativo el que le pasa los datos ingresados de alguna manera.

El primer aspecto ya alcanza para descartar scanf como herramienta válida para leer la entrada del teclado. La función está hecha para leer con un formato determinado y es imposible estar seguros de que un humano va a respetar ese formato sin ningún error. A pesar de esto, miles y miles de ejemplos de código que lee de teclado en C usan scanf aportando confusión y haciendo caer en el error a quienes recién se inician.

El otro aspecto es más técnico y requiere más explicación. Cuando se dice que programa lee del teclado en realidad lo que está haciendo es leer de un área de memoria llamada “buffer del teclado”. Por la naturaleza de los teclados su buffer es un “buffer de línea”, eso significa que los datos que vienen del teclado se incorporan de a líneas enteras y no carácter por carácter. Es por eso que cuando leemos de ahí, por más que ingresemos muchos caracteres, hasta que no ingresamos el fin de línea (enter) el programa no lee nada y se queda trabado esperando que haya algo en el buffer.

La forma en la que se llena el buffer del teclado nos permite estar seguros de 2 cosas:

  1. si el buffer no está vacío, sí o sí hay al menos un fin de línea;
  2. el buffer siempre tiene un carácter de fin de línea al final.

El primer problema que presenta scanf es que no siempre lee el carácter de fin de línea. Eso depende del formato que se usa. Cuando se usan “%d” o “%s” el fin de línea queda en el buffer y será leído en la próxima invocación a una función que lea del teclado. Esto genera incertidumbre en cuanto a qué hay en el buffer al momento de ir a leer. La primera vez está vacío, pero en las veces subsiguientes puede ser que tenga algunos caracteres que quedaron sin leer.

El segundo problema es que scanf no lee hasta el fin de línea (aunque lo deje sin leerlo) sino que lee hasta que se acaba el formato especificado. Si se le pide que lea un entero (“%d”) y se ingresan 2 enteros separados por un espacio o por una letra lee el primero y deja el resto en el buffer. Nuevamente esto genera incertidumbre en las lecturas subsiguientes no sabemos si quedaron caracteres ingresados anteriormente.

El tercer problema es más grave porque puede ocasionar que se cuelgue el programa. Scanf no controla que la cantidad de caracteres leídos sea menor o igual al tamaño del arreglo o string que se le pasa por parámetro. Si se esperan leer 10 caracteres es típico leerlos en un arreglo de 11 (uno para el fin del string), pero si el usuario ingresa 11 caracteres, scanf va a escribir más allá del fin de arreglo con el consiguiente error de memoria. Si usamos un formato que especifique la longitud máxima entonces el problema no será pasarse del tamaño del arreglo, sino que nuevamente dejaremos sin leer caracteres en el buffer.

El último problema es que scanf no permite detectar la mayoría de los errores de ingreso de datos del usuario. Cuando mandamos a leer un entero (“%d”) y el usuario ingresa letras, scanf nos devuelve un entero fruto de una conversión de las letras a número. ¿Cómo determinar si el usuario ingresó letras o realmente ese número? Imposible.

Para leer en forma correcta del teclado, es necesario asegurarse de 3 cosas.

  1. Dejar el buffer vacío luego de la lectura (leer siempre hasta el fin de línea). Es la única manera de evitar que las lecturas posteriores lean caracteres ingresados anteriormente y que no esperan.
  2. Leer carácter por carácter almacenándolos en un arreglo (o donde se quiera) teniendo en cuenta el tamaño disponible. Los caracteres que no entren igual deben ser leídos (para cumplir el punto anterior). No hay forma de predecir cuantos caracteres se van a leer.
  3. Ver de qué manera se puede informar que lo leído no es lo esperado de manera de que no se confundan un valor válido con uno erróneo. Por ejemplo, si se están leyendo números y se ingresan letras, devolver -1 puede interpretarse como un valor que señala el error o como un valor leído correctamente. Hay que elegir un valor que no pueda ser confundido con uno válido.

A modo de ejemplo del punto 3, la función getchar a pesar de que que lee caracteres el tipo de retorno no es char sino int. De esa manera cuando devuelve un valor -1 no puede confundirse con un carácter ya que estos son todos positivos.

Alguno se preguntará para qué existe scanf que es una función que lee de teclado, si no hay que usarla para ello. La respuesta radica en un detalle: scanf y algunas otras funciones como getchar no leen de teclado, sino que leen de la entrada estándar. Si no se indica lo contrario, la entrada estándar proviene del teclado, pero es muy simple hacer que provenga de un archivo cualquiera. Hay muchos programas en C que jamás se invocan sin establecer que la entrada estándar provenga de algún archivo. En esos casos, scanf puede ser útil, sobre todo si el archivo fue generado con printf.

Actualizado: Hace un tiempo publiqué cómo leer un entero de teclado. Ahora le corregí un bug y le mejoré el uso de memoria.

Leer un entero de teclado en C

Sin usar scanf, porque no sirve para leer de teclado.

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

#define MAX_DIGITS 9
#define BASE 10

static int* BUFFER = NULL;

int* readInt()
{
int i = 0,
j = 0,
base = 1,
factor = 1,
*result,
ch = getchar();
if(BUFFER == NULL){
BUFFER = malloc(MAX_DIGITS * sizeof(int));
}
if(ch == 45){
factor = -1;
ch = getchar();
}
while(i < MAX_DIGITS && ch != 10 && ch != -1 && 48 <= ch && ch <= 57){
BUFFER[i++] = ch;
ch = getchar();
}

if(i == 0 || (ch != 10 && ch != -1)){
while(ch != 10 && ch != -1){
/* consumo lo que haya en el buffer para vaciarlo.*/
ch = getchar();
}
return NULL;
}
/* transformo el string en un entero según la base.*/
result = (int*) malloc(sizeof(int));
*result = 0;
/* De atrás para adelante, multiplicando por la base y sumando.*/
for(j = i – 1; j >= 0; j = j-1){
*result += ( BUFFER[j] – 48 ) * base;
base *= BASE;
}
/* aplico el signo*/
*result *= factor;
return result;
}