assembly - Optimice para una multiplicación rápida pero una adición lenta: FMA y doble doble

floating-point x86

Question

Cuando obtuve un procesador Haswell por primera vez, intenté implementar FMA para determinar el conjunto de Mandelbrot. El algoritmo principal es este:

intn = 0;
for(int32_t i=0; i<maxiter; i++) {
floatn x2 = square(x), y2 = square(y); //square(x) = x*x
floatn r2 = x2 + y2;
booln mask = r2<cut; //booln is in the float domain non integer domain
if(!horizontal_or(mask)) break; //_mm256_testz_pd(mask)
n -= mask
floatn t = x*y; mul2(t); //mul2(t): t*=2
x = x2 - y2 + cx;
y = t + cy;
}

Esto determina si hay n píxeles en el conjunto de Mandelbrot. Entonces, para doble punto flotante se ejecuta sobre 4 píxeles ( floatn = __m256d , intn = __m256i ). Esto requiere 4 multiplicaciones de coma flotante SIMD y cuatro adiciones de coma flotante SIMD.

Luego modifiqué esto para trabajar con FMA como este

intn n = 0;
for(int32_t i=0; i<maxiter; i++) {
floatn r2 = mul_add(x,x,y*y);
booln mask = r2<cut;
if(!horizontal_or(mask)) break;
add_mask(n,mask);
floatn t = x*y;
x = mul_sub(x,x, mul_sub(y,y,cx));
y = mul_add(2.0f,t,cy);
}

donde mul_add llama a _mm256_fmad_pd y mul_sub llama a _mm256_fmsub_pd . Este método utiliza 4 operaciones SIMD de FMA y dos multiplicaciones SIMD, que son dos operaciones aritméticas menos que sin FMA. Además, FMA y la multiplicación pueden usar dos puertos y además solo uno.

Para hacer mis pruebas menos sesgadas, amplié una región que está completamente en el conjunto de Mandelbrot, por lo que todos los valores son maxiter . En este caso, el método que usa FMA es aproximadamente un 27% más rápido. Eso es ciertamente una mejora, pero pasar de SSE a AVX duplicó mi rendimiento, así que esperaba quizás otro factor de dos con FMA.

Pero luego encontré this respuesta con respecto a FMA donde dice

El aspecto importante de la instrucción de fusión múltiple multiplicada es la precisión (prácticamente) infinita del resultado intermedio. Esto ayuda con el rendimiento, pero no tanto porque dos operaciones están codificadas en una sola instrucción. Ayuda con el rendimiento porque la precisión prácticamente infinita del resultado intermedio es a veces importante y muy costosa de recuperar con la multiplicación y la suma ordinarias cuando este nivel de la precisión es realmente lo que busca el programador.

y luego da un ejemplo de multiplicación doble * doble a double-double

high = a * b; /* double-precision approximation of the real product */
low = fma(a, b, -high); /* remainder of the real product */

A partir de esto, llegué a la conclusión de que estaba implementando FMA de manera no óptima y decidí implementar SIMD doble-doble. Implementé doble doble basado en los números de coma flotante de precisión extendida en papel para computación GPU . El papel es para doble flotación, así que lo modifiqué para doble doble. Además, en lugar de empaquetar un valor doble doble en un registro SIMD, empaque 4 valores doble doble en un registro alto AVX y un registro bajo AVX.

Para el conjunto de Mandelbrot, lo que realmente necesito es multiplicación y suma doble-doble. En ese documento, estas son las funciones df64_add y df64_mult . La imagen a continuación muestra el ensamblaje de mi función df64_mult para el software FMA (izquierda) y el hardware FMA (derecha). Esto muestra claramente que el hardware FMA es una gran mejora para la multiplicación doble-doble.

