[Numpy-discussion] Rookie problems - Why is C-code much faster?

[Mads Ipsen]

On Tue, 21 Feb 2006, Tim Hochberg wrote:

> Mads Ipsen wrote:
>
> >On Tue, 21 Feb 2006, Tim Hochberg wrote:
> >
> >
> >
> >>This all makes perfect sense, but what happended to box? In your
> >>original code there was a step where you did some mumbo jumbo and box
> >>and rint. Namely:
> >>
> >>
> >
> >It's a minor detail, but the reason for this is the following
> >
> >Suppose you have a line with length of box = 10 with periodic boundary
> >conditions (basically this is a circle). Now consider two points x0 = 1
> >and x1 = 9 on this line. The shortest distance dx between the points x0
> >and x1 is dx = -2 and not 8. The calculation
> >
> >  dx  = x1 - x0                 ( = +8)
> >  dx -= box*rint(dx/box)        ( = -2)
> >
> >will give you the desired result, namely dx = -2. Hope this makes better
> >sense. Note that fmod() won't work since
> >
> >  fmod(dx,box) = 8
> >
> >
> I think you could use some variation like "fmod(dx+box/2, box) - box/2"
> but rint seems better.
>
> >Part of my original post was concerned with the fact, that I initially was
> >using around() from numpy for this step. This was terribly slow, so I made
> >some custom changes and added rint() from the C-math library to the numpy
> >module, giving a speedup factor about 4 for this particular line in the
> >code.
> >
> >Best regards // Mads
> >
> >
> >
> OK, that all makes sense. You might want to try the following, which
> factors out all the divisions and half the multiplies by box and
> produces  several fewer temporaries. Note I replaced x**2 with x*x,
> which for the moment is much faster (I don't know if you've been
> following the endless yacking about optimizing x**n, but x**2 will get
> fast eventually). Depending on what you're doing with r2, you may be
> able to avoid the last multiple by box as well.
>
>
>      # Loop over all particles
>      xbox = x/box
>      ybox = y/box
>      for i in range(n-1):
>         dx = xbox[i+1:] - xbox[i]
>         dy = ybox[i+1:] - ybox[i]
>         dx -= rint(dx)
>         dy -= rint(dy)
>         r2 = (dx*dx + dy*dy)
>         r2 *= box
>
>
>
> Regards,
>
> -tim
>

Thanks Tim,

I am only a factor 2.5 slower than the C loop now, thanks to your
suggestions.

// Mads