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.