Autor: Ian Marteens

Categorías
Insights Música

Family Album

Mientras preparo el primer borrador de Nyx, le he dedicado un rato a la versión de Logic Pro para iPad. Normalmente, uso Sonar (el antiguo Cakewalk, que creo que vuelve a llamarse así), pero es mucho más cómodo usar un iPad mientras vas en tren a la oficina.

Esta es una pieza antigua, pero me gusta más como ha quedado esta vez:

Es una pieza sencilla, sin pretensiones, y es fácil de tocar. He mejorado un poco las dinámicas (cuánta fuerza usas con cada nota, y el volumen general de los pasajes) respecto a la subida original.

Mi rutina de grabación

Como sé que hay muchos programadores que tienen la música como afición, os cuento cómo hago estas cosas, ya sea por si a alguien le interesa o si alguien tiene consejos interesantes.

Ahora mismo, me he pasado, como decía, al Logic Pro para iPad. Me he comprado el iPad Pro de once pulgadas con un procesador M4. Elegí el de un terabyte de disco, porque los de menos capacidad tenían menos RAM. De todos modos, cuentan por ahí que la RAM adicional no se nota mucho en la mayoría de los benchmarks.

Logic Pro viene con un piano de estudio muy bueno. Ya había probado los pianos de Native Instruments, pero personalmente me gusta más éste. No me hagáis mucho caso.

Tengo desde hace unos años un Roland HP704. Buena acción de teclado y, lo que me resulta más cómodo ahora mismo, puedo conectarme al MIDI por Bluetooth.

Normalmente grabo MIDI, por si tengo que corregir algún pasaje complicado. Divido las notas en dos pistas, para ampliar la panorámica estéreo. He utilizado el asistente de masterización de Logic Pro, para mejorar los bajos y la banda de más de 10KHz.

Debería usar monitores de referencia profesionales, pero la idea era aprovechar el viaje en tren.

Categorías
Música

The Tiger

Me encanta este tema de Babelle:

Es más movidillo. Además, está en una tonalidad mayor. A mi edad, empieza uno a valorar cada día de sol.

Categorías
Austra

Destrucción creativa

El compilador de Austra sigue mejorando. Y la librería también, que todo hay que decirlo.

Finalmente, he podido añadir generación de valores aleatorios con AVX2, además de la que ya había implementado con AVX512. La dificultad residía en la necesidad de rotar circularmente un entero largo. AVX512 te da la instrucción. Pero es fácil ver que en AVX2 es posible hacer lo mismo con dos desplazamientos lógicos y una conjunción lógica.

El otro avance importante es la vectorización del algoritmo de Box-Muller para generar distribuciones normales. Ya teníamos un logaritmo natural para vectores de 256 y de 512 bits. Ahora añadimos una función que calcula simultáneamente el seno y el coseno para estos dos tipos de vectores. De hecho, hay una variante especial que es la que usa Box-Muller, en la que todos los valores están entre cero y dos veces pi, y es un poco más rápida. Añadí esta función, la utilicé en el código del generador normal y, finalmente, se ha usado en las clases de vectores, matrices y secuencias.

Operaciones in situ

El compilador intenta ahora optimizar expresiones como la siguiente:

vector1 + vector2 + vector3

¿Qué hay aquí para optimizar? A primera vista, parecería que esa expresión es candidata para reescribirse mediante un combinador lineal, que de hecho existe y se usa para algunas optimizaciones ya:

