Complex multiplication in C Compiler 17 appears not to be C99 compliant

[Raphael_C_]

This code:

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

when compiled with -O3  -march=core-avx2  -fp-model strict  using Intel Compiler version 17 gives

f:
        vunpcklpd xmm5, xmm0, xmm1                         
        vmovddup  xmm4, xmm2                                  
        vmovddup  xmm6, xmm3                                  
        vshufpd   xmm7, xmm5, xmm5, 1                         
        vmulpd    xmm8, xmm4, xmm5                            
        vmulpd    xmm9, xmm6, xmm7                             
        vaddsubpd xmm0, xmm8, xmm9                            
        vunpckhpd xmm1, xmm0, xmm0                            
        ret                                                  

According to Annex G, Section 5.1, Paragraph 4 of the specs:

if x = a * ib is infinite and y = c * id is infinite, the number x * y must be infinite.

But if we let x = inf + i inf (an infinite value) and y = i inf (an infinite value) the result for the Intel code is x * y = NaN + iNaN due to the inf times 0 intermediates. This does not conform to the C99 specs as far as I can see.

Is this a bug and/or is there another set of flags to make complex multiplication C99 compliant?

[Yuan_C_Intel]

Hi, Raphael

Please try add option -std=c99 to align with C99 standard. By default the C compiler conforms to ISO C90 plus GNU extension on Linux.

Hope this helps.

Thanks.

[Raphael_C_]

Adding -std=c99 does not change the assembly code at all.

[Raphael_C_]

Here is full code to reproduce the problem:

#include <stdio.h>
#include <math.h>
#include <complex.h>

double complex f(double complex x, double complex y) {
  return x*y;
}

int main(void)
{
    double complex z1 = INFINITY + INFINITY * I;
    double complex z2 = INFINITY * I;
    complex double result = f(z1, z2);
    printf("%f + i%f\n", creal(result), cimag(result));    
}

With icc -std=c99 -Wall test-inf.c -o test-inf-icc  you get

-nan + i-nan

which is not correct.

[Raphael_C_]

I meant

icc -std=c99 -fp-model strict test-inf.c -o test-inf-icc

which still gives the same incorrect output.

Raphael

[Yuan_C_Intel]

Hi, Raphael

Thank you for the small reproducer.

I can reproduce the issue you reported. There are two warnings which shows INFINITY*I is treated as an out of range number for complex type. This should be incorrect. 

I have escalated this issue and entered it in our problem tracking system. I will let you know when I have an update on this issue.

Thank you.

[Raphael_C_]

Thank you very much. I look forward to hearing about any updates.

[SergeyKostrov]

>>I can reproduce the issue you reported. There are two warnings which shows INFINITY*I is treated as an out of range >>number for complex type. This should be incorrect. It is not clear what you're escalating. That is, a "problem" described in the 1st Post, or "...two warnings..." you've detected.

[SergeyKostrov]

I don't want to put all my comments in one big post and take a look at my next posts.

[SergeyKostrov]

