[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
>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 



More information about the Numpy-discussion mailing list