x = a + i*b
y = c + i*d

z = x*y = (a*c - b*d) + i*(a*d + b*c)

#include <algorithm>
#include <chrono>
#include <complex>
#include <iostream>
#include <random>
#include <vector>


template<typename T>
const char*
prettyType();


template<>
const char*
prettyType<float>()
{
    return "float";
}


template<>
const char*
prettyType<double>()
{
    return "double";
}


#ifdef direct
#warning "use C complex numbers in sum"

extern "C" {
    __complex__ float __mulsc3(float, float, float, float);
    __complex__ double __muldc3(double, double, double, double);
}


template<typename T>
inline
std::complex<T>
sum(const std::vector<std::complex<T>>&);


template<>
inline
std::complex<float>
sum<float>(const std::vector<std::complex<float>>& v)
{
    const size_t length = v.size();

    __complex__ float sum = 0.0;
    for (unsigned int i = 0; i < length; ++i)
        for (unsigned int j = 0; j < length; ++j)
            sum += __mulsc3(v[i].real(), v[i].imag(), v[j].real(), v[j].imag());

    return std::complex<float>(__real__ sum, __imag__ sum);
}


template<>
inline
std::complex<double>
sum<double>(const std::vector<std::complex<double>>& v)
{
    const size_t length = v.size();

    __complex__ double sum = 0.0;
    for (unsigned int i = 0; i < length; ++i)
        for (unsigned int j = 0; j < length; ++j)
            sum += __muldc3(v[i].real(), v[i].imag(), v[j].real(), v[j].imag());

    return std::complex<double>(__real__ sum, __imag__ sum);
}

#else

template<typename T>
inline
std::complex<T>
sum(const std::vector<std::complex<T>>& v)
{
    const size_t length = v.size();

    std::complex<T> sum;
    for (unsigned int i = 0; i < length; ++i)
        for (unsigned int j = 0; j < length; ++j)
            sum += v[i] * v[j];

    return sum;
}

#endif


template<typename T, const size_t length = 100>
void
test()
{
    std::mt19937 rndgen(123456);
    std::uniform_real_distribution<T> dist(0.0, 12345.6789);

    std::vector<std::complex<T>> v(length);
    std::generate(v.begin(), v.end(), [&dist,&rndgen]() { const T r = dist(rndgen); const T i = dist(rndgen); return std::complex<T>(r, i); });

    auto start = std::chrono::system_clock::now();
    std::complex<T> s = sum(v);
    auto end = std::chrono::system_clock::now();

    const std::chrono::duration<double> time(end-start);
    std::cout << "sum for " << prettyType<T>() << " took " << time.count() << " s: " << s << std::endl;
}


int
main(int argc, char* argv[])
{
    test<float>();
    test<double>();
}

$ g++ -std=c++11 -O3 bla.cc -o bla
$ ./bla
sum for float took 3.15104 s: (4.42991e+12,4.53069e+15)
sum for double took 0.994963 s: (3.35809e+12,3.08329e+16)

$ g++ -std=c++11 -O3 bla.cc -o bla -ffast-math
$ ./bla
sum for float took 0.431998 s: (4.42991e+12,4.53069e+15)
sum for double took 0.445863 s: (3.35809e+12,3.08329e+16)

$ clang++ -O3 -std=c++11 bla.cc -o bla
$ ./bla
sum for float took 1.79012 s: (4.42991e+12,4.53069e+15)
sum for double took 2.4145 s: (3.35809e+12,3.08329e+16)

z=x*y:

value_type e = x.real()*y.real - x.imag()*y.imag()
value_type f = x.real()*y.imag + x.imag()*y.real()
if (isnan(e) || isnan(f))
z = __mul?c3(x.real(), x.imag(), y.real(), y.imag());
else
z = e + i*f

$ clang++ -O3 -std=c++11 bla.cc -o bla -ffast-math
$ ./bla
sum for float took 1.36594 s: (4.55313e+12,3.77658e+15)
sum for double took 0.579735 s: (3.35809e+12,3.08329e+16)

$ clang++ -O3 -std=c++11 bla.cc -o bla -Ddirect
$ ./bla
sum for float took 2.18363 s: (4.42991e+12,4.53069e+15)
sum for double took 2.34664 s: (3.35809e+12,3.08329e+16)

$ g++ -std=c++11 -O3 bla.cc -o bla -Ddirect
$ ./bla
sum for float took 3.33212 s: (4.42991e+12,4.53069e+15)
sum for double took 0.945766 s: (3.35809e+12,3.08329e+16)

    std::generate(v.begin(), v.end(), [&dist,&rndgen]() { return std::complex<T>(dist(rndgen), dist(rndgen)); });

    std::generate(v.begin(), v.end(), [&dist,&rndgen]() { const T r = dist(rndgen); const T i = dist(rndgen); return std::complex<T>(r, i); });