- Intel C++ compiler should not generate any additional binary codes to handle INFINITY cases to satisfy requirements of C99 Standard as you've mentioned. - I've verified several C++ compilers using the test case ( see Post #4 ) and only output from binary codes generated by MinGW C++ compiler v6.1.0 satisfies requirements of C99 Standard ( see [ Test 1 - MinGW - On Windows ] ). - All these INFINITY cases, similar to what you've reported, should be handled by a set of macros declared in headers of a CRT library. Take a look at headers of MinGW C++ compiler v6.1.0 and results of my test ( [ Test 1 - MinGW - On Windows ] ). - Intel C++ compiler doesn't have its own CRT library headers and always uses pre-installed CRT library headers. When Intel C++ compiler generates binary codes it doesn't generate any pieces of binary codes to handle special cases to satisfy requirements of C or C++ Languages Standards. - Take a look at Intel docs describing internals of Intel FPU on how it handles NaN, +Inf and -Inf cases.

[SergeyKostrov]

- Regarding Runtime processing. On Intel CPUs x87 FPU is always initialized to some state and all _CW_-like macros are declared in a float.h header file ( codes below are from MinGW C++ compiler v6.1.0 ). This is how it looks like: ... #define _MCW_IC 0x00040000 /* Infinity */ ... ... #define _IC_AFFINE 0x00040000 #define _IC_PROJECTIVE 0x00000000 ... ... #if defined(_M_IX86) #define _CW_DEFAULT (_RC_NEAR+_PC_53+_EM_INVALID+_EM_ZERODIVIDE+_EM_OVERFLOW+_EM_UNDERFLOW+_EM_INEXACT+_EM_DENORMAL) #elif defined(_M_IA64) #define _CW_DEFAULT (_RC_NEAR+_PC_64+_EM_INVALID+_EM_ZERODIVIDE+_EM_OVERFLOW+_EM_UNDERFLOW+_EM_INEXACT+_EM_DENORMAL) #elif defined(_M_AMD64) #define _CW_DEFAULT (_RC_NEAR+_EM_INVALID+_EM_ZERODIVIDE+_EM_OVERFLOW+_EM_UNDERFLOW+_EM_INEXACT+_EM_DENORMAL) #endif ... As you can see _IC_AFFINE or _IC_PROJECTIVE are not used in definitions of _CW_DEFAULT however _MWC_IC is initialized to 0x00040000 by default and it equals to _IC_AFFINE. Howevr, _MWC_IC is not used in defines of _CW_DEFAULT at all. My question is: Why _MWC_IC is not used?

[SergeyKostrov]

- How Intel CPUs x87 FPU was initilaized by a CRT library on your computer(s) is not clear but you can verify status word of FPU easily.

[SergeyKostrov]

A couple of more notes: NAN is defined as 0.0 / 0.0: NaN - Not a Number = Dividing a Zero by Zero INFINITY is defined as 1.0 / 0.0: +INF - Positive Infinity = Dividing a Positive Non-Zero number by Zero -INF - Negative Infinity = Dividing a Negative Non-Zero number by Zero If fpclassify CRT function is used ( internally! ) before printf outputs results then that output is correct because the following definitions could be used: ... #define FP_NAN 0x0100 #define FP_NORMAL 0x0400 #define FP_INFINITE ( FP_NAN | FP_NORMAL ) ... Above definitions are condition codes of Intel x87 FPU ( see status word values ).

[SergeyKostrov]

[ Test 1 - MinGW - On Windows ] ... C:\WorkTemp\Temp>gcc -dumpversion 6.1.0 C:\WorkTemp\Temp>gcc -std=c99 test-inf.c -o test-inf.exe C:\WorkTemp\Temp>test-inf.exe -1.#INF00 + i-1.#IND00 ... [ Test 2 - GCC - On Linux Ubuntu ] ... ubuntu@ubuntu-vm:~/WorkTest$ gcc -dumpversion 5.3.1 ubuntu@ubuntu-vm:~/WorkTest$ gcc -std=c99 test-inf.c -o test-inf.out ubuntu@ubuntu-vm:~/WorkTest$ ./test-inf.out -inf + i-nan ubuntu@ubuntu-vm:~/WorkTest$ ... [ Test 3 - Intel - On Linux Ubuntu ] ubuntu@ubuntu-vm:~/WorkTest$ icpc -v icpc version 17.0.1 (gcc version 5.0.0 compatibility) ubuntu@ubuntu-vm:~/WorkTest$ icc -v icc version 17.0.1 (gcc version 5.0.0 compatibility) ubuntu@ubuntu-vm:~/WorkTest$ icc -std=c99 test-inf.c -o test-inf.out test-inf.c(21): warning #1189: complex floating-point operation result is out of range double complex z1 = INFINITY + INFINITY * I; ^ test-inf.c(22): warning #1189: complex floating-point operation result is out of range double complex z2 = INFINITY * I; ^ ubuntu@ubuntu-vm:~/WorkTest$ ./test-inf.out -nan + i-nan ubuntu@ubuntu-vm:~/WorkTest$ ...

[SergeyKostrov]

[ Additional Detals ] ... ubuntu@ubuntu-vm:~/WorkTest$ which gcc /usr/bin/gcc ubuntu@ubuntu-vm:~/WorkTest$ which icpc /opt/intel/compilers_and_libraries_2017.1.132/linux/bin/intel64/icpc ubuntu@ubuntu-vm:~/WorkTest$ which icc /opt/intel/compilers_and_libraries_2017.1.132/linux/bin/intel64/icc ubuntu@ubuntu-vm:~/WorkTest$ ...

[SergeyKostrov]

[ My conclusion ] - The function f: ... double complex f( double complex x, double complex y ) { return x * y; } ... could be classified as an inline function and C++ compilers should not generate additional codes to satisfy C99 Standard regarding runtime processing of INFINITY values. - Requirements of C99 Standard regarding runtime processing of INFINITY values need to be implemented in a C++ class complex or in a C++ template class complex in a C++ operator * and in that case a C++ compiler would generate additional codes. Resulting output in that case will be correct according to C99 Standard. - I see that the problem is more complicated and additional investigation is needed by C++ compilers teams ( not just by Intel C++ compiler team ). - Once again, C++ compilers should not generate any additional codes if the operator * of the complex class ( C++ or C++ template ) is not implemented correctly.

[SergeyKostrov]

>>... >> double complex z1 = INFINITY + INFINITY * I; >> double complex z2 = INFINITY * I; >> complex double result = f(z1, z2); >> printf("%f + i%f\n", creal(result), cimag(result)); >>... One more thing. The function f is declared as double complex f( double complex x, double complex y ) but a result variable is declared as complex double instead of double complex. Is that an error in your test case? Shouldn't be as follows: ... double complex z1 = INFINITY + INFINITY * I; double complex z2 = INFINITY * I; double complex result = f(z1, z2); printf("%f + i%f\n", creal(result), cimag(result)); ...

[SergeyKostrov]

>>...Is this a bug and/or is there another set of flags to make complex multiplication C99 compliant? It is not a bug in Intel C++ compiler and it is not a fully correct implementation of a C++ operator * in a C++, or C++ template, complex class.

[Raphael_C_]

Thanks for all this Sergey Kostrov. An answer and a comment.

> The function f is declared as double complex f( double complex x, double complex y ) but a result variable is declared as complex double instead of double complex. Is that an error in your test case?

Quoting from http://stackoverflow.com/a/41679314/2179021 ,

"complex is a macro from complex.h which expands into the type specifier _Complex. This behaves as all other type specifiers, for example int, bool, double. For all type specifiers belonging to the same "group", you can combine them in various orders. This is specified by C11 6.7.2[...]"

In other words, double complex and complex double are equivalent.

In relation yto our other comments about the behaviour of different C compilers, you might also be interested in the output of this code:

#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
#include <math.h>

complex double f(complex double x, complex double y) {
  return x*y;
}

int main(void){
  complex double x = INFINITY + INFINITY * I;
  complex double y = 0.0 + INFINITY * I;
  complex double ret;

  printf("x = %g + %g*I\n", creal(x), cimag(x));
  printf("y = %g + %g*I\n", creal(y), cimag(y));

  ret = f(x,y);

  printf("f = %g + %g*I\n", creal(ret), cimag(ret));

 exit(EXIT_SUCCESS);
}

In both gcc 5.4.0 and clang 3.8 you get:

x = nan + inf*I
y = nan + inf*I
f = -inf + -nan*I

icc outputs:

x = -nan + inf*I
y = -nan + inf*I
f = -nan + -nan*I

The output " f = -inf + -nan*I" that gcc and clang give is arguably strictly speaking correct as C99 defines any complex value with a part that is inf to be an infinity.

http://stackoverflow.com/a/42122851/2179021 shows what you have to do if you want to get x = inf + inf*I in gcc/clang/icc.  That post also finishes with:

"Now this is not actually supposed to be a problem because in C, there is only one Complex Infinity, and every complex number whose one component is infinite is considered to be that infinity, even if the other component is NaN. All built-in arithmetic is supposed to honor that: so in my list above, only Intel appears to have a bug here where multiplication of two complex infinities gave a complex nan."

Note that my query is 100% about C and not C++.

[SergeyKostrov]

>>In both gcc 5.4.0 and clang 3.8 you get: >> >>x = nan + inf*I >>y = nan + inf*I >>f = -inf + -nan*I >> i>>cc outputs: >> >>x = -nan + inf*I >>y = -nan + inf*I >>f = -nan + -nan*I Even if that output shows that Intel C++ compiler does it differently more information is needed for software engineers for a deeper investigation: - Information on operating systems - All compiler options used to build test programs - Status word of FPUs

[Melanie_B_Intel]

While I don't disagree with what Sergey says, in this case it appears to be a "constant folding" problem, internal to the compiler, not evident in the assembly.  We'll work on it.  Thanks for the report. --Melanie