miércoles, 26 de septiembre de 2012

El número de Reynolds (II): difusión superlumínica

En la entrada anterior hablé de los procesos de advección y difusión de un contaminante, la diferente importancia de cada uno y de cómo el número de Péclet (Reynolds) muestra la importancia de cada uno.

Al final había comenzado a hablar de cómo las ecuaciones de Navier-Stokes (NS) no son más que ecuaciones de conservación del momento con sus términos de advección y difusión, cuya importancia relativa está dada por el número de Reynolds.

En esta entrada voy a volver un momento atrás y comentar algo que se me pasó por alto en la anterior. No es que no me acordase de contarlo sino que no se me ocurrió hasta tiempo después.

Difusión superlumínica

Voy a empezar con un gráfico similar al de la entrada anterior, pero en esta ocasión he añadido una línea que representa un punto que se aleja del punto de interés a la velocidad de la luz.


Un físico no se encuentra cómodo hasta que no mete la velocidad de la luz en algún sitio.

Como dije anteriormente, para cualquier velocidad, existe un radio en las cercanías de cualquier punto, dentro del cual la difusión es la forma más rápida de "transportar información". Evidentemente, esto no puede ser correcto físicamente, es sólo un resultado de nuestro modelo de la difusión/viscosidad.

Esto también explica algo que afirmé en la anterior entrada, pero que no me molesté en argumentar: que en escalas pequeñas la difusión deja de ser una aproximación válida.

¿Cómo de pequeñas? Veámoslo. Vamos a suponer que hablamos de agua y de viscosidad. ¿A que distancia máxima la difusión es más rápida que la velocidad de la luz?

Una forma de hacerlo es calcular la distancia L que hace que el número de Reynolds sea 1 para la velocidad de la luz:

\begin{aligned} Re = 1 = \frac {Lc}{\nu} \end{aligned}

Usando valores estándar:

\begin{aligned} L = \frac {\nu}{c} = \frac {10^{-6} (m^2/s)}{3~10^5 m/s}\approx 10^{-11} m \end{aligned}

Es decir, más o menos la décima parte de un átomo de Helio. Parece que esto no nos va a dar problemas, al menos para el caso del agua. La brea tiene una viscosidad que es casualmente \(10^{11}\) veces mayor que la del agua, así que en ese caso, el radio es del orden de un metro. ¿Cómo tendrán en cuenta eso en las ecuaciones?

Me habrán dado un Ig-Nobel por parsimonioso, pero la difusión va que se las pela

En realidad, la difusión no puede ser más rápida que el proceso físico que la crea, que no es otro que la agitación térmica de las moléculas del agua. Nos quedamos sin ver radiación Cherenkov asociada a la difusión.

Sin embargo, la razón por la que decidí escribir esta entrada no es sólo porque me pareció curioso. Hay un motivo más fundamental.

Una de las cosas que he oído repetir como un mantra desde que llegué a Cornell es: "una vez que escojas un modelo de la realidad, olvídate de la realidad". La idea es que una vez que has hecho una serie de hipótesis simplificadoras no es lógico usar argumentos que caen fuera de esas simplificaciones.

En este caso, una de las hipótesis simplificadoras es la hipótesis del continuo: todas las magnitudes consideradas son continuas a todos los efectos. Por lo tanto, debemos olvidarnos de moléculas, agitación térmica y demás.

Es decir, la hipótesis del continuo implica la transferencia de información a velocidades mayores que la de la luz. Es parte de nuestro modelo y tenemos que vivir con ello. Al fin y al cabo, el suponer que un fluído es incompresible también tiene la misma implicación, y nadie se queja. Qué le vamos a hacer.







jueves, 20 de septiembre de 2012

El número de Reynolds (I): Introducción.

Aunque tengo por ahí empezada una serie sobre las ecuaciones de Navier-Stokes, quiero hablar un poco sobre el número de Reynolds sin tener que entrar a fondo en las ecuaciones, así que empiezo una serie nueva.

En esta primera entrada voy a hablar de conceptos básicos para luego ir entrando en materia. De hecho, antes de hablar del número de Reynolds en sí, voy a tratar un tema similar pero más sencillo.

Difusión y advección de sustancias


