Calculating complex exponential numbers using OpenMP in Cython

Question

I am trying to create a function that can either beat numexpr or perform comparably for the vectorized mathematical operation k * exp(ee), where k and ee are 1D arrays of complex numbers. I had tried using OpenMP in Cython directly, but I have resorted to creating a DLL in C++ and then wrapping with Cython. Here are the components (of current interest is the function multi_mlt_exp_c4). I am using Windows 10, 64 bit, Intel processor, and Python 3.x.

C++ DLL, Header:

#define EXTERN_DLL_EXPORT extern "C" __declspec(dllexport)
#include <complex>
#include <stdlib.h>

EXTERN_DLL_EXPORT void cexp_scale_c4(std::complex<float> k, std::complex<float> ee,
                                     std::complex<float>*res);
EXTERN_DLL_EXPORT void multi_mlt_exp_c4(std::complex<float>* k, std::complex<float>* ee,
                                        int sz, std::complex<float>* res, int threads);
EXTERN_DLL_EXPORT void multi_mlt_exp_c8(std::complex<double>* k, std::complex<double>* ee,
                                        int sz, std::complex<double>* res, int threads);
EXTERN_DLL_EXPORT void serial_mlt_exp_c4(std::complex<float>* k, std::complex<float>* ee,
                                         int sz, std::complex<float>* res);

C++ DLL, Source

#include "stdafx.h"
#include <stdio.h>
#include <omp.h>
#include "C_Fun_Lib.h"

void cexp_scale_c4(std::complex<float> k, std::complex<float> ee, std::complex<float>*res)
{
    *res = k*exp(ee);
}

void multi_mlt_exp_c4(std::complex<float>* k, std::complex<float>* ee, int sz, std::complex<float>* res, int threads)
{
    #pragma omp parallel num_threads(threads)
    {
        #pragma omp for
        for (int i = 0; i < sz; i++) res[i] = k[i] * exp(ee[i]);
    }
}

void multi_mlt_exp_c8(std::complex<double>* k, std::complex<double>* ee, int sz, std::complex<double>* res, int threads)
{
    #pragma omp parallel num_threads(threads)
    {
        #pragma omp for
        for (int i = 0; i < sz; i++) res[i] = k[i] * exp(ee[i]);
    }
}

void serial_mlt_exp_c4(std::complex<float>* k, std::complex<float>* ee, int sz, std::complex<float>* res)
{
    for (int i = 0; i < sz; i++) res[i] = k[i] * exp(ee[i]);
}

Here are the relevant Python / Cython components:

PXD

cdef extern from "C_Fun_Lib.h":
    void multi_mlt_exp_c4(float complex* k, float complex* ee, int sz,
                          float complex* res, int threads);
    void serial_mlt_exp_c4(float complex* k, float complex* ee, int sz,
                           float complex* res);
    void multi_mlt_exp_c8(double complex* k, double complex* ee, int sz,
                          double complex* res, int threads);

PYX

cimport csample
import numpy as np
cimport numpy as np

# Import some functionality from Python and the C stdlib
from cpython.pycapsule cimport *

from libc.stdlib cimport malloc, free

def mlt_exp_c4(np.ndarray[np.complex64_t, ndim=1] k,
               np.ndarray[np.complex64_t, ndim=1] ee, int sz,
               np.ndarray[np.complex64_t, ndim=1] res, int threads):
    csample.multi_mlt_exp_c4(&k[0], &ee[0], sz, &res[0], threads)

def mlt_exp_c4_serial(np.ndarray[np.complex64_t, ndim=1] k,
               np.ndarray[np.complex64_t, ndim=1] ee, int sz,
               np.ndarray[np.complex64_t, ndim=1] res):
    csample.serial_mlt_exp_c4(&k[0], &ee[0], sz, &res[0])

def mlt_exp_c8(np.ndarray[np.complex128_t, ndim=1] k,
               np.ndarray[np.complex128_t, ndim=1] ee, int sz,
               np.ndarray[np.complex128_t, ndim=1] res, int threads):
    csample.multi_mlt_exp_c8(&k[0], &ee[0], sz, &res[0], threads)

Setup.py

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
import numpy as np

ext_modules = [
    Extension('sample',
              ['sample.pyx'],
              language="c++",
              libraries=['C_Fun_Lib'],
              library_dirs=['.'])
    ]

setup(
    name = 'sample',
    cmdclass = {'build_ext': build_ext},
    ext_modules = ext_modules,
    include_dirs=[np.get_include()]
)

And then finally here is my timing code and results. I am specifically trying to beat Numexpr because, currently, that application up-casts complex64 arrays to complex128 arrays, and I am trying to save memory. Writing my own code gives me control over the memory, but I can't quite match the speed of Numexpr.

import sample
import numpy as np
import numexpr as ne
import time

k = np.ones(int(2.5e8), dtype='complex64')*1.1234 + 1.1234j
ee = k
sz = k.sz
res = np.zeros(int(2.5e8), dtype='complex64')  # Results

tic = time.time()
sample.mlt_exp_c4(k, ee, sz, res, 8)
toc = time.time()

print('Elasped MP: %f' % (toc - tic))

# Reset results
res = np.zeros(int(2.5e8), dtype='complex64')

tic = time.time()
res = ne.evaluate('k*exp(ee)')
toc = time.time()

print('Elasped NE: %f' % (toc - tic))

>>>>>>>
Elasped MP: 2.292317
Elasped NE: 1.184981

Things I have tried

First, out of all the things I've tried, this is my current best. I have experimented to see if using complex128s makes a difference (since I'm on a 64 bit machine). It does not. I have tried doing the complex exponential and multiplication directly and using Euler's formula instead of using the C++ built in complex type and functions... that is slower. I have tried using Cython code directly with OpenMP support... that is slower.

It may very well be that trying to match the highly optimized Numexpr code is a fool's errand. I just wanted you all to review this code to see if there were any obvious speed ups that I am missing.