Entonces, ¿cómo funciona el hardware FMA en el cálculo doble-doble de Mandelbrot? La respuesta es que solo es un 15% más rápido que con el software FMA. Eso es mucho menos de lo que esperaba. El cálculo de Mandelbrot doble-doble necesita 4 adiciones doble-doble y cuatro multiplicaciones doble-doble ( x*x , y*y , x*y , y 2*(x*y) ). Sin embargo, la multiplicación 2*(x*y) es trivial para doble-doble, por lo que esta multiplicación puede ignorarse en el costo. Por lo tanto, la razón por la que creo que la mejora con el uso de hardware FMA es tan pequeña es que el cálculo está dominado por la lenta adición doble-doble (ver el ensamblaje a continuación).

Solía ​​ser que la multiplicación era más lenta que la suma (y los programadores usaron varios trucos para evitar la multiplicación), pero con Haswell parece que es al revés. No solo debido a FMA sino también porque la multiplicación puede usar dos puertos, pero además solo uno.

Entonces mis preguntas (finalmente) son:

1. ¿Cómo se optimiza cuando la suma es lenta en comparación con la multiplicación?
2. ¿Hay alguna forma algebraica de cambiar mi algoritmo para usar más multiplicaciones y menos adiciones? Sé que hay un método para hacer lo contrario, por ejemplo (x+y)*(x+y) - (x*x+y*y) = 2*x*y que usan dos adiciones más para una multiplicación menos.
3. ¿Hay alguna manera de simplemente la función df64_add (por ejemplo, usando FMA)?

En caso de que alguien se pregunte, el método doble-doble es aproximadamente diez veces más lento que el doble. Eso no es tan malo, creo que si hubiera un tipo de hardware de precisión cuádruple, probablemente sería al menos dos veces más lento que el doble, por lo que mi método de software es aproximadamente cinco veces más lento de lo que esperaría para el hardware si existiera.

ensamblaje df64_add

vmovapd 8(%rsp), %ymm0
movq %rdi, %rax
vmovapd 72(%rsp), %ymm1
vmovapd 40(%rsp), %ymm3
vaddpd %ymm1, %ymm0, %ymm4
vmovapd 104(%rsp), %ymm5
vsubpd %ymm0, %ymm4, %ymm2
vsubpd %ymm2, %ymm1, %ymm1
vsubpd %ymm2, %ymm4, %ymm2
vsubpd %ymm2, %ymm0, %ymm0
vaddpd %ymm1, %ymm0, %ymm2
vaddpd %ymm5, %ymm3, %ymm1
vsubpd %ymm3, %ymm1, %ymm6
vsubpd %ymm6, %ymm5, %ymm5
vsubpd %ymm6, %ymm1, %ymm6
vaddpd %ymm1, %ymm2, %ymm1
vsubpd %ymm6, %ymm3, %ymm3
vaddpd %ymm1, %ymm4, %ymm2
vaddpd %ymm5, %ymm3, %ymm3
vsubpd %ymm4, %ymm2, %ymm4
vsubpd %ymm4, %ymm1, %ymm1
vaddpd %ymm3, %ymm1, %ymm0
vaddpd %ymm0, %ymm2, %ymm1
vsubpd %ymm2, %ymm1, %ymm2
vmovapd %ymm1, (%rdi)
vsubpd %ymm2, %ymm0, %ymm0
vmovapd %ymm0, 32(%rdi)
vzeroupper
ret

Answer 1

Mencionas el siguiente código:

vsubpd %ymm0, %ymm4, %ymm2
vsubpd %ymm2, %ymm1, %ymm1 <-- immediate dependency ymm2
vsubpd %ymm2, %ymm4, %ymm2
vsubpd %ymm2, %ymm0, %ymm0 <-- immediate dependency ymm2
vaddpd %ymm1, %ymm0, %ymm2 <-- immediate dependency ymm0
vaddpd %ymm5, %ymm3, %ymm1
vsubpd %ymm3, %ymm1, %ymm6 <-- immediate dependency ymm1
vsubpd %ymm6, %ymm5, %ymm5 <-- immediate dependency ymm6
vsubpd %ymm6, %ymm1, %ymm6 <-- dependency ymm1, ymm6
vaddpd %ymm1, %ymm2, %ymm1
vsubpd %ymm6, %ymm3, %ymm3 <-- dependency ymm6
vaddpd %ymm1, %ymm4, %ymm2
vaddpd %ymm5, %ymm3, %ymm3 <-- dependency ymm3
vsubpd %ymm4, %ymm2, %ymm4
vsubpd %ymm4, %ymm1, %ymm1 <-- immediate dependency ymm4
vaddpd %ymm3, %ymm1, %ymm0 <-- immediate dependency ymm1, ymm3
vaddpd %ymm0, %ymm2, %ymm1 <-- immediate dependency ymm0
vsubpd %ymm2, %ymm1, %ymm2 <-- immediate dependency ymm1