Supongamos que tenemos dos situaciones diferentes, en una tenemos agua inmovil, en otra agua moviéndose a velocidad constante \(u\). En el instante t=0 ponemos una gota de sustancia en un punto y esperamos a que pase un cierto tiempo. Supongamos además que la sustancia del primer caso tiene un coeficiente de difusión en el agua de \(\nu\), mientras que la del segundo tiene un coeficiente nulo (no se difunde).

Para los físicos, advección = convección. Para los ingenieros, sustancia = contaminante. Para todos 2D=3D.

La figura muesta el comportamiento de la sustancia en ambos casos. En el primero se diluye debido a la agitación térmica de las moléculas de agua, en el segundo avanza sin cambiar de forma a la misma velocidad que el agua. La ecuación que rige ambos casos (aunque no es necesaria para entender lo que sigue) es:

\begin{aligned} \frac {\partial c}{\partial t} + \nabla \cdot (c\mathbb{u}) = {\nu}\nabla^2c \end{aligned}

Una forma de comparar ambas situaciones es estudiar cuanto avanza la sustancia en un tiempo \(t\) desde el punto en el que fue depositada. En el caso de la difusión es necesario hacer una suposición adicional, así que definiré el avance como la distancia a la que la concentración es igual a la mitad de la concentración en el centro.

En el caso de la advección, claramente la sustancia avanza a una velocidad  \(u\) o, en términos de la distancia recorrida \(L_{adv}=ut\).

No voy a entrar en detalles de cómo se calcula, pero en el caso de la difusión la distancia recorrida es \(L_{dif}=\sqrt{4\nu t}\). Intuitivamente, es fácil ver que la difusión es menor cuanto más aumenta el tamaño de la mancha y disminuyen los gradientes, así que una dependencia con la raíz cuadrada del tiempo parece razonable.

La siguiente gráfica muestra la comparación entre ambos avances:



Como puede verse, una consecuencia interesante de que la difusión avance con la raíz cuadrada del tiempo es que existe un radio \(L_0(u,\nu)\) en el cual la difusión es más rápida que la advección. Esto es cierto independientemente de que velocidades o coeficientes de difusión tengamos; siempre habrá un valor de \(L_0\) dentro del cual la difusión avanza más rápido que la advección. Si pensamos en términos de información, todo punto dentro de ese radio recibirá noticia de la sustancia depositada más rapido debido a la difusión que a la advección.

Para tener en cuenta este efecto, se suele usar el número de Peclet;

\begin{aligned} Pe_L=\frac {LU} {\nu} \end{aligned}

Que no es más que un número de Reynolds aplicado a la difusión de sustancias en lugar de la difusión del momento (viscosidad). Como yo he usado la misma letra \(\nu\) que se suele usar para la viscosidad, ambos números tendrán la misma expresión.

También lo he hecho porque nunca me ha gustado dar nombres diferentes a cosas similares (¿series de Taylor y McLaurin?). Pienso que no hace más que confundir, el número de Peclet no es más que una versión del número de Reynolds.

Siguiendo con lo que estábamos, en la expresión anterior \(L\) es una longitud, digamos que entre el punto donde se ha vertido una sustancia y un punto de interés, \(U\) una velocidad característica del flujo de agua y \(\nu\), como ya he dicho, el coeficiente de difusión. Para entenderlo mejor podemos usar las expresiones anteriores para calcular el tiempo que tarda en avanzar la concentración desde el punto inicial a un punto a una distancia \(L\):

\begin{aligned} t_1 = \frac {L}{U},~t_2 = \frac{L^2}{4\nu} \end{aligned}

Dividiendo \(t_2\) entre \(t_1\):

\begin{aligned} \frac {t_2}{t_1} = \frac {LU}{4\nu} \end{aligned}
Es decir, el número de Peclet es una relación entre los tiempos que tarda en llegar la información. Si es mucho menor que uno significa que en nuestro problema, la difusión es más importante que la advección y su es mucho mayor que uno, todo lo contrario.

Si utilizamos \(L\) y \(U\) para adimensionalizar la ecuación de advección difusión obtendremos:

\begin{aligned} \frac {\partial c}{\partial t} + \nabla \cdot (c\mathbb{u}) = \frac{1}{Pe_L}\nabla^2c \end{aligned}

Es decir, el número de Peclet indica la importancia del término difusivo respecto al resto de términos.

