Pyrex - Numeric demo doesn't build in 0.7.2

[Greg Ewing (using news.cis.dfn.de)]

Fernando Perez wrote:
>             mat[k,l] = sqrt_N_inv * \
>                        exp(I*(alpha*(k*k-k*l+l*l) + kap_al*sin(alpha*l)))

The array indexing and the calls to exp and sin are all
being done as Python operations.

For maximum speed, you need to remove all operations
involving Python objects from the inner loop. This will
involve declaring appropriate functions from the
C math library (e.g. sin) so they will be called directly,
and using the techniques from the Numeric demo to poke
data directly into the Numeric array.

You'll also want to avoid dealing with Python complex
number objects inside the loop. This will mean using
some suitable C library for doing complex exp() etc.,
or handling the real and imaginary components yourself.

Something like this (warning, utterly untested):

   cdef extern from "math.h":
     double sin(double)
     double cos(double)
     double sqrt(double)

   cdef extern from "Numeric/arrayobject.h":
     # stuff from numeric_demo.py

   ctypedef struct ccomplex:
     double r, i

   cdef ccomplex cexp(ccomplex x):
     # something suitable here

   cdef ccomplex cmul(ccomplex x, ccomplex y):
     # something suitable here

   ...
   cdef ccomplex x, y
   ...
       c.r = kap_al*sin(alpha*l)
       c.i = alpha*(k*k-k*l+l*l)
       (<ccomplex*>mat.data)[N * k + l] = cmul(sqrt_N_inv, cexp(c))


-- 
Greg Ewing, Computer Science Dept,
University of Canterbury,	
Christchurch, New Zealand
http://www.cosc.canterbury.ac.nz/~greg