resueltos - ¿ICC cumple con las especificaciones de C99 para la multiplicación de números complejos?

multiplicacion de numeros complejos ejercicios resueltos

Question

Considere este código simple:

#include <complex.h>
complex float f(complex float x) {
  return x*x;
}

Si lo compila con -O3 -march=core-avx2 -fp-model strict utilizando el compilador de Intel, obtendrá:

f:
        vmovsldup xmm1, xmm0                                    #3.12
        vmovshdup xmm2, xmm0                                    #3.12
        vshufps   xmm3, xmm0, xmm0, 177                         #3.12
        vmulps    xmm4, xmm1, xmm0                              #3.12
        vmulps    xmm5, xmm2, xmm3                              #3.12
        vaddsubps xmm0, xmm4, xmm5                              #3.12
        ret 

Este es un código mucho más simple que el que obtienes de gcc y clang y también mucho más simple que el código que encontrarás en línea para multiplicar números complejos. No aparece, por ejemplo, explícitamente para tratar con NaN complejos o infinitos.

¿Cumple este conjunto las especificaciones para la multiplicación compleja de C99?

Answer 1

El código no es conforme.

El Anexo G, Sección 5.1, Párrafo 4 dice

Los operadores * y / satisfacen las siguientes propiedades infinitas para todos los operandos reales, imaginarios y complejos:

- si un operando es un infinito y el otro es un número finito distinto de cero o un infinito, entonces el resultado del operador * es un infinito;

Entonces, si z = a * i b es infinito y w = c * i d es infinito, el número z * w debe ser infinito.

El mismo anexo, Sección 3, Párrafo 1 define lo que significa que un número complejo sea infinito:

Un valor complejo o imaginario con al menos una parte infinita se considera un infinito (incluso si su otra parte es una NaN).

Entonces z es infinito si a o b son.
Esta es de hecho una elección sensata ya que refleja el marco matemático ¹ .

Sin embargo, si permitimos que z = ∞ + i ∞ (un valor infinito) yw = i ∞ (y un valor infinito) el resultado para el código de Intel es z * w = NaN + i NaN debido a los intermedios ∞ · 0 ² .

Esto basta para etiquetarlo como no conforme.

Podemos confirmar esto echando un vistazo a la nota al pie de la primera cita (la nota al pie no se informó aquí), menciona la directiva pragma CX_LIMITED_RANGE .

En la sección 7.3.4, se lee el párrafo 1.

Las fórmulas matemáticas habituales para la multiplicación compleja, la división y el valor absoluto son problemáticas debido a su tratamiento de los infinitos y debido a un exceso de desbordamiento y desbordamiento. El pragma CX_LIMITED_RANGE se puede usar para informar a la implementación que (donde el estado está "encendido") las fórmulas matemáticas habituales [que producen NaN] son ​​aceptables.

Aquí, el comité estándar está tratando de aliviar la enorme cantidad de trabajo para la compleja multiplicación (y división).
De hecho, GCC tiene una bandera para controlar este comportamiento :

-fcx-limited-range
Cuando está habilitada, esta opción indica que no se necesita un paso de reducción de rango al realizar una división compleja.

Además, no se verifica si el resultado de una multiplicación o división compleja es NaN + I * NaN, con un intento de rescatar la situación en ese caso.

El valor predeterminado es -fno-cx-limited-range , pero está habilitado por -ffast-math .
Esta opción controla la configuración predeterminada del pragma ISO C99 CX_LIMITED_RANGE .

Es solo esta opción la que hace que GCC genere código lento y verificaciones adicionales , sin él, el código que genera tiene los mismos defectos que el de Intel (traduje la fuente a C ++)

