Bresenham's algorithm would be a good way to quickly generate each
square. My DADS entry has links to implementations and explanations
that might help.
http://www.nist.gov/dads/HTML/bresenham.html

The problem is that Bresenham's algorithm is NOT a good way as far as it
would not generate all the squares crossed. If you have a look at the
picture there:
http://en.wikipedia.org/wiki/Image:Bresenham.svg
you see that there are some squares where the line pass, but they are
not grayed. And as this is for collision detection purpose, I can not
miss any square!
Thanks anyway.

Are you using Open-GL? It has built-ins for this. gluUnProject()
sounds like exactly what you're doing, here.
Hi

I have a 2D grid of n*n squares of the same size s. The first one starts
at (x0,y0). So the square nx,ny spans from (x0+nx*s,y0+ny*s) to
(x0+(nx+1)*s,y0+(ny+1)*s).
Then I have a half-line whose equation is ax+by+c=0 starting at (x1,y1).
What would be a quick way to know the squares crossed from (x1,y1)
further?

When I did this in the past, I just inverted the matrix that I was using
to draw the grid and multiplied it into the mouse position. This gave
me the square index in floating point, where the whole portion was the
square being clicked on and the fraction was how far across the square I
was clicking.

Granted, this was for a flat plane, but it worked well for any
projection matrix. If you're trying to click on squares that have been
manipulated in a heightfield, you might need to do something trickier.