Antes de pasar a otra cosa, quiero hacer un último comentario. La definición del número de Peclet (y del número de Reynolds) se basa en una cantidad objetiva \(\nu\), otra más o menos objetiva \(U\) y finalmente una que es subjetiva (en el sentido de que tenemos libertad para elegirla) \(L\). Esto tiene utilidad desde el punto de vista práctico, ya que en la realidad solemos estar interesados en calcular la difusión y advección entre el punto de vertido y la zona de interés.

Pero a nivel de procesos físicos, no existe una \(L\) única. La sustancia está, en cada punto y a cada instante, sujeta a procesos de advección y difusión locales, sin que la distancia \(L\) tenga ningún significado físico.

Más aún, podríamos argumentar que en cada punto, la difusión es infinitamente más importante que la advección (ver la gráfica anterior). ¿Cómo es posible que exista advección si localmente la difusión es infinitamente más importante?

Parte de la respuesta es que la difusión es siempre una aproximación (no es válida en un entorno infinitesimal), otra parte es que la ecuación es lineal (así que podemos ver ambos procesos como independientes) y luego está la parte de la respuesta que desconozco.

Pero no hablemos de cosas tristes. Quedémonos con la idea de que en cada punto existe un radio \(L_0\) dentro del cual la difusión es más importante que la advección.

Las ecuaciones de Navier-Stokes

La ecuación de conservación del momento, en forma no dimensional, puede escribirse así:

\begin{aligned} \frac {\partial \mathbb{u}}{\partial t} + \nabla \cdot(\mathbb{u}\mathbb{u}^T) = -\frac{1}{\rho}\nabla p +  \frac{1}{Re} \nabla^2(\mathbb{u}) + \mathbb{g} \end{aligned}

Donde \(Re=\frac {LU} {\nu}\)

Como es sabido, esta ecuación no es más que una ecuación de advección-difusión para el momento, más un par de términos adicionales: presión y gravedad. En esta ocasión, el número de Reynolds sustituye al número de Peclet, que como ya había dicho, son más o menos la misma cosa.

Continuara...








viernes, 14 de septiembre de 2012

En el principio fue el código máquina (II) : números y constantes

Como ya dije en un artículo anterior, me encuentro actualizando un código para incluir en él transporte de sedimentos. Se trata de un programa que originalmente "sólo" era capaz de simular multiples fluídos inmiscibles, incluyendo el efecto de la temperatura, electromagnetismo y otras cosas que nunca he tenido tiempo de mirar con detalle. Afortunadamente, el código está bastante bien compartimentado.
El programa inicialmente se desarrolló en Los Álamos. Un  insituto especializado en modelado numérico y que a la mayoría os sonará porque se creó para desarrollar la bomba atómica. Si en algún sitio deberían tener buenos programadores -pensaba yo- sería allí.

Stanisław Ulam, Richard Feynman y John Von Neumann en Los Alamos
Sin embargo, la realidad es que el programa es un caos, como todos los modelos numéricos que he tenido la desgracia de tener que modificar. Hay que decir que el diseño inicial no es del todo malo (aunque claramente no pensaban en eficiencia cuando lo hicieron); pero el resultado final es horrible: Trozos de código repetidos, funciones que hacen lo mismo que otras, sin explicación de por qué. Códigos que cláramente no hacen nada o hacen cosas erróneas... La lista sigue. El hecho de que el programa funcione aún así, da idea de la cantidad de redundancia que tiene.

En cualquier caso, no empecé este artículo para quejarme del programa en general (tal vez en otro artículo), sino para comentar un detalle que me parece curioso. Echemos un vistazo al módulo de constantes:
! Real whole numbers
  real(real_kind), parameter :: zero          =   0.0
  real(real_kind), parameter :: one           =   1.0
  real(real_kind), parameter :: two           =   2.0
[...]
  ! Real exponential numbers
  real(real_kind), parameter :: ten_tominus14 = 1.0e-14
  real(real_kind), parameter :: ten_tominus12 = 1.0e-12

[...]
  ! Real fractions
  real(real_kind), parameter :: one_168th     = one/168.0
  real(real_kind), parameter :: one_105th     = one/105.0
[...]
  ! Real constants
  real(real_kind), parameter :: pi  = 3.1415926535