f(std::complex<float>):
        movq    QWORD PTR [rsp-8], xmm0
        movss   xmm0, DWORD PTR [rsp-8]
        movss   xmm2, DWORD PTR [rsp-4]
        movaps  xmm1, xmm0
        movaps  xmm3, xmm2
        mulss   xmm1, xmm0
        mulss   xmm3, xmm2
        mulss   xmm0, xmm2
        subss   xmm1, xmm3
        addss   xmm0, xmm0
        movss   DWORD PTR [rsp-16], xmm1
        movss   DWORD PTR [rsp-12], xmm0
        movq    xmm0, QWORD PTR [rsp-16]
        ret

Sin él el código es

f(std::complex<float>):
        sub     rsp, 40
        movq    QWORD PTR [rsp+24], xmm0
        movss   xmm3, DWORD PTR [rsp+28]
        movss   xmm2, DWORD PTR [rsp+24]
        movaps  xmm1, xmm3
        movaps  xmm0, xmm2
        call    __mulsc3
        movq    QWORD PTR [rsp+16], xmm0
        movss   xmm0, DWORD PTR [rsp+16]
        movss   DWORD PTR [rsp+8], xmm0
        movss   xmm0, DWORD PTR [rsp+20]
        movss   DWORD PTR [rsp+12], xmm0
        movq    xmm0, QWORD PTR [rsp+8]
        add     rsp, 40
        ret

y la función __mulsc3 es prácticamente la misma que recomienda el estándar C99 para la multiplicación compleja.
Incluye los controles mencionados anteriormente.

¹ Donde el módulo de un número se extiende desde el caso real | z | al complejo ‖z‖, manteniendo la definición de infinito como resultado de límites ilimitados. En pocas palabras, en el plano complejo hay toda una circunferencia de valores infinitos y solo se necesita una "coordenada" para ser infinito para obtener un módulo infinito.

² La situación empeora si recordamos que z = NaN + i ∞ o z = ∞ + i NaN son valores infinitos válidos

Answer 2

Obtengo un código similar, pero no idéntico, de -O2 -march=core-avx2 -ffast-math 3.8 en -O2 -march=core-avx2 -ffast-math : no estoy muy familiarizado con las características de punto flotante de x86 de los últimos días, pero creo que está haciendo el mismo cálculo pero usando diferentes instrucciones para mezclar valores alrededor de los registros.

f:
        vmovshdup       %xmm0, %xmm1    # xmm1 = xmm0[1,1,3,3]
        vaddss  %xmm0, %xmm0, %xmm2
        vmulss  %xmm2, %xmm1, %xmm2
        vmulss  %xmm1, %xmm1, %xmm1
        vfmsub231ss     %xmm0, %xmm0, %xmm1
        vinsertps       $16, %xmm2, %xmm1, %xmm0 # xmm0 = xmm1[0],xmm2[0],xmm1[2,3]
        retq

GCC 6.3, con las mismas opciones, nuevamente parece hacer el mismo cálculo pero baraja los valores de una tercera manera:

f:
        vmovq   %xmm0, -8(%rsp)
        vmovss  -4(%rsp), %xmm2
        vmovss  -8(%rsp), %xmm0
        vmulss  %xmm2, %xmm2, %xmm1
        vfmsub231ss     %xmm0, %xmm0, %xmm1
        vmulss  %xmm2, %xmm0, %xmm0
        vmovss  %xmm1, -16(%rsp)
        vaddss  %xmm0, %xmm0, %xmm0
        vmovss  %xmm0, -12(%rsp)
        vmovq   -16(%rsp), %xmm0
        ret

Sin -ffast-math , ambos compiladores generan códigos sustancialmente diferentes que parecen verificar al menos NaN, al menos.

Concluyo de esto que el compilador de Intel no está generando una multiplicación compleja totalmente compatible con IEEE incluso con -fp-model strict . Puede haber algún otro interruptor de línea de comandos que lo haga generar código totalmente compatible con IEEE.

Si esto califica o no como una violación de C99 depende de si el compilador de Intel está documentado para cumplir con los Anexos F y G (que especifican qué significa para una implementación de C proporcionar aritmética compleja y real compatible con IEEE), y si es , qué opciones de línea de comandos debe dar para obtener el modo de conformidad.