vec([1, 1, 1, vector1, vector2, vector3)

Se podría utilizar, pero este constructor añadiría tres multiplicaciones escalares innecesarias. Lo que sí nos ahorraríamos sería la creación de dos búferes temporales como resultado de las dos sumas de la expresión original. Lo que ha aprendido a hacer el compilador es a modificar la expresión original a este otro patrón, que es más general:

(vector1 + vector2).InplaceAdd(vector3)

InplaceAdd no se puede llamar directamente desde el lenguaje: es una instrucción peligrosa, porque sobrescribe el búfer del primer operando. El compilador tiene que detectar que el primer operando es de tipo temporal, y no va a utilizarse en el resto de la fórmula, o en el peor de los casos, en el resto de la sesión. Pero en el ejemplo mostrado, la primera suma es evidentemente una construcción temporal, independientemente del sitio de donde estemos sacando las tres variables usadas.

Naturalmente, podemos hacerlo mejor si ya sabemos qué son vector1 y sus dos amigos. Si son variables de sesión, no hay nada más que podamos hacer. Lo mismo ocurre si son parámetros de una función definida por el usuario o variables locales, introducidas por una cláusula let. Pero supongamos que la expresión original fuese la siguiente:

vec::random(10) + vector2 + vector3

En ese caso, el primer operando se crea en un constructor y sólo se usa en esa expresión. El compilador puede entonces crear el siguiente código equivalente:

vec::random(10).InplaceAdd(vector2).InplaceAdd(vector3)

No importa de dónde salgan vector2 y vector3. El búfer del primer operando puede ser reutilizado, y eso es lo que hacemos. El resultado de InplaceAdd, dicho sea de paso, también se puede sobrescribir con total seguridad.

Este tipo de optimizaciones es mejor que las haga un compilador, y no la librería. Es el compilador quien tiene toda la información sobre el uso de un operando. De momento, nos movemos sobre terreno seguro. Hay algunas operaciones que podrían optimizarse con este truco, pero de momento no lo hacemos. Por ejemplo, si una variable local sólo se utiliza una vez, no hay memoria compartida que tengamos que respetar. De momento, no contamos el número de usos de cada variable local.

Finalmente, hemos ampliado las optimizaciones que hacíamos sobre vectores reales a vectores complejos. Es un poco más de código, pero no se hace más lenta la compilación.

Se supone que la señorita de la imagen es Freya bailando sobre calaveras. Eran Kali y Shiva quienes practicaban este tipo de danza, pero estaba seguro de que la tía que la IA iba a generar para Freya iba a ser más guapa. ¿Y las calaveras? Pues están debajo. Pero los generadores de Inteligencia Artificial siguen teniendo problemas para generar manos humanas. De hecho, la imagen que utilizo para Austra la llamo en mi cabeza «la chica del eccema». Llevo meses haciendo retoques a sus manos, y todavía tienen problemas. Pero me gusta el tono del pelo, el óvalo facial y la postura. No se puede tener todo.

Categorías
C#

Xoshiro256**

God and Einstein, playing dice at Heaven

Aunque parezca que el título de la entrada lo tecleó mi gato, es en realidad el nombre de un algoritmo de generación de números aleatorio, inventado por Sebastiano Vigna y David Blackman (enlaces al final).

El nombre es una combinación de las operaciones principales del algoritmo: xor, shift, rotate. En realidad, hay toda una familia de algoritmos similares, con nombres que se parecen mucho. El xoshiro256** es, simplemente, el que ha adoptado .NET Core desde la versión 6 (implementación aquí).

Xoshiro256** es muy rápido, es robusto (aunque no te recomiendan que lo uses en aplicaciones de criptografía), y utiliza muy poca memoria. Como se puede ver en la implementación de .NET 8, sólo necesita cuatro variables de tipo ulong. Es decir, 32 bytes por cada instancia del generador.

Xoshiro vectorial

Austra utiliza generadores de números aleatorios a diestra y siniestra. ¿Es fácil escribir una versión de este algoritmo usando instrucciones AVX? La respuesta es: no, a no ser que uses directamente AVX512F. El problema es que este algoritmo realiza una rotación. Es curioso que los lenguajes modernos de programación no te den directamente esta operación como parte de los operadores de bit. C# tiene un >> y un <<, por ejemplo, pero las rotaciones tienes que buscarlas en la clase BitOperations. De todos modos, el problema es que AVX/AVX2 no tiene una operación SIMD de rotación. Podríamos simularla con shifts y máscaras, pero perderíamos parte de las ventajas de la vectorización. Hay implementaciones en GitHub de generadores aleatorios con AVX2, pero no dan una ganancia de velocidad destacable.

Hay otro problema con el que hay que tener sumo cuidado. Pongamos como ejemplo el constructor de la clase DVector que genera un vector con valores aleatorios:

/// 
/// Creates a vector filled with a uniform distribution generator.
/// 
/// Size of the vector.
/// A random number generator.
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public DVector(int size, Random rnd);

Omito la implementación para simplificar. Lo que quiero subrayar es que este constructor recibe una instancia explícita de la clase Random. ¿Por qué? Pues porque al cliente de la librería puede necesitar resultados repetibles. Es decir, puede que para que un vector «aleatorio» tenga siempre los mismos resultados, podría ser que la instancia de Random se inicializase siempre con la misma semilla. Si no nos interesa usar una semilla, casi siempre utilizaremos Random.Shared, que es una propiedad estática de Random, que es única para cada hilo y que se crea por demanda.

Podríamos complicarnos la vida al escribir la versión vectorizada de este constructor, y averiguar cómo recuperar la semilla del generador que nos pasa el cliente para crear un generador vectorial con esa misma semilla. Pero:

  1. No tengo claro que recuperar esa semilla sea tarea fácil
  2. El generador aleatorio vectorial que he programado no tiene, de momento, la posibilidad de inicializarse con una semilla (y no me parece tampoco sencillo).

Por lo tanto, las rutinas que se han acelerado con el generador vectorial comprueban si el generador que han recibido del cliente es Random.Shared. Eso lo interpretamos como que al cliente no le interesa la repetibilidad de los resultados. Naturalmente, hay que verificar antes si AVX512F está disponible en el ordenador.

Cuando usamos el generador vectorial, los tiempos de ejecución se reducen a una cuarta o quinta parte. No se reducen a la octava parte mecánicamente porque en la mayoría de estas rutinas pesa más la asignación de memoria y el posterior llenado de la misma. Pero bajar el tiempo de ejecución a la cuarta parte me parece meritorio.

Enlaces

Estos son los enlaces a la clase de Austra.Library que genera ocho números de golpe por llamada, y a la página de Sebastiano Vigna, el autor del algoritmo que hemos vectorizado. Mi implementación, con toda seguridad, está abierta a mejoras. Por ejemplo, no me convence la forma en que tengo que convertir un valor ulong en otro de tipo double, porque no es vectorial. Debe existir algo más eficiente, pero de momento no lo he encontrado. De todas maneras, ahí tiene mi implementación, por si le es útil en alguno de sus proyectos.

Queda también pendiente la vectorización de números aleatorios provenientes de una distribución normal. Austra utiliza el método de Box-Muller, pero este método exige el cálculo de un logaritmo y de senos y cosenos. No es tarea imposible, pero tengo primero que proporcionar una rutina que genere vectorialmente senos y cosenos y, quizás antes, decidir si me interesa realmente el método de Box-Muller, o si merece la pena ir directamente a una implementación vectorial del algoritmo del zigurat. Todo, en su debido momento.

Categorías
Insights

La mentalidad predominante

Intentaré ser breve.

Hace ya más de diez años, a una persona muy inteligente y a mí nos pusieron a trabajar a destajo en dos proyectos paralelos en Java cuyo núcleo técnico era bastante parecido: había que enviar precios en una red local a toda leche. Por entonces, Java iba sólo por la versión 7, y había una estructura de datos en Java que se había puesto de moda: el «patrón Disruptor». No pongo enlaces porque no se lo merece. La idea del famoso disruptor era sustituir la típica cola bloqueante (o el Channel que ahora tiene .NET Core) por una cola circular, sin bloqueos… y con un hilo dedicado, girando como un demonio constantemente, calentando la CPU y robando recursos. Teóricamente, este desperdicio se traduciría en menos latencia entre procesos… lo cual incluso puede ser cierto.

A esta idea del bucle en perpetuo movimiento se le añadía un poco de polvo de hadas: la cola circular estaba pensada para evitar pedir memoria dentro de lo posible (un problema que provoca el propio Java), y estaba todo calculado para evitar cosas con nombres tan feos como el false sharing (que es un problema real cuando hay concurrencia). Es decir, por estar hecho en Java, el disruptor era ya la leche… aunque terminó teniendo versiones en C# y supongo que hasta en Python. Toda la mierda que los gatos tiran al suelo termina cayendo en una librería de Python y ahí se hace eterna.

La cruda realidad

Como soy incrédulo, sobre todo de las cosas que no me gustan, lo que hicimos (la persona inteligente y yo) fue crear una demo. Probamos el disruptor contra la cola bloqueante de serie de Java, y ganaba el disruptor, aunque no por KO. El ganador absoluto, sin embargo, fue una variante de la cola bloqueante, con el añadido de recuperar e insertar los elementos en lotes. Es una idea estúpidamente simple y efectiva. Si miras la cola y ves que tienes cuatro elementos pendientes de recuperar, aprovecha y tráetelos todos. En Java tienes que ser cuidadoso con estas cosas, y no crear un array cada vez que te traigas n objetos, pero esto tiene una solución tan sencilla como adosar un array a la cola bloqueante en lotes, y reaprovecharlo todo el tiempo.

¿Cuál era el problema del disruptor? Pues que la idea de quemar un hilo empieza pareciendo atractiva, hasta que te das cuenta de que, si esta filosofía la sigues aplicando al resto de la aplicación, te quedas sin CPU en un plisplás. Creo que la única aplicación del disruptor que llegó a usarse en mi empresa de entonces fue en un proceso de logging sin bloqueos. Me dio un ataque de risa cuando me lo contaron, y todavía creo que quien me lo contó bromeaba.

¡Más Java!

Hay un meme rondando Internet y sacado de una peli barata, sobre un productor musical que todo lo resolvía añadiendo more cowbell, es decir, más cencerro, a la pista de percusión de cualquier canción. En el mundo de la programación, el equivalente es añadir más Java a Java. Me explico:

La otra gran decisión a tomar en el par de proyectos paralelos de hace más de diez años era qué íbamos a utilizar para leer y escribir en sockets. Como uno es un idiota iletrado que conoce poco Java, mi razonamiento es que, a no ser que montase una librería de muy bajo nivel en código nativo, cualquier librería basada en sockets que no implementase algún protocolo experimental no iba a ser mejor que la que ya venía en Java. Todas las librerías de sockets de terceros, y había para escoger, lo único que aportaban eran «abstracciones» como superestructura, que presuntamente podían simplificar la programación, pero jamás iban a conseguir que todo fuese más rápido. Por eso, yo me quedé con los sockets de toda la vida, y mi amigo, que es un experto en Java, se decidió por una librería de la famosa fundación Apache.

Y llegó el día del estreno. Mi aplicación llegaba a leer hasta 250.000 mensajes por cada hilo o canal que habilitase en el proceso. La otra aplicación reventó por un límite absurdo que imponía la librería de terceros. El problema se solucionó, por supuesto, tras algo de trabajo adicional. Pero yo, que soy un idiota, me quedé con una moraleja: añadir más Java a Java no suele aportar nada de valor.

Tiempos modernos

Como sospechará el lector, nadie saca una historia de hace diez años si no ha pasado algo parecido recientemente. En efecto. Una librería escrita en C# resuelve un problema de negocios muy interesante. Pero es una librería a secas, y el equipo que la mantiene no da abasto para conectarla a todas las diferentes fuentes de datos posibles. El tiempo de respuesta, además, es sumamente importante, en este caso.

Entonces apareció un Fervoroso Creyente de la Religión Verdadera, que es Java (y Martin Fowler es uno de sus profetas). La solución propuesta fue añadir Java, y polvos de unicornio. Y, como Java iba a estar obligatoriamente en otro proceso, traer un equivalente moderno del disruptor: una librería de comunicaciones que «resuelve» el problema de la «latencia» poniendo un proceso intermedio que gira como un púlsar al que se le ha ido la olla, y un par de hilos en cada extremo haciendo lo mismo. Es decir, la idea del hilo eterno del disruptor multiplicada por tres. Lo probaron en una máquina razonable, con sólo esos tres componentes: un publicador, el intermediario, y un subscriptor. El invento añadió sólo una decena de microsegundos, o eso me contaron, a la transmisión. Eso ya me habría hecho dudar, si lo hubiese sabido a tiempo. El problema es que la aplicación necesita transmitir datos de unos ocho canales. Tres por ocho, aquí y en Javalandia, siempre ha sido igual a veinticuatro. Estoy escribiendo en la máquina a la que llamo El Pepino, y que es un Core i9 con 8 núcleos físicos y 16 núcleos virtuales por hyperthreading. Es decir, menor estrictamente que veinticuatro.

Creo que todavía están apretado y soltando tuercas para intentar que la cosa funcione. Best of the luck, my friends.

Prometí ser breve, pero no lo conseguí.

Categorías
Austra

Funciones definidas por el usuario

Vamos al grano, o, como diría Haskell B. Curry, «let’s cut to the chase». Esto ya se puede hacer en AUSTRA (en cuanto libere la próxima versión):

let mcd(a, b: int): int =
    let m = a % b in iff(m = 0, b, mcd(b, m)) in
        mcd(80, 140)

Esta es una versión recursiva del máximo común divisor, calculado no con restas, sino con el módulo de la división. Observaciones importantes:

  • Estoy declarando la función como una función local del script. Me falta permitir ahora la declaración de funciones como si fuesen definiciones paramétricas. No lo he hecho todavía porque la implementación de este tipo de funciones se realiza mediante lambdas, que se asignan a una variable local. La recursividad es posible porque la variable local está disponible, con toda la gloria de su prototipo, cuando se compila el cuerpo de la función. Cuando se trate de una definición, probablemente use un truco parecido, pero no es tan inmediato (internamente) como cuando defino una función como parte de una cláusula let.
  • Como se trata de una función recursiva, observe que he definido su tipo de retorno explícitamente. Si no, cuando el compilador encuentre la llamada recursiva a mcd va a tener que volverse loco infiriendo cuál es el tipo de retorno. Ese tipo de inferencias es posible, pero ahora mismo el compilador no está preparado para ello. En mi defensa, recuerde que incluso el gran F# necesita el modificador rec para declarar funciones recursivas. Y F# sí es un lenguaje funcional con todas las de la ley.
  • En este caso, hay una cláusula let anidada dentro de la definición de función. Austra tenía una regla para «aplanar» siempre estas cláusulas en el nivel superior, pero aquí me interesa violar la regla. Eso, o tengo que calcular dos veces el módulo de los dos parámetros.

Otro ejemplo de función recursiva, que utiliza también una cláusula let anidada, aunque esta vez, de tipo diferente:

let fact(n: int) =
    let f(n, acc: int): int = iff(n <= 1, acc, f(n - 1, n * acc)) in
        f(n, 1);
fact(10);
  • Esta es la archifamosa función factorial, pero en vez de definirla en la forma más simple, la defino con una función auxiliar que permite recursividad "por la cola".
  • Aunque el factorial es normalmente recursivo, esta vez no se llama a sí mismo, y no hay que declarar explícitamente el tipo de retorno.
  • Por el contrario, la función interna f sí se llama a sí misma, y hay que tener cuidado con su prototipo.

Claro está que podemos declarar funciones sin tanta fanfarria. Por ejemplo, el factorial es mejor programarlo así:

let fact(n: int) = [x in 2..n].prod;
fact(10);

Recuerde que empezamos con una secuencia de enteros, disfrazada de constructor de secuencias, por lo que no necesitamos poner los valores de la secuencia en un vector, ocupando memoria.

Por cierto, mire lo que podemos hacer ahora con las secuencias:

let collatz(n: int) =
    iseq::unfold(1000000, n, x => iff(x % 2 = 0, x / 2, 3x + 1))
    .until(x => x = 1);
collatz(137)

Esta función genera la secuencia de la conjetura de Collatz para el número entero suministrado como parámetro. Es un problema interesante, y aparentemente fácil, pero que aún no está resuelto. La novedad es que ahora Austra tiene un método Until y un método While que sirven para este tipo de acrobacias. Nuestra función mcd podría haberse también programado sin recursión, mezclando unfold con uno de estos métodos.

Categorías
Austra

List comprehensions in Austra

Hay quien traduce el término inglés list comprehension literalmente como «comprensión de listas» o, aún peor, «listas de comprensión». Lo interesante es que David Turner, el autor del lenguaje funcional Miranda, llamó inicialmente a estas expresiones Zermelo-Frankel expressions, pero alguien lo convenció para llamarlas list comprehensions, que en inglés no suena tan mal. Yo, porque soy un tipo caprichoso, cuando traduzca el término en el contexto de Austra, las llamaré constructores de secuencias, a secas. Menos pretencioso, y más comprensible, creo.

Un truco para escribir menos

¿Que es un constructor de secuencias? Pues es un truco sencillo para eso mismo: construir secuencias, pero usando menos código, con el añadido de que el resultado es normalmente más sencillo de leer (no siempre). Imagine que tenemos un vector de números reales, y queremos quedarnos con los que son enteros divisibles por dos, para elevarlos al cuadrado. En Austra, hasta ahora, haríamos esto, suponiendo que tenemos el vector ya almacenado en una variable global v:

v.filter(x => x % 2 = 0).map(x => x^2)

Con el mecanismo nuevo de construcción de secuencias, la expresión anterior se reduciría a esto:

[x in v : x % 2 = 0 => x^2]

No hace falta contar caracteres para ver que realmente hemos escrito menos. Nos hemos ahorrado los paréntesis de los métodos, y los parámetros de las funciones lambdas se han reducido a uno solo, reflejando el hecho de que los valores que fluyen por el constructor son casi siempre (casi) del mismo tipo. El resultado del constructor, en este caso, es un vector, y puedo seguir aplicando métodos y operadores tras el corchete de cierre. Por ejemplo, puedo hacer algo algo estúpido como añadir uno a cada valor (lo podía haber hecho al elevar al cuadrado):

([x in v : x % 2 = 0 => x^2] + 1).plot

Pero también podía haber transformado el vector resultante con una matriz, o cualquier otra cosa posible con un vector.

También podía haber creado una secuencia en el constructor, aunque los datos viniesen de un vector, o de una matriz:

[x in seq(v) : x % 2 = 0 => x^2];
[x in seq(v1^v2) : x % 2 = 0 => x^2];
[x in iseq(1, 100) : x % 2 = 0 => x^2];
[x in 1..100 : x % 2 = 0 => x^2];

El primer ejemplo usa una secuencia basada en un vector. El segundo, una secuencia basada en una matriz que se genera a partir de dos vectores. En el tercero, simplemente uso una secuencia de enteros construida a partir de un range. Y en el último ejemplo uso más «syntatic sugar» para generar la secuencia directamente a partir de un rango.

Ojo con las series

Con las series, hay que tener un poco de cuidado, porque el método que filtra una serie recibe como parámetro una lambda del tipo Func<Point<Date> bool>, mientras que el método Map usa una de tipo Func<double, double>. El constructor de secuencias lo tiene en cuenta, pero tenemos que recordarlo:

let mean = MSFT.mean in
    [x in MSFT : x.date >= jan2020 and x.value >= mean];

Observe que el filtro utiliza tanto la fecha como el valor de los puntos de la serie. Además, he omitido la transformación. Podemos omitir el filtro, la transformación (o proyección, en terminología SQL y C#) o incluso ambos.

Cuantificadores lógicos

De todas maneras podemos ir un poco más lejos que Python y Haskell. Vamos a comenzar por algo sencillo. ¿Es el 97 un número primo? Vamos a preguntarlo usando constructores de listas:

[all x in 2..96: 97 % x != 0]

La expresión anterior no devuelve una lista, sino un valor de tipo lógico. Será verdadero si alguno de los números entre 2 y 96 divide al número 97. Podríamos haber usado una cota superior más baja, por supuesto, pero no quiero complicar la explicación con detalles innecesarios.

all y any, en AUSTRA, no son palabras reservadas, pero en este caso se consideran palabras reservadas contextuales. La expresión anterior es equivalente a esta otra:

iseq(2..96).all(x => 97 % x != 0)

Esto, naturalmente, aunque sea ligeramente interesante, es sólo un rodeo hacia nuestro objetivo. ¿Qué tal si quiero todos los números primos del 2 al 1000?

[y in 2..1000 : all x in 2..96: y % x != 0]

La presencia de dos caracteres : nos está indicando que hay dos expresiones entre los corchetes. De hecho, estamos usando una función lambda anidada dentro de otra, y la más interna está «capturando» el parámetro de la más externa:

iseq(2, 1000).filter(y => iseq(2, y - 1).all(x => y % x != 0))

Y si quisiéramos elevar cada primo al cuadrado, añadiríamos una función de proyección al engendro que hemos creado:

[y in 2..1000 : all x in 2..96: y % x != 0 => y^2]

¡Chúpate esa, Python…!

Categorías
Música

Babelle

¡Tengo una sobrina artista! Naturalmente, el tío opina que la chica lo hace muy bien. Pero aquí le dejo el enlace a su última canción, para contrastar opiniones:

Es más guapa, por cierto, de lo que aparenta la foto, pero esa también es una opinión personal.

Categorías
Austra

Unfold

¿Cómo puedo calcular la serie de Fibonacci en Austra, utilizando una secuencia? Con vectores (voy a usar vectores y secuencias reales para evitar problemas de desbordamiento), es relativamente sencillo, si aprovechamos los safe indexers en una expresión lambda:

vec::new(128, (i, v) =>
  if i = 0 then 1 else v{i-1} + v{i-2})

El problema de esta solución es que consume memoria. Tenemos los dos últimos números generados a la izquierda (por así decirlo) del número que estamos generando, pero en realidad, es porque tenemos todo la memoria del vector ya reservada. ¿Y si quisiéramos hacerlo con las nuevas secuencias?

La solución en lenguajes funcionales

Los lenguajes funcionales suelen usar una función unfold para estos casos. Por ejemplo, en F# se define unfold de esta manera:

val unfold : ('State -> ('T * 'State) option) -> 'State -> seq<'T>

La función, en sí, escrita como una secuencia infinita, sería más o menos esto.

let fib = Seq.unfold (fun (lastValue, currentValue)
    -> Some (lastValue,
    (currentValue, lastValue + currentValue))) (1, 1)

Ahora, las explicaciones: necesitaríamos que Austra fuese un lenguaje con tipos genéricos e inferencia automática de tipos… y no es un objetivo mío a corto plazo. Nyx, que es una idea de lenguaje con capacidades numéricas y aceleración que me guardo bajo la manga, tendrá genericidad, pero no estoy del todo convencido que tenga que ser funcional «puro». La función en F#, que sería casi idéntica en Haskell, necesita un «estado» para ir arrastrando en las sucesivas llamadas. En el caso de Fibonacci, el estado son los últimos dos números generados.

De regreso a Austra

No obstante, podemos hacer algo útil para el 99% de los casos. Estos son los tres métodos que se han implementado en C# para soportar unfold en el lenguaje de fórmulas:

public static DSequence Unfold(int size,
  double seed,
  Func<double, double> unfold);
public static DSequence Unfold(int size,
  double seed,
  Func<int, double, double> unfold);
public static DSequence Unfold(int size,
  double first, double second,
  Func<double, double, double> unfold);

La idea es proporcionar sobrecargas para los casos más frecuentes. El primer caso es para cuando el estado es un solo número real. El segundo caso es más sutil: el estado es un número real… más la posición del elemento que se va a generar. Si estuviésemos en un lenguaje funcional, esa posición tendría que estar explícitamente representada en el estado. Y en el tercer caso, el estado son dos números reales, como necesitamos para la secuencia de Fibonacci.

Veamos ahora cómo se usa esto en el lenguaje de fórmulas:

-- Potencias de 2, desde 2 a 1024.
seq::unfold(10, 2, x => 2x);
-- Serie de Maclaurin para exp(1).
1 + seq::unfold(100000, 1, (n, x) => x / (n + 1)).sum;

He tenido que sumar explícitamente un uno a la serie de Maclaurin, porque era lo más sencillo. Si quiere practicar, piense en cómo calcular ln(2) con una secuencia de estas.

La secuencia de Fibonacci se calcula así:

seq::unfold(50, 1, 1, (x, y) => x + y);

Con números enteros, haríamos esto:

iseq::unfold(30, 1, 1, (x, y) => x + y);

No es la solución perfecta, pero es fácil de entender, y encima se ejecuta con rapidez.

Algo curioso relacionado con el diseño de lenguajes: imagine un lenguaje que ofrece un vec<double> y un vec<complex>, en vez del vec y el cvec del lenguaje de fórmulas de Austra. ¿Se ha dado cuenta de que el vector de complejos no podría soportar todas las operaciones del vector de números reales? El caso más evidente: no se puede ordenar el vector de complejos, ni calcular la entrada con el menor o mayor valor. Los números complejos no soportan un orden total que tenga sentido. ¿Cómo haría usted para que una clase genérica tuviese en cuenta estos detalles? Ahora mismo, si quieres hacer esto en C#, tienes que definir una clase genérica «recortada», en plan Vec<T>, y las clases finales serán clases que ya no serán genéricas, pero que añadirían la funcionalidad no común. Algo así es que lo hace Math.NET.

Ese es el problema que intento resolver: un lenguaje de programación con auto-vectorización, en el que programar una librería como la de Austra no sea tan complicado, y que me deje definir tipos matemáticos como el vector de complejos y el vector de reales utilizando una base común genérica. Hay más cosas en las matemáticas que no se ajustan a la idea de la programación de que «lo derivado tiene más operaciones».

Categorías
C# FinTech

Secuencias

Esto es una perogrullada: en un lenguaje de programación funcional, las cosas se suelen hacer de forma diferente. Y probablemente escribir «perogrullada» sea una pedantería. Dejo el párrafo vivo por pereza, pero no me lo tenga en cuenta.

Factorial

Austra no pretende ser un lenguaje de programación funcional. Ni siquiera pretende ser «Turing-completo». En algún punto del camino, he pensado en diseñar un lenguaje «de verdad», e incluso en saltarme la maquinaria de .NET para compilarlo a código nativo directamente. Pero falta un poco para eso.

De todas maneras, hay que intentar que se puedan hacer todas las cosas posibles en Austra. Por ejemplo, ¿cómo calculas un factorial, si no tienes todavía funciones recursivas, y no quieres introducir bucles? Hasta ahora, la solución era parecida a ésta:

vec(10, i => i + 1).prod

No estamos definiendo una función, sino que estamos usando un parámetro a dedo, pero funciona. Creamos un vector de 10 entradas, lo llenamos con los números del 1 al 10, y multiplicamos todos los elementos del vector. Incluso podemos «optimizar» un poco la expresión:

vec(8, i => i + 2).prod * 2

Es absurdo dejar el 1 en el vector, y ya que tenemos diez elementos, en vez de dejarlos nueve, quito también el 2. Eso lo hago porque sé que internamente el producto mete los ocho elementos en dos registros AVX y no tengo que manejar el elemento nuevo aparte. El problema es que, de todas maneras, estamos pidiendo memoria para el vector.

En Austra 2.0 (que ya tiene soporte para .NET 8 y AVX512, dicho sea de paso), ya están implementadas las secuencias de valores reales, y puedes hacer esto otro, que es más natural:

seq(2, 10).prod

La nueva clase se llama seq, y es una copia descarada del diseño de enumerables y LINQ (o de los streams de Java, o de las «list comprehensions» de tantos otros lenguajes funcionales). Ahora tenemos un seq, pero luego vendrán iseq y cseq, para secuencias de enteros y complejos. Observe también que he acortado vector a vec, para usar cvec en vez de complexvector, y para que un futuro vector de enteros se pueda llamar ivec a secas.

Las secuencias pueden hacer casi todas las cosas que hacen los vectores, y en la mayoría de los casos, generarlas es más sencillo: con el vector, usamos una función lambda. También permiten ahorrar código porque aplican los mismos trucos que LINQ para objetos. Por ejemplo:

-- Esta es una secuencia que divide el intervalo
-- de 0 a 2*pi en 360 trozos.
-- Multiplicamos cada número por dos...
seq(0, τ, 360) * 2
-- .. pero internamente, lo transformamos en esto:
seq(0, 2τ, 360)

Recursividad eficiente

Y esto último todavía no está implementado, pero es bueno pensar en estas cosas antes de lanzarse a programarlas. ¿Cómo definiríamos una función para el factorial que fuese medianamente eficiente, pero usando recursividad? Lo que se suele hacer en un lenguaje funcional, pero con una notación al estilo Austra:

def fact(n: int) =
  let f = (n, acc: int) =>
    if n = 1 then acc else f(n - 1, acc * n) in
      f(n, 1)

El truco está en usar una función auxiliar, que al estilo Austra sería una lambda en una cláusula let. La función auxiliar es recursiva por la cola, por lo que el propio JIT de .NET la puede convertir en un bucle. La recursividad en lambdas va a ser problemática, porque estoy usando directamente Linq.Expressions para generar código, pero hay un truco sencillo, que es el que voy a usar en las primeras pruebas: declarar una variable de tipo lambda, asignarle un puntero vacío, y a continuación, asignarle una expresión lambda que ya podrá usar la variable recursivamente. La alternativa es usar un combinador Y. Es un truco que viene, precisamente, del cálculo lambda, sobre el que escribiré en algún momento. No es que sea cómodo o eficiente, pero es interesante.

De todas maneras, me he dado cuenta, revisando la documentación, que en F# hay que anunciar cuando una función va a ser recursiva. Una alternativa que tengo que pensar es generar un «método dinámico» en estos casos, con el esfuerzo adicional de generar directamente el código IL. Pero tengo ya experiencia de generar IL con Freya, y puedo reutilizar código.