Los [...] son míos, el código define todas las constantes que nos podamos imaginar. La última que he puesto, es la única que parece tener algo de sentido. La razón por la que la puse es porque me recordó el famoso comentario en un manual de Fortran de Xerox:

The primary purpose of the DATA statement is to give names to constants; instead of referring to pi as 3.141592653589793 at every appearance, the variable PI can be given that value with a DATA statement and used instead of the longer form of the constant. This also simplifies modifying the program, should the value of pi change.
Supongo que los programadores de Los Alamos, aún más cautelosos, tuvieron en cuenta la posibilidad de que números cómo 1 ó 2 cambien un día de estos. Nunca se es lo suficientemente precavido.

Lo cierto es que no se me ocurre ninguna razón válida para hacer algo así. Quizás la razón sea que así se aseguran de que one es un double (real(real_kind)) y no un float (real(4)), aunque hay pocas situaciones en las que eso sea un problema y siempre es posible escribir: 1.0d0.

En cualquier caso, una cosa me llevó a pensar otra y de ahí a una idea para un artículo. ¿Cómo manejan los compiladores las constantes reales?

A diferencia de las constantes enteras, que se pueden incorporar directamente en el código máquina; las reales se suelen guardar como datos (aunque en el segmento de código). Es decir, cada vez que en nuestro código hacemos algo como:


x=2.3*y;

El compilador crea una entrada para esa constante. Si le echamos un vistazo al ensamblador que genera el compilador, veremos una lista con todas los números (sin repetir) usados en el programa. Lo mismo pasa con los parameters listados antes.

Es decir, los números reales no se incluyen en el programa cada vez que usan (como sucede con los enteros), sino se guardan en una única dirección de memoria... memoria que se puede modificar.

Dicho de otra forma, es posible escribir código que modifique todos los 1.0 de nuestro programa (o todos los ones) y los cambie por otro número.

Aprovecho el impacto que la frase anterior os habrá provocado, para cambiar momentáneamente de tema.

Parte de lo que he hecho con el programa es escribir todo el código nuevo, y reescribir parte del viejo, en c++. Esto tiene sus ventajas e inconvenientes. Las ventajas las voy a dar por sabidas, las desventajas están relacionadas con la interacción entre ambos lenguajes, ya que no existe un estándar. Lo que significa que hay que tener en cuenta diferentes compiladores.

En el caso de gcc+gfortran, las funciones Fortran siempre pasan los parámetros por referencia. Incluso cuando el parámetro es un número, la función pasa el puntero a la dirección en memoria donde ese número está guardado. Si uno no se sale de Fortran no pasa nada, pero cuando una función c++ recibe el parámetro, recibe el puntero.

Afortunadamente, como dije antes, esos datos están en el segmento de código, que es de sólo lectura, así que normalmente el mayor peligro que tenemos es que obtengamos un segmentation fault que no sabemos de donde viene.

Pero claro, nadie dijo que no se podían cambiar las propiedades del segmento de código. Eso me trae de vuelta al tema original. Veamos el siguiente pedazo de código:

program hello
    interface
      subroutine func(d)
          real(8), intent(in)   :: d
      end subroutine func
    end interface


    print *, '2 + 2 = ',4.0D0
    call func(4.0D0);
    print *, '2 + 2 = ',4.0D0

end program

El código en rojo define la interfaz a una función en c++ hecha por mi. En particular,  intent(in) significa que el parámetro es sólo de entrada y que no debemos esperar efectos secundarios.

Antes de poner el código c++, voy a mostrar el resultado de ejecutar el código anterior:
 2 + 2 =    4.0000000000000000    


 2 + 2 =    5.0000000000000000    

Process returned 0 (0x0)   execution time : 0.003 s

No hace falta saber Fortran para darse cuenta de que algo ha ido muy mal. El número 4 ha pasado a ser 5 y así seguirá durante todo el programa. Si esto no es un regalo para un desarrollador de librerías psicópata, no sé lo que es.


Empezamos cambiando constantes y luego pasa lo que pasa

