[Numpy-discussion] Rookie problems - Why is C-code much faster?
Tim Hochberg
tim.hochberg at cox.net
Tue Feb 21 14:24:07 CST 2006
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
More information about the Numpy-discussion
mailing list