Si insisto tanto en estos temas de instrucciones vectoriales, es parar preparar a mis lectores potenciales para que usen estas técnicas cuando lo crean necesario. Naturalmente, las primeras explicaciones van a ser del tipo heroico: todo hay que hacerlo a mano, razonando desde los primeros principios. Pero, como decía Molly Jones, la mujer de Desmond Jones, la la, how their life goes on: la vida sigue su curso, y lo normal en la vida de un programador es usar las técnicas que todos conocemos y apreciamos para ahorrarnos esfuerzo.
Lo primero para que todos nos ahorremos trabajo, por supuesto, consiste en usar estas cosas a través de una librería bien probada. Ya existen unas cuantas de estas, pero estoy creando Austra porque hay un nicho muy concreto, y porque la librería tiene méritos propios. La idea es hacer de Austra un proyecto open source, por supuesto. En este momento está en GitHub, pero no es por ahora un repositorio público porque la aplicación de «pruebas» está basada en WPF y DevExpress. O me compro personalmente una licencia comercial, o cambio los controles por algo gratuito. No es sencillo. Hay cosas como Avalonia o Uno Platform, que además, son multiplataformas, pero son bastante pobres en componentes. Y no quiero ni oír hablar de .NET MAUI: es un fracaso, de momento. Antes que usar .NET MAUI, prefiero «degradar» la aplicación a Windows Forms (tengo un editor de código bastante bueno) y currarme los gráficos y las cosas que falten. Pero a eso vamos: a que Austra esté disponible gratuitamente como open source, que cualquier que quiera pueda colaborar, que haya un conjunto de tests y benchmarks exhaustivo, etc, etc.
Pero no se escribe un post para explicar lo obvio. Mi verdadero objetivo ahora es explicar un par de técnicas muy tontas que bajan el listón del heroísmo al escribir este tipo de código, y que todos conocemos, pero que nos puede dar miedo usar a primeras, porque no sabes cómo va a interferir en la generación de código de .NET 7 y .NET 8 (que tiene sus cagadas, como todo).
Primer ejemplo: resulta que mucho de estos algoritmos algebraicos calculan el producto escalar de dos zonas de memoria en muchos casos. Hay una clase Vector con su correspondiente producto escalar, e incluso un ComplexVector, pero es mejor tener una rutina que maneje directamente zonas de memoria con un tamaño predeterminado. Austra tiene una clase estática CommonMatrix que precisamente implementa estas cosas. Ésta es la última versión del método en cuestión:
/// <summary>
/// Computes the dot product of two double arrays.
/// </summary>
/// <param name="p">Pointer to the first array.</param>
/// <param name="q">Pointer to the second array.</param>
/// <param name="size">Number of items in each array.</param>
/// <returns>A sum of products.</returns>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public unsafe static double DotProduct(
double* p, double* q, int size)
{
double sum;
int i = 0;
if (Avx.IsSupported)
{
Vector256<double> acc = Vector256<double>.Zero;
for (int top = size & AVX_MASK; i < top; i += 4)
acc = acc.MultiplyAdd(p + i, q + i);
sum = acc.Sum();
}
else
sum = 0;
for (; i < size; i++)
sum += p[i] * q[i];
return sum;
}
A primera vista, lo único raro aquí es el atributo que fuerza un inlining agresivo. Pero hay un par de cosillas más que no son métodos habituales de AVX. Los muestro aparte aquí, porque son métodos de extensión de la misma clase CommonMatrix:
/// <summary>Sums all the elements in a vector.</summary>
/// <param name="v">
/// A intrinsics vector with four doubles.
/// </param>
/// <returns>The total value.</returns>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal static double Sum(this Vector256<double> v)
{
v = Avx.HorizontalAdd(v, v);
return v.ToScalar() + v.GetElement(2);
}
Resulta que sumar los cuatro elementos de un vector de cuatro dobles no es moco de pavo. Hay teorías que pululan por la Internet sobre cuál es la forma eficiente de lograrlo. Yo no estoy seguro aún de que la mía sea la más eficiente, pero gracias al encapsulamiento en un método separado, si descubro algo mejor, tendré que tocar el código en un único lugar. Obvio, ¿no? Pero hasta comprobar que el JIT de .NET Core no se volvía loco, no me atreví a dar este paso. Ya he comprobado que no pasan cosas raras.
Es curioso, sin embargo, que el mínimo y el máximo de un vector se calcule más eficientemente con un método como el siguiente:
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal static double Max(this Vector256<double> v)
{
var x = Sse2.Max(v.GetLower(), v.GetUpper());
return Math.Max(x.ToScalar(), x.GetElement(1));
}
En este caso, no parece que la transición temporal a una instrucción SSE haga mucho daño al código que la rodea. De todas maneras, observe que para la reducción final al escalar hay que ir a por el segundo elemento directamente. Cosas de Intel.
Éste es otro método que usa el código original, omitiendo esta vez los comentarios XML:
internal unsafe static Vector256<double> MultiplyAdd(
this Vector256<double> summand,
double* multiplicand,
double* multiplier) =>
Fma.IsSupported
? Fma.MultiplyAdd(
Avx.LoadVector256(multiplicand),
Avx.LoadVector256(multiplier),
summand)
: Avx.Add(summand, Avx.Multiply(
Avx.LoadVector256(multiplicand),
Avx.LoadVector256(multiplier)));
Hay un par de variantes más que usan vectores directamente en vez de direcciones, y variantes adicionales para MultiplySubtract y MultiplyAddNegated, que sons métodos FMA. Estas aparentes tonterías ahorran líneas y líneas de código. En mi caso, como tengo memoria fotográfica, prefiero que el código fuente sea lo más pequeño posible para poder recordarlo más adelante.
CNOT significa «conditional NOT», una negación condicional, y el parentesco con la negación a secas se evidencia en la simbología asociada: la parte inferior del símbolo es la cruz dentro del círculo que se suele usar para la negación. Se trata de una compuerta que trabaja siempre con dos qubits. Todas las transformaciones que hemos visto hasta ahora sobre dos qubits se creaban componiendo transformaciones sobre un qubit. Con esta compuerta, ya saltamos a transformaciones más complejas.
En notación algebraica, este circuito se puede escribir así:
Efectivamente, nos advierte que hay dos posibles lecturas, las antes mencionadas, y que cada una se obtiene en un 50% de los casos.
Las compuertas $H$, naturalmente, corresponden a transformadas de Hadamard aplicadas a qubits individuales.
El objetivo de esta entrada no es ponernos ya a diseñar circuitos: cuando llegue el momento, es preferible mostrar el código en Python. Ahora nos interesa la parte matemática, sobre cómo definir transformaciones que afecten a más de un qubit. Utilizo Qiskit para no tener que dibujar los circuitos a mano.
Esta vez queremos invertir solamente el qubit menos significativo. Aparentemente, además, sólo tenemos una matriz, pero eso no es del todo cierto. Lo que queremos calcular ahora es la matriz combinada $I\otimes X$, donde $I$ es la matriz identidad.