El código de la función en c++ es el siguiente:

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <unistd.h>
#include <sys/mman.h>
#include <stdint.h>
#include <math.h>
extern "C" void func_(double &d)
{
 void *page =
     (void *) ((unsigned long) (&&checkpoint) &
        ~(getpagesize() - 1));

 /* mark the code section we are going to overwrite
  * as writable.
  */
 mprotect(page, getpagesize(), PROT_READ | PROT_WRITE | PROT_EXEC);

 /* Use the labels to avoid having GCC
  * optimize them out */
 switch ((int)(d)) {
   case 33:
     goto checkpoint;
   case 44:
     goto checkpoint2;

 }

 /* Replace code in checkpoint with code from
  * newcode.
  */
 //memcpy(&&checkpoint, &&newcode, &&newcode_end - &&newcode);

checkpoint:
 printf("\n");

checkpoint2:
 printf("\n");


  d=5.0;
}

La parte en rojo la copié de un ejemplo sin intentar averigüar cómo funciona. Seguro que se puede hacer mucho mejor.

lunes, 10 de septiembre de 2012

Expression templates y lazy evaluation

Estos días estoy dedicandome a escribir mi programa (simulación de transporte de sedimentos por oleaje, para los interesados) usando técnicas de template metaprogramming. Más concretamente, los expression templates. Ambas cosas tienen fama de complicadas (y hay que reconocer que lo son), pero también son muy interesantes. En este artículo voy a intentar explicar un poco cómo funcionan y por qué merecen la pena.

El artículo me ha salido largo y un tanto árido; pero quiero insistir en que creo que merece la pena ;-)

Templates

Lo primero que voy a hacer es hablar un poco de los templates según c++. Como casi todo en c++, tienen una sintaxis compleja y su funcionamiento y reglas lo son más, de hecho hay muchas cosas que aún no entiendo.

Otros lenguajes han optado por sistemas más sencillos, que en este artículo voy a englobar con el nombre de generics.

El problema es que los generics sirven para lo que hace años, la mayoría pensábamos que servían los templates; es decir, para cosas como:

vector<double> v;

O dicho de otra forma: para especificar clases que dependían de otros tipos con la intención de optimizar el código. Esto es lo que se llama instanciación explícita de los templates.

Por otro lado, c++ hablaba de la instanciación implícita o automática. El tipo de cosas de c++ que uno suele pasar por alto, porque las explicaciones de lo que significa suelen tener la luminosidad poética del BOE. Afortunadamente, otros no lo pasaron por alto. La idea básica es que si uno tiene una función template del tipo:

template<class T> sqr(T val) { return val*val; }

Cuando el código hagamos sqr(a), el compilador obtendrá el tipo de a y instanciará la versión adecuada de sqr, sin que tengamos que decirle el tipo explícitamente. Las reglas sobre cómo se deduce el tipo y cómo se escoge la versión adecuada del template son complicadas y no voy a entrar en ellas (por que no me las sé y lo que hago es procurar evitar ambigüedades).

¿Cuál es la ventaja de esto? Es difícil de explicar. Para hacerlo voy a usar el ejemplo que conozco mejor, y para ello abro otro tema.

La sobrecarga de operadores: ayer y hoy

Una de las novedades que trajo el c++ al mundo de la programación (al menos para el gran público) fue la sobrecarga de operadores. Por fin uno podía definir una clase de números complejos y operar con ella como si fueran números reales. Lo mismo podía decirse de vectores, matrices, etc.

Pero había un problema, cada vez que se realiza una operación hay una llamada a una función y se devuelve un valor. Algo que no es muy importante en el caso de objetos pequeños como números complejos, pero que es un problema grande en el caso de clases de matrices y vectores (pensemos en vectores de cientos de miles de elementos). Supongamos la siguiente operación:

\begin{aligned} a\times(x+2\times(y+z)) \end{aligned},

Donde todas las variables son vectores y se supone que las operaciones ocurren elemento a elemento. En este caso, la suma y+z genera un nuevo vector, la multiplicación por 2 otro más, la suma con x otro y, finalmente, el producto por a, otro más. Ni siquiera es seguro que esos temporales no se destruyan hasta el final de la operación, con lo cual tendremos la memoria ocupada con temporales innecesarios.

Otro problema asociado es que este tipo de expresiones no son adecuadas para la optimización. En una expresión como la anterior pero para números reales, el compilador puede encontrar todo tipo de optimizaciones. Por ejemplo, el producto por dos puede ponerse como una suma del valor consigo mismo. En el caso de vectores, no es posible.