Si verifica cuidadosamente, estas son operaciones en su mayoría dependientes, y no se cumple una regla básica sobre la eficacia de latencia / rendimiento. La mayoría de las instrucciones dependen del resultado del anterior, o de 2 instrucciones anteriores. Esta secuencia contiene una ruta crítica de 30 ciclos (aproximadamente 9 o 10 instrucciones sobre "latencia de 3 ciclos" / "rendimiento de 1 ciclo").

Su IACA informa "CP" => instrucción en la ruta crítica, y el costo evaluado es de 20 ciclos de rendimiento. Debería obtener el informe de latencia porque es el que importa si está interesado en la velocidad de ejecución.

Para eliminar el costo de esta ruta crítica, debe intercalar unas 20 instrucciones más similares si el compilador no puede hacerlo (por ejemplo, porque su código doble-doble está en una biblioteca separada compilada sin optimizaciones -flto y vzeroupper en todas partes en la entrada y salida de la función , vectorizer solo funciona bien con código en línea).

Una posibilidad es ejecutar 2 cálculos en paralelo (ver acerca de la unión de código en una publicación anterior para mejorar la canalización)

Si supongo que su código doble doble se parece a esta implementación "estándar"

// (r,e) = x + y
#define two_sum(x, y, r, e)
do { double t; r = x + y; t = r - x; e = (x - (r - t)) + (y - t); } while (0)
#define two_difference(x, y, r, e) /
do { double t; r = x - y; t = r - x; e = (x - (r - t)) - (y + t); } while (0)
.....

Luego debe considerar el siguiente código, donde las instrucciones se entrelazan con un grano bastante fino.

// (r1, e1) = x1 + y1, (r2, e2) x2 + y2
#define two_sum(x1, y1, x2, y2, r1, e1, r2, e2)
do { double t1, t2 /
r1 = x1 + y1; r2 = x2 + y2; /
t1 = r1 - x1; t2 = r2 - x2; /
e1 = (x1 - (r1 - t1)) + (y1 - t1); e2 = (x2 - (r2 - t2)) + (y2 - t2); /
} while (0)
....

Luego, esto crea un código como el siguiente (aproximadamente la misma ruta crítica en un informe de latencia y aproximadamente 35 instrucciones). Para obtener detalles sobre el tiempo de ejecución, la ejecución fuera de orden debe pasar por alto sin detenerse.

