Numeric: fromfunction() performance

[Pearu Peterson]

On 15 Oct 2002, jsaul wrote:

> Hi there,
> 
> I have a code in which an array is created 'on the fly' using
> Numeric's fromfunction() command:
> 
>     import Numeric
> 
>     x1,x2,nf = 1.08E-8, 6.28E-3, 100000
>     q = Numeric.fromfunction (lambda i,k: - x1*i*i - 1j*i*x2 , (nf,1))
>     q.shape = (nf,) # make it a one-dimensional array (arhhh)
>     q = Numeric.exp(q)
> 
> Surprisingly, the fromfunction() statement takes quite long to
> execute, while exp() is very quick, because it's coded in C.
> Is there any way to improve the speed of the fromfunction() call?
> 
> And, is there a way to avoid that 'shape' statement?
> fromfunction() apparently refuses to produce a one-dimensional
> array of length nf, so I currently achieve that by computing a
> nfX1 matrix, which I then re-shape to a 1-D array. Not quite
> an elegant way...
> 
> What I actually want is a (quicker) construct corresponding to
> 
>     for i in range(0, nf):
>         q[i] = cmath.exp(-x1*i*i -1j*i*x2)
> 
> Which is much too slow.

It seems that you are trying to evaluate genus 1 Riemann theta function
that, by definition, reads (I hope you can read Latex)

  \theta(\phi) = \sum_{n=-\infty}^{\infty} \exp(-b*n*n+sqrt(-1)*n*\phi)

It can be rewritten in the following form (after using Poisson
formula)

  \theta(\phi) = \sqrt{2\pi/b} \sum_{n=-\infty}^{\infty} 
                     \exp(-(\phi-2\pi n)^2/(4 b))

Now, if you want to evaluate \theta(\phi) numerically, then
the first formula is most efficient if b<2\pi and use the second formula
if b>2\pi. Also note that when truncating the series, it is enough to
use only terms that have
  |n| <= ceil(2.2\sqrt(p/b))  for the first formula
  |n| <= ceil(0.2\sqrt(p b))  for the second formula

where p is the number of digits you want to be correct, take p=16 if you
are working on 32-bit computer. If you are using the right formula for
the given b, then the max(|n|) never exceeds 4 (c.f. your nf=100000).

Another hint is that when you evaluate series numerically then you should
sum starting from the terms with smallest value.

I have a Fortran program that evaluates \theta function and its
derivatives upto 20th order. It can be called from Python using f2py. 
If you are interested, I can send the sources to you.

HTH,
	Pearu