El efecto de todo esto es que en c++, si queríamos hacer un código "bonito" teníamos que renunciar a la optimización, y no hablamos de porcentajes pequeños. En las pruebas que he hecho he visto reducciones de velocidad de menos de un 50%.

Así estaba la situación hasta que a alguien se le ocurrieron los expression templates. La idea es incluir en nuestras clases de vectores un método double Eval(index). Este método lo único que hace es devolver el valor en la posición index del vector.

El siguiente paso es crear un operador template como este:

template<class L, class R>
ExprAdd<L,R> operator+(const L & op1, const R & op2)

{
  return
ExprAdd<L,R>(op1,op2);
}
Basicamente, lo que hace es devolver una clase que construye sobre la marcha basándose en los operadores. En el resto de la explicación voy a hablar sólo de sumas, pero es claramente extensible a otros operadores.

Claramente, este operador es aplicable a cualquier suma, tengan los operandos el tipo que tengan. Afortunadamente, si hay otros operadores definidos sin template, el compilador escogerá esas versiones.

La parte interesante consiste en que si hacemos algo como:

\begin{aligned} a+b+c+d+e \end{aligned},

El compilador coge (d+e) y crea una nueva clase ExprAdd<T,T> donde T es el tipo de las variables. (c+d+e) tendrá tipo: ExprAdd<T, ExprAdd<T,T> >...y así sucesivamente. La expresión final tendrá tipo:

ExprAdd<T, ExprAdd<T, ExprAdd<T, ExprAdd<T, ExprAdd<T,T>>>>>



Afortunadamente para nosotros, la instanciación implícita hace que no sea necesario que conozcamos el tipo de esa clase.

Lo que hemos hecho hasta ahora es construir clases, pero hasta ahora no se ha hecho ninguna operación. Para eso tenemos que ver cómo es la clase ExprAdd. Una versión sencilla sería algo como:

template<class L, class R>
struct ExprAdd
{
  const L &Lop;
  const R &Rop;

  Expr(L &ALop, R &ARop) : Lop(ALop), Rop(ARop) {};

  double Eval(int idx) {return Lop.Eval(idx)+Rop.Eval(idx); };
};

En pocas palabras, es una clase que guarda referencias a los operandos Lop y Rop y, lo más importante, tiene el mismo método Eval(index) que tenían los vectores. Este método llama a los Eval() de los operandos y suma el resultado.

El caso es que los operandos pueden ser vectores (con el Eval definido en la clase vector) o clases del tipo ExprAdd. La función Eval final de la expresión llamará a las funciones Eval de subexpresiones que a su vez llamarán a otras. Todas estas llamadas son conocidas en tiempo de compilación. Aquí viene otra de las ideas importantes, al ser llamadas conocidas en tiempo de compilación, el compilador las inlineará, de forma que, al final, la función Eval de la expresión contendrá un código similar a: a+b+c+d+e, pero donde las variables son ahora reales.

El detalle final es que los cálculos se realizan en el operador igual. Cuando tenemos una expresión del tipo: z = a+b+c+d+e, donde todas las variables son vectores del mismo tamaño,  en a+b+c+d+e se construye una clase que incluye un método Eval y en el operador = se lo llama y se realizan las operaciones, elemento a elemento.

Es decir, no se construyen vectores temporales y el compilador es capaz de hacer todo tipo de optimizaciones. Además, nos permite escribir código claro y bonito. Todos contentos.

Lazy evaluation

De hecho es aún mejor. Esta técnica nos permite una expresividad nueva, que con los métodos antiguos no hubiésemos usado. La técnica está basada en la lazy evaluation.

En pocas palabras, lazy evaluation significa que los cálculos no se hagan hasta que no hagan falta. Para ver cómo va, es mejor poner un ejemplo.

Supongamos que tenemos un vector v y que en ciertas partes del código necesitamos usar su cuadrado. Una posibilidad es hacer algo asi:

v2 = v*v;
a = v2+b;
c = v2-e;
etc...


Que requiere más memoria para v2, a lo que hay que sumar los consabidos problemas de rendimiento. Otra posibilidad es construir la clase:

struct Tsqr
{
  const Vector &V;

  Tsqr(Vector &AV) : V(AV) {};

  double Eval(int idx) { return V.Eval(idx)*V.Eval(idx); };
};