vsubsd %xmm2, %xmm0, %xmm8
vsubsd %xmm3, %xmm1, %xmm1
vaddsd %xmm4, %xmm4, %xmm4
vaddsd %xmm5, %xmm5, %xmm5
vsubsd %xmm0, %xmm8, %xmm9
vsubsd %xmm9, %xmm8, %xmm10
vaddsd %xmm2, %xmm9, %xmm2
vsubsd %xmm10, %xmm0, %xmm0
vsubsd %xmm2, %xmm0, %xmm11
vaddsd %xmm14, %xmm4, %xmm2
vaddsd %xmm11, %xmm1, %xmm12
vsubsd %xmm4, %xmm2, %xmm0
vaddsd %xmm12, %xmm8, %xmm13
vsubsd %xmm0, %xmm2, %xmm11
vsubsd %xmm0, %xmm14, %xmm1
vaddsd %xmm6, %xmm13, %xmm3
vsubsd %xmm8, %xmm13, %xmm8
vsubsd %xmm11, %xmm4, %xmm4
vsubsd %xmm13, %xmm3, %xmm15
vsubsd %xmm8, %xmm12, %xmm12
vaddsd %xmm1, %xmm4, %xmm14
vsubsd %xmm15, %xmm3, %xmm9
vsubsd %xmm15, %xmm6, %xmm6
vaddsd %xmm7, %xmm12, %xmm7
vsubsd %xmm9, %xmm13, %xmm10
vaddsd 16(%rsp), %xmm5, %xmm9
vaddsd %xmm6, %xmm10, %xmm15
vaddsd %xmm14, %xmm9, %xmm10
vaddsd %xmm15, %xmm7, %xmm13
vaddsd %xmm10, %xmm2, %xmm15
vaddsd %xmm13, %xmm3, %xmm6
vsubsd %xmm2, %xmm15, %xmm2
vsubsd %xmm3, %xmm6, %xmm3
vsubsd %xmm2, %xmm10, %xmm11
vsubsd %xmm3, %xmm13, %xmm0

Resumen:

• incorpore su código fuente doble-doble: el compilador y el vectorizador no pueden optimizar las llamadas a funciones debido a restricciones ABI, y a través de accesos a la memoria por temor a alias.

• puntada de código para equilibrar el rendimiento y la latencia y maximizar el uso de los puertos de la CPU (y también maximizar las instrucciones por ciclo), siempre que el compilador no derrame demasiados registros en la memoria.

Puede rastrear los impactos de optimización con la utilidad perf (paquetes linux-tools-generic y linux-cloud-tools-generic) para obtener la cantidad de instrucciones ejecutadas y la cantidad de instrucciones por ciclo.

Answer 2

Para acelerar el algoritmo, utilizo una versión simplificada basada en 2 fma, 1 mul y 2 add. Proceso 8 iteraciones de esa manera. Luego calcule el radio de escape y revierta las últimas 8 iteraciones si es necesario.

El siguiente bucle crítico X = X ^ 2 + C escrito con intrínsecos x86 está muy bien desenrollado por el compilador, y verá después de desenrollar que las 2 operaciones FMA no dependen mucho unas de otras.

// IACA_START;
for (j = 0; j < 8; j++) {
Xrm = _mm256_mul_ps(Xre, Xim);
Xtt = _mm256_fmsub_ps(Xim, Xim, Cre);
Xrm = _mm256_add_ps(Xrm, Xrm);
Xim = _mm256_add_ps(Cim, Xrm);
Xre = _mm256_fmsub_ps(Xre, Xre, Xtt);
} // for
// IACA_END;

Y luego calculo el radio de escape (| X | <umbral), que cuesta otra fma y otra multiplicación, solo cada 8 iteraciones.

cmp = _mm256_mul_ps(Xre, Xre);
cmp = _mm256_fmadd_ps(Xim, Xim, cmp);
cmp = _mm256_cmp_ps(cmp, vec_threshold, _CMP_LE_OS);
if (_mm256_testc_si256((__m256i) cmp, vec_one)) {
i += 8;
continue;
}

Usted menciona "la suma es lenta", esto no es exactamente cierto, pero tiene razón, el rendimiento de la multiplicación aumenta con el tiempo en arquitecturas recientes.

Las latencias y dependencias de multiplicación son la clave. FMA tiene un rendimiento de 1 ciclo y una latencia de 5 ciclos. La ejecución de instrucciones independientes de FMA puede superponerse.

Las adiciones basadas en el resultado de una multiplicación obtienen el impacto de latencia completa.

Por lo tanto, debe romper estas dependencias inmediatas haciendo "unión de código" y calcular 2 puntos en el mismo bucle, y simplemente intercalar el código antes de verificar con IACA lo que sucederá. El siguiente código tiene 2 conjuntos de variables (con el sufijo 0 y 1 para X0 = X0 ^ 2 + C0, X1 = X1 ^ 2 + C1) y comienza a llenar los agujeros de FMA

