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

[Tim H. <tim...@co...>]

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