Y la función:

Tsqr sqr(Vector &v) { return Tsqr(v); }

La función sólo sirve para construir la clase Tsqr, que es la que realiza el trabajo. Puesto que el único requerimiento de una clase para poder ser usada, es que incluya el método Eval, la clase funcionará como cualquier otro vector. Ahora podemos usarla en expresiones como:

a = sqr(v)+b;
c = sqr(v)-e;
etc...

Donde, ahora no se guarda nada en memoria. Los cálculos del cuadrado de v se realizan de forma lazy, cuando son necesarios.

De hecho, usando la nueva instrucción auto  de c++11, podemos hacer cosas como:


auto X = a+(b*c);

Donde X es ahora la expressión a+(b*c), sin que se realice ningún cálculo en esa línea (auto permite definir variables deduciendo el tipo de la expresión tras el igual, lo cual añade un nuevo uso a la instanciación automática). Luego podemos usar X más adelante en otras expresiones y será entonces cuando se realizen los cálculos, de nuevo, de forma lazy.

Esto es todo. Como no he usado ninguna imagen hasta ahora, pondré una como premio a los que hayáis llegado hasta aquí.

Los 80 también tuvieron cosas buenas


SVG y asteroides (duros)



Hasta la semana pasada, para mí el SVG era un formato vectorial de gráficos que los navegadores de internet deberían aceptar (y que no todos lo hacían), además de ser el formato usado por Inkscape.

Mi opinión sobre él era que era el típico resultado de un diseño por comité: lo justo para contentar a todos pero insuficiente para gustar realmente a nadie.

El caso es que siempre me había llamado la atención el hecho de que Inkscape permite asignar eventos a los elementos del dibujo. Intentando averiguar qué era eso me encontré con que es posible hacer SVG dinámicos; aunque, cómo siempre, la posibilidad no está implementada en todos los navegadores.

Sin embargo me picó la curiosidad. Más que nada porque desconozco casi todo sobre la programación web, y cuando digo casi todo quiero decir que el otro día descubrí lo que era el DOM.

Para probarlo (y en el proceso, "probar" que soy un friki  que vivió en los 70), se me ocurrió hacer una versión del juego Asteroides más o menos fiel al original, con algún toque de color, pero sin convertirse en uno de esos engendros de remake que se ven por ahí.
Cuando digo "engendro" me refiero a algo así
La cosa resultó sorprendentemente fácil. En un par de horas tenía una nave girando y un asteroide. Después de tres días de baja dedicación ya tenía casi todo lo importante funcionando. El movimiento es bastante suave, incluso para tamaños grandes de pantalla. Eso sí, el código dejó de funcionarme en firefox cuando añadí la detección de choques entre disparos y asteroides; pero en Opera y Chrome sí funciona (aunque aún hay problemas con el zoom). Si queréis ver cómo va podéis probarlo aquí.

Y si no lo queréis probar, esta es una imagen
En este ejemplo, el SVG se abre directamente desde la URL, pero también es posible incrustarlo en una página HTML normal. Pienso que gran parte de las cosas que se hacen en Flash (salvo vídeo, que quizás sea, ahora, la más importante) pueden hacerse mejor en SVG.

Como bonus, aquí tenéis lo que, creo, es el diseño original de los asteroides que encontré en esta página. El que hayan puesto la raya en la O es lo que me convence de que es un documento original de la época, aunque ahora que lo pienso, la raya es precisamente para distinguir ceros de oes...


domingo, 9 de septiembre de 2012

El canal de Garonne

Viendo el puente siguiente, se me dió por pensar brevemente la siguiente pregunta:
¿Como calcular el peso máximo del barco/s que pueden pasar por el puente?.
Os lo dejo para pensar (poco pero algo).

lunes, 3 de septiembre de 2012

En el principio fue el código máquina (I)