for (j = 0; j < 8; j++) {
Xrm0 = _mm256_mul_ps(Xre0, Xim0);
Xrm1 = _mm256_mul_ps(Xre1, Xim1);
Xtt0 = _mm256_fmsub_ps(Xim0, Xim0, Cre);
Xtt1 = _mm256_fmsub_ps(Xim1, Xim1, Cre);
Xrm0 = _mm256_add_ps(Xrm0, Xrm0);
Xrm1 = _mm256_add_ps(Xrm1, Xrm1);
Xim0 = _mm256_add_ps(Cim0, Xrm0);
Xim1 = _mm256_add_ps(Cim1, Xrm1);
Xre0 = _mm256_fmsub_ps(Xre0, Xre0, Xtt0);
Xre1 = _mm256_fmsub_ps(Xre1, Xre1, Xtt1);
} // for

Para resumir,

• puedes reducir a la mitad la cantidad de instrucciones en tu ciclo crítico
• puede agregar instrucciones más independientes y aprovechar el alto rendimiento frente a la baja latencia de las multiplicaciones y multiplicar y sumar fusionados.

Answer 3

Para responder a mi tercera pregunta, encontré una solución más rápida para la suma doble-doble. Encontré una definición alternativa en el documento Implementación de operadores float-float en hardware de gráficos .

Theorem 5 (Add22 theorem) Let be ah+al and bh+bl the float-float arguments of the following
algorithm:
Add22 (ah ,al ,bh ,bl)
1 r = ah ⊕ bh
2 if | ah | ≥ | bh | then
3 s = ((( ah ⊖ r ) ⊕ bh ) ⊕ b l ) ⊕ a l
4 e l s e
5 s = ((( bh ⊖ r ) ⊕ ah ) ⊕ a l ) ⊕ b l
6 ( rh , r l ) = add12 ( r , s )
7 return (rh , r l)

Así es como implementé esto (pseudocódigo):

static inline doubledoublen add22(doubledoublen const &a, doubledouble const &b) {
doublen aa,ab,ah,bh,al,bl;
booln mask;
aa = abs(a.hi); //_mm256_and_pd
ab = abs(b.hi);
mask = aa >= ab; //_mm256_cmple_pd
// z = select(cut,x,y) is a SIMD version of z = cut ? x : y;
ah = select(mask,a.hi,b.hi); //_mm256_blendv_pd
bh = select(mask,b.hi,a.hi);
al = select(mask,a.lo,b.lo);
bl = select(mask,b.lo,a.lo);
doublen r, s;
r = ah + bh;
s = (((ah - r) + bh) + bl ) + al;
return two_sum(r,s);
}

Esta definición de Add22 utiliza 11 adiciones en lugar de 20, pero requiere un código adicional para determinar si |ah| >= |bh| |ah| >= |bh| . Aquí hay una discusión sobre cómo implementar las funciones SIMD minmag y maxmag . Afortunadamente, la mayoría del código adicional no usa el puerto 1. Ahora solo 12 instrucciones van al puerto 1 en lugar de 20.

Aquí hay un formulario de análisis de rendimiento IACA para el nuevo Add22