Por petición popular (de una sóla persona) inicio aquí una serie de entradas sobre diferentes aspectos de la optimización de programas: ¿Qué cosas funcionan y cuáles no? ¿cuándo es necesario optimizar a mano y cuándo es mejor dejar que el compilador lo haga por tí? etc. Me da que no es este un tema popular, entre otras cosas porque la historia nos dice que, debido a los cambios en el hardware, lo que era una optimización ayer, ahora resulta que es más lento que un código menos trabajado. La sensación es que el esfuerzo es fútil. De hecho ese es uno de los temas que pienso tocar en el futuro: ¿Por qué el paralelismo y las jerarquías de memoria (registros->caches->RAM->Disco) no son meros caprichos de la tecnología actual, sino que hay leyes físicas fundamentales que garantizan que vamos a tenerlas en el futuro? El tema de hoy comenzó, precisamente, como una pregunta referente a la memoria: ¿cuándo es preferible guardar unos cálculos en memoria para no tener que rehacerlos? La respuesta es muy compleja y totalmente caso-dependiente. Existen reglas generales, pero cada caso concreto tiene sus diferentes problemas. Por lo tanto, vamos a considerar el caso concreto, relacionado con un programa de ajedrez:

const int UP = 8;

inline int Sign(int color)
{ return (1 - (color * 2)); }  // White sign is +1, Black sign is -1

inline int Ahead(int color)
{ return (Sign(color) * UP); }  // UP is -8, DOWN is +8


Supongamos que la función Ahead se llama muchas veces siendo color un valor que sólo admite los valores 0 y 1. ¿Es más rápido precalcular un array ahead de dos valores?

En este caso creo que la respuesta es clara. Los cálculos son extremadamente sencillos y pueden hacerse sobre registros, así que ni hablar de precálculos. Aunque para dos únicos valores es seguro que la cache se utilizaría eficientemente, acceder a registros es más rápido que acceder a cache.

Sin embargo aún quedan algunas preguntas por responder:
  • ¿Los inline hacen algo?
  • ¿Se conseguiría alguna ganancia optimizando a mano Ahead y sign?
Para responder a esas preguntas lo mejor es recurrir al código máquina.

La solución está aquí


Para ver el código que genera el compilador escribí el siguiente programa:

int main()
{
    int b=rand() & 1;  // Me aseguro de que b=0,1

    int a=Ahead(b) ;

    printf("Hello world %i !\n", a);
    return 0;
}


El rand() es necesario porque si ponemos una constante literal el compilador no se molesta en incluir las funciones y llama a printf con el valor calculado en tiempo de compilación. El printf es necesario porque si no usamos el valor de a, el compilador no se molesta en incluir nada.

El código máquina generado por gcc con la máxima optimización (-O3), pero sin más florituras es:

sub    rsp,0x8              
call   0x400470 <rand@plt>  // Llama a rand()

and    eax,0x1              // b = rand() & 1  

    
mov    esi,0x4006bc         // N.P.I.

mov    edi,0x1              // N.P.I.

neg    eax                  // eax = -b  
add    eax,eax              // eax = eax+eax = -2*b

lea    edx,[rax*8+0x8]      // edx = 8*eax+8 = 8 - 16*b

 
xor    eax,eax              // eax = 0
call   0x400460 <__printf_chk@plt> // Llama a printf

xor    eax,eax
add    rsp,0x8
ret    
                             


La parte del programa donde se hace la "llamada" a Ahead son simplemente tres instrucciones en código máquina que realizan la operación a = 8 - 16*b.

Se pueden hacer muchos comentarios respecto al código y la cantidad de trucos que usa:

1. Se da cuenta de que la llamada combinada a las dos funciones se puede simplificar como:

Ahead(color) = (1 - (color * 2)) * UP = 8 - 16*color

2. Aprovecha que teníamos color en eax y lo suma consigo mismo para obtener el doble de una forma rápida y sin necesidad de parámetros adicionales.

3. Usa la instrucción lea (load effective address) como un truco para obtener un cálculo aritmético. Esta función suele usarse para calcular punteros y moverse por un array (índice + offset), pero aquí se usa para hacer un cálculo. Desconozco si es posible escoger cualquier multiplicador o si hay sólo unos cuantos a elegir (lo que explicaría por qué no hace simplemente [rax*16+0x8]).
 
En resumen, reto a cualquiera a que escriba un código máquina más eficiente que el que genera el compilador para este caso.

Un detalle interesante es que, en este caso, si hubiéramos usado un array, los "cálculos" necesarios para acceder a los elementos son casi igual de complicados que el cálculo en sí.

Respecto a la otra pregunta: pues no, el gcc genera el mismo código. Da igual si hemos puesto la directiva inline o no.

(Continuará)