Throughput Analysis Report
--------------------------
Block Throughput: 12.05 Cycles Throughput Bottleneck: Port1
Port Binding In Cycles Per Iteration:
---------------------------------------------------------------------------------------
| Port | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 |
---------------------------------------------------------------------------------------
| Cycles | 0.0 0.0 | 12.0 | 2.5 2.5 | 2.5 2.5 | 2.0 | 10.0 | 0.0 | 2.0 |
---------------------------------------------------------------------------------------
| Num Of | Ports pressure in cycles | |
| Uops | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 | |
---------------------------------------------------------------------------------
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | | vmovapd ymm3, ymmword ptr [rip]
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | | vmovapd ymm0, ymmword ptr [rdx]
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | | vmovapd ymm4, ymmword ptr [rsi]
| 1 | | | | | | 1.0 | | | | vandpd ymm2, ymm4, ymm3
| 1 | | | | | | 1.0 | | | | vandpd ymm3, ymm0, ymm3
| 1 | | 1.0 | | | | | | | CP | vcmppd ymm2, ymm3, ymm2, 0x2
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | | vmovapd ymm3, ymmword ptr [rsi+0x20]
| 2 | | | | | | 2.0 | | | | vblendvpd ymm1, ymm0, ymm4, ymm2
| 2 | | | | | | 2.0 | | | | vblendvpd ymm4, ymm4, ymm0, ymm2
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | | vmovapd ymm0, ymmword ptr [rdx+0x20]
| 2 | | | | | | 2.0 | | | | vblendvpd ymm5, ymm0, ymm3, ymm2
| 2 | | | | | | 2.0 | | | | vblendvpd ymm0, ymm3, ymm0, ymm2
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm3, ymm1, ymm4
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm2, ymm1, ymm3
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm1, ymm2, ymm4
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm1, ymm1, ymm0
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm0, ymm1, ymm5
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm2, ymm3, ymm0
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm1, ymm2, ymm3
| 2^ | | | | | 1.0 | | | 1.0 | | vmovapd ymmword ptr [rdi], ymm2
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm0, ymm0, ymm1
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm1, ymm2, ymm1
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm3, ymm3, ymm1
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm0, ymm3, ymm0
| 2^ | | | | | 1.0 | | | 1.0 | | vmovapd ymmword ptr [rdi+0x20], ymm0

y aquí está el análisis de rendimiento del antiguo

Throughput Analysis Report
--------------------------
Block Throughput: 20.00 Cycles Throughput Bottleneck: Port1
Port Binding In Cycles Per Iteration:
---------------------------------------------------------------------------------------
| Port | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 |
---------------------------------------------------------------------------------------
| Cycles | 0.0 0.0 | 20.0 | 2.0 2.0 | 2.0 2.0 | 2.0 | 0.0 | 0.0 | 2.0 |
---------------------------------------------------------------------------------------
| Num Of | Ports pressure in cycles | |
| Uops | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 | |
---------------------------------------------------------------------------------
| 1 | | | 1.0 1.0 | | | | | | | vmovapd ymm0, ymmword ptr [rsi]
| 1 | | | | 1.0 1.0 | | | | | | vmovapd ymm1, ymmword ptr [rdx]
| 1 | | | 1.0 1.0 | | | | | | | vmovapd ymm3, ymmword ptr [rsi+0x20]
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm4, ymm0, ymm1
| 1 | | | | 1.0 1.0 | | | | | | vmovapd ymm5, ymmword ptr [rdx+0x20]
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm2, ymm4, ymm0
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm1, ymm1, ymm2
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm2, ymm4, ymm2
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm0, ymm0, ymm2
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm2, ymm0, ymm1
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm1, ymm3, ymm5
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm6, ymm1, ymm3
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm5, ymm5, ymm6
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm6, ymm1, ymm6
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm1, ymm2, ymm1
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm3, ymm3, ymm6
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm2, ymm4, ymm1
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm3, ymm3, ymm5
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm4, ymm2, ymm4
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm1, ymm1, ymm4
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm0, ymm1, ymm3
| 1 | | 1.0 | | | | | | | CP | vaddpd ymm1, ymm2, ymm0
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm2, ymm1, ymm2
| 2^ | | | | | 1.0 | | | 1.0 | | vmovapd ymmword ptr [rdi], ymm1
| 1 | | 1.0 | | | | | | | CP | vsubpd ymm0, ymm0, ymm2
| 2^ | | | | | 1.0 | | | 1.0 | | vmovapd ymmword ptr [rdi+0x20], ymm0

Una mejor solución sería si hubiera tres instrucciones de modo de redondeo simple de operando además de FMA. Me parece que debería haber instrucciones de modo de redondeo único para

a + b + c
a * b + c //FMA - this is the only one in x86 so far
a * b * c